Dear community.
Given a MWE of the HDG Poisson example:
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import fem, mesh
from dolfinx.cpp.mesh import cell_num_entities
from petsc4py import PETSc
def compute_cell_boundary_facets(msh):
"""Compute the integration entities for integrals around the
boundaries of all cells in msh.
Parameters:
msh: The mesh.
Returns:
Facets to integrate over, identified by ``(cell, local facet
index)`` pairs.
"""
tdim = msh.topology.dim
fdim = tdim - 1
n_f = cell_num_entities(msh.topology.cell_type, fdim)
n_c = msh.topology.index_map(tdim).size_local
return np.vstack((np.repeat(np.arange(n_c), n_f), np.tile(np.arange(n_f), n_c))).T.flatten()
comm = MPI.COMM_WORLD
rank = comm.rank
dtype = PETSc.ScalarType
# Number of elements in each direction
n = 8
# Create the mesh
msh = mesh.create_unit_cube(comm, n, n, n, ghost_mode=mesh.GhostMode.none)
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_entities(fdim)
facet_imap = msh.topology.index_map(fdim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
facets = np.arange(num_facets, dtype=np.int32)
facet_mesh, facet_mesh_to_mesh = mesh.create_submesh(msh, fdim, facets)[:2]
# Define function spaces
k = 3 # Polynomial order
V = fem.functionspace(msh, ("Discontinuous Lagrange", k))
Vbar = fem.functionspace(facet_mesh, ("Discontinuous Lagrange", k))
# Trial and test functions in mixed space
W = ufl.MixedFunctionSpace(V, Vbar)
u, ubar = ufl.TrialFunctions(W)
v, vbar = ufl.TestFunctions(W)
# Define integration measures
# Cell
dx_c = ufl.Measure("dx", domain=msh)
# Cell boundaries
# We need to define an integration measure to integrate around the
# boundary of each cell. The integration entities can be computed
# using the following convenience function.
cell_boundary_facets = compute_cell_boundary_facets(msh)
cell_boundaries = 1 # A tag
# Create the measure
ds_c = ufl.Measure("ds", subdomain_data=[(cell_boundaries, cell_boundary_facets)], domain=msh)
# Create a cell integral measure over the facet mesh
dx_f = ufl.Measure("dx", domain=facet_mesh)
# We write the mixed domain forms as integrals over msh. Hence, we must
# provide a map from facets in msh to cells in facet_mesh. This is the
# 'inverse' of facet_mesh_to_mesh, which we compute as follows:
mesh_to_facet_mesh = np.full(num_facets, -1)
mesh_to_facet_mesh[facet_mesh_to_mesh] = np.arange(len(facet_mesh_to_mesh))
entity_maps = {facet_mesh: mesh_to_facet_mesh}
How can we modify ds_c
, and hence subdomain_data
in order to exclude the cell facets that coincide with the domain boundary?.
Thanks in advance.