BoundaryMesh

[c-abird]

Is there currently a way to extract a boundary mesh from a given Mesh object (similar to the BoundaryMesh functionality in FEniCS)? In the best case this would also be possible for the boundary of a subdomain of the mesh.

I want to implement FEM/BEM coupling, so I need both, the coordinates and the connectivity of the boundary mesh. Moreover an efficient way to map the values of functions living on either mesh would be nice. I guess this is somehow connected to #723.

[wence-]

Some of this stuff is available via the exterior facet maps (but it's a bit of work to pull it out). When you say you want the "connectivity" of the boundary mesh, do you want something that behaves like a normal Mesh object, but knows a relationship to its "parent"?

[c-abird]

By connectivity I just mean, how the vertices are connected to triangles. So, something that behaves like a mesh should be alright. I don't need a relationship to the parent mesh. However, I would need some method to pull the data from a Function defined on the full Mesh to a Function defined on the boundary mesh and vice versa (ideally without breaking parallel execution). Could you give me some pointers on how to use exterior facet maps?

[c-abird]

Hi Lawrence,

I managed to put together a basic code for the extraction of a boundary mesh from a given bulk mesh:

import firedrake as fd         
import numpy as np             
from pyop2.datatypes import IntType

mesh = fd.UnitCubeMesh(3,3,3)  
V = fd.FunctionSpace(mesh, "CG", 1) 

# this is adapted from firedrake/dmplex.pyx (boundary_nodes)
facet_dim = mesh.facet_dimension()
mesh_dim = mesh.geometric_dimension()

boundary_dofs = V.finat_element.entity_support_dofs()[facet_dim]

local_nodes = np.empty((len(boundary_dofs), 
                        len(boundary_dofs[0])),         
                       dtype=IntType)                  
for k, v in boundary_dofs.items():
    local_nodes[k, :] = v      

facets = V.mesh().exterior_facets
local_facets = facets.local_facet_dat.data_ro_with_halos

nlocal = local_nodes.shape[1]  
nfacet = facets.set.total_size 

facet_node_list = V.exterior_facet_node_map().values_with_halo

# collect facet data and set up mesh object
facet_data = np.zeros((nfacet, nlocal), dtype=IntType)
for facet in range(nfacet):    
    l = local_nodes[local_facets[facet]]
    facet_data[facet,:] = facet_node_list[facet].take(l)

node_map = np.unique(facet_data.flatten()) 
boundary_facets = np.searchsorted(node_map, facet_data)
boundary_coords = mesh.coordinates.vector().dat.data_ro[node_map]

plex = fd.mesh._from_cell_list(facet_dim, boundary_facets, boundary_coords, mesh.comm)
boundary_mesh = fd.Mesh(plex, dim = mesh_dim, reorder = False)

# create some function on boundary mesh and save it
VB = fd.FunctionSpace(boundary_mesh, "CG", 1)
f = fd.Function(VB).interpolate(fd.Constant(1.))
fd.File("f.pvd").write(f)

However, there are some problems:

• The orientation of the resulting boundary mesh facets is rather arbitrary which is bad for BEM. Is there a simple way to get facet vertex ordering leading to outward pointing facet normals?
• The generation of the node mapping and its application will most certainly break parallel execution

Looking at the API of DMPLEX it should be possible to extract a boundary mesh by using DMPlexMarkBoundaryFaces and DMPlexCreateSubmesh. This looks to me like the best way in order to preserve parallel features of the mesh. However, DMPlexCreateSubmesh is apparently deactivated in petsc4py. Are you involved with the DMPLEX development or do you know why this method is not available in Python?

Thanks a lot and best wishes, Claas

[dham]

I'm not sure about the parallel issues. However the situation with orientation is as follows. We use the FEniCS convention that the triangles are oriented by the indices of the vertices, not geometrically. That is baked in to how we do assembly so isn't easily changeable. However for meshes of immersed manifolds, the mesh does carry orientation information which enables the construction of a consistent normal field. I don't believe it is guaranteed to be outward (after all, the immersed mesh might not be the surface of a closed region).

[c-abird]

OK. I understand that the boundary mesh has an attribute cell_orientations, which can be initialized with init_cell_orientations(expr). However, I did not find a way to initialize this properly. Ideally, I would use the FacetNormal of the full mesh to intialize the cell orientations of the boundary mesh, but this doesn't work as an input for init_cell_orientations. Any hints? Thanks, Claas

[wence-]

Hi Claas,

I think you might be able to do something like this:

V = VectorFunctionSpace(mesh, "DGT", 0)
u = TrialFunction(V)
v = TestFunction(V)
fn = Function(V)
solve(u*v*ds == FacetNormal(mesh)*v*ds, fn, solver_parameters={"pc_type": "jacobi"})

Now you have a field on the trace that contains the facet normal on each exterior facet of the original mesh.

You would then need to extract the correct pieces for your boundary mesh (doing something like what you do above to identify the correct facets).

Then you should be able to do

V = VectorFunctionSpace(boundary_mesh, "DG", 0)
fn_boundary = Function(V)
fn_boundary.dat.data[:] = fn.dat.data_ro[magic_permutation]

boundary_mesh.init_cell_orientations(fn_boundary)

It's trivial to wrap up DMPlexCreateSubmesh in petsc4py. This might be helpful, but I suspect you still need to do a lot of this index-permutation dance.

We'd love to do this properly (as you note #723 is relevant), but it's not completely trivial.

[wence-]

I am hopeful that something along these lines will arise out of @ksagiyam's work on submesh support.

[ksagiyam]

Yes. Boundary mesh extraction related to DMPlexCreateSubmesh is under active development. We are planning to then enable Firedrake to solve couplings between domain and boundary PDEs. Extraction of boundary of a boundary should also arise somewhere in this line.

[c-abird]

Thanks for the update. Is there any issue/branch I can follow in order to stay updated on this feature set? Is there a rough timeframe for the development?

[wence-]

Hi @c-abird, @ksagiyam can probably provide more information. It's presently actively being worked on. We have a project tracking various bits here: https://github.com/firedrakeproject/firedrake/projects/1 (if you have particular requirements, you could perhaps add some notes: it would be useful to know what people who want this stuff need).

The relevant branch in Firedrake is, I think, https://github.com/firedrakeproject/firedrake/tree/dmplex-submesh

[ksagiyam]

Hi @c-abird, ability to define functions on arbitrary subdomains (or subdomains of subdomains, etc) itself will probably be available in a day or two on the firedrake/dmplex-submesh branch that Lawrence mentioned. I think I need a few more weeks to enable us to use full-domain function data on a subdomain. We will keep you posted.

Although it might not be very useful at this moment, you can install this version as:

python3 firedrake-install --package-branch firedrake dmplex-submesh --package-branch petsc dmplex-submesh1 --package-branch petsc4py dmplex-submesh

[c-abird]

Hi @ksagiyam, I saw that you worked quite a bit on the dmplex-submesh branch. So I started trying out a few things. When I try to extract a boundary mesh from a given mesh with

mesh = UnitCubeMesh(2,2,2)
boundary_mesh = SubMesh(mesh, "exterior_facets", 1, "facet")

I get

RuntimeError: Facet must have 1 or 2 supports.  Something is wrong.

I this due to work in progress or am I doing something wrong?