Geometry and meshes

In the getting started tutorial, we saw how to solve a simple Helmholtz scattering problem in 2D. We will now dig deeper into how to create and manipulate more complex geometrical shapes, as well the associated meshes.

Overview

Inti.jl provides a flexible way to define geometrical entities and their associated meshes. Simply put, the GeometricEntity type is the atomic building block of geometries: they can represent points, curves, surfaces, or volumes. Geometrical entities of the same dimension can be combined to form Domain, and domains can be manipulated using basic set operations such union and intersection. Meshes on the other hand are collections of (simple) elements that approximate the geometrical entities. A mesh element is a just a function that maps points from a ReferenceShape to the physical space.

In most applications involving complex three-dimensional surfaces, an external meshing software is used to generate a mesh, and the mesh is imported using the import_mesh function (which relies on Gmsh). The entities can then be extracted from the mesh based on e.g. their dimension or label. Here is an example of how to import a mesh from a file:

using Inti
using Gmsh
using LinearAlgebra
filename = joinpath(Inti.PROJECT_ROOT,"docs", "assets", "piece.msh")
msh = Inti.import_mesh(filename)

Inti.Mesh{3, Float64} containing:
	 126 elements of type StaticArraysCore.SVector{3, Float64}
	 7156 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{2}, 3, StaticArraysCore.SVector{3, Float64}}
	 14035 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{3}, 4, StaticArraysCore.SVector{3, Float64}}
	 852 elements of type Inti.LagrangeElement{Inti.ReferenceHyperCube{1}, 2, StaticArraysCore.SVector{3, Float64}}

The imported mesh contains elements of several types, used to represent the segments, triangles, and tetras used to approximate the geometry:

Inti.element_types(msh)

