# Differences between dolfin and dolfinx

Hi everyone,

I was trying to solve a simple hyperelastic problem with prescribed facet displacements and was noticing discrepancies between the results from legacy dolfin and dolfinx.

I am running the following versions:

• dolfin: version 2019.2.0.dev0
• commit 0d37780440ad11e15f8aeec1057001761fdd2ba7 if relevant,
• dolfinx: version 0.8.0
• commit ef9b7c15e6bb063acba3a568678881bc2f6e3a20.

The physical problem is a unit cube with the left (x=0) face fixed and prescribed displacements in the x direction on the top/bottom (y=1, y=0) and front/back (z=1, z=0) faces. The prescribed displacement is u(x,y,z) = x \cdot u_0, where u_0 = 0.05.

Legacy dolfin code:

from dolfin import (
UnitCubeMesh, NonlinearVariationalProblem, NonlinearVariationalSolver, MPI,
parameters, set_log_level,
FunctionSpace, Function, LocalSolver,
grad, inner, det, Identity, VectorElement, FiniteElement, MixedElement, split,
TrialFunction, TestFunction, derivative, Constant,
dx, CompiledSubDomain, Expression,
DirichletBC, XDMFFile
)

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["linear_algebra_backend"] = "PETSc"
ffc_options = {"optimize": True, "eliminate_zeros": True, "precompute_basis_const": True,"precompute_ip_const": True}
set_log_level(10)
comm = MPI.comm_world

mesh = UnitCubeMesh(2,2,2)

Wu = VectorElement("CG", mesh.ufl_cell(), 2)
Wp = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(
mesh, MixedElement([Wu, Wp])
)
u_p = Function(V)
v_q = TestFunction(V)
du_dp = TrialFunction(V)
_u, _p = split(u_p)
I1 = inner(_F, _F)
J = det(_F)

mu = Constant(0.1)
energy = _p * (J - 1) + mu * (I1 - 3)
energy *= dx
residual = derivative(energy, u_p, v_q)
jacobian = derivative(residual, u_p, du_dp)

# Mark boundary subdomains
bottom_fc=CompiledSubDomain("near(x[1],0.0) && on_boundary")                                 # bottom face
top_fc=CompiledSubDomain("near(x[1],1) && on_boundary")                                      # top face
right_fc=CompiledSubDomain("near(x[0],1) && on_boundary")                                    # right face
left_fc=CompiledSubDomain("near(x[0],0.0) && on_boundary")                                   # left face
front_fc=CompiledSubDomain("near(x[2],1) && on_boundary")                                  # front face
back_fc=CompiledSubDomain("near(x[2],0.0) && on_boundary")                                   # back face

strtch = Constant(0.05)
applied_disp = Expression("x[0] * strtch", strtch=float(strtch.values()[0]), degree=1)
# applied_disp.compute_vertex_values()
bc_left = DirichletBC(V.sub(0), Constant((0., 0., 0.)), left_fc)
bc_top_x = DirichletBC(V.sub(0).sub(0), applied_disp, top_fc)
bc_bottom_x = DirichletBC(V.sub(0).sub(0), applied_disp, bottom_fc)
bc_back= DirichletBC(V.sub(0).sub(0), applied_disp, back_fc)
bc_front_x = DirichletBC(V.sub(0).sub(0), applied_disp, front_fc)

bcs = [
bc_left, bc_top_x, bc_bottom_x, bc_back, bc_front_x
]
problem = NonlinearVariationalProblem(F=residual, u=u_p, bcs=bcs, J=jacobian, form_compiler_parameters=ffc_options)
solver = NonlinearVariationalSolver(problem)
solver.solve()

u, p = u_p.split(True)

with XDMFFile("dolfin_results.xdmf") as f:
f.parameters["flush_output"] = True
f.parameters["functions_share_mesh"]=True
f.parameters["rewrite_function_mesh"] = False
f.write(u, 0.)
f.write(p, 0.)


gives plausible results for the displacement field (x component in the contour below):

whereas the corresponding dolfinx code

import numpy as np
from mpi4py import MPI
from dolfinx.fem import (
Constant,
Function,
functionspace,
dirichletbc,
locate_dofs_topological,
petsc,
)
from dolfinx.mesh import (
CellType,
locate_entities_boundary,
create_box,
)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
from ufl import (
Measure,
inner,
Identity,
det,
TrialFunction,
TestFunction,
derivative,
split,
)
from dolfinx.log import set_log_level, LogLevel

set_log_level(LogLevel.INFO)
comm = MPI.COMM_WORLD
mesh = create_box(
comm,
[np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])],
n=[2, 2, 2],
cell_type=CellType.tetrahedron,
)
ord = 2
elemu = element("P", mesh.basix_cell(), ord, shape=(3,))
elem_u_project = element("P", mesh.basix_cell(), ord - 1, shape=(3,))
elemp = element("P", mesh.basix_cell(), ord - 1)
elem = mixed_element([elemu, elemp])

