Tutorial 1.2: Mesh Construction¶

This tutorial introduces the construction and manipulation of curvilinear meshes.

A PlanarMesh object is fundamentally a list of curvilinear Edge objects. Each edge consists of a pair of vertices and enough information to parameterize the edge as a curve in the plane. Each edge also is assigned to two mesh cells to which the edge is part of the boundary. One of these cells is taken to be "positive" and the other "negative." For the positive cell, the edge is oriented counterclockwise if it lies on the outer boundary of the cell (i.e. the edge is not the boundary of a hole in the cell), and clockwise otherwise. For the negative cell, the opposite is true. In this fashion, the entire mesh can be constructed from a list of edges.

We begin by importing the puncturedfem package, as well as numpy and matplotlib for the sake of this example.

In [1]:

import puncturedfem as pf
import numpy as np
import matplotlib.pyplot as plt

import puncturedfem as pf import numpy as np import matplotlib.pyplot as plt

Define the vertices¶

First we begin by defining the vertices of the mesh, which we will later join with edges.

In [2]:

verts: list[pf.Vert] = []

# rectangle corners
verts.append(pf.Vert(x=0.0, y=0.0))  # 0
verts.append(pf.Vert(x=1.0, y=0.0))  # 1
verts.append(pf.Vert(x=3.0, y=0.0))  # 2
verts.append(pf.Vert(x=4.0, y=0.0))  # 3
verts.append(pf.Vert(x=4.0, y=1.0))  # 4
verts.append(pf.Vert(x=3.0, y=1.0))  # 5
verts.append(pf.Vert(x=1.0, y=1.0))  # 6
verts.append(pf.Vert(x=0.0, y=1.0))  # 7

# pacman
pacman_scale = 0.4
verts.append(pf.Vert(x=0.5, y=0.5))  # 8
verts.append(
    pf.Vert(x=0.5 + pacman_scale * (np.sqrt(3) / 2), y=0.5 + pacman_scale * 0.5)
)  # 9
verts.append(
    pf.Vert(x=0.5 + pacman_scale * (np.sqrt(3) / 2), y=0.5 - pacman_scale * 0.5)
)  # 10
verts.append(
    pf.Vert(x=0.5 + pacman_scale * -0.1, y=0.5 + pacman_scale * 0.5)
)  # 11

# dots
verts.append(pf.Vert(x=1.5, y=0.5))  # 12
verts.append(pf.Vert(x=2.0, y=0.5))  # 13
verts.append(pf.Vert(x=2.5, y=0.5))  # 14

# ghost
ghost_scale = 0.6
ghost_x_shift = 3.5
ghost_y_shift = 0.5
verts.append(
    pf.Vert(
        x=ghost_x_shift + ghost_scale * (-0.5),
        y=ghost_y_shift + ghost_scale * (-0.6),
    )
)  # 15
verts.append(
    pf.Vert(
        x=ghost_x_shift + ghost_scale * (0.5),
        y=ghost_y_shift + ghost_scale * (-0.6),
    )
)  # 16
verts.append(
    pf.Vert(
        x=ghost_x_shift + ghost_scale * (0.5),
        y=ghost_y_shift + ghost_scale * (0.2),
    )
)  # 17
verts.append(
    pf.Vert(
        x=ghost_x_shift + ghost_scale * (-0.5),
        y=ghost_y_shift + ghost_scale * (0.2),
    )
)  # 18
verts.append(
    pf.Vert(
        x=ghost_x_shift + ghost_scale * (-0.25),
        y=ghost_y_shift + ghost_scale * (0.1),
    )
)  # 19
verts.append(
    pf.Vert(
        x=ghost_x_shift + ghost_scale * (0.25),
        y=ghost_y_shift + ghost_scale * (0.1),
    )
)  # 20

