R/mesh.R

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

#' Create a Triangular Mesh for SPDE Spatial Models
#'
#' 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.
#'
#' @param coords A matrix or data.frame with columns x and y, or a formula
#'   like `~ x + y` evaluated in `data`.
#' @param data Optional data.frame for formula evaluation.
#' @param boundary Optional boundary specification: a matrix of boundary
#'   vertex coordinates (N x 2), an sf polygon, or NULL for convex hull.
#' @param max_edge Maximum edge length. A single value or a vector of two
#'   values `c(inner, outer)` where `inner` sets the lattice spacing inside
#'   the study area (the boundary, or the convex hull of `coords`) and
#'   `outer` sets the lattice spacing and boundary segment length in the
#'   extension region. A single value is used for both.
#' @param cutoff Minimum distance between mesh vertices. Points closer than
#'   this are merged. Default 0 (no merging).
#' @param extend Width of the extension region around the convex hull of
#'   `coords`, as a fraction of the domain diameter (the largest distance
#'   between two points of `coords`). Default 0.1 (10% of domain diameter).
#'   The extended boundary is the outward buffer of the convex hull, with
#'   rounded corners. Set to 0 for no extension. Ignored when `boundary` is
#'   supplied.
#' @param min_angle Minimum angle (degrees) for Ruppert refinement. If
#'   specified, Steiner points are inserted at circumcenters of triangles
#'   with angles below this threshold. Theoretical maximum is ~20.7 degrees;
#'   values up to 30 usually work. Default `NULL` (no refinement).
#' @param max_area Maximum triangle area for refinement. Triangles larger
#'   than this are refined regardless of angle quality. Default `NULL`
#'   (no area constraint).
#' @param max_steiner Maximum number of Steiner points to insert during
#'   Ruppert refinement. Default 10000.
#'
#' @return A `tulpa_mesh` object with components:
#'   \itemize{
#'     \item `vertices`: N x 2 matrix of vertex coordinates. Every vertex
#'       belongs to at least one triangle.
#'     \item `triangles`: M x 3 integer matrix of vertex indices (1-based)
#'     \item `edges`: K x 2 integer matrix of edge vertex indices (1-based)
#'     \item `n_vertices`: number of mesh vertices
#'     \item `n_triangles`: number of triangles
#'     \item `n_edges`: number of edges
#'     \item `n_input_points`: number of original input points
#'   }
#'
#' @export
#'
#' @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)
tulpa_mesh <- function(coords, data = NULL, boundary = NULL,
                       max_edge = NULL, cutoff = 0, extend = 0.1,
                       min_angle = NULL, max_area = NULL,
                       max_steiner = 10000L) {
  # Resolve coordinates
  if (inherits(coords, "formula")) {
    if (is.null(data)) stop("data must be provided when coords is a formula")
    vars <- all.vars(coords)
    if (length(vars) != 2) stop("formula must have exactly 2 variables (x and y)")
    coord_mat <- cbind(data[[vars[1]]], data[[vars[2]]])
  } else if (is.data.frame(coords)) {
    coord_mat <- as.matrix(coords[, 1:2])
  } else {
    coord_mat <- as.matrix(coords)
  }

  if (ncol(coord_mat) != 2) stop("coords must have exactly 2 columns")
  storage.mode(coord_mat) <- "double"

  # Extract CRS from sf boundary if available
  input_crs <- NULL
  if (!is.null(boundary) && inherits(boundary, c("sf", "sfc")) &&
      requireNamespace("sf", quietly = TRUE)) {
    input_crs <- sf::st_crs(boundary)
    if (is.na(input_crs)) input_crs <- NULL
  }

  # Build boundary edges if provided. `inner_ring`/`inner_holes` describe the
  # study area that receives `max_edge[1]` lattice points; `outer_ring` is the
  # extended boundary whose ring outside the study area receives `max_edge[2]`.
  boundary_edges <- NULL
  all_points <- coord_mat
  holes <- NULL
  inner_ring <- NULL
  outer_ring <- NULL
  outer_me <- if (!is.null(max_edge)) max_edge[length(max_edge)] else Inf

  if (!is.null(boundary)) {
    if (inherits(boundary, c("sf", "sfc", "sfg"))) {
      sf_result <- sf_to_boundary(boundary)
      boundary <- sf_result$boundary
      holes <- sf_result$holes
    }
    boundary <- as.matrix(boundary)
    if (ncol(boundary) != 2) stop("boundary must have 2 columns (x, y)")
    storage.mode(boundary) <- "double"
    inner_ring <- boundary

    # Add boundary vertices to point set
    n_orig <- nrow(coord_mat)
    n_bnd <- nrow(boundary)
    all_points <- rbind(coord_mat, boundary)

    # Create boundary edge loop (connecting consecutive boundary vertices)
    bnd_start <- n_orig + 1L
    boundary_edges <- make_edge_loop(bnd_start, n_bnd)

    # Add hole constraint edges
    if (length(holes) > 0) {
      for (hole in holes) {
        hole <- as.matrix(hole)
        n_hole <- nrow(hole)
        hole_start <- nrow(all_points) + 1L
        all_points <- rbind(all_points, hole)
        boundary_edges <- rbind(boundary_edges,
                                make_edge_loop(hole_start, n_hole))
      }
    }
  } else {
    hull_pts <- convex_hull_ring(coord_mat)
    inner_ring <- hull_pts
    diam <- sqrt(max(outer(hull_pts[, 1], hull_pts[, 1], "-")^2 +
                       outer(hull_pts[, 2], hull_pts[, 2], "-")^2))
    if (extend > 0 && diam > 0) {
      # Extended boundary: every point within extend * diameter of the
      # convex hull of the data (outward buffer with rounded corners)
      outer_ring <- buffer_convex_ring(hull_pts, extend * diam, outer_me)
      bnd_start <- nrow(coord_mat) + 1L
      all_points <- rbind(coord_mat, outer_ring)
      boundary_edges <- make_edge_loop(bnd_start, nrow(outer_ring))
    }
  }

  # Add refinement points if max_edge is specified. Track the boundary +
  # hole range so dedup never drops a constraint-edge endpoint: dropping
  # one collapses the constraint loop and CDT's eraseOuterTrianglesAndHoles
  # then wipes every triangle as "outer".
  n_coord <- nrow(coord_mat)
  n_protected <- nrow(all_points)   # coords + boundary + holes (pre-hex)
  if (!is.null(max_edge)) {
    constraint_rings <- c(list(if (is.null(outer_ring)) inner_ring else outer_ring),
                          holes)
    # Hexagonal lattice: produces more equilateral triangles than rectangular.
    # Lattice points are generated only where the mesh will have triangles.
    grid_pts <- domain_lattice(inner_ring, holes = holes,
                               spacing = max_edge[1],
                               constraint_rings = constraint_rings)
    if (!is.null(outer_ring)) {
      grid_pts <- rbind(grid_pts,
                        domain_lattice(outer_ring, holes = list(inner_ring),
                                       spacing = outer_me,
                                       constraint_rings = c(constraint_rings,
                                                            list(inner_ring))))
    }
    all_points <- rbind(all_points, grid_pts)
  }

  # Dedup pass: boundary/hole vertices are protected; coords and hex points
  # may be merged. `cutoff` sets the user-facing merge tolerance for coords
  # (and tightens hex dedup when larger). When cutoff == 0 and max_edge is
  # supplied, the hex tail is still deduplicated at max_edge[1] * 0.3 to
  # avoid hex/boundary coincidence.
  hex_tol_default <- if (!is.null(max_edge)) max_edge[1] * 0.3 else 0
  need_dedup <- cutoff > 0 || hex_tol_default > 0
  if (need_dedup) {
    dedup <- dedup_with_protected(
      all_points,
      n_coord = n_coord,
      n_protected = n_protected,
      tol_coord = cutoff,
      tol_hex = max(cutoff, hex_tol_default),
      boundary_edges = boundary_edges
    )
    all_points <- dedup$points
    boundary_edges <- dedup$boundary_edges
  }

  # Triangulate
  use_ruppert <- !is.null(min_angle) || !is.null(max_area)

  if (use_ruppert) {
    # Scale tolerance to mesh: 1e-4 of domain diameter
    xr <- range(all_points[, 1]); yr <- range(all_points[, 2])
    diam <- sqrt(diff(xr)^2 + diff(yr)^2)
    dist_tol <- max(cutoff, diam * 1e-8, 1e-12)

    result <- cpp_ruppert_refine(
      points = all_points,
      edges_nullable = boundary_edges,
      min_angle_deg = if (!is.null(min_angle)) min_angle else 0,
      max_area = if (!is.null(max_area)) max_area else 0,
      max_steiner = as.integer(max_steiner),
      min_dist_tolerance = dist_tol
    )
  } else {
    result <- cpp_triangulate(
      points = all_points,
      edges_nullable = boundary_edges,
      min_dist_tolerance = max(cutoff, 1e-10)
    )
  }

  # A constrained triangulation that erases every triangle as "outer" leaves
  # an empty mesh. Downstream FEM assembly would then produce all-zero C and G
  # and a theta-independent SPDE precision. Fail loudly here instead of
  # returning a degenerate mesh.
  if (is.null(result$triangles) || nrow(result$triangles) == 0L) {
    stop("tulpa_mesh() produced 0 triangles. The constrained triangulation ",
         "erased every triangle, which usually means the boundary/constraint ",
         "loop collapsed. Check the boundary, cutoff, and max_edge arguments.")
  }

  result$n_input_points <- n_coord
  result <- structure(result, class = "tulpa_mesh")
  if (!is.null(input_crs)) result$crs <- input_crs
  result
}