V = functionspace(mesh, elem)
Vu = functionspace(mesh, elem_u_project)
Vp = functionspace(mesh, elemp)
u_p = Function(V)
v_q = TestFunction(V)
du_dp = TrialFunction(V)
u_h = Function(Vu, name="u")

u, p = split(u_p)

mu = Constant(mesh, default_scalar_type(1.0))
lmbda = Constant(mesh, default_scalar_type(1.0e3))

J = det(F)
I1 = inner(F, F)

energy = p * (J - 1) + mu * (I1 - 3)
energy *= dx
residual_form = derivative(energy, u_p, v_q)
jacobian_form = derivative(residual_form, u_p, du_dp)

def left(x):
return np.isclose(x[0], 0.0)

def right(x):
return np.isclose(x[0], 1)

def front(x):
return np.isclose(x[2], 1)

def back(x):
return np.isclose(x[2], 0.0)

def bottom(x):
return np.isclose(x[1], 0.0)

def top(x):
return np.isclose(x[1], 1)

# facets
left_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
front_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, front)
back_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, back)
bottom_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, bottom)
top_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, top)

boundary_dofs_bottom_x = locate_dofs_topological(
(V.sub(0).sub(0), V), mesh.topology.dim - 1, bottom_boundary_facets
)
boundary_dofs_top_x = locate_dofs_topological(
(V.sub(0).sub(0), V), mesh.topology.dim - 1, top_boundary_facets
)
boundary_dofs_front_x = locate_dofs_topological(
(V.sub(0).sub(0), V), mesh.topology.dim - 1, front_boundary_facets
)
boundary_dofs_back_x = locate_dofs_topological(
(V.sub(0).sub(0), V), mesh.topology.dim - 1, back_boundary_facets
)
boundary_dofs_left = locate_dofs_topological(
(V.sub(0), V), mesh.topology.dim - 1, left_boundary_facets
)

_V, _ = V.sub(0).collapse()
_V0, _ = V.sub(0).sub(0).collapse()
u_zero = Function(_V)
u_zero.x.array[:] = 0.0
fixed_bc = dirichletbc(u_zero, boundary_dofs_left, V.sub(0))

u_bdry = Function(_V0, name="u_bdry")
disp_x = Constant(mesh, default_scalar_type(0.05))

def disp_x_facets(x):
return x[0] * disp_x.value

u_bdry.interpolate(disp_x_facets)

# u_bdry.sub(0).sub(0).interpolate(disp_x_facets)

top_facet_bc = dirichletbc(u_bdry, boundary_dofs_top_x, V.sub(0).sub(0))
bottom_facet_bc = dirichletbc(u_bdry, boundary_dofs_bottom_x, V.sub(0).sub(0))
front_facet_bc = dirichletbc(u_bdry, boundary_dofs_front_x, V.sub(0).sub(0))
back_facet_bc = dirichletbc(u_bdry, boundary_dofs_back_x, V.sub(0).sub(0))
right_bc = dirichletbc(
Constant(mesh, default_scalar_type(0.01)),
locate_dofs_topological(V.sub(0).sub(0), 2, right_boundary_facets),
V.sub(0).sub(0),
)

bcs = [fixed_bc, front_facet_bc, back_facet_bc, bottom_facet_bc, top_facet_bc]

problem = petsc.NonlinearProblem(F=residual_form, u=u_p, bcs=bcs, J=jacobian_form)
solver = NewtonSolver(comm, problem=problem)

solver.solve(u=u_p)
u, p = u_p.split()

u_h.interpolate(u)

import os

trial_dir = os.path.join(os.getcwd(), "dolfinx_results")
os.makedirs(trial_dir, exist_ok=True)
with XDMFFile(comm, os.path.join(trial_dir, "dolfinx_edgeDisp.xdmf"), "w") as f:
f.write_mesh(mesh)
f.write_function(u_h)
# f.write_function(p)


crashes with

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0


Stepping through a debugger, it seems to be happening in apply_lifting inside the solve call, which could suggest a problem with the bcs but its not obvious to me what that is.

Would anyone know what could be going wrong?
Thanks

I think it is due to passing the wrong function spaces in locate_dofs_topological. Although it is still not intuitive to me as to what function spaces should be passed…

Here is a slightly modified working code:

import numpy as np
from mpi4py import MPI
from dolfinx.fem import (
Constant,
Function,
functionspace,
dirichletbc,
locate_dofs_topological,
petsc,
locate_dofs_geometrical,
)
from dolfinx.mesh import (
CellType,
locate_entities_boundary,
create_box,
)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
from ufl import (
Measure,
inner,
Identity,
det,
TrialFunction,
TestFunction,
derivative,
split,
)
from dolfinx.log import set_log_level, LogLevel
from petsc4py import PETSc

