I am trying to solve a simple linear-elasticity problem (as a starting point for a more complex model) and I am encountering an error that I cannot sort out. When I specify the stress in a DG0-tensor-valued function space with a CG1-vector-valued function space for displacement, everything works as expected. However, if I use a (DG1,CG2) pairing, I am getting the following error (with no other output):
double free or corruption (!prev)
Aborted (core dumped)
which suggests some sort of memory/allocation issue in the linear solver I am using.
I am currently running Dolfinx through a Conda environment. The version is 0.5.2
based on GIT commit: e7e3a2c1e7a26c251669d747686973628dc53593 of GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment
I am running this on Linux Mint 20.3 Cinnamon.
Below is the code which I believe is an MWE for this problem:
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx import fem
from dolfinx.fem import (Function, FunctionSpace, dirichletbc,
form, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx import mesh
from dolfinx.mesh import locate_entities_boundary, meshtags, create_rectangle
#Create rectangular geometry with a triangular meshing
msh = create_rectangle(comm=MPI.COMM_WORLD,
points=((-10.0, -10.0), (10.0, 10.0)), n=(32, 32),
cell_type=mesh.CellType.triangle,)
# Stress (Se) and displacement (Ue) elements
# ------------------------------------------------------------------------ #
# setting degree=1 and degree=2 for Se and Ue respectively leads
# to a core dump/apparent memory allocation issue when attempting to solve
# However, with degree=0 and degree=1, everything works fine, and if one
# writes the solutions to xdmf files, they look reasonable.
# ------------------------------------------------------------------------ #
Se = ufl.TensorElement("DG", msh.ufl_cell(), degree=1, shape=(2,2))
Ue = ufl.VectorElement("Lagrange", msh.ufl_cell(), degree=2)
#Create mixed stress-displacement function space, V
Te = ufl.MixedElement([Ue,Se])
V = FunctionSpace(msh,Te)
#Extract displacement space from V for use in setting boundary conditions
U0, _ = V.sub(0).collapse()
# Locate all facets at the free end and assign them value 1. Sort the
# facet indices (requirement for constructing MeshTags)
free_end_facets = np.sort(locate_entities_boundary(msh, 1, lambda x: np.isclose(x[1], 10.0)))
mt = meshtags(msh, 1, free_end_facets, 1)
ds = ufl.Measure("ds", subdomain_data=mt)
x = ufl.SpatialCoordinate(msh)
# Homogeneous boundary condition in displacement
u_bc1 = Function(U0)
u_expr = fem.Expression(ufl.as_vector((0.001*x[1],0.0)), U0.element.interpolation_points())
u_bc1.interpolate(u_expr)
# Displacement BC is applied to the left side
left_facets = locate_entities_boundary(msh, 1, lambda x: np.isclose(x[0], -10.0))
bdofs = locate_dofs_topological((V.sub(0),U0), 1, left_facets)
bcL = dirichletbc(u_bc1, bdofs,V.sub(0))
bc = [bcL]
# Elastic stiffness tensor and Poisson ratio
E, nu = 10.0, 1.0 / 3.0
def sigma_u(u):
"""Constitutive relation for stress-strain. Assuming plane-stress in XY"""
eps = 0.5 * (ufl.grad(u) + ufl.grad(u).T)
Sigma = E / (1. - nu ** 2) * ((1. - nu) * eps + nu * ufl.Identity(2) * ufl.tr(eps))
return Sigma
#Extract test and trial functions from the mixed space
v, tau = ufl.TestFunctions(V)
u, sigma = ufl.TrialFunctions(V)
#Create ufl representation of variational forms
a00 = ufl.inner(sigma, tau) * ufl.dx
a10 = - ufl.inner(sigma, ufl.grad(v)) * ufl.dx
a01 = - ufl.inner(sigma_u(u), tau) * ufl.dx
#Construct linear forms (e.g. right-hand side)
F = fem.Constant(msh, ScalarType((0,1)))
L = (-ufl.inner(F, v) * ds(1))
#Combine the different blocks
a = (a00+a10+a01)
#Solve monolithic (combined displacement-stress) problem
problem = fem.petsc.LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
Thank you in advance for you help with this issue.