Non-axisymmetric cylindrical mesh

[MartinEssink]

I want to do some simulations in polar coordinates (or cylindrical to be more precise), which are not axisymmetric. To generate the mesh, I transform a rectangular mesh into a cylinder. However, that leaves me with points at the origin that are defined in multiple elements, but not 'connected'. Is it possible to somehow merge these points?

[joostvanzwieten]

Since bases are generated on a topology, and the topology doesn't know about your coordinate transformation, you have to modify a basis by hand. I'll explain the procedure for a circle using tensorial topologies.

Assume we have a 1D topology $R$ and a 1D topology $Θ$. Using the removedofs parameter of Topology.basis we can create basis for $R Θ$ with basis functions that are nonzero at $r = 0$ removed:

basis0 = (R * Θ).basis(btype, degree, removedofs=[[0], []])

Then we select the single basis function for $R$ that is nonzero at $r = 0$ and concatenate these two:

basis1 = R.basis(btype, degree)[:1]
basis = numpy.concatenate([basis0, basis1])

Below you can find a full example of a Laplace problem in a cylinder using the above defined basis. The chosen problem is a bit boring though, because the exact solution is independent of $θ$.

import functools
import numpy
from nutils import cli, export, function, mesh, solver
from nutils.expression_v2 import Namespace
import scipy.special
import treelog

class jv(function.Custom):
    'Bessel function of the first kind.'

    def __init__(self, v: float, z: function.IntoArray) -> None:
        z = function.Array.cast(z)
        super().__init__(args=(v, z), shape=(), dtype=z.dtype, npointwise=z.ndim)

    evalf = staticmethod(scipy.special.jv)

def main(
        nelems: int = 4, # number of elements per dimension
        btype: str = 'spline', # type of basis
        degree: int = 3, # degree of basis
        rmax: float = 1., # upper bound of topology in r
        zmax: float = 1., # upper bound of topology in z
        mode: int = 1): # mode of exact solution: 1, 2, ...

    ns = Namespace()
    ns.j0 = functools.partial(jv, 0)

    # We create a topology per dimension in cylinder space and define cartesian
    # coordinates.

    ns.rmax = rmax
    ns.zmax = zmax
    R, ns.r = mesh.line(numpy.linspace(0, rmax, nelems + 1), space='R', bnames=('z-axis', 'wall'))
    Θ, ns.θ = mesh.line(numpy.linspace(0, 2 * numpy.pi, nelems + 1), space='Θ', periodic=True)
    Z, ns.z = mesh.line(numpy.linspace(0, zmax, nelems + 1), space='Z', bnames=('bottom', 'top'))
    topo = R * Θ * Z
    ns.x = 'r cos(θ)'
    ns.y = 'r sin(θ)'
    ns.xyz = numpy.stack([ns.x, ns.y, ns.z])
    ns.define_for('xyz', gradient='∇', jacobians=('dV', 'dS'))

    # We create a basis in the R-Θ plane with a single dof at the r = 0 axis.
    # The basis is a concatenation of a regular basis for R-Θ, with all basis
    # functions that are nonzero at r = 0 removed, together with the single
    # basis function fo R that is nonzero at r = 0.

    RΘbasis = numpy.concatenate([
        (R * Θ).basis(btype, degree, removedofs=[[0], []]),
        R.basis(btype, degree)[:1],
    ])

    # We create a regular basis for Z and define fields u and v.

    Zbasis = Z.basis(btype, degree)
    ns.u = function.dotarg('u', RΘbasis, Zbasis)
    ns.v = function.dotarg('v', RΘbasis, Zbasis)

    # Weak form of the Laplace problem.

    res = topo.integral('∇_i(v) ∇_i(u) dV' @ ns, degree = 2 * degree)

    # Exact solution.

    ns.α = scipy.special.jn_zeros(0, mode)[-1] / rmax
    ns.uexact = 'j0(α r) sinh(α z) / sinh(α zmax)'

    # We constrain u to the exact solution at the boundary of the domain. Note
    # that we explicitly select boundary groups, excluding `z-axis`.

    sqr = topo.boundary['wall,top,bottom'].integral('(u - uexact)^2 dS' @ ns, degree = 2 * degree)
    cons = solver.optimize(('u',), sqr, droptol=1e-10)

    # Solve and plot.

    args = solver.solve_linear(('u',), (res.derivative('v'),), constrain = cons)

    err = numpy.sqrt(topo.integral('(u - uexact)^2 dV' @ ns, degree = 2 * degree)).eval(**args)
    treelog.info(f'L2 err: {err:.3e}')

    smpl = (R * Z).sample('bezier', 5) * Θ[0:].locate(ns.θ, [0], tol = 1e-10)
    r, z, u = smpl.eval(['r', 'z', 'u'] @ ns, **args)
    export.triplot('u_θ=0.png', numpy.stack([r, z], axis = 1), u, tri = smpl.tri, hull = smpl.hull)

if __name__ == '__main__':
    cli.run(main)

[MartinEssink]

It works like a charm, but I am still having one issue. When I try to sample the full mesh to export to a VTK file, I encounter a NotImplementedError. It seems either triangulation a 3d space is unsupported, or sampling along the three dimensions is not supported in this case only.

smpl2 = (R * Z * Θ).sample('bezier', 5)
print(smpl2.tri)

[joostvanzwieten]

This should work:

smpl2 = (R * (Θ * Z)).sample('bezier', 5)

I've opened #711 to track this problem.