// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package integrate

import (
	"fmt"
	"math"
	"slices"

	"sourcedock.dev/petrbalvin/tensor/internal/base"
	"sourcedock.dev/petrbalvin/tensor/internal/core"
	linalg "sourcedock.dev/petrbalvin/tensor/linalg"
)

// The finite element groundwork for second-order problems in three
// dimensions, the volumetric sibling of the triangular surface in
// fem2d.go: piecewise-linear (P1) elements on a conforming
// tetrahedral mesh, the stiffness matrix assembled per tetrahedron
// from the gradient-of-basis formula over the element's edge vectors,
// the load integrated per element with a collapsed Gauss rule,
// Dirichlet values eliminated by lifting, Neumann fluxes integrated
// on prescribed boundary faces, and the reduced system handed to the
// same sparse Cholesky factorisation the two-dimensional path uses.

// TetraMesh3D carries a conforming tetrahedral mesh: vertex
// coordinates as x,y,z triples and tetrahedra as quadruples of vertex
// indices in positive orientation, meaning the signed volume
// (b−a)·((c−a)×(d−a)) of every stored tetrahedron is positive. A
// tetrahedron with zero volume or negative orientation does matter
// and is refused at construction.
type TetraMesh3D struct {
	// Vertices holds x,y,z for every vertex: three entries per vertex.
	Vertices []float64
	// Tetrahedra holds four vertex indices per tetrahedron.
	Tetrahedra []int64
}

// NewTetraMesh3D builds a mesh from a vertex table with three columns
// and a tetrahedron table with four columns of vertex indices.
// Indices must lie in range, every coordinate must be finite, and a
// degenerate (zero-volume) or inverted (negative-orientation)
// tetrahedron is an error naming the element and its vertices: its
// stiffness contribution is undefined.
func NewTetraMesh3D(vertices *core.Array, tetrahedra *core.Array) (*TetraMesh3D, error) {
	const name = "NewTetraMesh3D"
	if vertices.Dtype() == core.Complex || tetrahedra.Dtype() == core.Complex {
		return nil, base.Errf("%s: complex mesh data is not supported", name)
	}
	if vertices.NDim() != 2 || vertices.Shape()[1] != 3 {
		return nil, base.Errf("%s: the vertex table must be rank 2 with three columns, got shape %s",
			name, base.ShapeText(vertices.Shape()))
	}
	if tetrahedra.Dtype() != core.Int {
		return nil, base.Errf("%s: the tetrahedron table must hold integers, got %s", name, tetrahedra.Dtype())
	}
	if tetrahedra.NDim() != 2 || tetrahedra.Shape()[1] != 4 {
		return nil, base.Errf("%s: the tetrahedron table must be rank 2 with four columns, got shape %s",
			name, base.ShapeText(tetrahedra.Shape()))
	}
	n := vertices.Shape()[0]
	m := tetrahedra.Shape()[0]
	if n < 4 {
		return nil, base.Errf("%s: a mesh needs at least four vertices, got %d", name, n)
	}
	if m == 0 {
		// An empty tetrahedron table would surface deep in the sparse
		// factorisation on the zero rows of the free nodes, far from
		// the mesh that caused it.
		return nil, base.Errf("%s: the tetrahedron table must not be empty", name)
	}
	mesh := &TetraMesh3D{Vertices: make([]float64, 3*n), Tetrahedra: make([]int64, 4*m)}
	for i := range 3 * n {
		v := vertices.FloatAt(i)
		if math.IsNaN(v) || math.IsInf(v, 0) {
			return nil, base.Errf("%s: vertex coordinate %d is not finite", name, i)
		}
		mesh.Vertices[i] = v
	}
	for q := range 4 * m {
		idx := tetrahedra.RawInts()[q]
		if idx < 0 || idx >= int64(n) {
			return nil, base.Errf("%s: tetrahedron vertex index %d out of range for %d vertices", name, idx, n)
		}
		mesh.Tetrahedra[q] = idx
	}
	// Orientation and volume are checked where the caller can name the
	// tetrahedron and its vertices, not mid-assembly. Both messages
	// carry the coordinates, so a mis-ordered table can be fixed
	// without reopening a mesh debugger.
	for t := range m {
		a, b, c, d := int(mesh.Tetrahedra[4*t]), int(mesh.Tetrahedra[4*t+1]), int(mesh.Tetrahedra[4*t+2]), int(mesh.Tetrahedra[4*t+3])
		ax, ay, az := mesh.Vertices[3*a], mesh.Vertices[3*a+1], mesh.Vertices[3*a+2]
		bx, by, bz := mesh.Vertices[3*b], mesh.Vertices[3*b+1], mesh.Vertices[3*b+2]
		cx, cy, cz := mesh.Vertices[3*c], mesh.Vertices[3*c+1], mesh.Vertices[3*c+2]
		dx, dy, dz := mesh.Vertices[3*d], mesh.Vertices[3*d+1], mesh.Vertices[3*d+2]
		signed6 := signedTetraVolume(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz)
		at := func(v int) string {
			return fmt.Sprintf("(%g, %g, %g)", mesh.Vertices[3*v], mesh.Vertices[3*v+1], mesh.Vertices[3*v+2])
		}
		verts := fmt.Sprintf("vertices %d %s, %d %s, %d %s, %d %s", a, at(a), b, at(b), c, at(c), d, at(d))
		if signed6 == 0 {
			return nil, base.Errf("%s: tetrahedron %d is degenerate (zero volume), %s", name, t, verts)
		}
		if signed6 < 0 {
			return nil, base.Errf("%s: tetrahedron %d is inverted (signed volume %g), %s", name, t, signed6/6, verts)
		}
	}
	return mesh, nil
}

