Meshes.jl

Computational geometry and meshing algorithms in Julia.

Build Status Coverage Status Stable Documentation Latest Documentation License File

Overview

Meshes.jl provides efficient implementations of concepts from computational geometry and finite element analysis. It promotes rigorous mathematical definitions of spatial discretizations (a.k.a. meshes) that are adequate for describing general manifolds embedded in $\R^n$, including surfaces described with spherical coordinates, and geometries described with multiple coordinate reference systems. Our ambitious goal is to provide all the features of the CGAL project in pure Julia.

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.

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, MeshViz
import WGLMakie as Mke

Points and vectors

A Point is defined by its coordinates in a global reference system. The type of the coordinates is determined automatically based on the specified literals, or is forced to a specific type using helper constructors (e.g. Point2, Point3, Point2f, Point3f). Integer coordinates are converted to Float64 to fulfill the requirements of most geometric processing algorithms, which would be undefined in a discrete scale.

A vector Vec follows the same pattern. It can be constructed with the generic Vec constructor or with the variants Vec2 and Vec3 for double precision and Vec2f and Vec3f for single precision.

# 2D points
A = Point(0.0, 1.0) # double precision as expected
B = Point(0f0, 1f0) # single precision as expected
C = Point(0, 0) # Integer is converted to Float64 by design
D = Point2(0, 1) # explicitly ask for double precision
E = Point2f(0, 1) # explicitly ask for single precision

# 3D points
F = Point(1.0, 2.0, 3.0) # double precision as expected
G = Point(1f0, 2f0, 3f0) # single precision as expected
H = Point(1, 2, 3) # Integer is converted to Float64 by design
I = Point3(1, 2, 3) # explicitly ask for double precision
J = Point3f(1, 2, 3) # explicitly ask for single precision

for P in (A,B,C,D,E,F,G,H,I,J)
  println("Coordinate type: ", coordtype(P))
  println("Embedding dimension: ", embeddim(P))
end
Coordinate type: Float64
Embedding dimension: 2
Coordinate type: Float32
Embedding dimension: 2
Coordinate type: Float64
Embedding dimension: 2
Coordinate type: Float64
Embedding dimension: 2
Coordinate type: Float32
Embedding dimension: 2
Coordinate type: Float64
Embedding dimension: 3
Coordinate type: Float32
Embedding dimension: 3
Coordinate type: Float64
Embedding dimension: 3
Coordinate type: Float64
Embedding dimension: 3
Coordinate type: Float32
Embedding dimension: 3

Points can be subtracted to produce a vector:

B - A
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 0.0
 0.0

They can't be added, but their coordinates can:

coordinates(F) + coordinates(H)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 2.0
 4.0
 6.0

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

F + Vec(1, 1, 1)
Point(2.0, 3.0, 4.0)

And finally, we can create points at random with:

ps = rand(Point2, 10)
10-element Vector{Point2}:
 Point(0.6863230190736924, 0.04066383813371344)
 Point(0.5204207695014049, 0.5627072171598344)
 Point(0.4954869933911188, 0.23725680183841336)
 Point(0.08332470696643046, 0.6488558041934822)
 Point(0.2620307037258214, 0.40957347350142936)
 Point(0.923629471928636, 0.781765786054097)
 Point(0.3815643518212286, 0.6118632524663106)
 Point(0.8656048215603356, 0.6005199808468455)
 Point(0.9829584818222821, 0.6523303425846145)
 Point(0.4106717228297063, 0.38120867529081126)

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(0.0, 0.0, 0.0), Point(1.0, 1.0, 1.0))
centroid(s), radius(s)
(Point(0.0, 0.0, 0.0), 1.0)

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

measure(b)
1.0

We can sample random points on primitives using different methods:

vs = sample(s, RegularSampling(10)) # 10 points over the sphere
IterTools.IVec{Base.Generator{Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Meshes.var"#175#177"{Meshes.var"#r⃗#176"{Float64, Float64}, Point3}}}(Base.Generator{Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Meshes.var"#175#177"{Meshes.var"#r⃗#176"{Float64, Float64}, Point3}}(Meshes.var"#175#177"{Meshes.var"#r⃗#176"{Float64, Float64}, Point3}(Meshes.var"#r⃗#176"{Float64, Float64}(1.0), Point(0.0, 0.0, 0.0)), Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}((0.28559933214452665:0.28559933214452665:2.8559933214452666, 0.0:0.6283185307179586:5.654866776461628))))

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) == area(t) == 1/2
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{2,Float64}
  └─Point(0.5, 0.0)
  └─Point(1.0, 0.0)

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), (0,0))

