Hello everyone !
I’m currently using Dolfinx in the context of linear elasticity and must evaluate quite precisely the average strain in my simulation. My problem is that even when the problem is symmetrical with respect to the horizontal axis (meaning geometry, BC and mesh are symmetrical) the result of my simulation does not respect exactly the symmetry. You can see below a picture of the shear strain on the simple case of a rectangle clamped on the right end and with an imposed horizontal force on the left end.
I use a structured mesh with quadrilaterals and therefore I don’t see any reason for this break of symmetry. As a result I have a non zero shear average strain in my sample. The magnitude is low (1e-05 times lower than the compression) but should be close to machine precision if the symmetry was not broken. I would at least expect this average shear strain to converge to zero when I decrease the size of the mesh but it’s not the case (it seems even to be the contrary up to mesh size =1/150).
You can find attached the corresponding code for the MWE of the compressed rectangle. Any idea where does this break of symmetry comes from ?
import dolfinx # FEM in python
import matplotlib.pyplot as plt
import ufl # variational formulations
import numpy as np
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import gmsh # Mesh generation
#Geometry parameters
cell_size = 1 # height
x_length=3*cell_size #length
# Material
E = 1.0 # Young's modulus
nu = 0.3 # Poisson's ratio
# geometry and mesh
mesh_size = cell_size/10
domain=dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,[np.array([0, 0]), np.array([x_length,cell_size])],[int(x_length/mesh_size),int(cell_size/mesh_size)],dolfinx.mesh.CellType.quadrilateral)
#facets_tags
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], x_length)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], cell_size))]
indexLeft=1
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = dolfinx.mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
#define measure and function space
dx = ufl.Measure("dx", domain=domain)
one = dolfinx.fem.Constant(domain,ScalarType(1.))
ds= ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
V = dolfinx.fem.VectorFunctionSpace(domain,("Lagrange", 1),dim=2)
#dirichlet BC
facets_dim=1;
def right_boundary(x):
return np.isclose(x[0], x_length)
R_f=dolfinx.mesh.locate_entities_boundary(domain, facets_dim, right_boundary)
R_dofs = dolfinx.fem.locate_dofs_topological(V, facets_dim, R_f)
uright = np.array([0, 0],dtype=ScalarType)
bc = dolfinx.fem.dirichletbc(uright, R_dofs,V)
#constitutive law
def epsilon(u):
return ufl.sym(ufl.grad(u))
I2 = ufl.Identity(2)
def sigma(eps, E, nu):
mu = E/2/(1+nu)
lamb = 2*mu*nu/(1-nu)
return lamb*ufl.tr(eps)*I2 + 2*mu*eps
#define elastic problem
#test and trial functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
#elastic bilinear form
bilinear_form = ufl.inner(sigma(epsilon(u),E,nu),epsilon(v))*dx
#surface force (applied on boundary x=0)
surface_force = dolfinx.fem.Constant(domain, ScalarType((0.1,0)))#no lineic force distribution
linear_form = ( ufl.dot(surface_force,v) ) *ds(indexLeft)
problem=dolfinx.fem.petsc.LinearProblem(bilinear_form, linear_form, bcs=[bc])
#solve
u_solution = problem.solve()
#compute strains
V_eps = dolfinx.fem.FunctionSpace(domain,("DG", 0))
eps=epsilon(u_solution)
#compression along x
eps_xx_expr = dolfinx.fem.Expression(eps[0,0], V_eps.element.interpolation_points())
eps_xx = dolfinx.fem.Function(V_eps)
eps_xx.interpolate(eps_xx_expr)
#shear
eps_xy_expr = dolfinx.fem.Expression(eps[0,1], V_eps.element.interpolation_points())
eps_xy = dolfinx.fem.Function(V_eps)
eps_xy.interpolate(eps_xy_expr)
#average shear
comm=eps_xy.function_space.mesh.comm
shear_strain = dolfinx.fem.form(eps_xy *dx)
Shear_Resultant = dolfinx.fem.assemble_scalar(shear_strain)
print("Average Exy= ", Shear_Resultant)
#average compression
comm=eps_xx.function_space.mesh.comm
compressive_strain = dolfinx.fem.form(eps_xx*dx)
Exx_Resultant = dolfinx.fem.assemble_scalar(compressive_strain)
print("Average Exx = ", Exx_Resultant)
#write mesh+solution in XDMF file
# To have a name in Paraview
eps_xx.name = "exx"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output/strain_compressive_x_.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(eps_xx)
#write mesh+solution in XDMF file
# To have a name in Paraview
eps_xy.name = "exy"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output/strain_shear_.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(eps_xy)
