How to split the outer boundary of a mesh in a structured way

Question

Section 1. An overview of the problem

I have a planar mesh my_mesh which incorporates an internal boundary and whose outer boundary is a regular polygon centered at the origin.

My current task is to split the boundary of my_mesh in six arcs using a dolfin.MeshFunction, where there should be three boundary arcs $P_0, P_2, P_4$ with an angular aperture of $\frac{\pi}{6}$ and another three boundary arcs $P_1, P_3, P_5$ with an angular aperture of $\frac{\pi}{2}$.

In Section 2, I have included my attempt to implement this; the code also computes the lengths of $P_0, P_1, P_2, P_3, P_4, P_5$. According to the precise geometric description of my_mesh in Section 3, it should turn out that the lengths of $P_1, P_3, P_5$ are equal. However, the code in Section 2 computes:

• the length of $P_1$ as 1.55291427062
• the length of $P_3$ as 1.55291427062
• the length of $P_5$ as 1.48820950934

My question: If you could help me understand how to re-design my code so that $P_1, P_3, P_5$ are computed to have the same length, I would be very grateful.

Section 2. The implementation

# ------------------------
# BEGIN: import statements

import mshr
import dolfin

import numpy

# END: import statements
# ----------------------



# ------------------------
# BEGIN: class definitions

class BoundaryArc(dolfin.SubDomain):
    """ Given a planar domain P which is star-shaped
        with respect to a point x_0 in its interior,
        this class represents and manipulates 
        connected subsets of the boundary of P,
        where each such connected subset is
        specified as an angular aperture.

        Angles are ASSUMED to be measured:

            1) in the CCW direction from (the parallel translate of)
            the positive x-axis originating at x_0

            2) and in radians."""

    def __init__(self, theta_initial, theta_final):

        self.theta_initial = theta_initial
        self.theta_final = theta_final
        dolfin.SubDomain.__init__(self)


    def inside(self, x, on_boundary):

        angle = my_arctan(x[1], x[0], angular_unit = TYPE_ANGULAR_UNITS)
        is_angle_in_arc_aperture = (angle >= self.theta_initial - DOLFIN_TOL) \
                                and (angle <= self.theta_final + DOLFIN_TOL)

        return (on_boundary and is_angle_in_arc_aperture)

# END: class definitions
# ----------------------



# ---------------------------
# BEGIN: function definitions

def my_arctan(y, x, angular_unit = 'radians'):

    angle = numpy.arctan2(y, x)
    if (angle < 0):
        angle += 2 * numpy.pi

    if (angular_unit == 'radians'):
        return angle
    elif (angular_unit == 'degrees'):
        return angle * (180/numpy.pi)

# END: function definitions
# -------------------------



