Different results for the same mesh created differently (Hyperelasticity)

mesh
2

I claim that the meshes are not the same, see attached picture:

cube_comparison
cube_comparison1523×1238 142 KB
}
I would also recommend you not to use xml to load meshes, as it has been deprecated for a while. The suggested way of converting a mesh from gmsh is by using meshio.

Also, to make sure that your visualization is P2 data is saved, use xdmf.write_checkpoint.
Here is a script doing all of the suggested actions:

from dolfin import *
import meshio
import numpy as np
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


def problem(gmsh=True):
    if gmsh:
        msh = meshio.read("mesh.msh")
        points = msh.points
        cells = np.vstack([cells.data for cells in msh.cells
                           if cells.type == "tetra"])
        cell_data = np.vstack([msh.cell_data_dict["gmsh:physical"][key]
                            for key in
                            msh.cell_data_dict["gmsh:physical"].keys()
                            if key =="tetra"])
        mesh = meshio.Mesh(points=points,
                           cells=[("tetra", cells)],
                           cell_data={"name_to_read": cell_data})
        meshio.write("mesh.xdmf", mesh)


        # Write line mesh
        facet_cells = np.vstack([cells.data for cells in msh.cells
                                 if cells.type == "triangle"])
        facet_data = np.vstack([msh.cell_data_dict["gmsh:physical"][key]
                             for key in
                             msh.cell_data_dict["gmsh:physical"].keys()
                                 if key == "triangle"])

        facet_mesh = meshio.Mesh(points=points,
                             cells=[("triangle", facet_cells)],
                             cell_data={"name_to_read": facet_data})
        meshio.write("boundaries.xdmf", facet_mesh)

        mesh = Mesh()
        with XDMFFile("mesh.xdmf") as infile:
            infile.read(mesh)
            mvcv = MeshValueCollection("size_t", mesh, mesh.topology().dim())
            infile.read(mvcv, "name_to_read")
        subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvcv)
        mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
        with XDMFFile("boundaries.xdmf") as infile:
            infile.read(mvc, "name_to_read")
        boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)

    else:
        mesh = UnitCubeMesh.create(1,1,1,CellType.Type.tetrahedron)
        left = Left()
        right = Right()

        boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
        subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

        boundaries.set_all(0)
        left.mark(boundaries, 1)
        right.mark(boundaries, 2)

    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()
    if gmsh:
         XDMFFile("gmsh.xdmf").write_checkpoint(u,"u",0, append=False)
    else:
         XDMFFile("built_in.xdmf").write_checkpoint(u,"u",0,append=False)

problem(True)
problem(False)
3 Likes