// Vertices3 returns the vertex count.
func (m *TetraMesh3D) Vertices3() int { return len(m.Vertices) / 3 }

// Tetrahedra4 returns the tetrahedron count.
func (m *TetraMesh3D) Tetrahedra4() int { return len(m.Tetrahedra) / 4 }

// signedTetraVolume returns six times the signed volume of the
// tetrahedron (a, b, c, d): positive for the orientation the mesh
// stores, negative when the last two vertices are swapped, zero when
// the four points are coplanar.
func signedTetraVolume(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz float64) float64 {
	u := [3]float64{bx - ax, by - ay, bz - az}
	v := [3]float64{cx - ax, cy - ay, cz - az}
	w := [3]float64{dx - ax, dy - ay, dz - az}
	cross := [3]float64{v[1]*w[2] - v[2]*w[1], v[2]*w[0] - v[0]*w[2], v[0]*w[1] - v[1]*w[0]}
	return u[0]*cross[0] + u[1]*cross[1] + u[2]*cross[2]
}

// BoundaryFaces returns the mesh's boundary faces as flat triples of
// vertex indices: a face belongs to the boundary when exactly one
// tetrahedron carries it. The triples are sorted lexicographically,
// so the result is a pure function of the mesh.
func (m *TetraMesh3D) BoundaryFaces() []int {
	count := make(map[[3]int]int, len(m.Tetrahedra))
	key := func(a, b, c int) [3]int {
		if a > b {
			a, b = b, a
		}
		if b > c {
			b, c = c, b
		}
		if a > b {
			a, b = b, a
		}
		return [3]int{a, b, c}
	}
	for t := 0; t < m.Tetrahedra4(); t++ {
		a, b, c, d := int(m.Tetrahedra[4*t]), int(m.Tetrahedra[4*t+1]), int(m.Tetrahedra[4*t+2]), int(m.Tetrahedra[4*t+3])
		count[key(a, b, c)]++
		count[key(a, b, d)]++
		count[key(a, c, d)]++
		count[key(b, c, d)]++
	}
	sets := make([][3]int, 0, len(count))
	for f, n := range count {
		if n == 1 {
			sets = append(sets, f)
		}
	}
	slices.SortFunc(sets, func(x, y [3]int) int {
		for k := range 3 {
			if x[k] != y[k] {
				return x[k] - y[k]
			}
		}
		return 0
	})
	faces := make([]int, 0, 3*len(sets))
	for _, f := range sets {
		faces = append(faces, f[0], f[1], f[2])
	}
	return faces
}

