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

## ----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-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

## ----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------------------------------------------------------------------
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")

## ----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-----------------------------------------------------------------
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")

## ----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

## ----components---------------------------------------------------------------
comps <- mesh_components(mesh)
cat("Number of components:", max(comps), "\n")

## ----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