#' Compute FEM Matrices from a Mesh
#'
#' 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.
#'
#' @param mesh A `tulpa_mesh` object.
#' @param obs_coords Observation coordinates (N x 2 matrix). If NULL, the
#'   projection matrix A is the identity (observations at mesh nodes).
#' @param barrier Optional logical vector of length `n_triangles`. Triangles
#'   marked `TRUE` are treated as physical barriers (coastlines, rivers):
#'   their stiffness contributions are zeroed so the spatial field cannot
#'   smooth across them. Based on Bakka et al. (2019).
#' @param parallel Logical. If `TRUE`, uses parallel FEM assembly via
#'   RcppParallel (chunk-local triplet accumulation merged under a mutex).
#'   Beneficial for meshes with >50K triangles. Default `FALSE`.
#' @param lumped Logical. If `TRUE`, returns a diagonal lumped mass matrix
#'   C0 (vertex areas) in addition to the consistent mass matrix C. The
#'   lumped mass inverse is trivial and needed for the SPDE Q-builder.
#'   Default `FALSE`.
#'
#' @return A list with sparse matrices (dgCMatrix class from Matrix package):
#'   \itemize{
#'     \item `C`: consistent mass matrix (n_vertices x n_vertices)
#'     \item `G`: stiffness matrix (n_vertices x n_vertices)
#'     \item `A`: projection matrix (n_obs x n_vertices)
#'     \item `n_mesh`: number of mesh vertices
#'     \item `C0`: (only if `lumped = TRUE`) diagonal lumped mass matrix
#'     \item `va`: (only if `lumped = TRUE`) vertex areas (numeric vector)
#'     \item `ta`: (only if `lumped = TRUE`) triangle areas (numeric vector)
#'   }
#'
#' @export
fem_matrices <- function(mesh, obs_coords = NULL, barrier = NULL,
                        parallel = FALSE, lumped = FALSE) {
  if (!inherits(mesh, "tulpa_mesh")) {
    stop("mesh must be a tulpa_mesh object")
  }

  # A mesh with zero triangles produces all-zero C and G, which silently
  # collapses downstream SPDE precision builds to a theta-independent ridge.
  # Fail loudly here instead.
  if (is.null(mesh$triangles) || nrow(mesh$triangles) == 0L) {
    stop("mesh has 0 triangles; cannot assemble FEM matrices. ",
         "This usually means the constraint loop collapsed during ",
         "triangulation. Check the boundary, cutoff, and max_edge ",
         "arguments to tulpa_mesh().")
  }

  # Compute C and G — dispatch based on mesh type and parallel flag
  is_3d <- ncol(mesh$vertices) == 3
  if (is_3d) {
    fem <- cpp_fem_matrices_3d(mesh$vertices, mesh$triangles)
  } else if (parallel) {
    fem <- cpp_fem_matrices_parallel(mesh$vertices, mesh$triangles)
  } else {
    fem <- cpp_fem_matrices(mesh$vertices, mesh$triangles)
  }

  # Assemble into sparse matrices (0-based triplets to 1-based for Matrix)
  n <- fem$n
  C <- Matrix::sparseMatrix(
    i = fem$C_i + 1L, j = fem$C_j + 1L, x = fem$C_x,
    dims = c(n, n), repr = "C"
  )
  G <- Matrix::sparseMatrix(
    i = fem$G_i + 1L, j = fem$G_j + 1L, x = fem$G_x,
    dims = c(n, n), repr = "C"
  )

  # Projection matrix A
  if (!is.null(obs_coords)) {
    obs_coords <- as.matrix(obs_coords)
    storage.mode(obs_coords) <- "double"
    if (is_3d) {
      # For spherical meshes, convert lon/lat obs_coords to xyz if needed
      if (ncol(obs_coords) == 2) {
        radius <- if (!is.null(mesh$sphere)) mesh$sphere$radius else 1
        lon_rad <- obs_coords[, 1] * pi / 180
        lat_rad <- obs_coords[, 2] * pi / 180
        obs_coords <- cbind(
          radius * cos(lat_rad) * cos(lon_rad),
          radius * cos(lat_rad) * sin(lon_rad),
          radius * sin(lat_rad)
        )
      }
      proj <- cpp_projection_matrix_3d(obs_coords, mesh$vertices,
                                        mesh$triangles)
    } else {
      proj <- cpp_projection_matrix(obs_coords, mesh$vertices, mesh$triangles)
    }
    A <- Matrix::sparseMatrix(
      i = proj$i + 1L, j = proj$j + 1L, x = proj$x,
      dims = c(proj$nrow, proj$ncol), repr = "C"
    )
  } else {
    A <- Matrix::Diagonal(n)
  }

  # Barrier model: zero out stiffness contributions from barrier triangles
  # This prevents the spatial field from smoothing across barriers.
  # See Bakka et al. (2019) for the barrier model formulation.
  if (!is.null(barrier)) {
    barrier <- as.logical(barrier)
    if (length(barrier) != nrow(mesh$triangles)) {
      stop(sprintf(
        "barrier must be a logical vector of length n_triangles (%d), got %d",
        nrow(mesh$triangles), length(barrier)
      ))
    }
    # Rebuild G excluding barrier triangles
    # Re-run C++ FEM but with a triangle mask
    fem_barrier <- cpp_fem_matrices(mesh$vertices, mesh$triangles)
    # Identify which triplet entries come from barrier triangles
    # Each triangle contributes exactly 9 entries to both C and G
    n_tri <- nrow(mesh$triangles)
    tri_mask <- rep(!barrier, each = 9)

    # Rebuild G with barrier triangles zeroed out
    G <- Matrix::sparseMatrix(
      i = fem_barrier$G_i[tri_mask] + 1L,
      j = fem_barrier$G_j[tri_mask] + 1L,
      x = fem_barrier$G_x[tri_mask],
      dims = c(n, n), repr = "C"
    )
  }

  result <- list(C = C, G = G, A = A, n_mesh = n)

  if (lumped) {
    va <- as.numeric(Matrix::rowSums(C))
    result$C0 <- Matrix::Diagonal(x = va)
    result$va <- va

    # Triangle areas
    v <- mesh$vertices
    tri <- mesh$triangles
    ta <- numeric(nrow(tri))
    for (t in seq_len(nrow(tri))) {
      p0 <- v[tri[t, 1], ]
      p1 <- v[tri[t, 2], ]
      p2 <- v[tri[t, 3], ]
      ta[t] <- abs((p1[1] - p0[1]) * (p2[2] - p0[2]) -
                    (p2[1] - p0[1]) * (p1[2] - p0[2])) / 2
    }
    result$ta <- ta
  }

  result
}

