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
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[]
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"#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₂"]
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.
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:
- Only a single boundary entity can be associated with a given curved volume entity
- 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:
- The boundary parametrization must be global. Thus, a torus domain is possible but not a sphere.
- 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.
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.41172958308427027
0.39152798515716003
0.13290053704846747
Likewise, we can compute the jacobian
of the element, or its normal
at a given parametric coordinate.