Hi all,
I’m trying to leverage the entity_maps introduced in dolfinx 0.8.0 to solve a poisson equation on a unit square but with three submeshes codim 0.
Ultimately, I would like to model interface discontinuities using a DG method but have CG elements within the domains to obtain this:
For now, I’m simply trying to understand how to solve a poisson equation (without these discontinuities) on three submeshes.
This is the code that I have:
from dolfinx import mesh, fem
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from mpi4py import MPI
import ufl
from ufl import grad, dot, jump, avg
import numpy as np
from dolfinx.nls.petsc import NewtonSolver
msh = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
V = fem.functionspace(msh, ("DG", 1))
tol = 1e-16
def left_domain(x):
return x[0] < 0.5 + tol
def top_right_domain(x):
return np.logical_and(x[0] > 0.5 - tol, x[1] > 0.5 - tol)
def bottom_right_domain(x):
return np.logical_and(x[0] > 0.5 - tol, x[1] < 0.5 + tol)
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
c_to_f = msh.topology.connectivity(tdim, tdim - 1)
f_to_c = msh.topology.connectivity(tdim - 1, tdim)
facet_map = msh.topology.index_map(tdim - 1)
cell_imap = msh.topology.index_map(tdim)
left_cells = mesh.locate_entities(msh, tdim, left_domain)
top_right_cells = mesh.locate_entities(msh, tdim, top_right_domain)
bottom_right_cells = mesh.locate_entities(msh, tdim, bottom_right_domain)
left_submesh, left_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
msh, msh.topology.dim, left_cells
)
top_right_submesh, top_right_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
msh, msh.topology.dim, top_right_cells
)
bottom_right_submesh, bottom_right_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
msh, msh.topology.dim, bottom_right_cells
)
# Create three function spaces on each submesh
V_left = fem.functionspace(left_submesh, ("CG", 1))
V_top_right = fem.functionspace(top_right_submesh, ("CG", 1))
V_bottom_right = fem.functionspace(bottom_right_submesh, ("CG", 1))
# create functions and test functions
u_left = fem.Function(V_left)
v_left = ufl.TestFunction(V_left)
u_top_right = fem.Function(V_top_right)
v_top_right = ufl.TestFunction(V_top_right)
u_bottom_right = fem.Function(V_bottom_right)
v_bottom_right = ufl.TestFunction(V_bottom_right)
# create dx measures
cell_map = msh.topology.index_map(msh.topology.dim)
num_cells = cell_map.size_local + cell_map.num_ghosts
marker = np.arange(num_cells, dtype=np.int32)
mt = dolfinx.mesh.meshtags(msh, msh.topology.dim, marker, marker)
dx = ufl.Measure("dx", domain=msh, subdomain_data=mt)
dx_left = ufl.Measure("dx", domain=left_submesh)
dx_top_right = ufl.Measure("dx", domain=top_right_submesh)
dx_bottom_right = ufl.Measure("dx", domain=bottom_right_submesh)
# Define variational problem
F = 0
F += dot(grad(v_left), grad(u_left)) * dx_left
F += dot(grad(v_top_right), grad(u_top_right)) * dx_top_right
F += dot(grad(v_bottom_right), grad(u_bottom_right)) * dx_bottom_right
# source
f = 1
F += -v_left * f * dx_left
F += -v_top_right * f * dx_top_right
F += -v_bottom_right * f * dx_bottom_right
# boundary conditions
boundary_facets = mesh.exterior_facet_indices(msh.topology)
boundary_dofs = fem.locate_dofs_topological(V, tdim - 1, boundary_facets)
uD = fem.Function(V)
uD.interpolate(lambda x: np.full(x[0].shape, 0.0))
bc = fem.dirichletbc(uD, boundary_dofs)
# TODO add dS coupling terms for interfaces
# make entity maps
parent_to_sub_left = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_left[left_sub_to_parent] = np.arange(
len(left_sub_to_parent), dtype=np.int32
)
entity_map = {left_submesh._cpp_object: parent_to_sub_left}
parent_to_sub_right = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_right[top_right_sub_to_parent] = np.arange(
len(top_right_sub_to_parent), dtype=np.int32
)
entity_map.update({top_right_submesh._cpp_object: parent_to_sub_right})
parent_to_sub_bottom = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_bottom[bottom_right_sub_to_parent] = np.arange(
len(bottom_right_sub_to_parent), dtype=np.int32
)
entity_map.update({bottom_right_submesh._cpp_object: parent_to_sub_bottom})
# NonLinearProblem
formulation = dolfinx.fem.form(F == 0, entity_maps=entity_map)
u = fem.Function(V)
problem = fem.petsc.NonlinearProblem(
formulation,
u,
bcs=bc,
)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(u)
At the moment the code doesn’t run and produces the following error:
Traceback (most recent call last):
File "/home/remidm/FESTIM/interface_submeshes.py", line 118, in <module>
problem = fem.petsc.NonlinearProblem(
^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 905, in __init__
J = ufl.derivative(F, u, du)
^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/ufl/formoperators.py", line 404, in derivative
raise ValueError(f"Invalid argument type {type(form)}.")
ValueError: Invalid argument type <class 'ufl.equation.Equation'>.
I couldn’t find a demo for this application case so I was wondering if the approach was at least the right one?
Thanks for any pointers!
Cheers
Rem
