Meshes.jl

Computational geometry and meshing algorithms in Julia.

Overview

Meshes.jl provides efficient implementations of concepts from computational geometry. It promotes rigorous mathematical definitions of spatial discretizations (a.k.a. meshes) that are adequate for describing general manifolds embedded in $\mathbb{R}^3$, including surfaces described with spherical coordinates, and geometries described with multiple coordinate reference systems.

Unlike other existing efforts in the Julia ecosystem, this project is being carefully designed to facilitate the use of meshes across different scientific domains. We follow a strict set of good software engineering practices, and are quite pedantic in our test suite to make sure that all our implementations are free of bugs in both single and double floating point precision. Additionally, we guarantee type stability.

The design of this project was motivated by various issues encountered with past attempts to represent geometry, which have been originally designed for visualization purposes (e.g. GeometryTypes.jl, GeometryBasics.jl) or specifically for finite element analysis (e.g. JuAFEM.jl, MeshCore.jl). We hope to provide a smoother experience with mesh representations that are adequate for finite finite element analysis, advanced geospatial modeling and visualization, not just one domain.

The Geospatial Data Science with Julia book is a great resource to learn more about Meshes.jl and geospatial data (i.e. tables over meshes):

If you have questions or would like to brainstorm ideas, don't hesitate to start a thread in our zulip channel.

Installation

Get the latest stable release with Julia's package manager:

] add Meshes

Quick example

Although we didn't have time to document the functionality of the package properly, we still would like to illustrate some important features. For more information on available functionality, please consult the Reference guide and the suite of tests in the package.

In all examples we assume the following packages are loaded:

using Meshes
import CairoMakie as Mke

Points and vectors

A Point is defined by its coordinates in a coordinate reference system from CoordRefSystems.jl. By default, a Cartesian coordinates with NoDatum are used.

Integer coordinates are converted to Float64 to fulfill the requirements of most geometric processing algorithms, which would be undefined in a discrete scale:

Point(0.0, 1.0) # double precision as expected

Point with Cartesian{NoDatum} coordinates
├─ x: 0.0 m
└─ y: 1.0 m

Point(0f0, 1f0) # single precision as expected

Point with Cartesian{NoDatum} coordinates
├─ x: 0.0f0 m
└─ y: 1.0f0 m

Point(0, 0) # Integer is converted to Float64 by design

Point with Cartesian{NoDatum} coordinates
├─ x: 0.0 m
└─ y: 0.0 m

Point(1.0, 2.0, 3.0) # double precision as expected

Point with Cartesian{NoDatum} coordinates
├─ x: 1.0 m
├─ y: 2.0 m
└─ z: 3.0 m

Point(1f0, 2f0, 3f0) # single precision as expected

Point with Cartesian{NoDatum} coordinates
├─ x: 1.0f0 m
├─ y: 2.0f0 m
└─ z: 3.0f0 m

Point(1, 2, 3) # Integer is converted to Float64 by design

Point with Cartesian{NoDatum} coordinates
├─ x: 1.0 m
├─ y: 2.0 m
└─ z: 3.0 m

Points can be subtracted to produce a vector:

A = Point(1.0, 1.0)
B = Point(3.0, 3.0)
B - A

Vec(2.0 m, 2.0 m)

They can't be added, but the vectors from the origin to the points can:

to(A) + to(B)

Vec(4.0 m, 4.0 m)

We can add a point to a vector though, and get a new point:

A + Vec(1, 1)

Point with Cartesian{NoDatum} coordinates
├─ x: 2.0 m
└─ y: 2.0 m

Every point and vector has well-defined coordinates:

coords(A)

Cartesian{NoDatum} coordinates
├─ x: 1.0 m
└─ y: 1.0 m

Primitives

Primitive geometries such as Box, Ball, Sphere, Cylinder are those geometries that can be efficiently represented in a computer without discretization. We can construct such geometries using clean syntax:

b = Box((0.0, 0.0, 0.0), (1.0, 1.0, 1.0))

viz(b)

s = Sphere((0.0, 0.0, 0.0), 1.0)

viz(s)

The parameters of these primitive geometries can be queried easily:

extrema(b)

(Point(x: 0.0 m, y: 0.0 m, z: 0.0 m), Point(x: 1.0 m, y: 1.0 m, z: 1.0 m))

centroid(s), radius(s)

(Point(x: 0.0 m, y: 0.0 m, z: 0.0 m), 1.0 m)

As well as their measure (e.g. area, volume) and other geometric properties:

measure(b)

1.0 m^3

We can sample random points on primitives using different methods:

vs = sample(s, RegularSampling(10)) # 10 points over the sphere