// BoxTetraMesh3D builds the structured tetrahedralisation of the
// axis-aligned box [x0, x0+width] × [y0, y0+height] × [z0, z0+depth]
// with m by n by p cells, six tetrahedra per cell (the Kuhn
// subdivision along the cell diagonal, oriented positively). m, n and
// p must all be positive. The subdivision is conforming across cell
// faces, which makes the mesher the first port of call for tests and
// for boxes in general.
func BoxTetraMesh3D(x0, y0, z0, width, height, depth float64, m, n, p int) (*TetraMesh3D, error) {
	const name = "BoxTetraMesh3D"
	if m <= 0 || n <= 0 || p <= 0 {
		return nil, base.Errf("%s: the cell counts must be positive, got %d by %d by %d", name, m, n, p)
	}
	// The same guard the triangle mesher applies: a non-finite extent
	// or origin would lay out vertices at NaN or Inf and only surface
	// mid-factorisation, far from the cause.
	if !(width > 0) || !(height > 0) || !(depth > 0) ||
		math.IsInf(width, 0) || math.IsInf(height, 0) || math.IsInf(depth, 0) ||
		math.IsNaN(x0) || math.IsInf(x0, 0) || math.IsNaN(y0) || math.IsInf(y0, 0) || math.IsNaN(z0) || math.IsInf(z0, 0) {
		return nil, base.Errf("%s: the extents must be finite and positive and the origin finite, got origin (%g, %g, %g), extents %g by %g by %g",
			name, x0, y0, z0, width, height, depth)
	}
	vertices := make([]float64, 3*(m+1)*(n+1)*(p+1))
	for k := range p + 1 {
		for j := range n + 1 {
			for i := range m + 1 {
				v := 3 * ((k*(n+1)+j)*(m+1) + i)
				vertices[v] = x0 + width*float64(i)/float64(m)
				vertices[v+1] = y0 + height*float64(j)/float64(n)
				vertices[v+2] = z0 + depth*float64(k)/float64(p)
			}
		}
	}
	at := func(i, j, k int) int64 { return int64((k*(n+1)+j)*(m+1) + i) }
	// The six Kuhn paths from one cell corner to the opposite one,
	// given as axis orders. An odd permutation reaches the far corner
	// with negative orientation, so its last two vertices swap.
	perms := [6][3]int{{0, 1, 2}, {0, 2, 1}, {1, 0, 2}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}}
	tetrahedra := make([]int64, 0, 6*m*n*p)
	for k := range p {
		for j := range n {
			for i := range m {
				for _, pm := range perms {
					// The path walks from the cell corner to the far
					// corner, each vertex one axis-step beyond the
					// previous one.
					ox := [4]int{i, i, i, i}
					oy := [4]int{j, j, j, j}
					oz := [4]int{k, k, k, k}
					for s := range 3 {
						ox[s+1], oy[s+1], oz[s+1] = ox[s], oy[s], oz[s]
						switch pm[s] {
						case 0:
							ox[s+1]++
						case 1:
							oy[s+1]++
						default:
							oz[s+1]++
						}
					}
					odd := 0
					for s1 := range 3 {
						for s2 := s1 + 1; s2 < 3; s2++ {
							if pm[s1] > pm[s2] {
								odd++
							}
						}
					}
					v := [4]int64{at(ox[0], oy[0], oz[0]), at(ox[1], oy[1], oz[1]), at(ox[2], oy[2], oz[2]), at(ox[3], oy[3], oz[3])}
					if odd%2 == 1 {
						v[2], v[3] = v[3], v[2]
					}
					tetrahedra = append(tetrahedra, v[0], v[1], v[2], v[3])
				}
			}
		}
	}
	return &TetraMesh3D{Vertices: vertices, Tetrahedra: tetrahedra}, nil
}

