Tutorial 1.1: Vertices and Edges¶
This tutorial introduces the notion of vertices and edges, which will be used to construct curvilinear meshes in the next tutorial.
We begin by importing the puncturedfem package,
as well as numpy and matplotlib for the sake of this example.
import puncturedfem as pf
import numpy as np
import matplotlib.pyplot as plt
Creating an edge¶
The simplest type of edge is a straight line segment, which is the default
when initializing an Edge object.
# define vertices
v1 = pf.Vert(x=0.0, y=0.0)
v2 = pf.Vert(x=1.0, y=2.0)
# define a straight edge
e1 = pf.Edge(v1, v2)
For curvilinear edges, we can supply additional arguments to define the edge.
The curvetype string is the name of a module located in
puncturedfem/mesh/edgelib
where the functions defining the parameterization of this curve are located.
We can also pass in any keyword arguments used to define the edge.
For instance, we can create a circular arc corresponding to a $120^\circ$ angle as follows:
# create a circular arc
e2 = pf.Edge(v1, v2, curve_type="circular_arc_deg", theta0=120)
Some of the common curvetypes used in these examples are
curvetype |
keyword argument(s) |
|---|---|
'circle' |
'radius' |
'circular_arc' |
'theta0' |
'ellipse' |
'a', 'b' |
'line' |
|
'sine_wave' |
'amp', 'freq' |
To define a custom curvetype, see the appendix at the end of this notebook.
Visualizing Edges¶
We can plot the edges using the MeshPlot class:
pf.plot.MeshPlot([e1, e2]).draw()
We can visualize the orientation of each edge by setting the
show_orientation keyword argument to True.
We can also introduce grid lines by setting the show_grid keyword argument
to True.
pf.plot.MeshPlot([e1, e2]).draw(show_orientation=True, show_grid=True)
Custom parameterizations¶
To create the points $0=t_0 < t_1 < \cdots < t_{2n}=2\pi$ where $x(t)$
will be sampled, we will create a QuadDict object using the get_quad_dict() function.
The QuadDict object is a dictionary containing Quad objects, which are used to sample the curve parameterization.
quad_dict = pf.get_quad_dict(n=32)
print(quad_dict)
{'trap': Quad object
type trap
n 32, 'kress': Quad object
type kress
n 32}
The points for the trapezoidal ("trap") quadrature scheme are,
of course, sampled at equidistant nodes
$t_k = hk$, where $h=\pi / n$ for a chosen natural number $n$.
plt.figure()
plt.plot(quad_dict["trap"].t, "k.")
plt.title("Trapezoid quadrature points")
plt.grid("on")
plt.show()
The Kress ("kress") quadrature should always be used to parameterized edges that terminate at a corner.
Since this is the most common case in practice, it is the default method to parameterize an edge.
We can see that the Kress scheme samples points more heavily near the endpoints:
plt.figure()
plt.plot(quad_dict["kress"].t, "k.")
plt.title("Kress quadrature points")
plt.grid("on")
plt.show()
Creating an Edge with a Cubic Spline¶
The curve_type="spline" can be used to construct an Edge by passing in the keyword argument pts, which is a list of two numpy.ndarrays, one with the $x$-coordinates, the other with the $y$-coordinates.
x = np.array([0.7, 0.3, 0.5, 0.2])
y = np.array([1.0, 0.6, 0.4, 0.4])
anchor = pf.Vert(x[0], y[0])
endpnt = pf.Vert(x[-1], y[-1])
cubic_spline_edge = pf.Edge(
anchor=anchor, endpnt=endpnt, curve_type="spline", pts=[x, y]
)
pf.plot.MeshPlot([cubic_spline_edge]).draw(
show_grid=True, show_orientation=True
)
Splitting an Edge¶
We can use the split_edge() function to split an Edge into two separate Edges.
anchor = pf.Vert(x=1, y=1)
endpnt = pf.Vert(x=3, y=2)
sinusoid_edge = pf.Edge(
anchor=anchor, endpnt=endpnt, curve_type="sine_wave", amp=0.2, freq=7
)
pf.plot.MeshPlot([sinusoid_edge]).draw(show_orientation=True, show_grid=True)
We provide the Edge object we wish to split, and t_split, the value of $t$ where we wish to split the edge parameterized by $x(t)$.
Curves defined in puncturedfem's edge library are by default defined from $0$ to $2\pi$.
The default value of t_split is $\pi$.
e1, e2 = pf.split_edge(sinusoid_edge, t_split=np.pi / 2)
Let's take a look at our new edges:
pf.plot.MeshPlot([e1, e2]).draw(show_orientation=True, show_grid=True)
pf.plot.MeshPlot([e1]).draw(show_orientation=True, show_grid=True)
pf.plot.MeshPlot([e2]).draw(show_orientation=True, show_grid=True)
Defining a custom curvetype¶
(!) Warning: This method of defining custom edges will be deprecated in a future release.
An edge $e$ is taken to be a $C^2$ smooth curve in $\mathbb{R}^2$ parameterized by $x(t)$ for $0\leq t\leq 2\pi$. We refer to $x(0)$ as the anchor point and $x(2\pi)$ as the terminal point, and $x(0),x(2\pi)$ are referred to collectively as the endpoints. We make the following assumptions:
- The edge $e$ is nontrivial: $e$ is not a single point.
- The edge $e$ is nonselfintersecting: $x(t_1)\neq x(t_2)$ for all $0<t_1<t_2<2\pi$.
- $x(\cdot)$ is regularizable: there is some fixed $\sigma>0$ such that $|x'(t)|\geq\sigma$ for all $0 < t < 2\pi$.
In the event that we need an edge that is not provided in the
puncturedfem/mesh/edgelib folder,
we can add to the edge library as follows.
- Create a file
puncturedfem/mesh/edgelib/mycurve.py, wheremycurvewill be the name of the curve that will be called during the initialization of the edge object. - Import the
numpypackage. - In
mycurve.py, define three functions calledX(),DX(), andDDX(). These will define $x(t)$, $x'(t)$, and $x''(t)$ respectively. - Each of these three functions will return a $2\times (2n+1)$ array,
where $2n+1$ is the number of sampled points specified by the chosen
Quadobject. - Row 0 of each array contains the $x_1$ component, and row 1 contains the $x_2$ component.
- Unpack any additional arguments from
**kwargs.
The contents of mycurve.py will look generically like the following:
"""
A short description of the curve.
A description of any parameters that are used.
"""
import numpy as np
def X(t, **kwargs):
my_parameter = kwargs["my_parameter"]
x = np.zeros((2,len(t)))
x[0,:] = # the x_1 component
x[1,:] = # the x_2 component
return x
def DX(t, **kwargs):
my_parameter = kwargs["my_parameter"]
dx = np.zeros((2,len(t)))
dx[0,:] = # the derivative of the x_1 component wrt t
dx[1,:] = # the derivative of the x_2 component wrt t
return dx
def DDX(t, **kwargs):
my_parameter = kwargs["my_parameter"]
ddx = np.zeros((2,len(t)))
ddx[0,:] = # the second derivative of the x_1 component wrt t
ddx[1,:] = # the second derivative of the x_2 component wrt t
return ddx