verts: list[pf.Vert] = [] # rectangle corners verts.append(pf.Vert(x=0.0, y=0.0)) # 0 verts.append(pf.Vert(x=1.0, y=0.0)) # 1 verts.append(pf.Vert(x=3.0, y=0.0)) # 2 verts.append(pf.Vert(x=4.0, y=0.0)) # 3 verts.append(pf.Vert(x=4.0, y=1.0)) # 4 verts.append(pf.Vert(x=3.0, y=1.0)) # 5 verts.append(pf.Vert(x=1.0, y=1.0)) # 6 verts.append(pf.Vert(x=0.0, y=1.0)) # 7 # pacman pacman_scale = 0.4 verts.append(pf.Vert(x=0.5, y=0.5)) # 8 verts.append( pf.Vert(x=0.5 + pacman_scale * (np.sqrt(3) / 2), y=0.5 + pacman_scale * 0.5) ) # 9 verts.append( pf.Vert(x=0.5 + pacman_scale * (np.sqrt(3) / 2), y=0.5 - pacman_scale * 0.5) ) # 10 verts.append( pf.Vert(x=0.5 + pacman_scale * -0.1, y=0.5 + pacman_scale * 0.5) ) # 11 # dots verts.append(pf.Vert(x=1.5, y=0.5)) # 12 verts.append(pf.Vert(x=2.0, y=0.5)) # 13 verts.append(pf.Vert(x=2.5, y=0.5)) # 14 # ghost ghost_scale = 0.6 ghost_x_shift = 3.5 ghost_y_shift = 0.5 verts.append( pf.Vert( x=ghost_x_shift + ghost_scale * (-0.5), y=ghost_y_shift + ghost_scale * (-0.6), ) ) # 15 verts.append( pf.Vert( x=ghost_x_shift + ghost_scale * (0.5), y=ghost_y_shift + ghost_scale * (-0.6), ) ) # 16 verts.append( pf.Vert( x=ghost_x_shift + ghost_scale * (0.5), y=ghost_y_shift + ghost_scale * (0.2), ) ) # 17 verts.append( pf.Vert( x=ghost_x_shift + ghost_scale * (-0.5), y=ghost_y_shift + ghost_scale * (0.2), ) ) # 18 verts.append( pf.Vert( x=ghost_x_shift + ghost_scale * (-0.25), y=ghost_y_shift + ghost_scale * (0.1), ) ) # 19 verts.append( pf.Vert( x=ghost_x_shift + ghost_scale * (0.25), y=ghost_y_shift + ghost_scale * (0.1), ) ) # 20

We need to label our vertices:

In [3]:

# TODO: future versions should do this automatically.
for k in range(len(verts)):
    verts[k].set_idx(k)

# TODO: future versions should do this automatically. for k in range(len(verts)): verts[k].set_idx(k)

Let's visualized these points:

In [4]:

plt.figure()
for v in verts:
    plt.plot(v.x, v.y, "ko")
plt.axis("equal")
plt.grid("on")
plt.show()

plt.figure() for v in verts: plt.plot(v.x, v.y, "ko") plt.axis("equal") plt.grid("on") plt.show()

Define the edges¶

Next, we define the edges. The parameter pos_cell_idx is the index of the cell with this edge oriented counterclockwise on the outer boundary (or clockwise if on a hole boundary). In the event that this edge lies on the boundary of the domain and there is no such cell, pos_cell_idx = -1 is taken as the default argument. The neg_cell_idx is the index of the cell where the opposite is true.

In [5]:

edges: list[pf.Edge] = []

# rectangles
edges.append(pf.Edge(verts[0], verts[1], pos_cell_idx=0))
edges.append(pf.Edge(verts[1], verts[2], pos_cell_idx=3))
edges.append(pf.Edge(verts[2], verts[3], pos_cell_idx=7))
edges.append(pf.Edge(verts[3], verts[4], pos_cell_idx=7))
edges.append(pf.Edge(verts[4], verts[5], pos_cell_idx=7))
edges.append(pf.Edge(verts[5], verts[6], pos_cell_idx=3))
edges.append(pf.Edge(verts[6], verts[7], pos_cell_idx=0))
edges.append(pf.Edge(verts[7], verts[0], pos_cell_idx=0))
edges.append(pf.Edge(verts[1], verts[6], pos_cell_idx=0, neg_cell_idx=3))
edges.append(pf.Edge(verts[2], verts[5], pos_cell_idx=3, neg_cell_idx=7))

# pacman
edges.append(pf.Edge(verts[8], verts[9], pos_cell_idx=1, neg_cell_idx=0))
edges.append(
    pf.Edge(
        verts[9],
        verts[10],
        pos_cell_idx=1,
        neg_cell_idx=0,
        curve_type="circular_arc_deg",
        theta0=300,
    )
)
edges.append(pf.Edge(verts[10], verts[8], pos_cell_idx=1, neg_cell_idx=0))
edges.append(
    pf.Edge(
        verts[11],
        verts[11],
        pos_cell_idx=2,
        neg_cell_idx=1,
        curve_type="circle",
        radius=0.25 * pacman_scale,
    )
)

# dots
edges.append(
    pf.Edge(
        verts[12],
        verts[12],
        pos_cell_idx=4,
        neg_cell_idx=3,
        curve_type="circle",
        radius=0.1,
    )
)
edges.append(
    pf.Edge(
        verts[13],
        verts[13],
        pos_cell_idx=5,
        neg_cell_idx=3,
        curve_type="circle",
        radius=0.1,
    )
)
edges.append(
    pf.Edge(
        verts[14],
        verts[14],
        pos_cell_idx=6,
        neg_cell_idx=3,
        curve_type="circle",
        radius=0.1,
    )
)

