GitHub

Geometry and meshes

Important points covered in this tutorial
  • Combine simple shapes to create domains
  • Import a mesh from a file
  • Iterative over mesh elements

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
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}}
	 852 elements of type Inti.LagrangeElement{Inti.ReferenceHyperCube{1}, 2, StaticArraysCore.SVector{3, Float64}}
	 14035 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{3}, 4, 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.ReferenceHyperCube{1}, 2, StaticArraysCore.SVector{…
  Inti.LagrangeElement{Inti.ReferenceSimplex{3}, 4, StaticArraysCore.SVector{3,…

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
Example block output
Mesh visualization

Note that although the mesh may be of high order and/or conforming, the visualization of a mesh is always performed on the underlying first order mesh, and therefore elements may look flat even if the problem is solved on a curved mesh.

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"#120#122"{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 = "Γ₂")
Example block output

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()
Example block output

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),))
Example block output

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
Example block output

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

Mesh quality

The quality of the generated mesh created through meshgen depends heavily on the quality of the underlying parametrization. For surfaces containing a degenerate parametrization, or for complex shapes, one is better off using a suitable CAD (Computer-Aided Design) software in conjunction with a mesh generator.

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)
Example block output

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)
Example block output
Limitations

At present only the transfinite interpolation for the logically quadrilateral domains is supported. In the future we hope to add support for three-dimensional transfinite interpolation, as well as transfinite formulas for simplices.

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)
Example block output

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.

Type-stable iteration over elements

Since a mesh in Inti.jl can contain elements of various types, the elements function above is not type-stable. For a type-stable iterator approach, one should first iterate over the element types using element_types, and then use elements(msh, E) to iterate over a specific element type E.

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.18054505519780645
 -0.14630186206938434
  0.09955887205829847

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