KeySet for a Dict{DataType, Matrix{Int64}} with 4 entries. Keys:
  StaticArraysCore.SVector{3, Float64}
  Inti.LagrangeElement{Inti.ReferenceSimplex{2}, 3, StaticArraysCore.SVector{3,…
  Inti.LagrangeElement{Inti.ReferenceSimplex{3}, 4, StaticArraysCore.SVector{3,…
  Inti.LagrangeElement{Inti.ReferenceHyperCube{1}, 2, StaticArraysCore.SVector{…

Note that the msh object contains all entities used to construct the mesh, usually defined in a .geo file, which can be extracted using the entities:

ents = Inti.entities(msh)

Filtering of entities satisfying a certain condition, e.g., entities of a given dimension or containing a certain label, can also be performed in order to construct a domain:

filter = e -> Inti.geometric_dimension(e) == 3
Ω = Inti.Domain(filter, ents)

Domain with 4 entities
 EntityKey: (3, 3) => Inti.GeometricEntity with labels String[]
 EntityKey: (3, 6) => Inti.GeometricEntity with labels String[]
 EntityKey: (3, 4) => Inti.GeometricEntity with labels String[]
 EntityKey: (3, 5) => Inti.GeometricEntity with labels String[]

Domains can be used to index the mesh, creating either a new object containing only the necessary elements:

Γ = Inti.boundary(Ω)
msh[Γ]

Inti.Mesh{3, Float64} containing:
	 6892 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{2}, 3, StaticArraysCore.SVector{3, Float64}}

or a SubMesh containing a view of the mesh:

Γ_msh = view(msh, Γ)

Inti.SubMesh{3, Float64} containing:
	 6892 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{2}, 3, StaticArraysCore.SVector{3, Float64}}

Finally, we can visualize the mesh using:

using Meshes, GLMakie
fig = Figure(; size = (800,400))
ax = Axis3(fig[1, 1]; aspect = :data)
viz!(Γ_msh; showsegments = true, alpha = 0.5)
fig

Parametric entities and meshgen

In the previous section we saw an example of how to import a mesh from a file, and how to extract the entities from the mesh. For simple geometries for which an explicit parametrization is available, Inti.jl provides a way to create and manipulate geometrical entities and their associated meshes.

Parametric curves

The simplest parametric shapes are parametric_curves, which are defined by a function that maps a scalar parameter t to a point in 2D or 3D space. Parametric curves are expected to return an SVector, and can be created as follows:

using StaticArrays
l1 = Inti.parametric_curve(x->SVector(x, 0.1 * sin(2π * x)), 0.0, 1.0, labels = ["l₁"])

EntityKey: (1, 199) => Inti.GeometricEntity with labels ["l₁"]

The object l1 represents a GeometricEntity with a known push-forward map:

Inti.pushforward(l1)

(domain = Inti.HyperRectangle{1, Float64}([0.0], [1.0]), parametrization = Inti.var"#168#170"{Main.var"Main".var"#3#4"}(Main.var"Main".var"#3#4"()))

For the sake of this example, let's create three more curves, and group them together to form a Domain:

l2 = Inti.parametric_curve(x->SVector(1 + 0.1 * sin(2π * x), x), 0.0, 1.0, labels = ["l₂"])
l3 = Inti.parametric_curve(x->SVector(1 - x, 1 - 0.1 * sin(2π * x)), 0.0, 1.0, labels = ["l₃"])
l4 = Inti.parametric_curve(x->SVector(0.1 * sin(2π * x), 1 - x), 0.0, 1.0, labels = ["l₄"])
Γ  = l1 ∪ l2 ∪ l3 ∪ l4

Domain with 4 entities
 EntityKey: (1, 201) => Inti.GeometricEntity with labels ["l₃"]
 EntityKey: (1, 202) => Inti.GeometricEntity with labels ["l₄"]
 EntityKey: (1, 199) => Inti.GeometricEntity with labels ["l₁"]
 EntityKey: (1, 200) => Inti.GeometricEntity with labels ["l₂"]

Domains for which a parametric representation is available can be passed to the meshgen function:

msh = Inti.meshgen(Γ; meshsize = 0.05)

We can use the Meshes.viz function to visualize the mesh, and use domains to index the mesh:

Γ₁ = l1 ∪ l3
Γ₂ = l2 ∪ l4
fig, ax, pl = viz(view(msh, Γ₁); segmentsize = 4,  label = "Γ₁")
viz!(view(msh, Γ₂); segmentsize = 4, color = :red, label = "Γ₂")

Note that the orientation of the curve determines the direction of the normal vector. The normal points to the right of the curve when moving in the direction of increasing parameter t:

pts, tangents, normals = Makie.Point2f[], Makie.Vec2f[], Makie.Vec2f[]
for l in [l1, l2, l3, l4]
      push!(pts, l(0.5)) # mid-point of the curve
      push!(tangents, vec(Inti.jacobian(l, 0.5)))
      push!(normals,Inti.normal(l, 0.5))
end
arrows!(pts, tangents, color = :blue, linewidth = 2, linestyle = :dash, lengthscale = 1/4, label = "tangent")
arrows!(pts, normals, color = :black, linewidth = 2, linestyle = :dash, lengthscale = 1/4, label = "normal")
axislegend()

Parametric surfaces

Like parametric curves, parametric surfaces are defined by a function that maps a reference domain $D \subset \mathbb{R}^2$ to a surface in 3D space. They can be constructed using the parametric_surface function:

# a patch of the unit sphere
lc = SVector(-1.0, -1.0)
hc = SVector(1.0, 1.0)
f = (u,v) -> begin
      x = SVector(1.0, u, v)   # a face of the cube
      x ./ sqrt(u^2 + v^2 + 1) # project to the sphere
end
patch = Inti.parametric_surface(f, lc, hc, labels = ["patch1"])
Γ  = Inti.Domain(patch)
msh = Inti.meshgen(Γ; meshsize = 0.1)
viz(msh[Γ]; showsegments = true, figure = (; size = (400,400),))

Since creating parametric surfaces that form a closed volume can be a bit more involved, Inti.jl provide a few helper functions to create simple shapes:

fig = Figure(; size = (600,400))
nshapes = Inti.length(Inti.PREDEFINED_SHAPES)
ncols = 3; nrows = ceil(Int, nshapes/ncols)
for (n,shape) in enumerate(Inti.PREDEFINED_SHAPES)
      Ω = Inti.GeometricEntity(shape) |> Inti.Domain
      Γ = Inti.boundary(Ω)
      msh = Inti.meshgen(Γ; meshsize = 0.1)
      i,j = (n-1) ÷ ncols + 1, (n-1) % ncols + 1
      ax = Axis3(fig[i,j]; aspect = :data, title = shape)
      hidedecorations!(ax)
      viz!(msh; showsegments = true)
end

See GeometricEntity(shape::String) for a list of predefined geometries.

Transfinite domains

It is possible to combine parametric curves/surfaces to form a transfinite domain where the parametrization is inherited from the curves/surfaces that form its boundary. At present, Inti.jl only supports transfinite squares, which are defined by four parametric curves:

l1 = Inti.parametric_curve(x->SVector(x, 0.1 * sin(2π * x)), 0.0, 1.0, labels = ["l₁"])
l2 = Inti.parametric_curve(x->SVector(1 + 0.1 * sin(2π * x), x), 0.0, 1.0, labels = ["l₂"])
l3 = Inti.parametric_curve(x->SVector(1 - x, 1 - 0.1 * sin(2π * x)), 0.0, 1.0, labels = ["l₃"])
l4 = Inti.parametric_curve(x->SVector(0.1 * sin(2π * x), 1 - x), 0.0, 1.0, labels = ["l₄"])
surf = Inti.transfinite_square(l1, l2, l3, l4; labels = ["Ω"])
Ω = Inti.Domain(surf)
msh = Inti.meshgen(Ω; meshsize = 0.05)
viz(msh; showsegments = true)

Note that the msh object contains all entities used to construct Ω, including the boundary segments:

Inti.entities(msh)

KeySet for a Dict{Inti.EntityKey, Dict{DataType, Vector{Int64}}} with 5 entries. Keys:
  EntityKey: (1, 309) => Inti.GeometricEntity with labels ["l₃"]
  EntityKey: (1, 308) => Inti.GeometricEntity with labels ["l₂"]
  EntityKey: (1, 310) => Inti.GeometricEntity with labels ["l₄"]
  EntityKey: (1, 307) => Inti.GeometricEntity with labels ["l₁"]
  EntityKey: (2, 102) => Inti.GeometricEntity with labels ["Ω"]

This allows us to probe the msh object to extract e.g. the boundary mesh:

viz(msh[Inti.boundary(Ω)]; color = :red)

Curving a given mesh

Inti.jl possesses some capability to create curved meshes from a given mesh, which can be useful when the mesh is not conforming to the geometry and the geometry's boundary is available in parametric form. Specifically the methods implement the work of C. Bernardi [1] which provides so-called 'exact' (sometimes called isogeometric) parametrizations of simplicial elements in arbitrary dimension.

The following example first creates a flat triangulation of a disk using splines through Gmsh:

gmsh.initialize()
meshsize = 2π / 32
gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize)
gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize)
# Two kites
f = (s) -> SVector(-1, 0.0) + SVector(cos(2π*s), sin(2π*s))
bnd1 = Inti.gmsh_curve(f, 0, 1; meshsize)
cl = gmsh.model.occ.addCurveLoop([bnd1])
disk = gmsh.model.occ.addPlaneSurface([cl])
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)
msh = Inti.import_mesh(; dim = 2)
gmsh.finalize()
Ω = Inti.Domain(Inti.entities(msh)) do ent
      return Inti.geometric_dimension(ent) == 2
end
viz(msh[Ω], showsegments=true)
Ω_quad = Inti.Quadrature(msh[Ω]; qorder = 10)
area = Inti.integrate(x->1.0, Ω_quad)
println("Error in area computation using P1 mesh: ", abs(area - π))

Info    : Meshing 1D...
Info    : Meshing curve 1 (BSpline)
Info    : Done meshing 1D (Wall 0.00118074s, CPU 0.001182s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0263588s, CPU 0.026356s)
Info    : 221 nodes 343 elements
Error in area computation using P1 mesh: 0.02014764218061682

As can be seen, despite the large quadrature order employed, the approximation error is still significant. To improve the accuracy, we can use the curve_mesh function to create a curved mesh based on the boundary of the domain:

θ = 5 # smoothness order of curved elements
crvmsh = Inti.curve_mesh(msh, f, θ)
Ω_crv_quad = Inti.Quadrature(crvmsh[Ω]; qorder = 10)
area = Inti.integrate(x->1.0, Ω_crv_quad)
println("Error in area computation using curved mesh: ", abs(area - π))

Error in area computation using curved mesh: 1.7763568394002505e-15

Multiple curved domains, subdomains, and curved surfaces

It may be desired to have multiple curved volumes / boundaries. Inti.jl supports this, associating a parametrization with each volumetric entity in a mesh. Note the delicate correspondence between the correct EntityKey and the parametrization in setting entity_parametrization, and note also the limitations listed below.

gmsh.initialize()
meshsize = 0.075
gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize)
gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize)

