Hi everyone in the community,
I’m quite a newbie in Fenics(2019.1.0)/Fenicsx(0.9), actually also a newbie in FEM simulation, so please forgive me if my questions are not so smart.
Recently I’m working on 2D contact problems using the third medium method. Here is a simple box example to model:
Initial state of the box:
Final state of the box:
A second-order mesh is supposed to be used, so I created the mesh using GMSH, here is the box.geo:
p0x = 0.0;
p0y = 0.0;
p1x = 2.0;
p1y = 0.5;
Nx = 40;
Ny = 10;
lc = 1.0;
Point(1) = {p0x, p0y, 0, lc};
Point(2) = {p0x, p1y, 0, lc};
Point(3) = {p1x, p1y, 0, lc};
Point(4) = {p1x, p0y, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Physical Surface("My surface") = {1};
Transfinite Curve {1, -3} = Ny+1;
Transfinite Curve {2, -4} = Nx+1;
Transfinite Surface{1};
//Recombine Surface{1};
Mesh.ElementOrder = 2; // Second-order elements
Mesh.Algorithm = 8; // Delaunay mesh algorithm
Mesh 2;
Then I import the mesh into Fenics:
mesh_name = "box"
msh = meshio.read(f"{mesh_name}.msh")
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
"name_to_read":[cell_data]})
return out_mesh
triangle_mesh = create_mesh(msh, "triangle6", prune_z=True)
meshio.write(f"{mesh_name}_mesh.xdmf", triangle_mesh)
mesh = Mesh()
with XDMFFile(f"{mesh_name}_mesh.xdmf") as infile:
infile.read(mesh)
But when running the script, the problem is not solved at all and the box looks very strange in the xdmf file:
n=: 1
Step: 1 ut.n = 1.0 Nsteps = 10.0
Prescribed top displacement = -0.1
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.776e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
n=: 2
Step: 2 ut.n = 2.0 Nsteps = 10.0
Prescribed top displacement = -0.2
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.776e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
…
If I modify the mesh order into 1, i.e. in box.geo, set Mesh.ElementOrder = 1, the problem can be solved successfully, and the results are shown in the first two images.
So I think I am not importing the second-order mesh correctly — there must be something missing. I would appreciate any help!
FYI, I have also succeeded in loading the box.msh(with second-order mesh) in a Fenicsx script using gmshio.read_from_msh, and the problem can be solved without any issues. I want to stick with Fenics because there are some arc_length solvers available here, while Fenicsx doesn’t have them yet.