# Internal: convert sf geometry to boundary coords + hole constraints.
# Returns a list with:
#   $boundary: N x 2 matrix of outer boundary coordinates
#   $holes: list of M x 2 matrices of hole boundary coordinates (may be empty)
sf_to_boundary <- function(geom) {
  if (!requireNamespace("sf", quietly = TRUE)) {
    stop("sf package required for polygon boundaries")
  }

  # Normalize to a single sfc_POLYGON or sfc_MULTIPOLYGON
  geom <- sf::st_geometry(geom)
  if (length(geom) > 1) {
    geom <- sf::st_union(geom)
  }
  geom <- geom[[1]]  # extract the raw sfg

  # Collect all polygon rings
  if (inherits(geom, "MULTIPOLYGON")) {
    # MULTIPOLYGON: list of polygons, each a list of rings
    polygons <- lapply(geom, identity)
  } else if (inherits(geom, "POLYGON")) {
    polygons <- list(geom)
  } else {
    stop(sprintf("unsupported geometry type: %s", class(geom)[1]))
  }

  # Extract outer rings and holes from all component polygons
  outer_rings <- list()
  hole_rings <- list()

  for (poly in polygons) {
    rings <- lapply(poly, identity)  # each ring is a matrix
    # First ring is outer boundary
    outer_rings <- c(outer_rings, list(close_ring(rings[[1]])))
    # Remaining rings are holes
    if (length(rings) > 1) {
      for (k in 2:length(rings)) {
        hole_rings <- c(hole_rings, list(close_ring(rings[[k]])))
      }
    }
  }

  # For the boundary, combine all outer rings
  # (for MULTIPOLYGON with separate components, the user should union first
  #  or we take the largest component as outer boundary)
  if (length(outer_rings) == 1) {
    boundary <- outer_rings[[1]]
  } else {
    # Use the ring with the largest area as the outer boundary,
    # treat other outer rings as additional constraints
    areas <- vapply(outer_rings, ring_area, numeric(1))
    main_idx <- which.max(areas)
    boundary <- outer_rings[[main_idx]]
    # Other outer rings become "island" constraints (treated like holes
    # but with reversed winding, so they remain filled)
    for (idx in seq_along(outer_rings)) {
      if (idx != main_idx) {
        hole_rings <- c(hole_rings, list(outer_rings[[idx]]))
      }
    }
  }

  list(boundary = boundary, holes = hole_rings)
}

