Title:	Constrained Delaunay Triangulation Meshes for Spatial 'SPDE' Models
Version:	0.1.1
Description:	Generate constrained Delaunay triangulation meshes for use with stochastic partial differential equation (SPDE) spatial models (Lindgren, Rue and Lindstroem 2011 <doi:10.1111/j.1467-9868.2011.00777.x>). 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.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
SystemRequirements:	C++17, GNU make
Depends:	R (≥ 4.1.0)
Imports:	Rcpp (≥ 1.0.12), Matrix (≥ 1.5-0), RcppParallel
Suggests:	testthat (≥ 3.0.0), sf, fmesher, ggplot2, svglite, knitr, rmarkdown
VignetteBuilder:	knitr
LinkingTo:	Rcpp, RcppParallel
Config/testthat/edition:	3
URL:	https://github.com/gcol33/tulpaMesh
BugReports:	https://github.com/gcol33/tulpaMesh/issues
NeedsCompilation:	yes
Packaged:	2026-03-30 19:16:03 UTC; Gilles Colling
Author:	Gilles Colling [aut, cre, cph], Artem Amirkhanov [ctb, cph] (CDT library (src/cdt/, MPL-2.0)), William C. Lenthe [ctb, cph] (Geometric predicates (src/cdt/predicates.h, BSD-3-Clause))
Maintainer:	Gilles Colling <gilles.colling051@gmail.com>
Repository:	CRAN
Date/Publication:	2026-04-03 08:20:02 UTC

tulpaMesh: Constrained Delaunay Triangulation Meshes for Spatial 'SPDE' Models

Description

Generate constrained Delaunay triangulation meshes for use with stochastic partial differential equation (SPDE) spatial models (Lindgren, Rue and Lindstroem 2011 doi:10.1111/j.1467-9868.2011.00777.x). 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.

Author(s)

Maintainer: Gilles Colling gilles.colling051@gmail.com (ORCID) [copyright holder]

Other contributors:

• Artem Amirkhanov (CDT library (src/cdt/, MPL-2.0)) [contributor, copyright holder]

• William C. Lenthe (Geometric predicates (src/cdt/predicates.h, BSD-3-Clause)) [contributor, copyright holder]

See Also

Useful links:

• https://github.com/gcol33/tulpaMesh

• Report bugs at https://github.com/gcol33/tulpaMesh/issues

Convert to a tulpa_mesh Object

Description

Generic function to convert mesh objects from other packages into tulpa_mesh objects. Currently supports fm_mesh_2d objects from the fmesher package.

Usage

as_tulpa_mesh(x, ...)

## S3 method for class 'fm_mesh_2d'
as_tulpa_mesh(x, ...)

## S3 method for class 'inla.mesh'
as_tulpa_mesh(x, ...)

Arguments

Value

A tulpa_mesh object.

Identify Barrier Triangles

Description

Determines which mesh triangles fall inside barrier regions (e.g., coastlines, rivers, lakes). A triangle is marked as a barrier if its centroid falls inside any barrier polygon.

Usage

barrier_triangles(mesh, barriers)

Arguments

Value

A logical vector of length n_triangles. TRUE indicates the triangle is inside a barrier region.

Compute FEM Matrices from a Mesh

Description

Computes the finite element mass (C) and stiffness (G) matrices from a triangular mesh, plus the projection matrix (A) that maps mesh vertices to observation locations.

Usage

fem_matrices(
  mesh,
  obs_coords = NULL,
  barrier = NULL,
  parallel = FALSE,
  lumped = FALSE
)

Arguments

Value

A list with sparse matrices (dgCMatrix class from Matrix package):

• C: consistent mass matrix (n_vertices x n_vertices)

• G: stiffness matrix (n_vertices x n_vertices)

• A: projection matrix (n_obs x n_vertices)

• n_mesh: number of mesh vertices

• C0: (only if lumped = TRUE) diagonal lumped mass matrix

• va: (only if lumped = TRUE) vertex areas (numeric vector)

• ta: (only if lumped = TRUE) triangle areas (numeric vector)

Non-stationary FEM Matrices with Spatially Varying Parameters

Description

Computes weighted FEM matrices for non-stationary SPDE models where the range and variance parameters vary spatially. The weights are per-vertex kappa(s) and tau(s) values, interpolated within each triangle using the element-average approximation.

Usage

fem_matrices_nonstationary(mesh, kappa, tau)

Arguments

Value