if __name__ == '__main__':

    MAX_DEGREES = 360.0
    TYPE_ANGULAR_UNITS = 'degrees'
    CONV_FCT_DEG_TO_RAD = numpy.pi / 180.0

    DOLFIN_TOL = 1E-6


    # construct a regular polygon P centered at the origin,
    # on whose interior we are to create an unstructured mesh my_mesh later
    number_of_vertices = 12
    ang_pos_of_vertices = \
        numpy.linspace(0, MAX_DEGREES, num = number_of_vertices, endpoint = False)

    x_coords_of_vertices = [numpy.cos(t * CONV_FCT_DEG_TO_RAD) for t in ang_pos_of_vertices]
    y_coords_of_vertices = [numpy.sin(t * CONV_FCT_DEG_TO_RAD) for t in ang_pos_of_vertices]

    coords_of_vertices = zip(x_coords_of_vertices, y_coords_of_vertices)
    Vertices_Of_P = [dolfin.Point(x, y) for x, y in coords_of_vertices]
    Outer_Bdry_Of_P = mshr.Polygon(Vertices_Of_P)

    # construct an internal boundary 
    # to be incorporated into the unstructured mesh my_mesh on the interior of P,
    # which we are to construct later
    num_vertices_internal_bdry = 32
    ang_pos_of_vertices_internal_bdry = \
        numpy.linspace(0, MAX_DEGREES, num = num_vertices_internal_bdry, endpoint = False)

    x_coords_of_vertices_internal_bdry = \
                            [0.35 + 0.25 * numpy.cos(t * CONV_FCT_DEG_TO_RAD) 
                                    for t in ang_pos_of_vertices_internal_bdry]
    y_coords_of_vertices_internal_bdry = \
                            [0.35 + 0.25 * numpy.sin(t * CONV_FCT_DEG_TO_RAD)
                                    for t in ang_pos_of_vertices_internal_bdry]

    coords_of_vertices_internal_bdry = \
        zip(x_coords_of_vertices_internal_bdry, y_coords_of_vertices_internal_bdry)
    Vertices_Of_Internal_Bdry = [dolfin.Point(x, y) for x, y in coords_of_vertices_internal_bdry]
    Internal_Bdry = mshr.Polygon(Vertices_Of_Internal_Bdry)

    # generate the mesh my_mesh
    mesh_resol_param = 15
    Outer_Bdry_Of_P.set_subdomain(1, Internal_Bdry)
    my_mesh = mshr.generate_mesh(Outer_Bdry_Of_P, mesh_resol_param)

    # mark boundary arcs on the outer boundary of my_mesh
    # by means of a dolfin.MeshFunction;
    # each boundary arc will be specified by its angular aperture
    # (measured in the counter-clockwise sense from the positive x-axis)
    angles = [0.0, 30.0, 120.0, 150.0, 240.0, 270.0]
    mesh_bdry_fn = dolfin.MeshFunction("size_t", my_mesh,
                                my_mesh.topology().dim()-1)

    for j in range(len(angles)-1):
        print 'angular aperture of boundary arc #{0}'.format(j), \
                ':', angles[j], angles[j+1]
        Current_Bdry_Arc = BoundaryArc(angles[j], angles[j+1])
        Current_Bdry_Arc.mark(mesh_bdry_fn, j)

    j = len(angles)-1
    print 'angular aperture of boundary arc #{0}'.format(j), \
                ':', angles[j], MAX_DEGREES
    Current_Bdry_Arc = BoundaryArc(angles[j], MAX_DEGREES)
    Current_Bdry_Arc.mark(mesh_bdry_fn, j)


    # compute the lengths of each boundary arc
    ds = dolfin.Measure("ds", domain = my_mesh, subdomain_data = mesh_bdry_fn)

    number_bdry_arcs = len(angles)

    print
    for j in range(number_bdry_arcs):
        length_of_curr_bdry_arc = dolfin.Constant(1) * ds(j)
        length_of_curr_bdry_arc = dolfin.assemble(length_of_curr_bdry_arc)

        print 'length of boundary arc {0}:'.format(j), length_of_curr_bdry_arc

Section 3. A precise geometric description of my_mesh

I would like to use the mshr module to:

• create a mesh my_mesh whose outer boundary is the boundary $Q$ of the regular polygon inscribed in the unit circle with 12 vertices and one vertex of $Q$ is at the point $(1, 0)$
• incorporate an internal polygonal boundary $Q_{internal}$ into my_mesh, where $Q_{internal}$ is the boundary of the regular polygon inscribed in the circle $\partial B( (0.35, 0.35); 0.25)$ and one vertex of $Q_{internal}$ is at the point $(0.35 + 0.25, 0.35)$
• use a dolfin.MeshFunction to mark parts of the outer boundary of my_mesh (i.e. parts of the boundary of $Q$) according to their angular aperture as measured in the counter-clockwise sense from the positive $x$-axis; specifically, I would like to divide the outer boundary of my_mesh into the following six angular apertures: [0, 30], [30, 120], [120, 150], [150, 240], [240, 270], [270, 360], where all angles are in degrees and all endpoints are included

By applying the law of sines to the mathematical situation described above, one sees that:

• each of the boundary parts corresponding to angular apertures [0, 30], [120, 150], [240, 270] has a length of $sin(\pi/6) / sin(5\pi / 12) \approx 0.517638$
• each of the boundary parts corresponding to angular apertures [30, 120], [150, 240], [270, 360] has a length of $3 * sin(\pi/6) / sin(5\pi / 12) \approx 1.552914$