# Remove closing point from a ring matrix (last row == first row)
close_ring <- function(ring) {
  ring <- ring[, 1:2, drop = FALSE]
  if (nrow(ring) > 1 && all(ring[1, ] == ring[nrow(ring), ])) {
    ring <- ring[-nrow(ring), , drop = FALSE]
  }
  ring
}

# Shoelace area of a ring (absolute value)
ring_area <- function(ring) {
  n <- nrow(ring)
  x <- ring[, 1]
  y <- ring[, 2]
  abs(sum(x * c(y[-1], y[1]) - c(x[-1], x[1]) * y)) / 2
}

# Deduplicate `points` so that:
#   - Boundary + hole rows (indices coord_n+1 .. protected_n) are never
#     dropped (they are constraint-edge endpoints).
#   - Coord rows (1 .. coord_n) may be dropped when within `tol_coord` of a
#     boundary point or an earlier kept coord.
#   - Hex rows (protected_n+1 .. nrow(points)) may be dropped when within
#     `tol_hex` of any earlier kept point.
# Iteration order is boundary -> coords -> hex so boundary always wins when
# a coord and a boundary point coincide (the coord is dropped, the
# constraint endpoint survives). Returns the dedupped points and the
# boundary_edges with indices remapped to the new row positions.
dedup_with_protected <- function(points, n_coord, n_protected,
                                 tol_coord, tol_hex, boundary_edges) {
  n <- nrow(points)
  boundary_idx <- if (n_protected > n_coord) seq.int(n_coord + 1L, n_protected) else integer(0)
  coord_idx <- if (n_coord >= 1L) seq_len(n_coord) else integer(0)
  hex_idx <- if (n > n_protected) seq.int(n_protected + 1L, n) else integer(0)

  order_idx <- c(boundary_idx, coord_idx, hex_idx)
  n_bnd <- length(boundary_idx)
  n_bnd_coord <- n_bnd + length(coord_idx)

  keep <- rep(TRUE, n)
  if (length(order_idx) >= 2L && (tol_coord > 0 || tol_hex > 0)) {
    for (k in seq.int(2L, length(order_idx))) {
      i <- order_idx[k]
      if (k <= n_bnd) next  # protected: constraint endpoint
      tol_i <- if (k <= n_bnd_coord) tol_coord else tol_hex
      if (tol_i <= 0) next
      earlier <- order_idx[seq_len(k - 1L)]
      kept_earlier <- earlier[keep[earlier]]
      if (!length(kept_earlier)) next
      diffs <- points[kept_earlier, , drop = FALSE] -
        matrix(points[i, ], nrow = length(kept_earlier), ncol = 2L, byrow = TRUE)
      if (any(sqrt(rowSums(diffs^2)) < tol_i)) keep[i] <- FALSE
    }
  }

  if (!is.null(boundary_edges) && any(!keep)) {
    remap <- cumsum(keep); remap[!keep] <- 0L
    boundary_edges[, 1] <- remap[boundary_edges[, 1]]
    boundary_edges[, 2] <- remap[boundary_edges[, 2]]
    valid <- boundary_edges[, 1] > 0 & boundary_edges[, 2] > 0
    boundary_edges <- boundary_edges[valid, , drop = FALSE]
  }

  list(points = points[keep, , drop = FALSE],
       boundary_edges = boundary_edges)
}

