Hello everyone,

I am running a linear elastic problem using Fenics. The mesh files are generated in Gmsh. I found that it works for smaller domain, while the results are all zeros for larger domain. It seems the displacement boundary conditions were not actually applied, but there are no errors during the running. The job finished normally, but the displacements in the output file are all zeros. Does anyone have any idea about this issue? Any comments and suggestions are appreciated.

```
from __future__ import print_function
from dolfin import *
from ufl import nabla_div
from mpi4py import MPI
import numpy as np
import meshio
import gmsh
import sys
if comm.rank == 0:
gmsh.initialize()
gmsh.merge('sample1_600.stl')
s = gmsh.model.getEntities(2)
l = gmsh.model.geo.addSurfaceLoop([e[1] for e in s])
V0 = gmsh.model.geo.addVolume([l])
gmsh.model.addPhysicalGroup(3, [V0], 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3)
#obtain the model dimension
dimensionBox = gmsh.model.getBoundingBox(-1, -1)
gmsh.write("S1_600.msh")
gmsh.finalize()
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
mesh_from_file = meshio.read("S1_600.msh")
tetra_mesh = create_mesh(mesh_from_file, "tetra", prune_z=False)
meshio.write("tetra_mesh.xdmf", tetra_mesh)
else:
dimensionBox = None
comm.Barrier()
dimensionBox = comm.bcast(dimensionBox, root=0)
xmin = dimensionBox[0]
ymin = dimensionBox[1]
zmin = dimensionBox[2]
xmax = dimensionBox[3]
ymax = dimensionBox[4]
zmax = dimensionBox[5]
mesh = Mesh()
with XDMFFile("tetra_mesh.xdmf") as infile:
infile.read(mesh)
V = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
tol = 1E-6
def boundary_Left(x, on_boundary):
return on_boundary and near(x[0], xmin, tol)
def boundary_Right(x, on_boundary):
return on_boundary and near(x[0], xmax, tol)
def boundary_Front(x, on_boundary):
return on_boundary and near(x[1], ymin, tol)
def boundary_Back(x, on_boundary):
return on_boundary and near(x[1], ymax, tol)
def boundary_Top(x, on_boundary):
return on_boundary and near(x[2], zmax, tol)
def boundary_Bottom(x, on_boundary):
return on_boundary and near(x[2], zmin, tol)
applied_strain = 1.0e-4
applied_disp = applied_strain * (xmax - xmin)
bc_l = DirichletBC(V.sub(0), Constant(0), boundary_Left)
bc_r = DirichletBC(V.sub(0), Constant(applied_disp), boundary_Right)
bc_f = DirichletBC(V.sub(1), Constant(0), boundary_Front)
bc_ba = DirichletBC(V.sub(1), Constant(0), boundary_Back)
bc_bo = DirichletBC(V.sub(2), Constant(0), boundary_Bottom)
bc_t = DirichletBC(V.sub(2), Constant(0), boundary_Top)
bcs = [bc_l, bc_r, bc_f, bc_ba, bc_bo, bc_t]
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx
u = Function(V)
solve(a == L, u, bcs, solver_parameters = {'linear_solver':'mumps'})
XDMFFile('disp_test1.xdmf').write(u)
```