Geometry and meshes
- 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:
852 elements of type Inti.LagrangeElement{Inti.ReferenceHyperCube{1}, 2, StaticArraysCore.SVector{3, Float64}}
126 elements of type StaticArraysCore.SVector{3, Float64}
14035 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{3}, 4, StaticArraysCore.SVector{3, Float64}}
7156 elements of type Inti.LagrangeElement{Inti.ReferenceSimplex{2}, 3, 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:
Inti.LagrangeElement{Inti.ReferenceHyperCube{1}, 2, StaticArraysCore.SVector{…
StaticArraysCore.SVector{3, Float64}
Inti.LagrangeElement{Inti.ReferenceSimplex{3}, 4, StaticArraysCore.SVector{3,…
Inti.LagrangeElement{Inti.ReferenceSimplex{2}, 3, 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[]
Domain
s 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

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_curve
s, 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₂"]
Domain
s 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.
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)

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)

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)

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.
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.