4.3.2. A helical line mesh & differential operators on manifolds
Until now the meshes have always be conforming in the number of dimensions, i.e. either one-dimensional meshes with one-dimensional elements or two-dimensional meshes with two-dimensional elements. However, you can also have a mesh with one-dimensional elements embedded in a two-dimensional or three-dimensional space. The same holds for a mesh consisting of two-dimensional elements embedded in a three-dimensional space. These meshes represent manifolds with a non-vanishing co-dimension. We will now create a line mesh that resembles a helical shape, i.e. has a co-dimension of 2:
from pyoomph import *
from pyoomph.equations.poisson import * # use the pre-defined Poisson equation
class HelicalLineMesh(MeshTemplate):
def __init__(self, N=100, radius=1, length=5, windings=4, domain_name="helix"):
super(HelicalLineMesh, self).__init__()
self.N = N
self.radius = radius
self.length = length
self.windings = windings
self.domain_name = domain_name
def define_geometry(self):
domain = self.new_domain(self.domain_name, 3) # Domain, but with 3d nodes
# function to get the node based on a parameter l from [0:1]
def node_at_parameter(l):
x = self.radius * cos(2 * pi * self.windings * l)
y = self.radius * sin(2 * pi * self.windings * l)
z = self.length * l
return self.add_node_unique(x, y, z)
# loop to generate the elements
for i in range(self.N):
n0 = node_at_parameter(i / self.N) # constructing nodes
n1 = node_at_parameter((i + 0.5) / self.N)
n2 = node_at_parameter((i + 1) / self.N)
domain.add_line_1d_C2(n0, n1, n2) # add a second order line element
if i == 0: # Marking the start boundary:
self.add_nodes_to_boundary("start", [n0])
elif i == self.N - 1: # Marking the end boundary:
self.add_nodes_to_boundary("end", [n2])
Note how we name the domain by default "helix", so that we must add equations to the domain "helix" in the add_equations() method of the problem class to restrict them on this helix. Furthermore, in the new_domain() calls, we add a second argument, 3, which sets the nodal dimension space of this domain to \(3\). The rest works analogous to the previous example, however, this time we create second order line elements instead of first order quadrilateral elements by the call of add_line_1d_C2(). Due to the second order, we must supply a third node so that we in total have a start, a center and an end node of each line element. pyoomph automatically converts first order elements to second order elements, if equations on the space "C2" are defined on this domain. Vice versa, if we only have "C1" fields defined on a domain, pyoomph will simplify all generated second order elements to first order elements.
A potential driver code could read
class MeshTestProblem(Problem):
def define_problem(self):
self.add_mesh(HelicalLineMesh())
eqs = MeshFileOutput()
x, y, z = var(["coordinate_x", "coordinate_y", "coordinate_z"])
source = x ** 2 + 5 * y * z + z
eqs += PoissonEquation(name="u", source=source, space="C2")
eqs += DirichletBC(u=0) @ "start"
eqs += DirichletBC(u=0) @ "end"
self.add_equations(eqs @ "helix")
if __name__ == "__main__":
with MeshTestProblem() as problem:
problem.solve(spatial_adapt=4)
problem.output_at_increased_time()
Fig. 4.13 Poisson equation solved on a helical manifold, where the differential operators have been implicitly replaced by their counterparts acting on manifolds.
At this stage, one might wonder, what the PoissonEquation actually does. It should solve
or likewise in weak formulation
However, what should \(\nabla^2 u\) or \(\nabla u\) mean here? If the helix is unfolded, it is just a one-dimensional function. We can parameterize \(u\) by a single parameter \(s\), i.e. \(u=u(s)\), where \(s\) could be e.g. the arc length along the helix. However, how should we apply a gradient or a Laplacian here? The conventional definition of \(\nabla=(\partial_x,\partial_y,\partial_z)\) cannot work, since we cannot derive \(u\) in all these directions, but only along the helix, i.e. with respect to the parameterization \(s\), i.e. \(\partial_s u\).
Whenever pyoomph notices that we have a co-dimension, all spatial derivatives like grad() or div() will be internally expanded in a different manner. Instead of the conventional gradient or divergence, the corresponding differential operator reasonable for the actual manifold is selected. These are defined as follows: Let \(n\) be the dimension of the embedding space, i.e. the dimension of the nodal coordinates, and \(e\) be the element dimension. The co-dimension is obviously \(k=n-e\). The manifold spanned by the elements can be parameterized (at least locally) by \(e\) parameters \(\xi^\alpha\) for \(\alpha=1,\ldots,e\). The \(n\)-dimensional position vector hence reads \(\vec{R}(\xi^\alpha)\). We construct the tangents, i.e. the covariant basis vectors by
And define the covariant metric tensor \(\mathbf{g}\) by \(g_{\alpha\beta}=\vec{g}_\alpha\cdot\vec{g}_\beta\). The inverse of \(\mathbf{g}\), the contravariant metric tensor, is denoted by the components \(g^{\alpha\beta}\). With that, we define the scalar gradient of a field \(\phi\) on the manifold as
where we sum over all \(\alpha\) and \(\beta\). The result is an \(n\)-dimensional vector in the local tangent space of the manifold, pointing in the direction of the largest increase of \(\phi\) along the surface. The choice of the particular parameterization of the manifold does not influence the result.
To illustrate that, let us consider the case of a 1d manifold embedded in a 3d space, as in our helix. Let \(\xi\) be our only parameter \(\xi^1\). We get the covariant basis vector, which is just a non-normalized tangent to the helix, by
is indeed the local normalized tangent on the helix. The metric tensor is just \(\mathbf{g}=[g_{11}]=[\vec{g}\cdot\vec{g}]=[(\partial_\xi\vec{R})^2]\). Hence, the contravariant metric tensor is just given by the single component \(g^{11}=1/g_{11}\). With that, \(\nabla_S\phi\) reads according to the definition (4.11)
When defining the normalized tangent to the helix as \(\vec{t}=\partial_\xi\vec{R}/\|\partial_\xi\vec{R}\|\), we get
Upon reparameterization with the arc length \(s=\xi\|\partial_\xi\vec{R}\|\), we arrive at
which is indeed the slope of \(\phi\) along the manifold pointing in tangential direction.
The divergence of a vector field \(\vec{\psi}\) defined on a manifold reads similar to (4.11), namely
Finally, when changing the coordinate system by set_coordinate_system(), the corresponding scale factors are also considered in the differential operators on manifolds.
As conclusion, we can just use grad() and div() in our equations. When any equation involving these differential operators is restricted to a manifold, the only reasonable differential operator is selected automatically. This allows to use the same PoissonEquation either in the bulk (with co-dimension 0) or at any manifold with co-dimension \(>0\). In the latter case, the Poisson equation reads
with the Laplace-Beltrami operator \(\nabla_S^2\) or likewise in weak formulation
Warning
Spatial adaptivity does not work yet on domains with a co-dimension, but it will be implemented soon. Hence, we cannot use at SpatialErrorEstimator here yet.
Tip
oomph-lib also has a tutorial on the definition and calculation of surface gradients and surface divergences at https://oomph-lib.github.io/oomph-lib/doc/navier_stokes/surface_theory/html/index.html.