# Build an edge loop connecting consecutive 1-based vertex indices
# start_idx: 1-based index of first vertex in the loop
# n: number of vertices in the loop
make_edge_loop <- function(start_idx, n) {
  edges <- matrix(0L, nrow = n, ncol = 2)
  for (i in seq_len(n)) {
    edges[i, 1] <- start_idx + i - 1L
    edges[i, 2] <- start_idx + (i %% n)
  }
  edges
}

# Generate hexagonal lattice points within a bounding box.
# Row spacing = spacing * sqrt(3)/2; odd rows offset by spacing/2.
# Produces more equilateral triangles than a rectangular grid.
hex_lattice <- function(x_min, x_max, y_min, y_max, spacing) {
  row_height <- spacing * sqrt(3) / 2
  ys <- seq(y_min, y_max, by = row_height)
  pts <- vector("list", length(ys))
  for (i in seq_along(ys)) {
    offset <- if (i %% 2 == 0) spacing / 2 else 0
    xs <- seq(x_min + offset, x_max, by = spacing)
    pts[[i]] <- cbind(xs, rep(ys[i], length(xs)))
  }
  do.call(rbind, pts)
}

# Convex hull of a point set as a counter-clockwise ring without repeated
# vertices.
convex_hull_ring <- function(points) {
  ring <- points[rev(grDevices::chull(points)), , drop = FALSE]
  if (nrow(ring) > 1L) {
    nxt <- c(seq.int(2L, nrow(ring)), 1L)
    ring <- ring[rowSums((ring - ring[nxt, , drop = FALSE])^2) > 0, ,
                 drop = FALSE]
  }
  ring
}