viz(p)

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

orientation(p)
:CCW

We can get the outer polygonal chain, and reverse it:

c = chains(p)[1]

reverse(c)
6-Chain{2,Float64}
  └─Point(0.0, 0.0)
  └─Point(0.0, 2.0)
  └─Point(1.0, 3.0)
  └─Point(2.0, 2.0)
  └─Point(2.0, 0.0)
  └─Point(0.0, 0.0)

A closed chain (a.k.a. ring) has circular vertices:

v = vertices(c)
5-element CircularVector(view(::Vector{Point2}, 1:5)):
 Point(0.0, 0.0)
 Point(2.0, 0.0)
 Point(2.0, 2.0)
 Point(1.0, 3.0)
 Point(0.0, 2.0)

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

v[1:10]
10-element CircularVector(::Vector{Point2}):
 Point(0.0, 0.0)
 Point(2.0, 0.0)
 Point(2.0, 2.0)
 Point(1.0, 3.0)
 Point(0.0, 2.0)
 Point(0.0, 0.0)
 Point(2.0, 0.0)
 Point(2.0, 2.0)
 Point(1.0, 3.0)
 Point(0.0, 2.0)

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

angles(c) * 180 / pi
5-element Vector{Float64}:
  -90.0
  -90.0
 -135.0
  -90.0
 -135.0

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

angles(reverse(c)) * 180 / pi
5-element Vector{Float64}:
  90.0
 135.0
  90.0
 135.0
  90.0

In case of closed chains, we can compute inner angles as well:

innerangles(c) * 180 / pi
5-element Vector{Float64}:
  90.0
  90.0
 135.0
  90.0
 135.0

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:

g = CartesianGrid(100, 100)

viz(g, showfacets = true)

No memory is allocated:

@allocated CartesianGrid(10000, 10000, 10000)
0

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

