constructing SubDomain from Mesh and MeshFunction

Question

In the documentation for SubDomain, it is alluded that a SubDomain can be defined with a Mesh and a MeshFunction.

However, I cannot figure out how to do this. Is there some constructor or factory method I should use?

(my underlying objective is to serialize/deserialize a pointwise DirichletBC, so ideas from that angle are also welcome)

Answer 1

I hope this python code example can help you - let's say you want to "submesh" the upper and lower half of the UnitCube (SubMesh only works in serial):

class Dom1(SubDomain):
   def inside(self, x, on_boundary):
       return x[2] < 0.5 + DOLFIN_EPS

class Dom2(SubDomain):
   def inside(self, x, on_boundary):
       return x[2] > 0.5 - DOLFIN_EPS

mesh = UnitCubeMesh(6, 6, 6)

dom1 = Dom1()
dom2 = Dom2()
domains = CellFunction("size_t", mesh)
domains.set_all(0)
dom1.mark(domains, 1)
dom2.mark(domains, 2)

submesh1 = SubMesh(mesh, domains, 1)
submesh2 = SubMesh(mesh, domains, 2)
df.File('submesh1.xml') << submesh1  # export if needed
df.File('submesh2.xml') << submesh2  #        -- "" --

Comments

Thanks. I modified the approach slightly for use with VertexFunction and individual points.

Pasting below in case it's useful to anyone:

import numpy as np

def pointwise_inside(vertex_fn, index):
    """
    Return a SubDomain-style inside() method based on the part of the specified VertexFunction matching the specified index.
    """

    mesh_coords = vertex_fn.mesh().coordinates()
    num_points = mesh_coords.shape[0]
    fn_values = vertex_fn.array()
    sq_tol = DOLFIN_EPS**2

    def inside(x, on_boundary):
        """
        Check if point at x belongs to this subdomain.

        If x is not sufficiently close to a mesh coordinate, this will always return False.
        """

        # difference between x and mesh coordinates
        difference_vector = np.outer(np.ones(num_points), x) - mesh_coords

        sq_dist = np.sum(difference_vector**2, axis=1)

        min_ind = np.argmin(sq_dist)

        if sq_dist[min_ind] > sq_tol:
            return False

        return fn_values[min_ind] == index

    return inside

# subdomain_points is a VertexFunction, marked with a different index for each subdomain

# e.g., specify boundary 1
boundary = pointwise_inside(subdomain_points, 1)
dbc = DirichletBC(V, value, boundary, method='pointwise')