# Outward buffer of a convex counter-clockwise ring by distance `d`: each
# edge is shifted along its outward normal and each corner is replaced by a
# circular arc of radius `d` around the hull vertex. Arcs use steps of at
# most pi/8; arcs and straight edges are subdivided so that no boundary
# segment is longer than `max_len`. Returns a counter-clockwise ring.
buffer_convex_ring <- function(ring, d, max_len = Inf) {
  n <- nrow(ring)
  nxt <- c(seq.int(2L, n), 1L)
  prv <- c(n, seq_len(n - 1L))
  edge <- ring[nxt, , drop = FALSE] - ring
  edge_len <- sqrt(rowSums(edge^2))
  normal <- cbind(edge[, 2], -edge[, 1]) / edge_len  # outward for CCW

  out <- vector("list", 2L * n)
  for (i in seq_len(n)) {
    n_in <- normal[prv[i], ]
    n_out <- normal[i, ]
    theta0 <- atan2(n_in[2], n_in[1])
    turn <- atan2(n_in[1] * n_out[2] - n_in[2] * n_out[1],
                  sum(n_in * n_out))
    turn <- if (turn < -1e-9) turn + 2 * pi else max(turn, 0)
    k <- max(1L, ceiling(turn / (pi / 8)), ceiling(turn * d / max_len))
    theta <- theta0 + turn * seq(0, 1, length.out = k + 1L)
    arc <- cbind(ring[i, 1] + d * cos(theta), ring[i, 2] + d * sin(theta))
    if (turn == 0) arc <- arc[1L, , drop = FALSE]
    out[[2L * i - 1L]] <- arc

    m <- ceiling(edge_len[i] / max_len)
    if (is.finite(m) && m > 1L) {
      t <- seq_len(m - 1L) / m
      a <- ring[i, ] + d * n_out
      out[[2L * i]] <- cbind(a[1] + t * edge[i, 1], a[2] + t * edge[i, 2])
    }
  }
  do.call(rbind, out)
}