collect(elements(g))
10000-element Vector{Quadrangle{2, Float64, Vector{Point2}}}:
 Quadrangle(Point(0.0, 0.0), Point(1.0, 0.0), Point(1.0, 1.0), Point(0.0, 1.0))
 Quadrangle(Point(1.0, 0.0), Point(2.0, 0.0), Point(2.0, 1.0), Point(1.0, 1.0))
 Quadrangle(Point(2.0, 0.0), Point(3.0, 0.0), Point(3.0, 1.0), Point(2.0, 1.0))
 Quadrangle(Point(3.0, 0.0), Point(4.0, 0.0), Point(4.0, 1.0), Point(3.0, 1.0))
 Quadrangle(Point(4.0, 0.0), Point(5.0, 0.0), Point(5.0, 1.0), Point(4.0, 1.0))
 Quadrangle(Point(5.0, 0.0), Point(6.0, 0.0), Point(6.0, 1.0), Point(5.0, 1.0))
 Quadrangle(Point(6.0, 0.0), Point(7.0, 0.0), Point(7.0, 1.0), Point(6.0, 1.0))
 Quadrangle(Point(7.0, 0.0), Point(8.0, 0.0), Point(8.0, 1.0), Point(7.0, 1.0))
 Quadrangle(Point(8.0, 0.0), Point(9.0, 0.0), Point(9.0, 1.0), Point(8.0, 1.0))
 Quadrangle(Point(9.0, 0.0), Point(10.0, 0.0), Point(10.0, 1.0), Point(9.0, 1.0))
 ⋮
 Quadrangle(Point(91.0, 99.0), Point(92.0, 99.0), Point(92.0, 100.0), Point(91.0, 100.0))
 Quadrangle(Point(92.0, 99.0), Point(93.0, 99.0), Point(93.0, 100.0), Point(92.0, 100.0))
 Quadrangle(Point(93.0, 99.0), Point(94.0, 99.0), Point(94.0, 100.0), Point(93.0, 100.0))
 Quadrangle(Point(94.0, 99.0), Point(95.0, 99.0), Point(95.0, 100.0), Point(94.0, 100.0))
 Quadrangle(Point(95.0, 99.0), Point(96.0, 99.0), Point(96.0, 100.0), Point(95.0, 100.0))
 Quadrangle(Point(96.0, 99.0), Point(97.0, 99.0), Point(97.0, 100.0), Point(96.0, 100.0))
 Quadrangle(Point(97.0, 99.0), Point(98.0, 99.0), Point(98.0, 100.0), Point(97.0, 100.0))
 Quadrangle(Point(98.0, 99.0), Point(99.0, 99.0), Point(99.0, 100.0), Point(98.0, 100.0))
 Quadrangle(Point(99.0, 99.0), Point(100.0, 99.0), Point(100.0, 100.0), Point(99.0, 100.0))

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 = Point2[(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{2,Float64}
  6 vertices
    └─Point(0.0, 0.0)
    └─Point(1.0, 0.0)
    └─Point(0.0, 1.0)
    └─Point(1.0, 1.0)
    └─Point(0.25, 0.5)
    └─Point(0.75, 0.5)
  4 elements
    └─Triangle(1, 5, 3)
    └─Triangle(4, 6, 2)
    └─Quadrangle(1, 2, 6, 5)
    └─Quadrangle(4, 3, 5, 6)
viz(mesh, showfacets = true)

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

collect(elements(mesh))
4-element Vector{Ngon{N, 2, Float64, V} where {N, V<:AbstractVector{Point2}}}:
 Triangle(Point(0.0, 0.0), Point(0.25, 0.5), Point(0.0, 1.0))
 Triangle(Point(1.0, 1.0), Point(0.75, 0.5), Point(1.0, 0.0))
 Quadrangle(Point(0.0, 0.0), Point(1.0, 0.0), Point(0.75, 0.5), Point(0.25, 0.5))
 Quadrangle(Point(1.0, 1.0), Point(0.0, 1.0), Point(0.25, 0.5), Point(0.75, 0.5))

Mesh data

To attach data to the geometries of a mesh, we can use the meshdata function, which combines a mesh object with a collection of Tables.jl tables. For example, it is common to attach a table vtable to the vertices and a table etable to the elements of the mesh:

d = meshdata(mesh,
  vtable = (temperature=rand(6), pressure=rand(6)),
  etable = (quality=["A","B","C","D"], state=[true,false,true,false])
)
4 MeshData{2,Float64}
  variables (rank 0)
    └─pressure (Float64)
    └─temperature (Float64)
  variables (rank 2)
    └─quality (String)
    └─state (Bool)
  domain: 4 SimpleMesh{2,Float64}

More generally, we can attach a table to any rank:

  • 0 (vertices)
  • 1 (segments)
  • 2 (triangles, quadrangles, ...)
  • 3 (tetrahedrons, hexahedrons, ...)

To retrieve the data table for a given rank we use the values function:

values(d, 0)
(temperature = [0.04739994939412351, 0.6047645273313935, 0.0022802414049802877, 0.13658628416817997, 0.24567571111872555, 0.35938745920066584], pressure = [0.15592623141552586, 0.3469999728191109, 0.661461532611197, 0.5553176618826767, 0.9127035477297065, 0.00798599156677171])
values(d, 2)
(quality = ["A", "B", "C", "D"], state = Bool[1, 0, 1, 0])

If we ommit the rank, the function will return the etable of the mesh:

values(d)
(quality = ["A", "B", "C", "D"], state = Bool[1, 0, 1, 0])

When a table is not available for a given rank, the value nothing is returned instead:

values(d, 1) === nothing
true

Finally, we can use the domain function to retrieve the underlying domain of the data, which in this case is a SimpleMesh:

domain(d)
4 SimpleMesh{2,Float64}
  6 vertices
    └─Point(0.0, 0.0)
    └─Point(1.0, 0.0)
    └─Point(0.0, 1.0)
    └─Point(1.0, 1.0)
    └─Point(0.25, 0.5)
    └─Point(0.75, 0.5)
  4 elements
    └─Triangle(1, 5, 3)
    └─Triangle(4, 6, 2)
    └─Quadrangle(1, 2, 6, 5)
    └─Quadrangle(4, 3, 5, 6)