// tetraGradients returns the gradients of the four P1 basis functions
// on the tetrahedron (a, b, c, d) and its volume. The gradients are
// the columns of the inverse of the edge matrix whose rows are the
// vectors from d to a, b and c, which is the standard
// gradient-of-basis formula over the element's edge vectors.
func tetraGradients(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz float64) (g [4][3]float64, volume float64) {
	// Rows of the edge matrix relative to d.
	r0 := [3]float64{ax - dx, ay - dy, az - dz}
	r1 := [3]float64{bx - dx, by - dy, bz - dz}
	r2 := [3]float64{cx - dx, cy - dy, cz - dz}
	// Cofactors of the edge matrix; the inverse is their transpose
	// over the determinant, so column j of the inverse is row j of the
	// cofactor matrix over det.
	c00 := r1[1]*r2[2] - r1[2]*r2[1]
	c01 := -(r1[0]*r2[2] - r1[2]*r2[0])
	c02 := r1[0]*r2[1] - r1[1]*r2[0]
	c10 := -(r0[1]*r2[2] - r0[2]*r2[1])
	c11 := r0[0]*r2[2] - r0[2]*r2[0]
	c12 := -(r0[0]*r2[1] - r0[1]*r2[0])
	c20 := r0[1]*r1[2] - r0[2]*r1[1]
	c21 := -(r0[0]*r1[2] - r0[2]*r1[0])
	c22 := r0[0]*r1[1] - r0[1]*r1[0]
	det := r0[0]*c00 + r0[1]*c01 + r0[2]*c02
	g[0] = [3]float64{c00 / det, c01 / det, c02 / det}
	g[1] = [3]float64{c10 / det, c11 / det, c12 / det}
	g[2] = [3]float64{c20 / det, c21 / det, c22 / det}
	for i := range 3 {
		for k := range 3 {
			g[3][k] -= g[i][k]
		}
	}
	volume = math.Abs(signedTetraVolume(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz)) / 6
	return g, volume
}

// tetraStiffness returns the P1 stiffness matrix of one tetrahedron:
// K[i][j] = κ·V·(∇λᵢ·∇λⱼ), the gradient-of-basis formula integrated
// over the element, where the gradients are constant on a linear
// element.
func tetraStiffness(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, kappa float64) [4][4]float64 {
	g, volume := tetraGradients(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz)
	var k [4][4]float64
	for i := range 4 {
		for j := range 4 {
			k[i][j] = kappa * volume * (g[i][0]*g[j][0] + g[i][1]*g[j][1] + g[i][2]*g[j][2])
		}
	}
	return k
}

// FEMPoisson3DOptions carries the data SolvePoissonFEM3D needs beside
// the mesh and the source: the conductivity, the prescribed boundary
// values, and the optional flux boundary.
type FEMPoisson3DOptions struct {
	// Kappa is the constant conductivity when KappaFunc is nil. It
	// must be positive.
	Kappa float64
	// KappaFunc, when set, gives the conductivity at a point. It is
	// evaluated at the tetrahedron centroids and must be positive
	// there for every element; a non-positive value names the element.
	KappaFunc func(x, y, z float64) float64
	// DirichletNodes lists the vertices with prescribed values and
	// DirichletValues the values in the same order. The nodes leave
	// the system with their rows and columns; at least one is
	// required, because a purely Neumann problem has no unique
	// solution.
	DirichletNodes  []int
	DirichletValues []float64
	// NeumannFaces lists boundary faces as flat triples of vertex
	// indices and NeumannFlux gives the flux κ∂u/∂n along each face's
	// outward normal: each face's integral is built from the degree-2
	// edge-midpoint rule, a third of area·flux at each edge midpoint
	// shared by that edge's two vertices. A nil flux means zero.
	NeumannFaces []int
	NeumannFlux  func(x, y, z float64) float64
	// Ordering selects the fill-reducing permutation for the sparse
	// Cholesky factorisation. The zero value is the natural order;
	// meshes usually want SparseOrderingReverseCuthillMcKee.
	Ordering linalg.SparseOrdering
}

