Constrained Delaunay Triangulation Meshes for Spatial 'SPDE' Models

Generate constrained Delaunay triangulation meshes for use with stochastic partial differential equation (SPDE) spatial models (Lindgren, Rue and Lindstroem 2011 ). Provides automatic mesh generation from point coordinates with boundary constraints, Ruppert refinement for mesh quality, finite element method (FEM) matrix assembly (mass, stiffness, projection), barrier models, spherical meshes via icosahedral subdivision, and metric graph meshes for network geometries. Built on the 'CDT' header-only C++ library (Amirkhanov 2024 < https://github.com/artem-ogre/CDT>). Designed as the mesh backend for the 'tulpa' Bayesian hierarchical modelling engine but usable standalone for any spatial triangulation task.

tulpaMesh

a mesh from scattered points

tulpaMesh builds triangular meshes from point coordinates, such as survey sites or sampling stations, and assembles the finite-element matrices that spatial models based on stochastic partial differential equations (SPDE) need.

The package is the mesh backend for the tulpa Bayesian modelling engine and runs standalone for any spatial triangulation task. It triangulates the points with the constrained Delaunay triangulation of the CDT C++ library (Amirkhanov 2024), which uses exact geometric predicates, and can refine the triangles until every angle exceeds a minimum. From the mesh it assembles the mass, stiffness and projection matrices of the SPDE approach (Lindgren, Rue and Lindström 2011) as sparse matrices. You set the mesh resolution, check the triangle quality and pass the matrices to the spatial model.

tulpaMesh supports:

• Meshes from points or boundaries. tulpa_mesh() takes a coordinate matrix, a formula with a data frame, or a boundary polygon from sf, with holes. A maximum edge length sets the density, and Ruppert refinement enforces a minimum angle or a maximum triangle area.
• Finite-element matrices. fem_matrices() returns the mass, stiffness and projection matrices, with a lumped mass option and parallel assembly. Quadratic elements and spatially varying range and variance have their own assembly functions.
• Barrier models. A field that must not cross a coastline or another region (Bakka et al. 2019) gets zero stiffness across the barrier triangles.
• Spheres, time and networks. Icosahedral meshes for global models, one-dimensional meshes for time, and metric-graph meshes for rivers and roads use the same matrix assembly.
• Conversion from fmesher. An existing mesh from fmesher converts with as_tulpa_mesh() and gives the same matrices.

The automated tests compare the mass and stiffness matrices with those fmesher assembles on the same mesh and require them to match.

Installation

Install the released version from CRAN:

install.packages("tulpaMesh")

CRAN provides binaries for Windows and macOS, which need no build tools.

The development version compiles C++17 code. Install the build tools for your operating system first:

Check that R can find the compiler:

pkgbuild::has_build_tools(debug = TRUE)

Then install tulpaMesh from GitHub:

# install.packages("pak")
pak::pak("gcol33/tulpaMesh")

The build needs GNU make. The triangulation library is bundled in the package sources and the threading library comes with RcppParallel, so no other system library is needed.

Optional: parallel assembly

parallel = TRUE in fem_matrices() assembles the matrices on several CPU threads through RcppParallel’s thread pool. No compiler setup is needed: RcppParallel ships its threading library, and on platforms where that library is unavailable it falls back to its own portable backend. The result is the same at any thread count.

Usage

1. Build a mesh from the coordinates with tulpa_mesh().
2. Check the triangles with mesh_summary() or plot().
3. Assemble the matrices with fem_matrices().

The input is a two-column matrix of coordinates in a planar projection, or a formula naming two columns of a data frame. Coordinates in longitude and latitude belong on a spherical mesh (Spherical and Temporal Meshes).

library(tulpaMesh)

coords <- cbind(runif(200), runif(200))
mesh <- tulpa_mesh(coords)
mesh
#> tulpa_mesh:
#>   Vertices:   236 
#>   Triangles:  434 
#>   Edges:      669 
#>   ( 200  input + 36  generated)

Without further arguments the mesh is the Delaunay triangulation of the points, extended by a margin around their hull. max_edge adds vertices on a hexagonal lattice so that no edge exceeds the given length, and min_angle or max_area run Ruppert refinement to the stated quality.

mesh <- tulpa_mesh(coords, max_edge = 0.1)
fem <- fem_matrices(mesh, obs_coords = coords)
dim(fem$A)
#> [1] 200 326

The result contains the mass matrix C, the stiffness matrix G and, when observation coordinates are given, the projection matrix A with one row per observation and one column per mesh vertex. Quick Start continues with mesh density, the formula interface and quality diagnostics.

Boundaries and barriers

A boundary polygon from sf restricts the mesh to a study area, with holes left untriangulated. A barrier model keeps a spatial field from crossing a region such as a coastline: barrier_triangles() identifies the triangles inside the barrier polygons, and the stiffness matrix then carries zero stiffness across them.

library(sf)

barrier <- st_polygon(list(rbind(c(3, 3), c(7, 3), c(7, 7), c(3, 7), c(3, 3))))
bt <- barrier_triangles(mesh, st_sfc(barrier))
fem <- fem_matrices(mesh, barrier = bt)

Spatial Workflows covers custom boundaries, polygons with holes, barrier models, subdivision and adaptive refinement.

Spheres, time and networks

tulpa_mesh_sphere() subdivides an icosahedron into a geodesic mesh for global models, with the matrices assembled on the surface. tulpa_mesh_1d() builds a mesh over time knots for space-time models. tulpa_mesh_graph() builds a mesh over the edges of a river or road network. Spherical and Temporal Meshes works through each, together with non-stationary and quadratic-element assembly.

Mesh quality and operations

mesh_quality() returns the angles, aspect ratio and area of every triangle, and mesh_summary() their minimum, median and maximum. The plot method shades the triangles by their minimum angle with color = "quality". subdivide_mesh() splits every triangle into four, subset_mesh() extracts a submesh, mesh_components() finds disconnected parts, and refine_mesh() refines the triangles with the largest error indicators.

Documentation

Start with Quick Start: one path from a coordinate matrix to the finite-element matrices, with mesh density, the formula interface and quality checks.

• Spatial Workflows: boundaries, polygons with holes, barrier models, subdivision, adaptive refinement and conversion from fmesher.
• Spherical and Temporal Meshes: spheres, time, metric graphs, non-stationary and quadratic elements, parallel assembly.
• Full reference: every function, grouped by job.

Getting help

Questions and bug reports go to the GitHub issues. A bug report is easiest to act on with a small reproducible example, for instance made with reprex.

Support

“Software is like sex: it’s better when it’s free.” — Linus Torvalds

I’m a PhD student who builds R packages in my free time because I believe good tools should be free and open. I started these projects for my own work and figured others might find them useful too.

If this package saved you some time, buying me a coffee is a nice way to say thanks. It helps with my coffee addiction.

Acknowledgements

The constrained Delaunay triangulation comes from the CDT header-only library by Artem Amirkhanov, with the exact geometric predicates by William C. Lenthe.

License

tulpaMesh is licensed under MIT. Two bundled components carry their own licences:

• The CDT library (src/cdt/) by Artem Amirkhanov, under MPL-2.0.
• The geometric predicates (src/cdt/predicates.h) by William C. Lenthe, under the 3-clause BSD licence.

Citation

Colling G (2026). _tulpaMesh: Constrained Delaunay Triangulation Meshes
for Spatial SPDE Models_. doi:10.32614/CRAN.package.tulpaMesh
<https://doi.org/10.32614/CRAN.package.tulpaMesh>.
<https://CRAN.R-project.org/package=tulpaMesh>.

BibTeX:

@Manual{,
  title = {tulpaMesh: Constrained Delaunay Triangulation Meshes for Spatial SPDE Models},
  author = {Gilles Colling},
  year = {2026},
  url = {https://CRAN.R-project.org/package=tulpaMesh},
  doi = {10.32614/CRAN.package.tulpaMesh},
}

Download the citation as BibTeX or RIS.