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.

using Meshes

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

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.

A = Point(0, 0) # point with integer coordinates
B = Point(1, 0) # another point in 2D space
C = Point(0.0, 1.0) # double precision
D = Point2(0, 1) # double precision from Int literal
E = Point(1, 2, 3) # a point in 3D space
F = Point3(1, 2, 3) # another point now with double precision
G = Point(1f0, 2f0, 3f0) # single precision
H = Point3f(1, 2, 3) # single precision from Int literal

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

Points can be subtracted to produce a vector:

B - A
2-element StaticArrays.SVector{2, Int64} with indices SOneTo(2):
 1
 0

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

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

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

G + Vec3f(1,1,1)
Point(2.0f0, 3.0f0, 4.0f0)

And finally, we can create points at random with:

ps = rand(Point2, 10)
10-element Vector{Point2}:
 Point(0.44810102760173476, 0.5555812005222929)
 Point(0.9410553298227295, 0.4123403433805599)
 Point(0.5087360257595488, 0.4462164757267175)
 Point(0.3477644048225135, 0.9987506574463716)
 Point(0.33548180273612993, 0.13057474549450587)
 Point(0.22244492290271567, 0.848779602420022)
 Point(0.009038419643565998, 0.139751929530054)
 Point(0.7978372957880202, 0.4081305637376098)
 Point(0.0609390665062961, 0.03146352962625376)
 Point(0.563049514290691, 0.7394020888833555)

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), (1,1))
Box{2, Int64}(Point(0, 0), Point(1, 1))
s = Sphere((0,0), 1)
Sphere{2, Int64}(Point(0, 0), 1)

The parameters of these primitive geometries can be queried easily:

extrema(b)
(Point(0, 0), Point(1, 1))
center(s), radius(s)
(Point(0, 0), 1)

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

measure(b) == 1
true

We can sample random points on primitives using different methods:

vs = sample(s, RegularSampling(10)) # 10 points over the sphere
IterTools.IVec{Base.Generator{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, Meshes.var"#120#122"{Meshes.var"#r⃗#121"{DataType, Int64}, Point{2, Int64}}}}(Base.Generator{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, Meshes.var"#120#122"{Meshes.var"#r⃗#121"{DataType, Int64}, Point{2, Int64}}}(Meshes.var"#120#122"{Meshes.var"#r⃗#121"{DataType, Int64}, Point{2, Int64}}(Meshes.var"#r⃗#121"{DataType, Int64}(Float64, 1), Point(0, 0)), 0.0:0.6283185307179586:5.654866776461628))

And collect the generator with:

collect(vs)
10-element Vector{Point2}:
 Point(1.0, 0.0)
 Point(0.8090169943749475, 0.5877852522924731)
 Point(0.30901699437494745, 0.9510565162951535)
 Point(-0.30901699437494734, 0.9510565162951536)
 Point(-0.8090169943749473, 0.5877852522924732)
 Point(-1.0, 1.2246467991473532e-16)
 Point(-0.8090169943749475, -0.587785252292473)
 Point(-0.30901699437494756, -0.9510565162951535)
 Point(0.30901699437494723, -0.9510565162951536)
 Point(0.8090169943749473, -0.5877852522924734)

Polytopes

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

t = Triangle((0,0), (1,0), (0,1))
Triangle{2,Int64}
  └─Point(0, 0)
  └─Point(1, 0)
  └─Point(0, 1)

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(Point2[(0,0), (2,0), (2,2), (1,3), (0,2), (0,0)])
PolyArea{2,Float64}
  outer
    └─6-chain

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

orientation(p)
1-element Vector{Symbol}:
 :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)
100×100 CartesianGrid{2,Float64}
  minimum: Point(0.0, 0.0)
  maximum: Point(100.0, 100.0)
  spacing: (1.0, 1.0)

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.([(3,1,5),(4,6,2)], Triangle)
quads = connect.([(1,2,5,6),(5,6,3,4)], 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 faces
    └─Triangle(3, 1, 5)
    └─Triangle(4, 6, 2)
    └─Quadrangle(1, 2, 5, 6)
    └─Quadrangle(5, 6, 3, 4)

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

collect(elements(mesh))
4-element Vector{Polygon{2, Float64}}:
 Triangle(Point(0.0, 1.0), Point(0.0, 0.0), Point(0.25, 0.5))
 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.25, 0.5), Point(0.75, 0.5))
 Quadrangle(Point(0.25, 0.5), Point(0.75, 0.5), Point(0.0, 1.0), Point(1.0, 1.0))