---
title: "Spherical and Temporal Meshes"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Spherical and Temporal Meshes}
  %\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)
```

## Spherical meshes

For global-scale species distribution models or climate fields, you need a mesh on the sphere. `tulpa_mesh_sphere()` generates one by recursive subdivision of an icosahedron.

```{r sphere-basic}
globe <- tulpa_mesh_sphere(subdivisions = 3)
globe
```

Each subdivision level quadruples the triangle count:

| Level | Vertices | Triangles |
|-------|----------|-----------|
| 0     | 12       | 20        |
| 1     | 42       | 80        |
| 2     | 162      | 320       |
| 3     | 642      | 1280      |
| 4     | 2562     | 5120      |

All vertices lie exactly on the sphere surface:

```{r sphere-radius}
radii <- sqrt(rowSums(globe$vertices^2))
range(radii)
```

### FEM on the sphere

The 3D FEM assembly uses the metric tensor on each triangle's tangent plane:

```{r sphere-fem}
fem <- fem_matrices(globe, lumped = TRUE)

# Total surface area should approximate 4*pi for unit sphere
cat("Total area:", sum(fem$va), "\n")
cat("4*pi:      ", 4 * pi, "\n")
cat("Error:     ", abs(sum(fem$va) - 4 * pi) / (4 * pi) * 100, "%\n")
```

### Projection with lon/lat coordinates

Pass observation coordinates as (longitude, latitude) in degrees:

```{r sphere-proj}
obs <- cbind(lon = c(0, 90, -45, 170), lat = c(0, 45, -30, 60))
fem <- fem_matrices(globe, obs_coords = obs)
dim(fem$A)
range(Matrix::rowSums(fem$A))  # row sums = 1
```

### Earth-radius meshes

```{r sphere-earth}
earth <- tulpa_mesh_sphere(subdivisions = 3, radius = 6371)
sqrt(sum(earth$vertices[1, ]^2))  # 6371 km
```

### Euler characteristic

A closed sphere satisfies V - E + T = 2:

```{r euler}
globe$n_vertices - globe$n_edges + globe$n_triangles
```

## 1D temporal meshes

For spatio-temporal SPDE models, the temporal dimension needs its own mesh. `tulpa_mesh_1d()` produces knots with tridiagonal FEM matrices that combine with the spatial mesh via Kronecker products.

```{r mesh1d}
# Monthly observations over 5 years
times <- seq(2020, 2025, by = 1/12)
m1d <- tulpa_mesh_1d(times)
m1d
```

The extension knots (controlled by `n_extend`) prevent boundary effects:

```{r mesh1d-range}
range(times)       # data range
range(m1d$knots)   # mesh range (extended)
```

### 1D FEM properties

Mass and stiffness matrices are tridiagonal and sparse:

```{r mesh1d-fem}
# No extension for cleaner demonstration
m <- tulpa_mesh_1d(seq(0, 1, by = 0.1), n_extend = 0)

# Symmetric, positive definite mass
Matrix::isSymmetric(m$C)
all(Matrix::diag(m$C) > 0)

# Stiffness row sums = 0
max(abs(Matrix::rowSums(m$G)))

# Total mass = domain length
sum(m$C)
```

### Irregular time spacing

The 1D mesh handles irregular intervals naturally:

```{r mesh1d-irregular}
# Denser sampling in summer, sparser in winter
times_irr <- c(
  seq(2020, 2020.25, by = 1/52),      # weekly Jan-Mar
  seq(2020.25, 2020.75, by = 1/365),   # daily Apr-Sep
  seq(2020.75, 2021, by = 1/52)        # weekly Oct-Dec
)
m_irr <- tulpa_mesh_1d(times_irr, n_extend = 0)
cat(m_irr$n, "knots from", length(times_irr), "unique time points\n")
```

## Metric graph meshes

For SPDE models along roads, rivers, or coastlines, `tulpa_mesh_graph()` builds a 1D FEM mesh along network edges.

```{r graph}
# Simple river network: main channel + tributary
edges <- list(
  cbind(x = seq(0, 10, by = 0.5), y = rep(0, 21)),     # main channel
  cbind(x = c(5, 5, 5, 5), y = c(0, 2, 4, 6))          # tributary
)

g <- tulpa_mesh_graph(edges)
g
```

### Junction detection

Nodes where edges meet are automatically detected:

```{r graph-junctions}
cat("Junctions (degree > 2):", sum(g$degree > 2), "\n")
cat("Endpoints (degree = 1):", sum(g$degree == 1), "\n")
```

### FEM on graphs

The graph FEM matrices are structurally identical to the 1D case but assembled across the full network:

```{r graph-fem}
Matrix::isSymmetric(g$C)
max(abs(Matrix::rowSums(g$G)))  # row sums ~ 0
```

### Subdividing long edges

```{r graph-refine}
g_fine <- tulpa_mesh_graph(edges, max_edge = 0.3)
cat("Coarse:", g$n_vertices, "vertices\n")
cat("Fine:  ", g_fine$n_vertices, "vertices\n")
```

### sf LINESTRING input

```{r graph-sf, eval = requireNamespace("sf", quietly = TRUE)}
library(sf)
line1 <- st_linestring(cbind(c(0, 5, 10), c(0, 3, 0)))
line2 <- st_linestring(cbind(c(5, 5), c(3, 8)))
g_sf <- tulpa_mesh_graph(st_sfc(line1, line2), max_edge = 1)
g_sf
```

## Non-stationary FEM

For fields where the range or variance changes across space, `fem_matrices_nonstationary()` computes weighted FEM matrices in C++:

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

# Range decreases from left to right
kappa <- sqrt(8) / (2 - mesh$vertices[, 1])  # shorter range on the right
tau <- rep(1, n)

ns <- fem_matrices_nonstationary(mesh, kappa, tau)
names(ns)
```

With constant kappa/tau, the weighted matrices are simply scaled versions of the standard ones:

```{r ns-constant}
ns_const <- fem_matrices_nonstationary(mesh, rep(2, n), rep(3, n))
fem <- fem_matrices(mesh)

max(abs(ns_const$Ck - 4 * fem$C))  # kappa^2 = 4
max(abs(ns_const$Ct - 9 * fem$C))  # tau^2 = 9
```

## P2 quadratic elements

For higher accuracy with fewer mesh nodes, use 6-node quadratic elements:

```{r p2}
set.seed(42)
mesh <- tulpa_mesh(cbind(runif(30), runif(30)))
p2 <- fem_matrices_p2(mesh)

cat("P1 nodes:", mesh$n_vertices, "\n")
cat("P2 nodes:", p2$n_mesh, "(", p2$n_vertices, "vertices +",
    p2$n_midpoints, "midpoints)\n")
```

The total area is preserved:

```{r p2-area}
fem_p1 <- fem_matrices(mesh)
cat("P1 total area:", sum(fem_p1$C), "\n")
cat("P2 total area:", sum(p2$C), "\n")
```

## Parallel FEM assembly

For large meshes (>50K triangles), enable parallel assembly:

```{r parallel}
set.seed(42)
mesh_large <- tulpa_mesh(cbind(runif(500), runif(500)), max_edge = 0.03)
cat(mesh_large$n_triangles, "triangles\n")

fem_seq <- fem_matrices(mesh_large, parallel = FALSE)
fem_par <- fem_matrices(mesh_large, parallel = TRUE)

# Results are identical
max(abs(fem_seq$C - fem_par$C))
```