# Even-odd point-in-polygon test, vectorised over points.
points_in_polygon <- function(px, py, poly) {
  poly <- as.matrix(poly)
  n <- nrow(poly)
  inside <- rep(FALSE, length(px))
  for (i in seq_len(n)) {
    j <- if (i == n) 1L else i + 1L
    x0 <- poly[i, 1]; y0 <- poly[i, 2]
    x1 <- poly[j, 1]; y1 <- poly[j, 2]
    crosses <- (y0 <= py & y1 > py) | (y1 <= py & y0 > py)
    x_cross <- x0 + (py - y0) / (y1 - y0) * (x1 - x0)
    flip <- crosses & px < x_cross
    inside[flip] <- !inside[flip]
  }
  inside
}

# Distance from each point to the nearest edge of a closed ring.
points_ring_distance <- function(px, py, ring) {
  ring <- as.matrix(ring)
  n <- nrow(ring)
  dmin <- rep(Inf, length(px))
  for (i in seq_len(n)) {
    j <- if (i == n) 1L else i + 1L
    ex <- ring[j, 1] - ring[i, 1]
    ey <- ring[j, 2] - ring[i, 2]
    len2 <- ex^2 + ey^2
    t <- if (len2 > 0) ((px - ring[i, 1]) * ex + (py - ring[i, 2]) * ey) / len2 else 0
    t <- pmin(pmax(t, 0), 1)
    dmin <- pmin(dmin, sqrt((px - ring[i, 1] - t * ex)^2 +
                              (py - ring[i, 2] - t * ey)^2))
  }
  dmin
}

# Hexagonal lattice with the given spacing, restricted to the region inside
# `ring` and outside every ring in `holes`. Points closer than half a spacing
# to any ring in `constraint_rings` are dropped so lattice points do not form
# slivers against constraint edges.
domain_lattice <- function(ring, holes = NULL, spacing, constraint_rings = list()) {
  ring <- as.matrix(ring)
  pts <- hex_lattice(min(ring[, 1]), max(ring[, 1]),
                     min(ring[, 2]), max(ring[, 2]), spacing = spacing)
  keep <- points_in_polygon(pts[, 1], pts[, 2], ring)
  for (h in holes) {
    keep <- keep & !points_in_polygon(pts[, 1], pts[, 2], h)
  }
  for (r in constraint_rings) {
    if (!any(keep)) break
    keep[keep] <- points_ring_distance(pts[keep, 1], pts[keep, 2], r) >=
      spacing / 2
  }
  pts[keep, , drop = FALSE]
}


# =====================================================================
# Mesh quality diagnostics
# =====================================================================