opts = PETSc.Options()

set_log_level(LogLevel.INFO)
comm = MPI.COMM_WORLD
Lx, Ly, Lz = 1., 1., 1.
mesh = create_box(
comm,
[np.array([0.0, 0.0, 0.0]), np.array([Lx, Ly, Lz])],
n=[3, 10, 5],
cell_type=CellType.tetrahedron,
)
ord = 2
elemu = element("P", mesh.basix_cell(), ord, shape=(3,))
elem_u_project = element("P", mesh.basix_cell(), ord - 1, shape=(3,))
elemp = element("P", mesh.basix_cell(), ord - 1)
elem = mixed_element([elemu, elemp])

V = functionspace(mesh, elem)
Vu = functionspace(mesh, elem_u_project)
Vp = functionspace(mesh, elemp)
u_p = Function(V)
v_q = TestFunction(V)
du_dp = TrialFunction(V)
u_h = Function(Vu, name="u")

u, p = split(u_p)

mu = Constant(mesh, default_scalar_type(1.0))
lmbda = Constant(mesh, default_scalar_type(1.0e3))

J = det(F)
I1 = inner(F, F)

energy = p * (J - 1) + mu * (I1 - 3)
energy *= dx
residual_form = derivative(energy, u_p, v_q)
jacobian_form = derivative(residual_form, u_p, du_dp)

def left(x):
return np.isclose(x[0], 0.0)

def right(x):
return np.isclose(x[0], Lx)

def front(x):
return np.isclose(x[2], Lz)

def back(x):
return np.isclose(x[2], 0.0)

def bottom(x):
return np.isclose(x[1], 0.0)

def top(x):
return np.isclose(x[1], Ly)

# facets
left_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
front_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, front)
back_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, back)
bottom_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, bottom)
top_boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, top)

_V, _ = V.sub(0).collapse()
_V0, _ = V.sub(0).sub(0).collapse()

boundary_dofs_bottom_x = locate_dofs_topological(
(V.sub(0).sub(0), _V0), mesh.topology.dim - 1, bottom_boundary_facets
)
boundary_dofs_top_x = locate_dofs_topological(
(V.sub(0).sub(0), _V0), mesh.topology.dim - 1, top_boundary_facets
)
boundary_dofs_front_x = locate_dofs_topological(
(V.sub(0).sub(0), _V0), mesh.topology.dim - 1, front_boundary_facets
)
boundary_dofs_back_x = locate_dofs_topological(
(V.sub(0).sub(0), _V0), mesh.topology.dim - 1, back_boundary_facets
)

boundary_dofs_left = locate_dofs_topological(
(V.sub(0), _V), mesh.topology.dim - 1, left_boundary_facets
)

u_zero = Function(_V)
u_zero.x.array[:] = 0.0
fixed_bc = dirichletbc(u_zero, boundary_dofs_left, V.sub(0))

u_bdry = Function(_V0, name="u_bdry")
disp_x = Constant(mesh, default_scalar_type(0.05))

def disp_x_facets(x):
return x[0] * disp_x.value

u_bdry.interpolate(disp_x_facets)

# u_bdry.sub(0).sub(0).interpolate(disp_x_facets)

top_facet_bc = dirichletbc(u_bdry, boundary_dofs_top_x, V.sub(0).sub(0))
bottom_facet_bc = dirichletbc(u_bdry, boundary_dofs_bottom_x, V.sub(0).sub(0))
front_facet_bc = dirichletbc(u_bdry, boundary_dofs_front_x, V.sub(0).sub(0))
back_facet_bc = dirichletbc(u_bdry, boundary_dofs_back_x, V.sub(0).sub(0))
right_bc = dirichletbc(
Constant(mesh, default_scalar_type(0.01)),
locate_dofs_topological(V.sub(0).sub(0), 2, right_boundary_facets),
V.sub(0).sub(0),
)

bcs = [fixed_bc, front_facet_bc, back_facet_bc, bottom_facet_bc, top_facet_bc]

problem = petsc.NonlinearProblem(F=residual_form, u=u_p, bcs=bcs, J=jacobian_form)
solver = NewtonSolver(comm, problem=problem)

solver.solve(u=u_p)
u, p = u_p.split()

u_h.interpolate(u)

import os

trial_dir = os.path.join(os.getcwd(), "trialDofs")
os.makedirs(trial_dir, exist_ok=True)
with XDMFFile(comm, os.path.join(trial_dir, "mixedFormulationDofs_Box.xdmf"), "w") as f:
f.write_mesh(mesh)
f.write_function(u_h)


A helpful place to look for me was the tests. Any docs elaborating what do the two function spaces in the V argument represent in locate_dofs_topological would be helpful!

dolfinx is an open source project. If you have ideas on how to improve the documentation feel free to open a pull request.