Exclude domain boundary facets from subdomain_data

dolfinx
1

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.

2

Not sure if I understand the question. ds integrates only on domain boundary facets, what would be point of using it if you want to exclude all of them?

3

He is using the custom integration measure, meaning that he is doing one-sided integrals over all facets from the view of every cell.
You can do this with:

from mpi4py import MPI
import numpy as np
import dolfinx

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 = dolfinx.cpp.mesh.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()


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.cpp.mesh.CellType.triangle)

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
all_facets = compute_cell_boundary_facets(mesh)

n_f = dolfinx.cpp.mesh.cell_num_entities(mesh.topology.cell_type, mesh.topology.dim - 1)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
exterior_entities = dolfinx.fem.compute_integration_domains(dolfinx.fem.IntegralType.exterior_facet, mesh.topology, exterior_facets)
# We can now exclude these entities from "all_facets"
remove_index = n_f*exterior_entities[::2] + exterior_entities[1::2]

remaining_facets = np.delete(all_facets.reshape(-1, 2), remove_index, axis=0).flatten()

# Verify that all remaining facets are interior
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim - 1)
for facet in remaining_facets.reshape(-1,2):
    f = c_to_f.links(facet[0])[facet[1]]
    assert not np.isin(f, exterior_facets)
1 Like
4

Thanks for the response, @dokken.

Do we need to modify/recompute the following?

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}

Also, I was thinking if we can define the interior skeleton mesh from the top, i.e.,

msh.topology.create_connectivity(fdim, tdim)
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)
exterior_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)
interior_facets= np.setdiff1d(facets, exterior_facets)

facet_mesh, facet_mesh_to_mesh = mesh.create_submesh(msh, fdim, interior_facets)[:2]

but then Iā€™m not sure about the computation of mesh_to_facet_mesh and entity_maps.

5

You do not need to modify the entity maps. Given the integration entities above:

it will look up the corresponding submesh cell in the facet mesh (using entity map).

1 Like
6

Then your function space would have no degrees of freedom on the top facets. Is that what you want? If yes, you could exclude it from the submesh creation. The mesh_to_facet_mesh and entity_maps constructions should really change. However, then you cannot involve a quantity from that submesh on an integral on the top boundary.

1 Like
7

Thanks, @dokken.

As a final extension of this solved post, I wonder how can we extend

to an external GMSH mesh with marked facets?. I gladly open a new post if required.

Cheers.

8

That would just be modifying these lines to only use selected/marked facets read in as a facet tag to dolfinx. For instance exterior entities could be facet_tag.find(marker_value)

1 Like