Hi
I am solving a hyperelasticity problem. The domain is a simple unit cube including 6 elements (1 element on each edge) which is fixed from the left face and a pressure is applied on the right face. I created the mesh using 2 different approaches. In the first approach I used FEniCS built-in function. In the second approach I made exactly the same mesh in GMSH and imported into the FEniCS. Here is the code when mesh was built in FEniCS:
from dolfin import *
class Left(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary
mesh = UnitCubeMesh.create(1,1,1,CellType.Type.tetrahedron)
#mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)
left = Left()
right = Right()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())
n = FacetNormal(mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
V = VectorFunctionSpace(mesh, 'CG', degree=2)
du = TrialFunction(V) # Trial function
v = TestFunction(V) # Test function
u = Function(V)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F
Ic = tr(C)
J = det(F)
Ic_bar = pow(J,-2./3.)*Ic
K = 1E6
Pressure = Constant((1.))
c_1 = 10.
psi = c_1 * (Ic_bar - 3.) * dx + 0.5 * K * pow(ln(J),2) * dx - inner(Pressure * n, u) * ds(2)
F1 = derivative(psi, u, v)
# Compute Jacobian of F
Jac = derivative(F1, u, du)
bc_1 = DirichletBC(V, Constant((0., 0., 0.)), boundaries, 1)
bcs = [bc_1]
problem = NonlinearVariationalProblem(F1, u, bcs, Jac)
solver = NonlinearVariationalSolver(problem)
solver.solve()
File("fenics.pvd") << u
Here is the geo file of the GMSH:
// Gmsh
Point(1) = {0, 0, 0,1};
Point(2) = {0, 1, 0,1};
Point(3) = {0, 0, 1,1};
Point(4) = {0, 1, 1,1};
Point(5) = {1, 0, 0,1};
Point(6) = {1, 1, 0,1};
Point(7) = {1, 0, 1,1};
Point(8) = {1, 1, 1,1};
Line(1) = {1, 3};
Line(2) = {7, 3};
Line(3) = {7, 5};
Line(4) = {5, 1};
Line(5) = {1, 2};
Line(6) = {2, 4};
Line(7) = {4, 3};
Line(8) = {2, 6};
Line(9) = {6, 8};
Line(10) = {8, 4};
Line(11) = {6, 5};
Line(12) = {8, 7};
Line Loop(13) = {8, 9, 10, -6};
Plane Surface(14) = {13};
Line Loop(15) = {8, 11, 4, 5};
Plane Surface(16) = {15};
Line Loop(17) = {6, 7, -1, 5};
Plane Surface(18) = {17};
Line Loop(19) = {4, 1, -2, 3};
Plane Surface(20) = {19};
Line Loop(21) = {11, -3, -12, -9};
Plane Surface(22) = {21};
Line Loop(23) = {7, -2, -12, 10};
Plane Surface(24) = {23};
Surface Loop(25) = {14, 16, 22, 20, 18, 24};
Volume(26) = {25};
Transfinite Surface {14,16,18,20,22,24};
Transfinite Line{1,2,3,4,5,6,7,8,9,10,11,12}=2;
Physical Volume(29) = {26};
// LEFT BOUNDARY
Physical Surface(1) = {18};
// RIGHT BOUNDARY
Physical Surface(2) = {22};
and this is the FEniCS code when the mesh is imported from GMSH:
from dolfin import *
mesh = Mesh("MESH.xml")
boundaries = MeshFunction("size_t", mesh, "MESH_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "MESH_physical_region.xml")
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, 'CG', degree=2)
du = TrialFunction(V) # Trial function
v = TestFunction(V) # Test function
u = Function(V)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F
Ic = tr(C)
J = det(F)
Ic_bar = pow(J,-2./3.)*Ic
K = 1E6
# Total potential energy
Pressure = Constant((1))
c_1 = 10.
psi = c_1 * (Ic_bar - 3.) * dx + 0.5 * K * pow(ln(J),2) * dx - inner(Pressure * n, u) * ds(2)
F1 = derivative(psi, u, v)
# Compute Jacobian of F
Jac = derivative(F1, u, du)
bc_1 = DirichletBC(V, Constant((0., 0., 0.)), boundaries, 1)
bcs = [bc_1]
problem = NonlinearVariationalProblem(F1, u, bcs, Jac)
solver = NonlinearVariationalSolver(problem)
solver.solve()
File("gmsh.pvd") << u
If we look into the results, we can see that the results are significantly different. For example this is displacement contour in Z direction:
As it was shown, the results are completely different while I expect the same results as everything are the same in 2 approaches. The only difference is the way the mesh has been created. Even if we increase the polynomial order from 1 to 2, the results are still different.
It is confusing why the results for such a simple mesh are not matched. Does anybody have any idea what is wrong in here?
