Hello everyone,
Does FEniCS update the coordinates of each node in each iteration ? For example, I am using an Expression as
normal = Expression(('x[0]', 'x[1]', '0.0'),degree=1)
but I want to use coordinates in Euler configuration but not Lagrangian. Does it update the coordinates in each iteration ? If not, how it can be possible ? Thanks
I am not sure what exactly you mean here. If I understand you question correctly, you ask whether fenics updates the corrdinates (of the dofs) when the mesh is moved.
The answer to that question is yes. If you move your mesh, then (almost) everything is updated automatically. In particular, if you have the above mentioned expression and want to use it, e.g., in a variational form then fenics evaluates this expression on the new coordinates when you perform a “solve” or “assemble”.
Hello plugged,
Thanks for your reply. Yes this is what I wanted to ask. However from my experience I created an Expression as:
normal = Expression(('x[0]', 'x[1]', '0.0'),degree=1)
and I applied torsion to the geometry. Normally x[0] and x[1] coordinates should change if it updates the coordinates in Expression. However, I checked the outputs of this “normal” vector by using projection on a vectorfunction space, there is nothing change during incremental solution.
Can you send a MWE? Also, I would not project the expression, but interpolate it (this gives the correct nodal values).
Hello plugged,
I converted hyperelasticity demo to my problem as following:
from dolfin import * import numpy as np # Optimization options for the form compiler parameters["form_compiler"]["cpp_optimize"] = True ffc_options = {"optimize": True, \ "eliminate_zeros": True, \ "precompute_basis_const": True, \ "precompute_ip_const": True} # Create mesh and define function space mesh = UnitCubeMesh(24, 16, 16) V = VectorFunctionSpace(mesh, "Lagrange", 1) # Mark boundary subdomians left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0) right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0) # Define Dirichlet boundary (x = 0 or x = 1) c = Expression(("0.0", "0.0", "0.0"),degree=1) r = Expression(("scale*0.0", "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])", "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"), scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3,degree=1) bcl = DirichletBC(V, c, left) bcr = DirichletBC(V, r, right) bcs = [bcl, bcr] load_start =0 load_stop =pi/3 load_step =15 def normal(a): normal = Expression(('x[0]', 'x[1]', 'x[2]'),degree=1) unit = a*normal / sqrt(dot(normal,normal)) asd=interpolate(normal,V) aaa11=asd.vector().array() print("lalalal1",aaa11[102:105]) print("lalalal2",aaa11[198:201]) print("lalalal3",aaa11[798:801]) #theta=project(theta1,P3) return normal inc_load = np.linspace(load_start,load_stop,load_step+1) # Define functions du = TrialFunction(V) # Incremental displacement v = TestFunction(V) # Test function u = Function(V) # Displacement from previous iteration B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary # Kinematics d = u.geometric_dimension() I = Identity(d) # Identity tensor F = I + grad(u) # Deformation gradient C = F.T*F # Right Cauchy-Green tensor # Invariants of deformation tensors Ic = tr(C) J = det(F) # Elasticity parameters E, nu = 10.0, 0.3 mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu))) # Stored strain energy density (compressible neo-Hookean model) psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 # Total potential energy Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds # Compute first variation of Pi (directional derivative about u in the direction of v) F = derivative(Pi, u, v) # Compute Jacobian of F J = derivative(F, u, du) for t in inc_load: r.theta=t # Solve variational problem solve(F == 0, u, bcs, J=J) print("t is",t) normal(1) # Save solution in VTK format file = File("displacement.pvd"); file << u; # Plot and hold solution plot(u, mode = "displacement", interactive = True)
If you run this code, it will print 3 of random coordinate vector which called as lalalal here. It does not change in each iteration. They are constant.
You have not moved the mesh here. The displacement mode for the plotting shows a deformed mesh, which is, however, not altered. You can check this by comparing mesh.coordinates() before and after the solve.
You can also move the mesh with, e.g.,
ALE.move(mesh, u)
and with this you will see that your coordinates actually do change. Note that it might be possible to update the bounding box tree of the mesh via
mesh.bounding_box_tree().build(mesh)
Hope this helps and fixes your problem.
Hello plugged,
ALE.move(mesh,u)
solved the problem perfectly for hyperelasticity demo. However when I try it on my exact Mixed Formulation it gives an error as following:
RuntimeError Traceback (most recent call last)
in ()
347 solver.solve()
348 (p,u,teta) = upteta.split()
→ 349 ALE.move(mesh, u)
350/usr/lib/python3/dist-packages/dolfin/cpp/function.py in move(*args)
2112
2113 “”"
→ 2114 return _function.ALE_move(*args)
2115
2116 move = staticmethod(move)RuntimeError:
Beside of that, actually I do not know anything about bounding box of the mesh, can you please explain that code briefly what is used for ?
First, try the following line
(p, u, theta) = upteta.split(True)
which creates a deepcopy of the function (assuming that u is the component of the mixed function space that is valid for a deformation, e.g., a CG1 vector function space of the physical dimension).
Second, the bounding box tree is important, e.g., when you want to evaluate a function via
value = u(0, 0, 0)
and it is not updated automatically when using ALE.move.
Thanks for your reply but I got still the same error:
RuntimeError Traceback (most recent call last)
in ()
347 solver.solve()
348 (p, u, theta) = upteta.split(True)
→ 349 ALE.move(mesh, u)
350
351/usr/lib/python3/dist-packages/dolfin/cpp/function.py in move(*args)
2112
2113 “”"
→ 2114 return _function.ALE_move(*args)
2115
2116 move = staticmethod(move)
I also tried
mesh1 = Mesh(mesh) X = mesh1.coordinates() X += np.vstack(map(u, X))
Actually it creates deformed mesh but It gives other error when I try to project the stress quantity for visualization.
I figured out that
ALE.move(mesh,u)
does not for CG=2 for hyperelasticity demo, but it works with CG=1. Is there any way to make it work for CG=2 ?
I am not sure. I think that the problem is the following: If you use CG1 Elements, then you map a linear mesh (consisting of, e.g., triangles) to a nonlinear one, where the points are not necessarily joint by straight lines. You could, of course, always interpolate your CG2 solution into a CG1 FunctionSpace, such that you use that nodal values of your CG2 solution, but keep the linear mesh.
Yes probably… It creates problems with CG2. Is there any alternative method rather than ALE.move ?
You can update the mesh.coordinates() manually. I’ve also tried this, and for a vector CG1 function it yields exactly the same results. So the interpolating would yield the same thing there. But otherwise you have a bit more freedom with that approach. I dont know about other things.
Ok I tried to mapped displacement on mesh and it seems like it works. However I am not able to create a copy or deepcopy of mesh. That mapping directly affects other variables (other expressions, boundary conditions etc.). Therefore I tried to create a copy and deepcopy of the mesh but it did not work. Is there any solution to use deformed mesh coordinates for a specific expression and use undeformed mesh coordinates for the rest of the expressions ? Thanks.
I think you should be able to do so via the domain keyword of expression. For example
mesh_copy = Mesh(mesh)
# transform the copied mesh
expr_1 = Expression('x[0]', degree=1, domain=mesh)
expr_2 = Expression('x[0]', degree=1, domain=mesh_copy)
I think this will do the trick.
I tried that method as following but it seems like does not work. I just try to keep update the normal vector in deformed mesh for radial stress calculation and keep the rest of the calculation in undeformed mesh. The problem might be about, I am solving displacement (u) in “mesh” but trying to map it by “copy_mesh”…But I am not sure about it.
> from dolfin import *
> import numpy as np
> import copy
> # Optimization options for the form compiler
> parameters["form_compiler"]["cpp_optimize"] = True
> ffc_options = {"optimize": True, \
> "eliminate_zeros": True, \
> "precompute_basis_const": True, \
> "precompute_ip_const": True}
>
> # Create mesh and define function space
> mesh = UnitCubeMesh(6, 4, 4)
> mesh_copy = Mesh(mesh)
> V = VectorFunctionSpace(mesh, "Lagrange", 2)
>
> # Mark boundary subdomians
> left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
> right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> c = Expression(("0.0", "0.0", "0.0"),degree=1)
> r = Expression(("scale*0.0",
> "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
> "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
> scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3,degree=1,domain=mesh)
>
> bcl = DirichletBC(V, c, left)
> bcr = DirichletBC(V, r, right)
> bcs = [bcl, bcr]
>
>
> load_start =0
> load_stop =pi/3
> load_step =15
>
> def normal_deformed(mesh_copy,u,V):
> mesh_copy.bounding_box_tree().build(mesh_copy)
> X = mesh_copy.coordinates()
> X += np.vstack(map(u, X))
> x = SpatialCoordinate(mesh_copy)
> normal = Expression(('x[0]', 'x[1]', 'x[2]'),degree=1,domain=mesh_copy)
> unit = normal / sqrt(dot(normal,normal))
> asd=project(unit,V)
> aaa11=asd.vector().array()
> print("Deformed1",aaa11[102:105])
> print("Deformed2",aaa11[198:201])
> print("Deformed3",aaa11[798:801])
> #theta=project(theta1,P3)
> #mesh.bounding_box_tree().build(mesh1)
> return dot(normal,normal)
>
> def normal_undeformed(mesh_copy,u,V):
> x = SpatialCoordinate(mesh)
> normal = Expression(('x[0]', 'x[1]', 'x[2]'),degree=1,domain=mesh)
> unit = normal / sqrt(dot(normal,normal))
> asd=project(unit,V)
> aaa11=asd.vector().array()
> print("Undeformed1",aaa11[102:105])
> print("Undeformed2",aaa11[198:201])
> print("Undeformed3",aaa11[798:801])
> #theta=project(theta1,P3)
> #mesh.bounding_box_tree().build(mesh1)
> return dot(normal,normal)
>
>
>
>
> inc_load = np.linspace(load_start,load_stop,load_step+1)
> # Define functions
> du = TrialFunction(V) # Incremental displacement
> v = TestFunction(V) # Test function
> u = Function(V) # Displacement from previous iteration
> B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume
> T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary
>
> # Kinematics
> d = u.geometric_dimension()
> I = Identity(d) # Identity tensor
> F = I + grad(u) # Deformation gradient
> C = F.T*F # Right Cauchy-Green tensor
>
> # Invariants of deformation tensors
> Ic = tr(C)
> J = det(F)
>
> # Elasticity parameters
> E, nu = 10.0, 0.3
> mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
>
> # Stored strain energy density (compressible neo-Hookean model)
> psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
>
> # Total potential energy
> Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
>
> # Compute first variation of Pi (directional derivative about u in the direction of v)
> F = derivative(Pi, u, v)
>
> # Compute Jacobian of F
> J = derivative(F, u, du)
>
>
>
> P3_Def_v= VectorFunctionSpace(mesh_copy, "CG", 1)
> P3_UnDef_v= VectorFunctionSpace(mesh, "CG", 1)
> P3_Def= FunctionSpace(mesh_copy, "CG", 1)
> P3_UnDef= FunctionSpace(mesh, "CG", 1)
> normal_def=Function (P3_Def, name="Deformed Normal")
> normal_undef=Function (P3_UnDef, name="Undeformed Normal")
>
> for t in inc_load:
> r.theta=t
> # Solve variational problem
> solve(F == 0, u, bcs, J=J)
> print("t is",t)
> normal_def.assign(project(normal_deformed(mesh_copy,u,P3_Def_v),P3_Def))
> normal_undef.assign(project(normal_undeformed(mesh,u,P3_UnDef_v),P3_UnDef))
>
>
> # Save solution in VTK format
> file = File("displacement.pvd");
> file << u;
>
> # Plot and hold solution
> plot(u, mode = "displacement", interactive = True)Is the following example helpful?
from dolfin import *
mesh = UnitCubeMesh(10,10,10)
p = 2
V = VectorFunctionSpace(mesh,"Lagrange",p)
# You can pass a Function into the Expression, and it will be evaluated at
# the coordinate x, i.e., inside the string of C++ code, u is a double [3],
# and you can index components directly:
u = Function(V)
d = mesh.geometry().dim()
x = Expression(tuple(("x["+str(i)+"]+u["+str(i)+"]" for i in range(0,d))),
u=u,degree=p)
# Verify:
u.interpolate(Expression(("x[0]*x[0]","0","0"),degree=2))
print(assemble(x[0]*dx(domain=mesh)))
# Simpler alternative, if you just need the spatial coordinates in UFL, and
# not specifically as an Expression:
X = SpatialCoordinate(mesh)
x = X + u
# Example from original post:
normal = as_vector([x[0],x[1],0.0])
Hello kamensky,
Thank you so much for your advice. Actually it seems like, it works without an error. But I am not sure about the correctness of the results. I applied your second simpler alternative to the my example code. It is a unit cube and I rotate it in the right face where x=1.0. The rotation angle is pi/2. So I expect that the tip coordinate of node which is (1.0,0.75,0) should become (0.75,1.0,0) with 90 degree rotation in deformed configuration. But it only becomes (0.986, 0.754,0.0). May it be about the projection that I perform on “mesh” which represent the undeformed configuration ?
> from dolfin import *
> import numpy as np
> import copy
> # Optimization options for the form compiler
> parameters["form_compiler"]["cpp_optimize"] = True
> ffc_options = {"optimize": True, \
> "eliminate_zeros": True, \
> "precompute_basis_const": True, \
> "precompute_ip_const": True}
>
> # Create mesh and define function space
> mesh = UnitCubeMesh(6, 4, 4)
> V = VectorFunctionSpace(mesh, "Lagrange", 2)
>
> # Mark boundary subdomians
> left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
> right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> c = Expression(("0.0", "0.0", "0.0"),degree=1)
> r = Expression(("scale*0.0",
> "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
> "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
> scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3,degree=1,domain=mesh)
>
> bcl = DirichletBC(V, c, left)
> bcr = DirichletBC(V, r, right)
> bcs = [bcl, bcr]
>
>
> load_start =0
> load_stop =pi/2
> load_step =15
>
> def normal_deformed(mesh,u,V):
> X = SpatialCoordinate(mesh)
> x = X + u
> normal = as_vector([x[0],x[1],0.0])
> unit = normal / sqrt(dot(normal,normal))
> asd=project(normal,V)
> aaa11=asd.vector().array()
> print("Deformed1",aaa11[6:9])
> print("Deformed2",aaa11[198:201])
> print("Deformed3",aaa11[498:501])
> #theta=project(theta1,P3)
> #mesh.bounding_box_tree().build(mesh1)
> return dot(normal,normal)
>
> def normal_undeformed(mesh,u,V):
> x = SpatialCoordinate(mesh)
> normal = as_vector([x[0],x[1],0.0])
> unit = normal / sqrt(dot(normal,normal))
> asd=project(normal,V)
> aaa11=asd.vector().array()
> print("Undeformed1",aaa11[6:9])
> print("Undeformed2",aaa11[198:201])
> print("Undeformed3",aaa11[498:501])
> #theta=project(theta1,P3)
> #mesh.bounding_box_tree().build(mesh1)
> return dot(normal,normal)
>
>
>
>
> inc_load = np.linspace(load_start,load_stop,load_step+1)
> # Define functions
> du = TrialFunction(V) # Incremental displacement
> v = TestFunction(V) # Test function
> u = Function(V) # Displacement from previous iteration
> B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume
> T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary
>
> # Kinematics
> d = u.geometric_dimension()
> I = Identity(d) # Identity tensor
> F = I + grad(u) # Deformation gradient
> C = F.T*F # Right Cauchy-Green tensor
>
> # Invariants of deformation tensors
> Ic = tr(C)
> J = det(F)
>
> # Elasticity parameters
> E, nu = 10.0, 0.3
> mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
>
> # Stored strain energy density (compressible neo-Hookean model)
> psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
>
> # Total potential energy
> Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
>
> # Compute first variation of Pi (directional derivative about u in the direction of v)
> F = derivative(Pi, u, v)
>
> # Compute Jacobian of F
> J = derivative(F, u, du)
>
>
>
> P3_v= VectorFunctionSpace(mesh, "CG", 1)
> P3= FunctionSpace(mesh, "CG", 1)
> normal_def=Function (P3, name="Deformed Normal")
> normal_undef=Function (P3, name="Undeformed Normal")
>
> for t in inc_load:
> r.theta=t
> # Solve variational problem
> solve(F == 0, u, bcs, J=J)
> print("t is",t)
> normal_def.assign(project(normal_deformed(mesh,u,P3_v),P3))
> normal_undef.assign(project(normal_undeformed(mesh,u,P3_v),P3))
>
>
> # Save solution in VTK format
> file = File("displacement.pvd");
> file << u;
>
> # Plot and hold solution
> plot(u, mode = "displacement", interactive = True)