How to preserve stress field when deforming the mesh using ALE.move()

I deformed my original mesh using ALE.move() . When I deform the mesh the function space of the stress field defined using that mesh changes. For my application. I need the stress field on the deformed shape. so I create a new mesh and deform that new mesh instead of the original mesh. since my tensor function space is changed I cannot get the stress field on the deformed object as the deformation of the object is somewhat large. can some one suggest an option.

import dolfin as df
import numpy as np
from Plot import *
#function definitions
def eps(v):
    return df.sym(df.grad(v))
def sigmaE(v, lmbda, mu):
    sigE = lmbda * df.tr(eps(v)) * df.Identity(2) + 2.0 * mu * eps(v)
    return sigE
def sigmaA(v, sigA, lmbda, mu):
    sig = sigA + lmbda * df.tr(eps(v)) * df.Identity(2) + 2.0 * mu * eps(v)
    return sig
def btm_Mid_pt(x,on_boundary,xval=5.0,yval=0.0,tolx=0.02,toly=0.03):
    return not_(on_boundary) and df.near(x[0],xval,tolx) and df.near(x[1],yval,toly)
def tp_mid_pt(x,on_boundary,xval=5.0,yval=5.0,tolx=0.02,toly=0.03):
    return not_(on_boundary) and df.near(x[0],xval,tolx) and df.near(x[1],yval,toly)

def sigA_init(V):
    # geometry
    #long = np.array(V.tabulate_dof_coordinates())[:, 0].max() - np.array(V.tabulate_dof_coordinates())[:, 0].min()
    ep = np.array(V.tabulate_dof_coordinates())[:, 1].max() - np.array(V.tabulate_dof_coordinates())[:, 1].min()
    # Residual stress
    depth1 = 1.5
    sigXX = [1000., 0.]
    sigYY = [0., 0.]
    sigXY = [0., 0.]
    depth = df.Expression('ep - x[1]', degree=1, ep=ep,d=depth1)
    
    sigA_xx = df.Expression('s0 + (s1 - s0) / d1 * d', degree=1, d=depth, d1=depth1, s0=sigXX[0], s1=sigXX[1])
    sigA_yy = df.Expression('s0 + (s1 - s0) / d1 * d', degree=1, d=depth, d1=depth1, s0=sigYY[0], s1=sigYY[1])
    #sigA_yy = df.Constant(sigYY[0])
    sigA_xy = df.Expression('s0 + (s1 - s0) / d1 * d', degree=1, d=depth, d1=depth1, s0=sigXY[0], s1=sigXY[1])
    sigA = df.Expression((('d <= d1 ? sigA_xx : 0.', 'd <= d1 ? sigA_xy : 0.'), \
                          ('d <= d1 ? sigA_xy : 0.', 'd <= d1 ? sigA_yy : 0.')), \
                         degree=1, ep=ep, d=depth, d1=depth1, sigA_xx=sigA_xx, sigA_yy=sigA_yy, sigA_xy=sigA_xy)
    return sigA

def resolution(V, sigA, lmbda, mu,bcs_1):
    # Problem defintion
    du = df.TrialFunction(V)
    u_ = df.TestFunction(V)
    a = df.inner(sigmaE(du, lmbda, mu), eps(u_)) * df.dx
    l = -df.inner(sigA, eps(u_)) * df.dx    
    # problem solving
    u = df.Function(V, name="Displacement")
    df.solve(a == l, u, bcs_1)
    
    return u

model = "plane_stress"
# Material
E = df.Constant(110000.)
nu = df.Constant(0.3)
mu = E / (2 * (1 + nu))
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
if model == "plane_stress":
    lmbda = 2 * mu * lmbda / (lmbda + 2 * mu)
mesh = df.RectangleMesh(df.Point(0,0),df.Point(10,5),100,50,diagonal="right")

#creating fuction spaces

V0 = df.VectorFunctionSpace(mesh, "Lagrange", degree=2)
V0sig = df.TensorFunctionSpace(mesh, 'DG', degree=1)

#creating boundqry conditions
bcs_0=df.DirichletBC(V0.sub(0), df.Constant(0.0),lambda x,on_boundary:tp_mid_pt(x,on_boundary,tolx=0.01,toly=0.01),"pointwise")
bcs_1=df.DirichletBC(V0, df.Constant((0.,0.)),lambda x,on_boundary:btm_Mid_pt(x,on_boundary,tolx=0.01,toly=0.01),"pointwise")
bcs = [bcs_0,bcs_1]
#creating new stress field
sig=sigA_init(V0)
u=resolution(V0, sig, lmbda, mu,bcs)
sig_after = sigmaA(u,sig,lmbda,mu)

