Hello everyone,
I am trying to interpolate a solution of a coarse mesh into a fine mesh. This does not seem to hard, however, somehow I keep getting the following error: “PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation”. I run my code on a dolfinx 0.8.0.
My MWE is the Poisson equation example from https://jsdokken.com/dolfinx-tutorial/chapter1/fundamentals_code.html where I try to interpolate the solution to a finer mesh. It looks as follows:
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import functionspace, function, Function, locate_dofs_topological,dirichletbc, Constant
import numpy
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = functionspace(domain, ("Lagrange", 1))
uD = Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = locate_dofs_topological(V, fdim, boundary_facets)
bc = dirichletbc(uD, boundary_dofs)
import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
from dolfinx import default_scalar_type
f = Constant(domain, default_scalar_type(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = Function(V)
uh.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
fine_domain = mesh.create_unit_square(MPI.COMM_WORLD, 16, 16, mesh.CellType.quadrilateral)
fine_V = functionspace(fine_domain, ("Lagrange",1))
fine_uh = Function(fine_V)
fine_uh.interpolate(uh)
Can someone please show me why this happens? I have been trying to solve it for the past days but I cannot find the reason for the error.
Many thanks in advance!
There are dedicated functionalities for interpolation on different meshes.
A starting point is for instance
V = functionspace(mesh, mixed_element([vlag_el, lagr_el]))
Eb_m = Function(V)
Eb_m.sub(1).interpolate(g)
diff = Eb_m[2].dx(1) - dgdy
assert assemble_scalar(form(ufl.dot(diff, diff) * ufl.dx)) == pytest.approx(error)
@pytest.mark.parametrize("xtype", [np.float64])
@pytest.mark.parametrize("cell_type0", [CellType.hexahedron, CellType.tetrahedron])
@pytest.mark.parametrize("cell_type1", [CellType.triangle, CellType.quadrilateral])
def test_nonmatching_mesh_interpolation(xtype, cell_type0, cell_type1):
mesh0 = create_unit_cube(MPI.COMM_WORLD, 5, 6, 7, cell_type=cell_type0, dtype=xtype)
mesh1 = create_unit_square(MPI.COMM_WORLD, 5, 4, cell_type=cell_type1, dtype=xtype)
def f(x):
return (7 * x[1], 3 * x[0], x[2] + 0.4)
el0 = element("Lagrange", mesh0.basix_cell(), 1, shape=(3,), dtype=xtype)
V0 = functionspace(mesh0, el0)
el1 = element("Lagrange", mesh1.basix_cell(), 1, shape=(3,), dtype=xtype)
V1 = functionspace(mesh1, el1)
if you are on the main branch
V = functionspace(mesh, mixed_element([vlag_el, lagr_el]))
Eb_m = Function(V)
Eb_m.sub(1).interpolate(g)
diff = Eb_m[2].dx(1) - dgdy
assert assemble_scalar(form(ufl.dot(diff, diff) * ufl.dx)) == pytest.approx(error)
@pytest.mark.parametrize("xtype", [np.float64])
@pytest.mark.parametrize("cell_type0", [CellType.hexahedron, CellType.tetrahedron])
@pytest.mark.parametrize("cell_type1", [CellType.triangle, CellType.quadrilateral])
def test_nonmatching_mesh_interpolation(xtype, cell_type0, cell_type1):
mesh0 = create_unit_cube(MPI.COMM_WORLD, 5, 6, 7, cell_type=cell_type0, dtype=xtype)
mesh1 = create_unit_square(MPI.COMM_WORLD, 5, 4, cell_type=cell_type1, dtype=xtype)
def f(x):
return (7 * x[1], 3 * x[0], x[2] + 0.4)
el0 = element("Lagrange", mesh0.basix_cell(), 1, shape=(3,))
V0 = functionspace(mesh0, el0)
el1 = element("Lagrange", mesh1.basix_cell(), 1, shape=(3,))
V1 = functionspace(mesh1, el1)
if you are on 0.8.0.
Once you read that, you’ll notice a few keywords (for instance, create_nonmatching_meshes_interpolation_data in 0.8.0) that may help you find several related posts in this forum which will contain other useful examples.