tests/testthat/test-extend.R

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

# Regression tests for issue #4: vertices outside every triangle and the
# reach of `extend`.

n_orphans <- function(mesh) {
  mesh$n_vertices - length(unique(as.vector(mesh$triangles)))
}

test_that("every vertex belongs to a triangle and carries mass", {
  set.seed(1)
  xy <- cbind(runif(300), runif(300))
  for (args in list(
    list(max_edge = c(0.05, 0.3), cutoff = 0.02),
    list(max_edge = c(0.05, 0.3), cutoff = 0),
    list(max_edge = 0.1, cutoff = 0),
    list(max_edge = 0.1, extend = 0),
    list(min_angle = 25, max_edge = c(0.1, 0.3)),
    list()
  )) {
    m <- do.call(tulpa_mesh, c(list(xy), args))
    expect_equal(n_orphans(m), 0L)
    expect_equal(m$n_vertices, nrow(m$vertices))
    expect_true(all(Matrix::rowSums(fem_matrices(m)$C) > 0))
  }
})

test_that("input points outside a user boundary are not returned as vertices", {
  set.seed(2)
  xy <- cbind(runif(100, -0.5, 1.5), runif(100, -0.5, 1.5))
  bnd <- cbind(c(0, 1, 1, 0), c(0, 0, 1, 1))
  m <- tulpa_mesh(xy, boundary = bnd, max_edge = 0.1)
  expect_equal(n_orphans(m), 0L)
  expect_true(all(m$vertices >= -1e-12 & m$vertices <= 1 + 1e-12))
})

test_that("extend reaches the documented fraction of the domain diameter", {
  set.seed(1)
  xy <- cbind(runif(300), runif(300))
  diam <- max(dist(xy))
  for (ext in c(0.1, 0.3)) {
    m <- tulpa_mesh(xy, max_edge = c(0.05, 0.3), extend = ext)
    reach <- c(min(xy[, 1]) - min(m$vertices[, 1]),
               max(m$vertices[, 1]) - max(xy[, 1]),
               min(xy[, 2]) - min(m$vertices[, 2]),
               max(m$vertices[, 2]) - max(xy[, 2]))
    # Arcs are polygonal with steps of at most pi/8, so the reach in an axis
    # direction lies between d cos(pi/16) and d.
    expect_true(all(reach <= ext * diam + 1e-9))
    expect_true(all(reach >= ext * diam * cos(pi / 16)))
  }
})

test_that("mesh area matches the buffered convex hull", {
  set.seed(3)
  xy <- cbind(runif(200), runif(200))
  hull <- xy[grDevices::chull(xy), ]
  perim <- sum(sqrt(rowSums((hull - hull[c(2:nrow(hull), 1), ])^2)))
  hull_area <- abs(sum(hull[, 1] * hull[c(2:nrow(hull), 1), 2] -
                         hull[c(2:nrow(hull), 1), 1] * hull[, 2])) / 2
  d <- 0.1 * max(dist(xy))
  m <- tulpa_mesh(xy, max_edge = c(0.1, 0.2))
  mesh_area <- sum(mesh_quality(m)$area)
  # Buffer area hull + perimeter * d + pi d^2, minus the arc-chord deficit
  expect_lt(mesh_area, hull_area + perim * d + pi * d^2)
  expect_gt(mesh_area, hull_area + perim * d + pi * d^2 * 0.97)
})

test_that("outer max_edge controls the extension boundary resolution", {
  set.seed(1)
  xy <- cbind(runif(300), runif(300))
  m <- tulpa_mesh(xy, max_edge = c(0.05, 0.1))
  tri <- m$triangles
  e <- rbind(tri[, 1:2], tri[, 2:3], tri[, c(3, 1)])
  e <- cbind(pmin(e[, 1], e[, 2]), pmax(e[, 1], e[, 2]))
  key <- paste(e[, 1], e[, 2])
  bnd <- e[!(key %in% key[duplicated(key)]), , drop = FALSE]
  len <- sqrt(rowSums((m$vertices[bnd[, 1], ] - m$vertices[bnd[, 2], ])^2))
  expect_gt(nrow(bnd), 0L)
  expect_true(all(len <= 0.1 + 1e-9))
})

test_that("n_input_points counts the supplied coordinates", {
  set.seed(1)
  xy <- cbind(runif(300), runif(300))
  m <- tulpa_mesh(xy, max_edge = c(0.05, 0.3))
  expect_equal(m$n_input_points, 300L)
})

test_that("collinear points with an extension give a buffered mesh", {
  collinear <- cbind(seq(0, 1, length.out = 20), seq(0, 1, length.out = 20))
  m <- tulpa_mesh(collinear)
  expect_gt(m$n_triangles, 0L)
  expect_equal(n_orphans(m), 0L)
})