Hi Everyone,
I’ve been using legacy dolfin 2019.1.0 for some time. I want to move some of my legacy code to (and build new code with) fenicsX. However, I’m having some reproducibility issues which manifest when using the built-in mesh generation tools for different elements. To illustrate this, I am considering the 2d Eshelby problem, which involves solving for the displacements in an isotropic and homogeneous elastic medium near to an artificially strained region.
To be explicit, the static elasticity problem is, in this case:
0 = \underline{\nabla}\cdot\underline{\underline{\sigma}}
With the stress \sigma defined as:
\underline{\underline{\sigma}} = 2\mu(\varepsilon - \varepsilon_{pl}) + \lambda \text{Tr}(\varepsilon-\varepsilon_{pl})
Where \varepsilon_{pl} represents the eigenstrain at a particular point in space. For purposes of testing, \varepsilon_{pl} is everywhere zero, except in a square of side length 1, at the center of the simulation box, with \varepsilon_{pl} = \begin{pmatrix} 0 & 1 \\ 1& 0\end{pmatrix}. Meanwhile, the strain is defined in the usual way: \varepsilon = \frac{1}{2}(\underline{\nabla} \underline{u} + (\underline{\nabla}\underline{u})^T)
By applying a simple shear to a small region of the material, we expect stress response outside of the strained region with a quadrupolar character. Within the artificially strained region, I expect a stress response in the opposite direction of the strain (i.e. if a strain with a positive _xy component is applied, then there should be a negative stress_xy component in the region of eigenstrain). This works well in legacy dolfin, regardless of the choice of triangular elements generated by the rectangular mesh (‘crossed’, ‘left/right’ or just ‘left’).
However, in dolfinX, I get very different results in the region of eigenstrain depending on the choice of mesh generation parameters. This doesn’t seem to happen in dolfin (where at least both ‘crossed’ and ‘left/right’ meshes give similar values, though dolfin didn’t support quadrilateral elements).
Consider the following MWE. To avoid redundancy, since I want to compare different meshes, I’ve split the solver part of the code into the second code snippet.
import numpy as np
import epm_phonic
import matplotlib.pyplot as plt
from dolfinx.mesh import (CellType, DiagonalType)
def compute_eshelby_stress(L,mesh_options):
epm = epm_phonic.epm_phonic(L,mesh_options = mesh_options)
#applying an eigenstrain into the center of the domain:
epm.add_inclusion_xy_instant(1/float(epm.mu.value),L//2,L//2)
#solving for the dipslacements
epm.solve_static()
#projecting out the 'xy' component of the stresses.
epm.solve_sig_xy()
#let's just look at the degrees of freedom for the central stress.
dof_coords = epm.U.tabulate_dof_coordinates()
delta = 0.5
for i,dof_coord in enumerate(dof_coords):
if(np.abs(dof_coord[0] - (L//2 + .5)) <= delta and np.abs(dof_coord[1] - (L//2+0.5)) <= delta):
print('stress_xy in inclusion: ',dof_coord,epm.sig_xy_U_vector.vector[i])
fig,ax = plt.subplots(1,1,dpi=100,figsize=[2,2])
sc = ax.scatter(dof_coords[:,0],dof_coords[:,1],c = (1e-99+(epm.sig_xy_U_vector.vector[:])))
plt.colorbar(sc)
plt.show()
mesh_options = {'cell_type':CellType.quadrilateral}
compute_eshelby_stress(11,mesh_options)
mesh_options = {'cell_type':CellType.triangle,'diagonal':DiagonalType.crossed}
compute_eshelby_stress(11,mesh_options)
mesh_options = {'cell_type':CellType.triangle,'diagonal':DiagonalType.left}
compute_eshelby_stress(11,mesh_options)
Here’s the stress distribution for quadrilateral elements:
And for the ‘left’ triangular elements:
Both show the expected quadrupolar kernel, but there’s a striking difference within the strained region.
Within the strained region, I find stress_xy components ranging from -0.04 (with crossed elements) to -0.5 (with the ‘left’ triangular elements).
Here is code for epm_phonic.py defining the forms and the solvers for u and for stress (which is imported as import epm_phonic)
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import dolfinx
from dolfinx import la
from dolfinx.fem import (Expression, Function, FunctionSpace,
VectorFunctionSpace, TensorFunctionSpace, dirichletbc, form,
locate_dofs_topological, Constant)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary,DiagonalType)
from ufl import dx, grad, inner
import numpy as np
#building the boundary conditions class
class simple_shear_bc():
def __init__(self,epsxy):
self.epsxy = epsxy
def __call__(self,x):
values = np.zeros(shape=(2,x.shape[1]))
#we have a simple shear condition, so:
values[0,:] = self.epsxy*x[1,:]
# values[1,:] = 0.0 #we omit this, because we are initializing with zeros anyways.
return values
# Strain function
def epsilon(u):
return 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
# Stress function
def sigma(u,eigenstrain,mu,lmbda):
return lmbda*(ufl.div(u)-ufl.tr(eigenstrain))*ufl.Identity(2) + 2*mu*(epsilon(u) - eigenstrain)
#pretty slow probably, but whatever.
class single_inclusion():
def __init__(self,epsxy,x0,y0):
self.epsxy = epsxy
self.x0=x0
self.y0=y0
def __call__(self,x):
values = np.zeros(shape=(4,x.shape[1]))
x0,y0 = self.x0,self.y0
eps = 1e-8
xfilt = np.logical_and( x[0,:] <= x0+1+eps, x[0,:] >= x0-eps)
yfilt = np.logical_and( x[1,:] <= y0+1+eps, x[1,:] >= y0-eps)
filt = np.logical_and(xfilt , yfilt )
values.T[filt] = [0, self.epsxy, self.epsxy, 0]
return values
class epm_phonic():
def __init__(self, gridN,mesh_options = {'cell_type':CellType.quadrilateral}):
self.gridN = N = gridN
#diagonal types: dolfinx.cpp.mesh.DiagonalType.left, left_right, crossed, etc.
self.mesh = mesh = create_rectangle(MPI.COMM_WORLD,[[0,0],[N,N]],[N,N],**mesh_options)#, CellType.quadrilateral)
# self.mesh = mesh = create_rectangle(MPI.COMM_WORLD,[[0,0],[N,N]],[N,N], CellType.quadrilateral)
#==material parameters
self.mu = mu = Constant(mesh,PETSc.ScalarType(1.)) #lame parameter
self.lmbda = lmbda = Constant(mesh,PETSc.ScalarType(2.)) #lame parameter
self.rho = rho = Constant(mesh,PETSc.ScalarType(1.)) #density
self.Gamma = Gamma = Constant(mesh,PETSc.ScalarType(10.0)) #damping
self.tau = tau = 1e-3 #rearrangement time
#==Integration parameters:
#newmark beta parameters:
self.beta,self.gamma = beta,gamma = Constant(mesh,PETSc.ScalarType(0.25)),Constant(mesh,PETSc.ScalarType(0.5))
#probably need to be careful about setting this, but we can make this adjustable.
self.dt = dt = Constant(mesh,PETSc.ScalarType(3e-4))
self.t = 0.0
#okay, we probably need to set up our appropriate forms, finite element spaces, etc.
self.U = U = FunctionSpace(mesh,('DG',0))
self.V = V = VectorFunctionSpace(mesh, ("Lagrange", 1))
self.W = W = TensorFunctionSpace(mesh,('DG',0),shape = (2,2))
#building the associated trial and test functions:
U_tr, U_test = ufl.TrialFunction(U), ufl.TestFunction(U)
V_tr,V_test = ufl.TrialFunction(V),ufl.TestFunction(V)
W_tr,W_test = ufl.TrialFunction(W),ufl.TestFunction(W)
#building the state variables for the EPm
self.u = u = Function(V) #displacements
self.ubar, self.du, self.ddu = ubar,du,ddu = Function(V),Function(V),Function(V) #displacement velocities, accel. etc.
#eigenstrain rules:
self.eigenstrain = eigenstrain = Function(W)
self.eigenstrain_inclusion_projector = Function(W) #temporary variable to project eigenstrain inclusions into
self.incident_eigenstrain = np.zeros(shape = eigenstrain.vector[:].shape)
#initializing the boundary conditions. Must be done before
self.set_bc_simple(0.0)
#=== building the forms corresponding to time-stepping ===
#=== building the static elasticity equations===
self.a_static = a_static = form((ufl.inner(mu*ufl.sym(ufl.grad(V_tr)) + lmbda*ufl.div(V_tr)*ufl.Identity(2),grad(V_test))) *dx)
self.L_static = L_static = form((ufl.inner(2*mu*eigenstrain + lmbda * ufl.Identity(2) * ufl.tr(eigenstrain), grad(V_test) )) *dx)
#perhaps we need ot set up the 'assemble A_static' as a function?
self.A_static = A_static = assemble_matrix(a_static,bcs=[self.bc])
A_static.assemble()
self.update_b_static()
#building the static solver:
opts = PETSc.Options()
opts['ksp_type'] = 'cg'
opts['ksp_atol'] = 1e-14
opts['ksp_rtol'] = 1e-8
opts['ksp_max_it'] = 1000
self.solver_static = solver_static = PETSc.KSP().create(mesh.comm)
solver_static.setFromOptions()
solver_static.setOperators(A_static)
#=== Building the sigma_xy solver ===
sigma_result = sigma(u,eigenstrain,mu,lmbda)
self.L_sig_xy = L_sig_xy = form(inner(sigma_result[0,1],U_test)*dx)
self.a_sig_xy = a_sig_xy = form(inner(U_tr,U_test)*dx)
self.A_sig_xy = A_sig_xy = assemble_matrix(a_sig_xy,bcs=[])
A_sig_xy.assemble()
self.update_b_sig_xy()
opts = PETSc.Options()
opts['ksp_type'] = 'cg'
opts['ksp_atol'] = 1e-14
opts['ksp_rtol'] = 1e-8
opts['ksp_max_it'] = 1000
self.U_sig_xy_solver = U_sig_xy_solver = PETSc.KSP().create(mesh.comm)
U_sig_xy_solver.setFromOptions()
U_sig_xy_solver.setOperators(A_sig_xy)
self.sig_xy_U_vector = sig_xy_U_vector = Function(U)
def add_inclusion_xy_instant(self,magn,x,y):
incl = single_inclusion(magn,x,y)
self.eigenstrain_inclusion_projector.interpolate(incl) #it looks like we can limit which cells we interpolate over. That might be a useful optimization later.
self.eigenstrain.vector[:] += self.eigenstrain_inclusion_projector.vector[:]
self.eigenstrain.vector.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
def perform_timestep(self):
pass
def set_bc_simple(self,epsxy):
self.u_bc_displacements = u_bc_displacements = Function(self.V)
u_bc_displacements.interpolate(simple_shear_bc(epsxy))
tdim = self.mesh.topology.dim #topological dimensions
fdim = tdim - 1 #facet dimension
self.mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(self.mesh.topology)
boundary_dofs = locate_dofs_topological(self.V, fdim, boundary_facets)
self.bc = dirichletbc(u_bc_displacements,boundary_dofs)
def update_b_static(self):
self.b_static = b_static = assemble_vector(self.L_static)
apply_lifting(self.b_static, [self.a_static], bcs=[[self.bc]])
b_static.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b_static, [self.bc])
def solve_static(self):
self.update_b_static()
self.solver_static.solve(self.b_static,self.u.vector)
def update_b_sig_xy(self):
self.b_sig_xy = b_sig_xy=assemble_vector(self.L_sig_xy)
b_sig_xy.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
def solve_sig_xy(self):
self.update_b_sig_xy()
self.U_sig_xy_solver.solve(self.b_sig_xy,self.sig_xy_U_vector.vector)