// SolvePoissonFEM3D solves −∇·(κ∇u) = f on the tetrahedral mesh with
// piecewise-linear elements: the stiffness matrix is assembled per
// tetrahedron (the conductivity evaluated at the centroids when it
// varies), the load is integrated per tetrahedron with the 3×3×3
// collapsed Gauss rule (exact through degree 5; the centroid lump
// does not hold the O(h²) rate on the structured Kuhn mesh), Neumann
// fluxes are integrated on their boundary faces with the degree-2
// edge-midpoint rule, and Dirichlet values are eliminated by lifting.
// f may be nil for the homogeneous equation. The error contract
// mirrors SolvePoissonFEM2D.
func SolvePoissonFEM3D(mesh *TetraMesh3D, f func(x, y, z float64) float64, opts FEMPoisson3DOptions) (*core.Array, error) {
	const name = "SolvePoissonFEM3D"
	if mesh == nil {
		return nil, base.Errf("%s: the mesh must not be nil", name)
	}
	// The same conductivity gate as the two-dimensional solve: with
	// KappaFunc nil the constant is the value used, so it must be
	// positive and finite; with the field set the constant is a
	// placeholder, but a non-finite one is still refused.
	if opts.KappaFunc == nil {
		if !(opts.Kappa > 0) || math.IsInf(opts.Kappa, 0) {
			return nil, base.Errf("%s: the conductivity must be positive, got %g", name, opts.Kappa)
		}
	} else if math.IsNaN(opts.Kappa) || math.IsInf(opts.Kappa, 0) {
		return nil, base.Errf("%s: the conductivity must be positive, got %g", name, opts.Kappa)
	}
	if len(opts.DirichletNodes) != len(opts.DirichletValues) {
		return nil, base.Errf("%s: %d Dirichlet nodes but %d values", name, len(opts.DirichletNodes), len(opts.DirichletValues))
	}
	if len(opts.DirichletNodes) == 0 {
		return nil, base.Errf("%s: a purely Neumann problem has no unique solution; prescribe at least one Dirichlet value", name)
	}
	n := mesh.Vertices3()
	// The Dirichlet nodes as a dense marker with their prescribed
	// values, exactly as the two-dimensional solve carries them: the
	// lifting and the unit rows each visit every assembled entry, and a
	// marker answers those visits in constant time where a set of nodes
	// answered with a hash. A node listed twice keeps its last value
	// and appears once, as it did in the set; the appended order does
	// not reach the assembled system, whose coordinate entries the
	// sparse conversion sorts and merges by coordinate.
	dirichletMark := make([]bool, n)
	dirichletVal := make([]float64, n)
	dirichletNodes := make([]int, 0, len(opts.DirichletNodes))
	for p, d := range opts.DirichletNodes {
		if d < 0 || d >= n {
			return nil, base.Errf("%s: Dirichlet node %d out of range for %d vertices", name, d, n)
		}
		v := opts.DirichletValues[p]
		if math.IsNaN(v) || math.IsInf(v, 0) {
			return nil, base.Errf("%s: Dirichlet value at node %d is not finite", name, d)
		}
		if !dirichletMark[d] {
			dirichletNodes = append(dirichletNodes, d)
		}
		dirichletMark[d] = true
		dirichletVal[d] = v
	}
	if len(opts.NeumannFaces)%3 != 0 {
		return nil, base.Errf("%s: %d Neumann face indices, want triples", name, len(opts.NeumannFaces))
	}
	for p := 0; p < len(opts.NeumannFaces); p += 3 {
		for _, v := range opts.NeumannFaces[p : p+3] {
			if v < 0 || v >= n {
				return nil, base.Errf("%s: Neumann face [%d %d %d] holds the out-of-range vertex %d",
					name, opts.NeumannFaces[p], opts.NeumannFaces[p+1], opts.NeumannFaces[p+2], v)
			}
		}
		if opts.NeumannFaces[p] == opts.NeumannFaces[p+1] ||
			opts.NeumannFaces[p] == opts.NeumannFaces[p+2] ||
			opts.NeumannFaces[p+1] == opts.NeumannFaces[p+2] {
			return nil, base.Errf("%s: Neumann face [%d %d %d] repeats a vertex",
				name, opts.NeumannFaces[p], opts.NeumannFaces[p+1], opts.NeumannFaces[p+2])
		}
	}
	// Assembly: sixteen entries per tetrahedron, symmetric by
	// construction, with the conductivity evaluated at the centroid
	// when it varies.
	entries := make([]float64, 0, 16*mesh.Tetrahedra4())
	rows := make([]int, 0, 16*mesh.Tetrahedra4())
	cols := make([]int, 0, 16*mesh.Tetrahedra4())
	load := make([]float64, n)
	// The collapsed Gauss rule's abscissae and weights are constants of
	// the scheme: built once here, not per tetrahedron.
	gl := [3]float64{(1 - math.Sqrt(3.0/5)) / 2, 0.5, (1 + math.Sqrt(3.0/5)) / 2}
	gw := [3]float64{5.0 / 18, 4.0 / 9, 5.0 / 18}
	for t := 0; t < mesh.Tetrahedra4(); t++ {
		a, b, c, d := int(mesh.Tetrahedra[4*t]), int(mesh.Tetrahedra[4*t+1]), int(mesh.Tetrahedra[4*t+2]), int(mesh.Tetrahedra[4*t+3])
		ax, ay, az := mesh.Vertices[3*a], mesh.Vertices[3*a+1], mesh.Vertices[3*a+2]
		bx, by, bz := mesh.Vertices[3*b], mesh.Vertices[3*b+1], mesh.Vertices[3*b+2]
		cx, cy, cz := mesh.Vertices[3*c], mesh.Vertices[3*c+1], mesh.Vertices[3*c+2]
		dx, dy, dz := mesh.Vertices[3*d], mesh.Vertices[3*d+1], mesh.Vertices[3*d+2]
		volume := math.Abs(signedTetraVolume(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz)) / 6
		if volume == 0 {
			return nil, base.Errf("%s: tetrahedron %d is degenerate (zero volume)", name, t)
		}
		kappa := opts.Kappa
		if opts.KappaFunc != nil {
			kappa = opts.KappaFunc((ax+bx+cx+dx)/4, (ay+by+cy+dy)/4, (az+bz+cz+dz)/4)
			if !(kappa > 0) || math.IsNaN(kappa) || math.IsInf(kappa, 0) {
				return nil, base.Errf("%s: the conductivity at tetrahedron %d is %g, want positive", name, t, kappa)
			}
		}
		k := tetraStiffness(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, kappa)
		nodes := [4]int{a, b, c, d}
		for i := range 4 {
			for j := range 4 {
				rows = append(rows, nodes[i])
				cols = append(cols, nodes[j])
				entries = append(entries, k[i][j])
			}
		}
		// The load on this element, integrated with the 3×3×3
		// collapsed Gauss rule: λ weights follow the Duffy collapse
		// toward vertex a, and the Jacobian of the map from the unit
		// cube is (1−r)²(1−s)·6V.
		if f != nil {
			for ir := range 3 {
				for is := range 3 {
					for it := range 3 {
						r, s, t := gl[ir], gl[is], gl[it]
						la := (1 - r) * (1 - s) * (1 - t)
						lb := (1 - r) * (1 - s) * t
						lc := (1 - r) * s
						ld := r
						x := la*ax + lb*bx + lc*cx + ld*dx
						y := la*ay + lb*by + lc*cy + ld*dy
						z := la*az + lb*bz + lc*cz + ld*dz
						w := gw[ir] * gw[is] * gw[it] * (1 - r) * (1 - r) * (1 - s) * 6 * volume
						fv := f(x, y, z)
						// A non-finite source value would flow into the
						// load and the solve would publish an all-NaN
						// solution with a nil error, the breach every
						// other integrator here refuses up front.
						if math.IsNaN(fv) || math.IsInf(fv, 0) {
							return nil, base.Errf("%s: the source returned the non-finite value %g at tetrahedron %d", name, fv, t)
						}
						load[a] += w * fv * la
						load[b] += w * fv * lb
						load[c] += w * fv * lc
						load[d] += w * fv * ld
					}
				}
			}
		}
	}
	// Neumann fluxes: the degree-2 edge-midpoint rule on every listed
	// face, a third of area·flux at each edge midpoint into that
	// edge's two vertices.
	if len(opts.NeumannFaces) > 0 && opts.NeumannFlux != nil {
		for p := 0; p < len(opts.NeumannFaces); p += 3 {
			a, b, c := opts.NeumannFaces[p], opts.NeumannFaces[p+1], opts.NeumannFaces[p+2]
			ax, ay, az := mesh.Vertices[3*a], mesh.Vertices[3*a+1], mesh.Vertices[3*a+2]
			bx, by, bz := mesh.Vertices[3*b], mesh.Vertices[3*b+1], mesh.Vertices[3*b+2]
			cx, cy, cz := mesh.Vertices[3*c], mesh.Vertices[3*c+1], mesh.Vertices[3*c+2]
			u := [3]float64{bx - ax, by - ay, bz - az}
			v := [3]float64{cx - ax, cy - ay, cz - az}
			cross := [3]float64{u[1]*v[2] - u[2]*v[1], u[2]*v[0] - u[0]*v[2], u[0]*v[1] - u[1]*v[0]}
			area := math.Sqrt(cross[0]*cross[0]+cross[1]*cross[1]+cross[2]*cross[2]) / 2
			w := area / 3
			// A non-finite flux lands in the load like a non-finite
			// source, so the same refusal answers it, naming the face.
			fab := w * opts.NeumannFlux((ax+bx)/2, (ay+by)/2, (az+bz)/2)
			fbc := w * opts.NeumannFlux((bx+cx)/2, (by+cy)/2, (bz+cz)/2)
			fca := w * opts.NeumannFlux((cx+ax)/2, (cy+ay)/2, (cz+az)/2)
			for _, fv := range []float64{fab, fbc, fca} {
				if math.IsNaN(fv) || math.IsInf(fv, 0) {
					return nil, base.Errf("%s: the Neumann flux returned a non-finite value on face [%d %d %d]", name, a, b, c)
				}
			}
			load[a] += fab/2 + fca/2
			load[b] += fab/2 + fbc/2
			load[c] += fbc/2 + fca/2
		}
	}
	// Dirichlet lifting: the known boundary values move to the right
	// hand side, then their rows and columns leave the system as unit
	// rows, exactly as in the two-dimensional solve.
	for p, i := range rows {
		if j := cols[p]; dirichletMark[j] {
			load[i] -= entries[p] * dirichletVal[j]
		}
	}
	keptRows := make([]int64, 0, len(rows))
	keptCols := make([]int64, 0, len(rows))
	keptVals := make([]float64, 0, len(rows))
	for p := range rows {
		i, j := rows[p], cols[p]
		if dirichletMark[i] || dirichletMark[j] {
			continue
		}
		keptRows = append(keptRows, int64(i))
		keptCols = append(keptCols, int64(j))
		keptVals = append(keptVals, entries[p])
	}
	for _, d := range dirichletNodes {
		keptRows = append(keptRows, int64(d))
		keptCols = append(keptCols, int64(d))
		keptVals = append(keptVals, 1)
		load[d] = dirichletVal[d]
	}
	indices, err := core.FromInts(pairInts(keptRows, keptCols), len(keptVals), 2)
	if err != nil {
		return nil, base.Errf("%s: %w", name, err)
	}
	coo, err := core.NewSparseCOO(indices, fromSlice(keptVals, len(keptVals)), []int{n, n})
	if err != nil {
		return nil, base.Errf("%s: %w", name, err)
	}
	factor, err := linalg.NewSparseCholesky(coo, opts.Ordering)
	if err != nil {
		return nil, base.Errf("%s: %w", name, err)
	}
	rhs := core.New(core.Float, []int{n}...)
	copy(rhs.RawFloats(), load)
	return factor.Solve(rhs)
}