# Three circles
c1 = gmsh.model.occ.addDisk(0, 0, 0, 1, 1)
c2 = gmsh.model.occ.addDisk(0, 3.0, 0, 1, 1)
c3 = gmsh.model.occ.addDisk(0, 8.0, 0, 2, 2)
gmsh.model.occ.synchronize()

# Add tags for stable identification of the entities
gmsh.model.addPhysicalGroup(2, [c1], -1, "c1")
gmsh.model.addPhysicalGroup(2, [c2], -1, "c2")
gmsh.model.addPhysicalGroup(2, [c3], -1, "c3")

gmsh.model.mesh.generate(2)
msh = Inti.import_mesh(; dim = 2)

Ω = Inti.Domain(Inti.entities(msh)) do ent
    return Inti.geometric_dimension(ent) == 2
end
gmsh.finalize()

Γ = Inti.external_boundary(Ω)
Ωₕ = view(msh, Ω)
Γₕ = view(msh, Γ)

# Three circles
ψ₁ = (t) -> [cos(2 * π * t), sin(2 * π * t)]
ψ₂ = (t) -> [cos(2 * π * t), 3.0 + sin(2 * π * t)]
ψ₃ = (t) -> [2 * cos(2 * π * t), 8.0 + 2 * sin(2 * π * t)]
entity_parametrizations = Dict{Inti.EntityKey,Function}()
for e in Inti.entities(Ω)
    l = Inti.labels(e)
    if "c1" in l
        entity_parametrizations[e] = ψ₁
    elseif "c2" in l
        entity_parametrizations[e] = ψ₂
    elseif "c3" in l
        entity_parametrizations[e] = ψ₃
    end
