Tag inner and outer boundary

Question

Hello,
I have demonstrative mesh bordered with outer and inner closed loop. Please, is it possible to apply different Dirichlet boundary conditions on each boundary in Fenics way?

Thank you,

Petr

Answer 1

The following code will do an automated marking based on a threshold angle between facets. It will also mark disconnected parts of the boundary (i.e. just what you need).

Edit: Bug fix.

from dolfin import *
import numpy as np

def angle_between(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'::

            >>> angle_between((1, 0, 0), (0, 1, 0))
            1.5707963267948966
            >>> angle_between((1, 0, 0), (1, 0, 0))
            0.0
            >>> angle_between((1, 0, 0), (-1, 0, 0))
            3.141592653589793
    """
    v1 = v1/np.linalg.norm(v1)
    v2 = v2/np.linalg.norm(v2)
    angle = np.arccos(np.dot(v1, v2))
    if np.isnan(angle):
        if (v1_u == v2_u).all():
            return 0.0
        else:
            return np.pi
    return angle

def _mark_loop(bmesh, normals, threshold_angles, domains, queue, visited, marker, threshold=30):
    queue_indices = np.where(queue.array())[0]

    num_marked = 0
    emap = bmesh.entity_map(2)
    while len(queue_indices) > 0:
        idx = queue_indices[0]
        facet = Cell(bmesh, idx)
        #FACET = Facet(mesh, emap[idx])
        visited[facet] = True
        queue[facet] = False
        domains[facet] = marker
        N = normals[idx]
        threshold_radians = np.pi*threshold_angles[facet]/180.

        neighbours = []
        for edge in edges(facet):
            for neighbour in cells(edge):
                if not visited[neighbour]:
                    neighbours.append(neighbour)

        for neighbour in neighbours:
            n = normals[neighbour.index()]

            angle = angle_between(N,n)
            if angle < threshold_radians:
                queue[neighbour] = not visited[neighbour]
        queue_indices = np.where(queue.array())[0]
        num_marked += 1

def mark_domains(mesh, threshold_angle=30):
    """Mark domains automatically by feature edges in mesh. Will connect
    facets to the same domain if angle between them is less than threshold_angle.
    Will also detect and mark disconnected parts of boundary.
    """
    assert MPI.size(mesh.mpi_comm()) == 1, "Currently only works on non-partitioned meshes."

    bmesh = BoundaryMesh(mesh, "exterior")
    visited_facets = CellFunction("bool", bmesh)
    visited_facets.set_all(False)
    bmesh.init(1,2)

    queue = CellFunction("bool", bmesh)
    queue.set_all(False)

    domains = CellFunction("size_t", bmesh)
    domains.set_all(0)

    normals = np.zeros((bmesh.num_cells(), 3), dtype=np.float_)
    for facet in cells(bmesh):
        N = Facet(mesh, bmesh.entity_map(2)[facet.index()]).normal()
        normals[facet.index()][0] = N.x()
        normals[facet.index()][1] = N.y()
        normals[facet.index()][2] = N.z()


    threshold_angles = CellFunction("size_t", bmesh)
    threshold_angles.set_all(threshold_angle)

    unvisited = np.where(visited_facets.array() == False)[0]
    I = 0
    while len(unvisited) > 0:
        I += 1
        queue[unvisited[0]] = True
        _mark_loop(bmesh, normals, threshold_angles, domains, queue, visited_facets, I)
        unvisited = np.where(visited_facets.array() == False)[0]

    domains_full = MeshFunction("size_t", mesh, 2)
    domains_full.set_all(0)
    for facet in cells(bmesh):
        domains_full[bmesh.entity_map(2)[facet.index()]] = domains[facet]


    return domains_full


if __name__ == '__main__':
    try:
        import mshr
        domain = mshr.Box(Point(0.0,0.0,0.0), Point(1.0,1.0,1.0)) - mshr.Sphere(Point(0.5,0.5,0.5), 0.2)
        mesh = mshr.generate_mesh(domain, 20)
    except:
        mesh = UnitSquareMesh(6,6,6)

    domains = mark_domains(mesh,30)
    File("domains.pvd") << domains

Comments

Hi Øyvind,
that's exactly what I was looking for. Thank you very much!

Answer 2

The easiest way is to mark boundaries using a FacetFunction which has different values for each boundary. You will need to fill the FacetFunction by iterating over the mesh, and deciding which facets belong to which boundary.

Comments

Thank you Chris,

Now, I understand the concept of FacetFunction and exterior() method. The method exterior() gives all boundaries, but is there any way how to decide whether boundary is inner or outer? Or should I know it explicitly (from meshing process)?

---

You need to have some way of describing which boundary is which. You may be able to do this from the geometry (i.e. by knowing the range of x[0], x[1], x[2] etc. for your boundaries). In a more complex case, you need to mark the boundaries in your mesh generation software.