A list with sparse matrices:

• Ck: kappa²-weighted mass matrix

• Gk: kappa²-weighted stiffness matrix

• Ct: tau²-weighted mass matrix

• C: unweighted consistent mass matrix

• G: unweighted stiffness matrix

• C0: lumped (diagonal) mass matrix

• n_mesh: number of mesh vertices

Compute P2 (Quadratic) FEM Matrices

Description

Computes finite element matrices using 6-node quadratic triangular elements. Each triangle edge gets a midpoint node, giving 6 basis functions per triangle instead of 3. Quadratic elements provide better approximation accuracy with fewer mesh nodes.

Usage

fem_matrices_p2(mesh)

Arguments

Value

A list with:

• C: consistent mass matrix (n_total x n_total sparse)

• G: stiffness matrix (n_total x n_total sparse)

• n_mesh: total number of nodes (vertices + midpoints)

• n_vertices: number of original mesh vertices

• n_midpoints: number of added midpoint nodes

• vertices: N x 2 matrix of all node coordinates

• triangles6: M x 6 connectivity matrix (columns: v0, v1, v2, mid01, mid12, mid20)

Identify Disconnected Mesh Components

Description

Finds connected components of a mesh via triangle adjacency (two triangles are connected if they share an edge).

Usage

mesh_components(mesh)

Arguments

Value

An integer vector of length n_triangles giving the component ID (1, 2, ...) for each triangle.

Get or Set the CRS of a Mesh

Description

Access or assign a coordinate reference system to a tulpa_mesh or tulpa_mesh_graph object. The CRS is stored as metadata and propagated through mesh operations.

Usage

mesh_crs(x)

set_crs(x, value)

Arguments

Value

mesh_crs() returns the CRS (an sf::crs object or NULL). set_crs() returns the mesh with CRS attached.

Per-Triangle Mesh Quality Metrics

Description

Computes quality metrics for each triangle in a mesh: minimum angle, maximum angle, aspect ratio, and area.

Computes per-triangle quality metrics: minimum angle, maximum angle, aspect ratio, and area. Useful for identifying degenerate triangles that may cause numerical issues in SPDE models.

Usage

mesh_quality(mesh)

mesh_quality(mesh)

Arguments

Value

A data.frame with one row per triangle and columns:

• min_angle: minimum interior angle (degrees)

• max_angle: maximum interior angle (degrees)

• aspect_ratio: ratio of circumradius to twice the inradius (1 for equilateral, larger for worse quality)

• area: triangle area

A data.frame with one row per triangle and columns:

• min_angle: minimum interior angle (degrees)

• max_angle: maximum interior angle (degrees)

• aspect_ratio: longest edge / shortest edge (1 = equilateral)

• area: triangle area

Examples

set.seed(42)
mesh <- tulpa_mesh(cbind(runif(50), runif(50)))
q <- mesh_quality(mesh)
summary(q)

Print Mesh Quality Summary

Description

Prints min/median/max of mesh quality metrics.

Prints a summary of mesh quality metrics.

Usage

mesh_summary(mesh)

mesh_summary(mesh)

Arguments

Value

Invisible data.frame of per-triangle quality metrics.

Invisibly returns the quality data.frame.

Plot a Triangular Mesh

Description

Draws the mesh using base R graphics: edges as line segments, vertices as points. Optionally colors triangles by a quality metric.

Usage

## S3 method for class 'tulpa_mesh'
plot(
  x,
  color = NULL,
  border = "grey50",
  vertex_col = NULL,
  vertex_cex = 0.5,
  palette = grDevices::hcl.colors,
  n_colors = 100L,
  main = "tulpa_mesh",
  ...
)

Arguments

Value

The tulpa_mesh object x, returned invisibly. Called for the side effect of producing a plot.

Adaptively Refine a Mesh Based on Error Indicators

Description