#exporting stress field to paraview
u.set_allow_extrapolation(True)
file_results = df.XDMFFile('Before_ALE_move' + ".xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

sig_ = df.Function(V0sig, name="Sig")
sig_.assign(df.project(sig_after, V0sig))#,mesh=stg_0_msh))
file_results.write(sig_, 0.)

# question - How to get deformed shape without changing the stress field and use it later in program
df.ALE.move(mesh,u)

# stress field after deformation 
u.set_allow_extrapolation(True)
file_results = df.XDMFFile('After_ALE_move' + ".xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

sig_ = df.Function(V0sig, name="Sig")
sig_.assign(df.project(sig_after, V0sig))#,mesh=stg_0_msh))
file_results.write(sig_, 0.)```

If you project your stress field before moving the mesh, the values at each degree of freedom is preserved (but moved along the mesh).
i.e.

# Program that gives you `u`
# ...
sig = project(sig_after, V0sig)
ALE.move(mesh, u)
# Do computations with `sig`
# ...

yes, But If I use that stress field to project on the deformed shape I get projection error as the points are way outside the initial undeformed domain.

I need the stress field on the deformed shape but not affected by the ALE.Move() function

Do you want me to send you the actual image of the problem?

As I said, you need to do the projection before moving the mesh. Then the stress field is preserved at the degrees of freedom (of course an integral over the field will give a different result after moving the mesh).

I understand that. Even If I preserve the stress field as you mentioned how do I project them on the deformed shape? since the deformed shape is outside the domain of the V0sig. it asks me to allow extrapolation which I already did.
BeforeDeformation


Before ALE

After ALE

I need the stress field from Before ALE integrable on the deformed domain since the deformed domain is outside the original domain

The following code shows how grad(u) can be transported with the mesh.

from dolfin import *

mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh, "CG", 2)
u = project(Expression(("pow(x[0],2)", "2*x[1]*x[0]"), degree=2), V)
Vgrad = TensorFunctionSpace(mesh, "DG", 1)
gradu = project(grad(u), Vgrad)
XDMFFile("grad_pre.xdmf").write_checkpoint(gradu, "grad(u)", 0)
ALE.move(mesh, u)
XDMFFile("grad_post.xdmf").write_checkpoint(gradu, "grad(u)", 0)

returning


As you observe here the stress field is the same in each node of the original mesh.

If you want to extrapolate the gradient field from the undeformed mesh to the deformed one you can do:

from dolfin import *

mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh, "CG", 2)
u = project(Expression(("pow(x[0],2)", "2*x[1]*x[0]"), degree=2), V)
Vgrad = TensorFunctionSpace(mesh, "DG", 1)
gradu = project(grad(u), Vgrad)

gradu.set_allow_extrapolation(True)
mesh_new = UnitSquareMesh(10,10)
V_new = VectorFunctionSpace(mesh_new, "CG", 2)
Vgrad_new = TensorFunctionSpace(mesh_new, "DG", 1)
u_new = Function(V_new)
u_new.vector()[:] = u.vector()[:]
print(u_new.vector()[:])
ALE.move(mesh_new, u_new)
gradu_new = interpolate(gradu, Vgrad_new)
XDMFFile("grad_post.xdmf").write_checkpoint(gradu_new, "grad(u)", 0)
2 Likes

It works, But I am having small difference in the stress field (not on the nodes). Is it possible that the error is due to re-projection while exporting to para-view?. When I plot (New_Sig - Old_Sig) for error, I get no error. So I presume that the stress on the nodal points are same.
Befor interpolation


After interpolation

This depends on how you save your files, and what code you are running. If you want any more help, please make a minimal code similar to:

to illustrate your issue

Hi Dokken,

I would like to convert this code to dolfinx. Could you please check and help me to complete the conversion?

from dolfinx import mesh
from dolfinx.fem import Function, FunctionSpace,TensorFunctionSpace
from dolfinx.mesh import CellType
from ufl import TestFunction,grad
from mpi4py import MPI
from dolfinx import fem

mesh= mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral)
V = fem.FunctionSpace(mesh, ("CG", 2))
u= fem.Function(V)
u.interpolate(lambda x: x[0]**2 + x[0]*x[1])
Vgrad = TensorFunctionSpace(mesh, ("CG", 1))
gradu= fem.Function(Vgrad)
gradu.interpolate(grad(u))
ALE.move(mesh, u)

This gives me an error for gradu.interpolate(grad(u)) as :
component, i = component[:-1], component[-1]
IndexError: tuple index out of range

Thanks!

You need to use dolfinx.fem.Expression in the interpolation, see

and to move the mesh use

deformation_array = u.x.array.reshape((-1, mesh.geometry.dim))
    mesh.geometry.x[:, :mesh.geometry.dim] += deformation_array

as done in

1 Like