end

θ = 6 # smoothness order of curved elements
crvmsh = Inti.curve_mesh(msh, entity_parametrizations, θ)

Γₕ = crvmsh[Γ]
Ωₕ = crvmsh[Ω]

qorder = 5
Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder)
Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Ellipse)
Info    : [ 40%] Meshing curve 2 (Ellipse)
Info    : [ 70%] Meshing curve 3 (Ellipse)
Info    : Done meshing 1D (Wall 0.000369461s, CPU 0.000369s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : [ 70%] Meshing surface 3 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.160406s, CPU 0.160391s)
Info    : 4143 nodes 8283 elements

We can verify once again that the correct area of the region is obtained.

area = Inti.integrate(x -> 1, Ωₕ_quad)
println("Error in computing area of three circles: ", abs(area - 6π))

Error in computing area of three circles: 3.197442310920451e-14

One can extract a subcomponent of the curved (volumetric) domain as usual:

Ω_sub = Inti.Domain(e -> "c3" in Inti.labels(e), Inti.entities(Ω))
Ωₕ_sub = crvmsh[Ω_sub]
Ωₕ_sub_quad = Inti.Quadrature(Ωₕ_sub; qorder = qorder)
area = Inti.integrate(x -> 1, Ωₕ_sub_quad)
println("Error in computing area of one (large) circle: ", abs(area - 4π))

Error in computing area of one (large) circle: 1.7763568394002505e-14

The curved mesh also contains surface elements which are, like their volume counterparts, 'exact' (or isogeometric). To demonstrate this we compute the area using Green's theorem:

F = (x) -> [1/2*x[1], 1/2*x[2]]
lineint = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad)
println("Error in computing area using line integral: ", abs(lineint - 6π))

Error in computing area using line integral: 2.1316282072803006e-14

Note the following restrictions that (currently) hold for 2D curved meshes:

1. Only a single boundary entity can be associated with a given curved volume entity
2. A curved boundary entity cannot be associated with multiple volume entities.

Curved 3D meshes with the same interface are also available with the following two (admittedly significant) restrictions:

1. The boundary parametrization must be global. Thus, a torus domain is possible but not a sphere.
2. Only a single curved domain is possible.

(The second item could be easily addressed in Inti.jl if there is user interest; the first is more difficult to address.)

Elements of a mesh

To iterate over the elements of a mesh, use the elements function:

filename = joinpath(Inti.PROJECT_ROOT,"docs", "assets", "piece.msh")
msh = Inti.import_mesh(filename)
ents = Inti.entities(msh)
Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 3, ents)
els = Inti.elements(view(msh, Ω))
centers = map(el -> Inti.center(el), els)
fig = Figure(; size = (800,400))
ax = Axis3(fig[1, 1]; aspect = :data)
scatter!([c[1] for c in centers], [c[2] for c in centers], [c[3] for c in centers], markersize = 5)

This example shows how to extract the centers of the tetrahedral elements in the mesh; and of course we can perform any desired computation on the elements.

Under the hood, each element is simply a functor which maps points x̂ from a ReferenceShape into the physical space:

el = first(els)
x̂ = SVector(1/3,1/3, 1/3)
el(x̂)

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 0.41172958308427027
 0.39152798515716003
 0.13290053704846747

Likewise, we can compute the jacobian of the element, or its normal at a given parametric coordinate.