#' Mesh Quality Diagnostics
#'
#' 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.
#'
#' @param mesh A `tulpa_mesh` object.
#'
#' @return A data.frame with one row per triangle and columns:
#'   \itemize{
#'     \item `min_angle`: minimum interior angle (degrees)
#'     \item `max_angle`: maximum interior angle (degrees)
#'     \item `aspect_ratio`: longest edge / shortest edge (1 = equilateral)
#'     \item `area`: triangle area
#'   }
#'
#' @export
#'
#' @examples
#' set.seed(42)
#' mesh <- tulpa_mesh(cbind(runif(50), runif(50)))
#' q <- mesh_quality(mesh)
#' summary(q)
mesh_quality <- function(mesh) {
  if (!inherits(mesh, "tulpa_mesh")) stop("mesh must be a tulpa_mesh object")

  v <- mesh$vertices
  tri <- mesh$triangles
  n_tri <- nrow(tri)

  min_angle <- numeric(n_tri)
  max_angle <- numeric(n_tri)
  aspect_ratio <- numeric(n_tri)
  area <- numeric(n_tri)

  for (t in seq_len(n_tri)) {
    p1 <- v[tri[t, 1], ]
    p2 <- v[tri[t, 2], ]
    p3 <- v[tri[t, 3], ]

    # Edge lengths
    a <- sqrt(sum((p2 - p3)^2))  # opposite p1
    b <- sqrt(sum((p1 - p3)^2))  # opposite p2
    c <- sqrt(sum((p1 - p2)^2))  # opposite p3

    # Angles via law of cosines
    angles <- numeric(3)
    angles[1] <- acos(pmin(1, pmax(-1, (b^2 + c^2 - a^2) / (2 * b * c))))
    angles[2] <- acos(pmin(1, pmax(-1, (a^2 + c^2 - b^2) / (2 * a * c))))
    angles[3] <- pi - angles[1] - angles[2]

    angles_deg <- angles * 180 / pi
    min_angle[t] <- min(angles_deg)
    max_angle[t] <- max(angles_deg)

    edges <- sort(c(a, b, c))
    aspect_ratio[t] <- edges[3] / max(edges[1], 1e-15)

    # Signed area
    area[t] <- abs((p2[1] - p1[1]) * (p3[2] - p1[2]) -
                   (p3[1] - p1[1]) * (p2[2] - p1[2])) / 2
  }

  data.frame(
    min_angle = min_angle,
    max_angle = max_angle,
    aspect_ratio = aspect_ratio,
    area = area
  )
}

#' Summarize Mesh Quality
#'
#' Prints a summary of mesh quality metrics.
#'
#' @param mesh A `tulpa_mesh` object.
#' @return Invisibly returns the quality data.frame.
#' @export
mesh_summary <- function(mesh) {
  q <- mesh_quality(mesh)
  cat("Mesh quality summary:\n")
  cat(sprintf("  Triangles:      %d\n", nrow(q)))
  cat(sprintf("  Min angle:      %.1f / %.1f / %.1f (min / median / max)\n",
              min(q$min_angle), median(q$min_angle), max(q$min_angle)))
  cat(sprintf("  Max angle:      %.1f / %.1f / %.1f\n",
              min(q$max_angle), median(q$max_angle), max(q$max_angle)))
  cat(sprintf("  Aspect ratio:   %.2f / %.2f / %.2f\n",
              min(q$aspect_ratio), median(q$aspect_ratio), max(q$aspect_ratio)))
  cat(sprintf("  Area:           %.4f / %.4f / %.4f\n",
              min(q$area), median(q$area), max(q$area)))

  # Flag bad triangles
  bad <- sum(q$min_angle < 10)
  if (bad > 0) {
    cat(sprintf("  WARNING: %d triangles with min angle < 10 degrees\n", bad))
  }
  very_bad <- sum(q$min_angle < 5)
  if (very_bad > 0) {
    cat(sprintf("  WARNING: %d triangles with min angle < 5 degrees (slivers)\n", very_bad))
  }

  invisible(q)
}