Refines triangles where error indicators exceed a threshold by inserting their centroids as new vertices and re-triangulating. Designed for iterative solve-refine-re-solve workflows where error indicators come from an SPDE solver (e.g., tulpa's posterior variance).

Usage

refine_mesh(
  mesh,
  indicators,
  threshold = NULL,
  fraction = NULL,
  max_iter = 1L,
  min_area = 0
)

Arguments

Value

A new tulpa_mesh object with refined triangulation.

Subdivide a Mesh

Description

Splits each triangle into 4 sub-triangles by inserting edge midpoints, then re-triangulates. Useful for multi-resolution workflows where a coarser mesh needs uniform refinement.

Usage

subdivide_mesh(mesh)

Arguments

Value

A new tulpa_mesh object with approximately 4x as many triangles.

Extract a Submesh from Triangle Indices

Description

Creates a new tulpa_mesh containing only the specified triangles, with vertex indices remapped.

Usage

subset_mesh(mesh, triangles)

Arguments

Value

A new tulpa_mesh object containing only the selected triangles and their vertices.

Create a Triangular Mesh for SPDE Spatial Models

Description

Generates a constrained Delaunay triangulation from point coordinates, optionally with boundary constraints. The resulting mesh can be used directly with tulpa's SPDE spatial fields.

Usage

tulpa_mesh(
  coords,
  data = NULL,
  boundary = NULL,
  max_edge = NULL,
  cutoff = 0,
  extend = 0.1,
  min_angle = NULL,
  max_area = NULL,
  max_steiner = 10000L
)

## S3 method for class 'tulpa_mesh'
print(x, ...)

Arguments

Value

A tulpa_mesh object with components:

• vertices: N x 2 matrix of vertex coordinates

• triangles: M x 3 integer matrix of vertex indices (1-based)

• edges: K x 2 integer matrix of edge vertex indices (1-based)

• n_vertices: number of mesh vertices

• n_triangles: number of triangles

• n_edges: number of edges

• n_input_points: number of original input points

The tulpa_mesh object x, returned invisibly.

Examples

# Simple mesh from random points
set.seed(42)
coords <- cbind(runif(50), runif(50))
mesh <- tulpa_mesh(coords)
print(mesh)

# Mesh with Ruppert refinement (min angle 25 degrees)
mesh_refined <- tulpa_mesh(coords, min_angle = 25)

Create a 1D Mesh for Temporal SPDE Models

Description

Generates a 1D mesh (interval partition) for temporal components of spatio-temporal SPDE models. Returns 1D FEM matrices (tridiagonal mass and stiffness) that can be combined with 2D spatial FEM matrices via Kronecker products in tulpa's space-time Q-builder.

Usage

tulpa_mesh_1d(knots, boundary = NULL, n_extend = 3L)

## S3 method for class 'tulpa_mesh_1d'
print(x, ...)

Arguments

Value

A tulpa_mesh_1d object with components:

• knots: sorted numeric vector of mesh knot locations

• n: number of knots

• C: consistent mass matrix (tridiagonal, n x n sparse)

• G: stiffness matrix (tridiagonal, n x n sparse)

• C0: lumped (diagonal) mass matrix

The tulpa_mesh_1d object x, returned invisibly.

Create a Metric Graph Mesh

Description

Builds a 1D FEM mesh along the edges of a spatial network (roads, rivers, coastlines). Each edge is discretized into segments with P1 linear elements. Junction nodes where edges meet are shared. Based on the Whittle-Matern SPDE formulation on metric graphs (Bolin, Simas & Wallin, 2024).

Usage

tulpa_mesh_graph(edges, max_edge = NULL, snap_tolerance = 1e-08)

## S3 method for class 'tulpa_mesh_graph'
print(x, ...)

Arguments

Value

A tulpa_mesh_graph object with components:

• vertices: N x 2 matrix of node coordinates

• segments: M x 2 integer matrix of segment connectivity (1-based)

• n_vertices, n_segments: counts

• C: consistent mass matrix (tridiagonal-block sparse)

• G: stiffness matrix (tridiagonal-block sparse)

• C0: lumped (diagonal) mass matrix

• degree: integer vector of node degrees (junction = degree > 2)

The tulpa_mesh_graph object x, returned invisibly.

Create a Triangular Mesh on the Sphere

Description

Generates a geodesic mesh on the unit sphere by recursive subdivision of an icosahedron. Optionally refines around data locations using stereographic projection and local Delaunay triangulation.

Usage

tulpa_mesh_sphere(subdivisions = 3L, coords = NULL, radius = 1)

Arguments

Value

A tulpa_mesh object with components:

• vertices: N x 3 matrix of (x, y, z) Cartesian coordinates on the sphere surface

• triangles: M x 3 integer matrix of vertex indices (1-based)

• edges: K x 2 integer matrix of edge vertex indices (1-based)

• lonlat: N x 2 matrix of (longitude, latitude) in degrees

• n_vertices, n_triangles, n_edges: counts

• n_input_points: number of original input points

• sphere: list with radius and subdivisions