---
title: "Spatial Workflows"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Spatial Workflows}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  dev = "svglite",
  fig.ext = "svg"
)
library(tulpaMesh)
```

## Custom boundaries

By default, `tulpa_mesh()` uses the convex hull extended by 10%. You can provide an explicit boundary as a coordinate matrix:

```{r boundary}
set.seed(42)
coords <- cbind(runif(80, 1, 9), runif(80, 1, 9))

# L-shaped boundary
bnd <- rbind(
  c(0, 0), c(10, 0), c(10, 5),
  c(5, 5), c(5, 10), c(0, 10)
)

mesh <- tulpa_mesh(coords, boundary = bnd, max_edge = 1, extend = 0)
plot(mesh, vertex_col = "black", main = "L-shaped domain")
```

## sf polygon boundaries

Pass an sf polygon directly. CRS is preserved automatically.

```{r sf-boundary, eval = requireNamespace("sf", quietly = TRUE)}
library(sf)

poly <- st_polygon(list(rbind(
  c(0, 0), c(10, 0), c(10, 10), c(0, 10), c(0, 0)
)))
sfc <- st_sfc(poly, crs = 32633)

mesh_sf <- tulpa_mesh(coords, boundary = sfc, max_edge = 1, extend = 0)
mesh_crs(mesh_sf)  # CRS attached
```

### Polygons with holes

Holes in sf polygons become constraint edges that the mesh respects:

```{r holes, eval = requireNamespace("sf", quietly = TRUE)}
outer <- rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 10), c(0, 0))
hole  <- rbind(c(3, 3), c(7, 3), c(7, 7), c(3, 7), c(3, 3))
poly_h <- st_polygon(list(outer, hole))

# Remove points inside the hole
inside <- coords[, 1] > 3 & coords[, 1] < 7 & coords[, 2] > 3 & coords[, 2] < 7
pts_outside <- coords[!inside, ]

mesh_h <- tulpa_mesh(pts_outside, boundary = st_sfc(poly_h), max_edge = 1, extend = 0)
plot(mesh_h, vertex_col = "black", main = "Mesh with hole")
```

## Barrier models

Physical barriers (coastlines, rivers) prevent the spatial field from smoothing across them. Mark barrier triangles and pass them to `fem_matrices()`:

```{r barrier}
set.seed(42)
coords <- cbind(runif(100, 0, 10), runif(100, 0, 10))
mesh <- tulpa_mesh(coords, max_edge = 0.8)

# A river running through the domain
river <- rbind(c(4, 0), c(6, 5), c(4, 10))
bt <- barrier_triangles(mesh, river)

cat(sum(bt), "barrier triangles out of", mesh$n_triangles, "\n")

# FEM with barrier: stiffness is zeroed for barrier triangles
fem_b <- fem_matrices(mesh, barrier = bt)

# Compare: barrier mesh has fewer stiffness nonzeros
fem_n <- fem_matrices(mesh)
cat("Nonzeros without barrier:", length(fem_n$G@x), "\n")
cat("Nonzeros with barrier:   ", length(fem_b$G@x), "\n")
```

## Mesh subdivision

Split every triangle into 4 for uniform refinement:

```{r subdivide}
set.seed(42)
mesh <- tulpa_mesh(cbind(runif(20), runif(20)))
sub <- subdivide_mesh(mesh)
cat("Original:", mesh$n_triangles, "triangles\n")
cat("Subdivided:", sub$n_triangles, "triangles\n")
```

## Adaptive refinement

Refine only where error indicators are high. This is the typical workflow after an initial SPDE solve: the solver returns posterior variance per triangle, and you refine where variance is large.

```{r adaptive}
set.seed(42)
mesh <- tulpa_mesh(cbind(runif(50), runif(50)), max_edge = 0.15)

# Simulate error indicators (high in one corner)
q <- mesh_quality(mesh)
centroids <- cbind(
  (mesh$vertices[mesh$triangles[,1], 1] + mesh$vertices[mesh$triangles[,2], 1] +
   mesh$vertices[mesh$triangles[,3], 1]) / 3,
  (mesh$vertices[mesh$triangles[,1], 2] + mesh$vertices[mesh$triangles[,2], 2] +
   mesh$vertices[mesh$triangles[,3], 2]) / 3
)
error <- exp(-3 * centroids[, 1])  # high error near x = 0

refined <- refine_mesh(mesh, error, fraction = 0.3)
cat("Before:", mesh$n_triangles, "triangles\n")
cat("After: ", refined$n_triangles, "triangles\n")
```

## Mesh operations

### Extract submesh

```{r subset}
set.seed(42)
mesh <- tulpa_mesh(cbind(runif(50, 0, 10), runif(50, 0, 10)), max_edge = 1)

# Keep only left half
q <- mesh_quality(mesh)
centroids_x <- (mesh$vertices[mesh$triangles[,1], 1] +
                mesh$vertices[mesh$triangles[,2], 1] +
                mesh$vertices[mesh$triangles[,3], 1]) / 3
left <- subset_mesh(mesh, centroids_x < 5)
left
```

### Find disconnected components

```{r components}
comps <- mesh_components(mesh)
cat("Number of components:", max(comps), "\n")
```

## Converting from fmesher

If you have an existing fmesher mesh, convert it directly:

```{r convert, eval = requireNamespace("fmesher", quietly = TRUE)}
library(fmesher)
fm <- fm_mesh_2d(loc = coords, max.edge = c(1, 3))
tm <- as_tulpa_mesh(fm)
tm
```

FEM matrices from the converted mesh match fmesher's output exactly.