Different results running Linear elasticity with Cube (gmsh vs fenics)

Question

Hi, I am running a same code for linear elasticity problem, once by creating a 3D cube mesh in Fenics itself and other by importing it from gmsh. However, the output I get for both is a little different. Either this is only because of the mesh or I think my understanding of the behaviour of DomainBoundary() is not correct. Can any one pls. shed some light on this. I can post my geo file if you want. Thanks.

I apologize if I am violating the rule for posting long code: both codes are same except for the mesh

Same code when imported from gmsh:

(Displays the following when run: *** Warning: Found no facets matching domain for boundary condition.)

from dolfin import *
mesh = Mesh("cube.xml")
boundaries = MeshFunction("size_t", mesh, "cube_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "cube_physical_region.xml")
V = VectorFunctionSpace(mesh, "CG", 2)
u, w = TrialFunction(V), TestFunction(V)
b = Constant((1.0, 0.0, 0.0))
E, nu = 10.0, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(w.cell().d)
F = inner(sigma, grad(w))*dx - dot(b, w)*dx
a, L = lhs(F), rhs(F)
c = Expression(("0.0", "0.0", "0.0"))
bc = DirichletBC(V, c, DomainBoundary())
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.parameters["symmetric"] = True
solver.solve()
plot(u, interactive=True)

Same code with Fenics mesh:

from dolfin import *
mesh = UnitCubeMesh(8, 8, 8)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
u, w = TrialFunction(V), TestFunction(V)
b = Constant((1.0, 0.0, 0.0))
E, nu = 10.0, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 -2.0*nu))
sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(w.cell().d)
F = inner(sigma, grad(w))*dx - dot(b, w)*dx
a, L = lhs(F), rhs(F)
c = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, c, DomainBoundary())
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.parameters["symmetric"] = True
solver.solve()
plot(u, interactive=True)

Answer 1

Hello,

At first sight, it seems that you used a finer mesh in the second example than in the first. Try setting the same mesh size for both meshes and then try again. To set the mesh size for the second example, you need to adjust the input arguments for UnitCubeMesh().

regards

Comments

Thanks. I tried doing that and the results come close. They definitely wont ever be absolutely similar in the graphics due to mesh differences but could you still try to explain what does DomainBoundary(), which I used in BC, return?

---

I think that DomainBoundary() is just a shortcut for specifying the domain boundary as an input argument for DirichletBC. According to this, the function DirichletBC expects (among others) an object of type SubDomain as input. And according to this, DomainBoundary returns an object of that type. It seems to be more or less equivalent to simply specifying your own boundary, and then return only on_boundary in the inside() function.