# ghost
edges.append(
    pf.Edge(
        verts[15],
        verts[16],
        pos_cell_idx=8,
        neg_cell_idx=7,
        curve_type="sine_wave",
        amp=0.1,
        freq=6,
    )
)
edges.append(pf.Edge(verts[16], verts[17], pos_cell_idx=8, neg_cell_idx=7))
edges.append(
    pf.Edge(
        verts[17],
        verts[18],
        pos_cell_idx=8,
        neg_cell_idx=7,
        curve_type="circular_arc_deg",
        theta0=180,
    )
)
edges.append(pf.Edge(verts[18], verts[15], pos_cell_idx=8, neg_cell_idx=7))
edges.append(
    pf.Edge(
        verts[19],
        verts[19],
        pos_cell_idx=9,
        neg_cell_idx=8,
        curve_type="ellipse",
        a=0.15 * ghost_scale,
        b=0.2 * ghost_scale,
    )
)
edges.append(
    pf.Edge(
        verts[20],
        verts[20],
        pos_cell_idx=10,
        neg_cell_idx=8,
        curve_type="ellipse",
        a=0.15 * ghost_scale,
        b=0.2 * ghost_scale,
    )
)

edges: list[pf.Edge] = [] # rectangles edges.append(pf.Edge(verts[0], verts[1], pos_cell_idx=0)) edges.append(pf.Edge(verts[1], verts[2], pos_cell_idx=3)) edges.append(pf.Edge(verts[2], verts[3], pos_cell_idx=7)) edges.append(pf.Edge(verts[3], verts[4], pos_cell_idx=7)) edges.append(pf.Edge(verts[4], verts[5], pos_cell_idx=7)) edges.append(pf.Edge(verts[5], verts[6], pos_cell_idx=3)) edges.append(pf.Edge(verts[6], verts[7], pos_cell_idx=0)) edges.append(pf.Edge(verts[7], verts[0], pos_cell_idx=0)) edges.append(pf.Edge(verts[1], verts[6], pos_cell_idx=0, neg_cell_idx=3)) edges.append(pf.Edge(verts[2], verts[5], pos_cell_idx=3, neg_cell_idx=7)) # pacman edges.append(pf.Edge(verts[8], verts[9], pos_cell_idx=1, neg_cell_idx=0)) edges.append( pf.Edge( verts[9], verts[10], pos_cell_idx=1, neg_cell_idx=0, curve_type="circular_arc_deg", theta0=300, ) ) edges.append(pf.Edge(verts[10], verts[8], pos_cell_idx=1, neg_cell_idx=0)) edges.append( pf.Edge( verts[11], verts[11], pos_cell_idx=2, neg_cell_idx=1, curve_type="circle", radius=0.25 * pacman_scale, ) ) # dots edges.append( pf.Edge( verts[12], verts[12], pos_cell_idx=4, neg_cell_idx=3, curve_type="circle", radius=0.1, ) ) edges.append( pf.Edge( verts[13], verts[13], pos_cell_idx=5, neg_cell_idx=3, curve_type="circle", radius=0.1, ) ) edges.append( pf.Edge( verts[14], verts[14], pos_cell_idx=6, neg_cell_idx=3, curve_type="circle", radius=0.1, ) ) # ghost edges.append( pf.Edge( verts[15], verts[16], pos_cell_idx=8, neg_cell_idx=7, curve_type="sine_wave", amp=0.1, freq=6, ) ) edges.append(pf.Edge(verts[16], verts[17], pos_cell_idx=8, neg_cell_idx=7)) edges.append( pf.Edge( verts[17], verts[18], pos_cell_idx=8, neg_cell_idx=7, curve_type="circular_arc_deg", theta0=180, ) ) edges.append(pf.Edge(verts[18], verts[15], pos_cell_idx=8, neg_cell_idx=7)) edges.append( pf.Edge( verts[19], verts[19], pos_cell_idx=9, neg_cell_idx=8, curve_type="ellipse", a=0.15 * ghost_scale, b=0.2 * ghost_scale, ) ) edges.append( pf.Edge( verts[20], verts[20], pos_cell_idx=10, neg_cell_idx=8, curve_type="ellipse", a=0.15 * ghost_scale, b=0.2 * ghost_scale, ) )

Create a mesh¶

With all of the edges of the mesh defined, we are prepared to define a planar_mesh object.

In [6]:

T = pf.PlanarMesh(edges)

T = pf.PlanarMesh(edges)

Building planar mesh...
PlanarMesh:
	num_verts: 15
	num_edges: 23
	num_cells: 11

Let's visualize the mesh skeleton, but first we should remember to parameterize the edges.

In [7]:

pf.plot.MeshPlot(T.edges).draw(show_axis=False)

pf.plot.MeshPlot(T.edges).draw(show_axis=False)

Moreover, we can visualize an individual cell of the mesh:

In [8]:

cell_idx = 8
K = T.get_cell(cell_idx)
pf.plot.MeshPlot(K.get_edges()).draw()

cell_idx = 8 K = T.get_cell(cell_idx) pf.plot.MeshPlot(K.get_edges()).draw()

Reference elements¶

This is a planned feature. Check back soon for updates!