Base.Iterators.Flatten{Base.Generator{Tuple{Base.Generator{Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, Meshes.var"#453#457"{Sphere{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, Base.Generator{Tuple{Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, typeof(identity)}}, typeof(identity)}}(Base.Generator{Tuple{Base.Generator{Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, Meshes.var"#453#457"{Sphere{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, Base.Generator{Tuple{Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, typeof(identity)}}, typeof(identity)}(identity, (Base.Generator{Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, Meshes.var"#453#457"{Sphere{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}(Meshes.var"#453#457"{Sphere{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}(Sphere(center: (x: 0.0 m, y: 0.0 m, z: 0.0 m), radius: 1.0 m)), Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}((0.09090909090909091:0.09090909090909091:0.9090909090909091, 0.0:0.1:0.9))), Base.Generator{Tuple{Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, typeof(identity)}(identity, (Point(x: 0.0 m, y: 0.0 m, z: 1.0 m), Point(x: 1.2246467991473532e-16 m, y: 0.0 m, z: -1.0 m))))))

And collect the generator with:

viz(collect(vs))

Polytopes

Polytopes are geometries with "flat" sides. They generalize polygons and polyhedra. Most commonly used polytopes are already defined in the project, including Segment, Ngon (e.g. Triangle, Quadrangle), Tetrahedron, Pyramid and Hexahedron.

t = Triangle((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))

viz(t)

Some of these geometries have additional functionality like the measure (or area):

measure(t)

0.5 m^2

measure(t) == area(t)

true

Or the ability to know whether or not a point is inside:

p = Point(0.5, 0.0)

p ∈ t

true

For line segments, for example, we have robust intersection algorithms:

s1 = Segment((0.0, 0.0), (1.0, 0.0))
s2 = Segment((0.5, 0.0), (2.0, 0.0))

s1 ∩ s2

Segment
├─ Point(x: 0.5 m, y: 0.0 m)
└─ Point(x: 1.0 m, y: 0.0 m)

Polytopes are widely used in GIS software under names such as "LineString" and "Polygon". We provide robust implementations of these concepts, which are formally known as polygonal Chain and PolyArea.

We can compute the orientation of a chain as clockwise or counter-clockwise, can open and close the chain, create bridges between the various inner rings with the outer ring, and other useful functionality:

p = PolyArea((0,0), (2,0), (2,2), (1,3), (0,2))

viz(p)

The orientation of the above polygonal area is counter-clockwise (CCW):

orientation(p)

CCW::OrientationType = 1

We can get the outer ring, and reverse it:

r = rings(p) |> first

reverse(r)

Ring
├─ Point(x: 0.0 m, y: 0.0 m)
├─ Point(x: 0.0 m, y: 2.0 m)
├─ Point(x: 1.0 m, y: 3.0 m)
├─ Point(x: 2.0 m, y: 2.0 m)
└─ Point(x: 2.0 m, y: 0.0 m)

A ring has circular vertices:

v = vertices(r)

5-element CircularVector(::Vector{Point{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}):
 Point(x: 0.0 m, y: 0.0 m)
 Point(x: 2.0 m, y: 0.0 m)
 Point(x: 2.0 m, y: 2.0 m)
 Point(x: 1.0 m, y: 3.0 m)
 Point(x: 0.0 m, y: 2.0 m)

This means that we can index the vertices with indices that go beyond the range 1:nvertices(r). This is very useful for writing algorithms:

v[1:10]

10-element CircularVector(::Vector{Point{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}):
 Point(x: 0.0 m, y: 0.0 m)
 Point(x: 2.0 m, y: 0.0 m)
 Point(x: 2.0 m, y: 2.0 m)
 Point(x: 1.0 m, y: 3.0 m)
 Point(x: 0.0 m, y: 2.0 m)
 Point(x: 0.0 m, y: 0.0 m)
 Point(x: 2.0 m, y: 0.0 m)
 Point(x: 2.0 m, y: 2.0 m)
 Point(x: 1.0 m, y: 3.0 m)
 Point(x: 0.0 m, y: 2.0 m)

We can also compute angles of any given chain, no matter if it is open or closed:

angles(r)

5-element Vector{Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(rad,), NoDims, nothing}}}:
 -1.5707963267948966 rad
 -1.5707963267948966 rad
  -2.356194490192345 rad
 -1.5707963267948966 rad
  -2.356194490192345 rad

The sign of these angles is a function of the orientation:

angles(reverse(r))

5-element Vector{Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(rad,), NoDims, nothing}}}:
 1.5707963267948966 rad
  2.356194490192345 rad
 1.5707963267948966 rad
  2.356194490192345 rad
 1.5707963267948966 rad

In case of rings (i.e. closed chains), we can compute inner angles as well:

innerangles(r)

5-element Vector{Unitful.Quantity{Float64, NoDims, Unitful.FreeUnits{(rad,), NoDims, nothing}}}:
 1.5707963267948966 rad
 1.5707963267948966 rad
  2.356194490192345 rad
 1.5707963267948966 rad
  2.356194490192345 rad

And there is a lot more functionality available like for instance determining whether or not a polygonal area or chain is simple:

issimple(p)

true

Meshes

Efficient (lazy) mesh representations are provided, including CartesianGrid and SimpleMesh, which are specific types of Domain:

Meshes.Domain — Type

Domain{M,CRS}

A domain is an indexable collection of geometries (e.g. mesh) in a given manifold M with point coordinates specified in a coordinate reference system CRS.

source

grid = CartesianGrid(100, 100)

viz(grid, showsegments = true)

No memory is allocated:

@allocated CartesianGrid(10000, 10000, 10000)

but we can still loop over the elements, which are quadrangles in 2D:

collect(grid)

10000-element Vector{Quadrangle{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Quadrangle((x: 0.0 m, y: 0.0 m), ..., (x: 0.0 m, y: 1.0 m))
 Quadrangle((x: 1.0 m, y: 0.0 m), ..., (x: 1.0 m, y: 1.0 m))
 Quadrangle((x: 2.0 m, y: 0.0 m), ..., (x: 2.0 m, y: 1.0 m))
 Quadrangle((x: 3.0 m, y: 0.0 m), ..., (x: 3.0 m, y: 1.0 m))
 Quadrangle((x: 4.0 m, y: 0.0 m), ..., (x: 4.0 m, y: 1.0 m))
 Quadrangle((x: 5.0 m, y: 0.0 m), ..., (x: 5.0 m, y: 1.0 m))
 Quadrangle((x: 6.0 m, y: 0.0 m), ..., (x: 6.0 m, y: 1.0 m))
 Quadrangle((x: 7.0 m, y: 0.0 m), ..., (x: 7.0 m, y: 1.0 m))
 Quadrangle((x: 8.0 m, y: 0.0 m), ..., (x: 8.0 m, y: 1.0 m))
 Quadrangle((x: 9.0 m, y: 0.0 m), ..., (x: 9.0 m, y: 1.0 m))
 ⋮
 Quadrangle((x: 91.0 m, y: 99.0 m), ..., (x: 91.0 m, y: 100.0 m))
 Quadrangle((x: 92.0 m, y: 99.0 m), ..., (x: 92.0 m, y: 100.0 m))
 Quadrangle((x: 93.0 m, y: 99.0 m), ..., (x: 93.0 m, y: 100.0 m))
 Quadrangle((x: 94.0 m, y: 99.0 m), ..., (x: 94.0 m, y: 100.0 m))
 Quadrangle((x: 95.0 m, y: 99.0 m), ..., (x: 95.0 m, y: 100.0 m))
 Quadrangle((x: 96.0 m, y: 99.0 m), ..., (x: 96.0 m, y: 100.0 m))
 Quadrangle((x: 97.0 m, y: 99.0 m), ..., (x: 97.0 m, y: 100.0 m))
 Quadrangle((x: 98.0 m, y: 99.0 m), ..., (x: 98.0 m, y: 100.0 m))
 Quadrangle((x: 99.0 m, y: 99.0 m), ..., (x: 99.0 m, y: 100.0 m))

We can construct a general unstructured mesh with a global vector of points and a collection of Connectivity that store the indices to the global vector of points:

points = [(0,0), (1,0), (0,1), (1,1), (0.25,0.5), (0.75,0.5)]
tris  = connect.([(1,5,3), (4,6,2)], Triangle)
quads = connect.([(1,2,6,5), (4,3,5,6)], Quadrangle)
mesh = SimpleMesh(points, [tris; quads])

4 SimpleMesh
  6 vertices
  ├─ Point(x: 0.0 m, y: 0.0 m)
  ├─ Point(x: 1.0 m, y: 0.0 m)
  ├─ Point(x: 0.0 m, y: 1.0 m)
  ├─ Point(x: 1.0 m, y: 1.0 m)
  ├─ Point(x: 0.25 m, y: 0.5 m)
  └─ Point(x: 0.75 m, y: 0.5 m)
  4 elements
  ├─ Triangle(1, 5, 3)
  ├─ Triangle(4, 6, 2)
  ├─ Quadrangle(1, 2, 6, 5)
  └─ Quadrangle(4, 3, 5, 6)

viz(mesh, showsegments = true)

The actual geometries of the elements are materialized in a lazy fashion like with the Cartesian grid:

collect(mesh)

4-element Vector{Ngon{N, 𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}} where N}:
 Triangle((x: 0.0 m, y: 0.0 m), (x: 0.25 m, y: 0.5 m), (x: 0.0 m, y: 1.0 m))
 Triangle((x: 1.0 m, y: 1.0 m), (x: 0.75 m, y: 0.5 m), (x: 1.0 m, y: 0.0 m))
 Quadrangle((x: 0.0 m, y: 0.0 m), ..., (x: 0.25 m, y: 0.5 m))
 Quadrangle((x: 1.0 m, y: 1.0 m), ..., (x: 0.75 m, y: 0.5 m))