506 lines
20 KiB
Go
506 lines
20 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package integrate
|
||
|
||
import (
|
||
"math"
|
||
"testing"
|
||
|
||
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
||
linalg "sourcedock.dev/petrbalvin/tensor/linalg"
|
||
)
|
||
|
||
// gridMesh builds the structured triangulation of the unit square
|
||
// with m cells per side, two triangles per cell, and returns the mesh
|
||
// plus the list of boundary vertices in the order (bottom row, top
|
||
// row, left column, right column), duplicates removed.
|
||
func gridMesh(t *testing.T, m int) (*TriangleMesh2D, []int) {
|
||
t.Helper()
|
||
mesh, err := GridTriangleMesh2D(0, 0, 1, 1, m, m)
|
||
if err != nil {
|
||
t.Fatalf("GridTriangleMesh2D: %v", err)
|
||
}
|
||
boundary := make([]int, 0, 4*m)
|
||
for i := range m + 1 {
|
||
boundary = append(boundary, i, m*(m+1)+i)
|
||
}
|
||
for j := 1; j < m; j++ {
|
||
boundary = append(boundary, j*(m+1), j*(m+1)+m)
|
||
}
|
||
return mesh, boundary
|
||
}
|
||
|
||
func TestSolvePoissonFEM2DConvergence(t *testing.T) {
|
||
// The manufactured solution u = sin(πx)·sin(πy) on the unit
|
||
// square drives f = 2π²·u; with the boundary lifted the P1 error
|
||
// must halve twice when the mesh is refined, the O(h²) the
|
||
// piecewise-linear theory promises.
|
||
solution := func(x, y float64) float64 { return math.Sin(math.Pi*x) * math.Sin(math.Pi*y) }
|
||
source := func(x, y float64) float64 { return 2 * math.Pi * math.Pi * solution(x, y) }
|
||
previous := 0.0
|
||
for _, m := range []int{8, 16, 32} {
|
||
mesh, boundary := gridMesh(t, m)
|
||
values := make([]float64, len(boundary))
|
||
for p, node := range boundary {
|
||
values[p] = solution(mesh.Vertices[2*node], mesh.Vertices[2*node+1])
|
||
}
|
||
u, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: values})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D(m=%d): %v", m, err)
|
||
}
|
||
worst := 0.0
|
||
for i := range mesh.Vertices2() {
|
||
if d := math.Abs(u.FloatAt(i) - solution(mesh.Vertices[2*i], mesh.Vertices[2*i+1])); d > worst {
|
||
worst = d
|
||
}
|
||
}
|
||
t.Logf("m=%2d: max nodal error %.3g", m, worst)
|
||
if previous > 0 && previous/worst < 2.5 {
|
||
t.Fatalf("m=%d: refinement ratio %.2f, want the O(h²) rate (previous %.3g, now %.3g)",
|
||
m, previous/worst, previous, worst)
|
||
}
|
||
if m == 32 && worst > 2e-3 {
|
||
t.Fatalf("m=32: error %.3g too large for the asymptotic range", worst)
|
||
}
|
||
previous = worst
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DLinearExactness is the patch test the P1
|
||
// elements must pass without compromise: a linear field lies in the
|
||
// approximation space, so with f = 0 and the boundary lifted the
|
||
// interior solution must equal the field to machine precision.
|
||
func TestSolvePoissonFEM2DLinearExactness(t *testing.T) {
|
||
mesh, boundary := gridMesh(t, 12)
|
||
field := func(x, y float64) float64 { return 1 + 2*x - 3*y }
|
||
values := make([]float64, len(boundary))
|
||
for p, node := range boundary {
|
||
values[p] = field(mesh.Vertices[2*node], mesh.Vertices[2*node+1])
|
||
}
|
||
u, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: values})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D: %v", err)
|
||
}
|
||
worst := 0.0
|
||
for i := range mesh.Vertices2() {
|
||
if d := math.Abs(u.FloatAt(i) - field(mesh.Vertices[2*i], mesh.Vertices[2*i+1])); d > worst {
|
||
worst = d
|
||
}
|
||
}
|
||
if worst > 1e-12 {
|
||
t.Fatalf("linear patch test error %.3g, want machine precision", worst)
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DNeumannNatural pins the natural boundary: a
|
||
// constant field with f = 0 satisfies the homogeneous Neumann
|
||
// condition everywhere, so pinning the constant at a single vertex
|
||
// must reproduce it across the whole mesh.
|
||
func TestSolvePoissonFEM2DNeumannNatural(t *testing.T) {
|
||
mesh, _ := gridMesh(t, 10)
|
||
u, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{0}, DirichletValues: []float64{5}})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D: %v", err)
|
||
}
|
||
for i := range mesh.Vertices2() {
|
||
if math.Abs(u.FloatAt(i)-5) > 1e-10 {
|
||
t.Fatalf("node %d: solution %.12g, want the constant 5", i, u.FloatAt(i))
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DOrderings runs the manufactured-solution solve
|
||
// under every ordering the factor offers: the ordering changes the
|
||
// fill, never the answer.
|
||
func TestSolvePoissonFEM2DOrderings(t *testing.T) {
|
||
solution := func(x, y float64) float64 { return math.Sin(math.Pi*x) * math.Sin(math.Pi*y) }
|
||
mesh, boundary := gridMesh(t, 10)
|
||
values := make([]float64, len(boundary))
|
||
for p, node := range boundary {
|
||
values[p] = solution(mesh.Vertices[2*node], mesh.Vertices[2*node+1])
|
||
}
|
||
source := func(x, y float64) float64 { return 2 * math.Pi * math.Pi * solution(x, y) }
|
||
reference, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: values})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D(natural): %v", err)
|
||
}
|
||
for _, ordering := range []linalg.SparseOrdering{
|
||
linalg.SparseOrderingReverseCuthillMcKee,
|
||
linalg.SparseOrderingMinimumDegree,
|
||
} {
|
||
u, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: values, Ordering: ordering})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D(%d): %v", ordering, err)
|
||
}
|
||
for i := range mesh.Vertices2() {
|
||
if math.Abs(u.FloatAt(i)-reference.FloatAt(i)) > 1e-9 {
|
||
t.Fatalf("ordering %d: node %d differs from the natural run", ordering, i)
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
func TestSolvePoissonFEM2DRefusals(t *testing.T) {
|
||
mesh, boundary := gridMesh(t, 5)
|
||
// Degenerate triangle: three collinear vertices.
|
||
if _, err := NewTriangleMesh2D(
|
||
floatsToArrayFEM(t, []float64{0, 0, 1, 0, 2, 0}, 3, 2),
|
||
intsToArrayFEM(t, []int64{0, 1, 2}, 1, 3)); err == nil || !stringsContains(err, "degenerate") {
|
||
t.Fatalf("a degenerate triangle: %v", err)
|
||
}
|
||
// Triangle index out of range.
|
||
if _, err := NewTriangleMesh2D(
|
||
floatsToArrayFEM(t, []float64{0, 0, 1, 0, 0, 1}, 3, 2),
|
||
intsToArrayFEM(t, []int64{0, 1, 3}, 1, 3)); err == nil || !stringsContains(err, "out of range") {
|
||
t.Fatalf("out of range index: %v", err)
|
||
}
|
||
// A float triangle table: the triangles must be integer indices.
|
||
if _, err := NewTriangleMesh2D(
|
||
floatsToArrayFEM(t, []float64{0, 0, 1, 0, 0, 1}, 3, 2),
|
||
floatsToArrayFEM(t, []float64{0, 1, 2}, 1, 3)); err == nil || !stringsContains(err, "integers") {
|
||
t.Fatalf("a float triangle table: %v", err)
|
||
}
|
||
// Dirichlet node out of range and a length mismatch.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{99}, DirichletValues: []float64{1}}); err == nil || !stringsContains(err, "out of range") {
|
||
t.Fatalf("an out of range Dirichlet node: %v", err)
|
||
}
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{0, 1}, DirichletValues: []float64{1}}); err == nil {
|
||
t.Fatal("a Dirichlet length mismatch was accepted")
|
||
}
|
||
// Non-positive conductivity.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 0, DirichletNodes: boundary, DirichletValues: make([]float64, len(boundary))}); err == nil || !stringsContains(err, "positive") {
|
||
t.Fatalf("zero conductivity: %v", err)
|
||
}
|
||
// Non-finite Dirichlet value.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{0}, DirichletValues: []float64{math.NaN()}}); err == nil || !stringsContains(err, "finite") {
|
||
t.Fatalf("a NaN Dirichlet value: %v", err)
|
||
}
|
||
// A NaN vertex coordinate in the mesh table.
|
||
if _, err := NewTriangleMesh2D(
|
||
floatsToArrayFEM(t, []float64{math.NaN(), 0, 1, 0, 0, 1}, 3, 2),
|
||
intsToArrayFEM(t, []int64{0, 1, 2}, 1, 3)); err == nil || !stringsContains(err, "not finite") {
|
||
t.Fatalf("a NaN vertex coordinate: %v", err)
|
||
}
|
||
// An odd number of Neumann edge indices: no complete pairs.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{0}, DirichletValues: []float64{0}, NeumannEdges: []int{0, 1, 2}}); err == nil || !stringsContains(err, "pairs") {
|
||
t.Fatalf("an odd Neumann edge count: %v", err)
|
||
}
|
||
// A degenerate Neumann edge a == b.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{0}, DirichletValues: []float64{0}, NeumannEdges: []int{3, 3}}); err == nil || !stringsContains(err, "valid vertex pair") {
|
||
t.Fatalf("a degenerate Neumann edge: %v", err)
|
||
}
|
||
// A Neumann edge index out of range.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: []int{0}, DirichletValues: []float64{0}, NeumannEdges: []int{0, 999}}); err == nil || !stringsContains(err, "valid vertex pair") {
|
||
t.Fatalf("an out of range Neumann edge: %v", err)
|
||
}
|
||
// A KappaFunc returning a non-positive conductivity names the
|
||
// triangle instead of assembling a singular stiffness matrix.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{
|
||
KappaFunc: func(float64, float64) float64 { return -1 },
|
||
DirichletNodes: []int{0},
|
||
DirichletValues: []float64{0},
|
||
}); err == nil || !stringsContains(err, "positive") {
|
||
t.Fatalf("a non-positive KappaFunc value: %v", err)
|
||
}
|
||
// A KappaFunc returning an infinite conductivity names the triangle
|
||
// the way the constant field's gate names itself, instead of
|
||
// surfacing as a factorisation failure far from the cause.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{
|
||
KappaFunc: func(float64, float64) float64 { return math.Inf(1) },
|
||
DirichletNodes: []int{0},
|
||
DirichletValues: []float64{0},
|
||
}); err == nil || !stringsContains(err, "positive") {
|
||
t.Fatalf("an infinite KappaFunc value: %v", err)
|
||
}
|
||
// A non-finite source value refuses the solve: it used to land in
|
||
// the load and publish an all-NaN solution with a nil error.
|
||
if _, err := SolvePoissonFEM2D(mesh, func(x, y float64) float64 { return math.NaN() },
|
||
FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: make([]float64, len(boundary))}); err == nil || !stringsContains(err, "non-finite") {
|
||
t.Fatalf("a NaN source value: %v", err)
|
||
}
|
||
// A non-finite Neumann flux refuses the solve the same way.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{
|
||
Kappa: 1,
|
||
DirichletNodes: boundary,
|
||
DirichletValues: make([]float64, len(boundary)),
|
||
NeumannEdges: []int{0, mesh.Vertices2() - 1},
|
||
NeumannFlux: func(x, y float64) float64 { return math.Inf(1) },
|
||
}); err == nil || !stringsContains(err, "non-finite") {
|
||
t.Fatalf("an infinite Neumann flux: %v", err)
|
||
}
|
||
// An ordering that does not exist.
|
||
if _, err := SolvePoissonFEM2D(mesh, nil, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: make([]float64, len(boundary)), Ordering: linalg.SparseOrdering(7)}); err == nil {
|
||
t.Fatal("an unknown ordering was accepted")
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DIsDeterministic solves the same problem twice
|
||
// and requires bit-identical nodal values, the contract every Tensor
|
||
// entry point carries.
|
||
func TestSolvePoissonFEM2DIsDeterministic(t *testing.T) {
|
||
solution := func(x, y float64) float64 { return math.Sin(math.Pi*x) * math.Sin(math.Pi*y) }
|
||
mesh, boundary := gridMesh(t, 10)
|
||
values := make([]float64, len(boundary))
|
||
for p, node := range boundary {
|
||
values[p] = solution(mesh.Vertices[2*node], mesh.Vertices[2*node+1])
|
||
}
|
||
source := func(x, y float64) float64 { return 2 * math.Pi * math.Pi * solution(x, y) }
|
||
u1, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: values})
|
||
if err != nil {
|
||
t.Fatalf("first solve: %v", err)
|
||
}
|
||
u2, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: values})
|
||
if err != nil {
|
||
t.Fatalf("second solve: %v", err)
|
||
}
|
||
for i := range mesh.Vertices2() {
|
||
if u1.FloatAt(i) != u2.FloatAt(i) {
|
||
t.Fatalf("node %d differs: %.17g vs %.17g", i, u1.FloatAt(i), u2.FloatAt(i))
|
||
}
|
||
}
|
||
}
|
||
|
||
func stringsContains(err error, fragment string) bool {
|
||
return err != nil && len(err.Error()) >= len(fragment) && indexOf(err.Error(), fragment) >= 0
|
||
}
|
||
|
||
func indexOf(s, fragment string) int {
|
||
for i := 0; i+len(fragment) <= len(s); i++ {
|
||
if s[i:i+len(fragment)] == fragment {
|
||
return i
|
||
}
|
||
}
|
||
return -1
|
||
}
|
||
|
||
func floatsToArrayFEM(t *testing.T, vals []float64, shape ...int) *core.Array {
|
||
t.Helper()
|
||
a, err := core.FromFloats(vals, shape...)
|
||
if err != nil {
|
||
t.Fatalf("FromFloats: %v", err)
|
||
}
|
||
return a
|
||
}
|
||
|
||
func intsToArrayFEM(t *testing.T, vals []int64, shape ...int) *core.Array {
|
||
t.Helper()
|
||
a, err := core.FromInts(vals, shape...)
|
||
if err != nil {
|
||
t.Fatalf("FromInts: %v", err)
|
||
}
|
||
return a
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DNeumannFlux pins the boundary-edge integrals:
|
||
// u = (x²+y²)/2 has −Δu = −2 and the flux κ∂u/∂n = 1 along the right
|
||
// and top edges' outward normals (0 along the bottom and left), so
|
||
// prescribing those fluxes with a single pinned vertex must
|
||
// reproduce the quadratic field. The midpoint edge rule is
|
||
// first-order consistent, so the error must halve with the mesh.
|
||
func TestSolvePoissonFEM2DNeumannFlux(t *testing.T) {
|
||
field := func(x, y float64) float64 { return (x*x + y*y) / 2 }
|
||
previous := 0.0
|
||
for _, m := range []int{10, 20} {
|
||
mesh, err := GridTriangleMesh2D(0, 0, 1, 1, m, m)
|
||
if err != nil {
|
||
t.Fatalf("GridTriangleMesh2D: %v", err)
|
||
}
|
||
// Boundary edges: pairs of neighbouring boundary vertices.
|
||
var edges []int
|
||
at := func(i, j int) int { return j*(m+1) + i }
|
||
for j := range m {
|
||
edges = append(edges, at(j, 0), at(j+1, 0)) // bottom: flux 0
|
||
edges = append(edges, at(j, m), at(j+1, m)) // top: flux 1
|
||
edges = append(edges, at(m, j), at(m, j+1)) // right: flux 1
|
||
edges = append(edges, at(0, j), at(0, j+1)) // left: flux 0
|
||
}
|
||
flux := func(x, y float64) float64 {
|
||
if x == 1 || y == 1 {
|
||
return 1
|
||
}
|
||
return 0
|
||
}
|
||
u, err := SolvePoissonFEM2D(mesh, func(float64, float64) float64 { return -2 },
|
||
FEMPoissonOptions{
|
||
Kappa: 1,
|
||
DirichletNodes: []int{at(0, 0)},
|
||
DirichletValues: []float64{0},
|
||
NeumannEdges: edges,
|
||
NeumannFlux: flux,
|
||
})
|
||
if err != nil {
|
||
t.Fatalf("m=%d: %v", m, err)
|
||
}
|
||
worst := 0.0
|
||
for i := range mesh.Vertices2() {
|
||
x := mesh.Vertices[2*i]
|
||
y := mesh.Vertices[2*i+1]
|
||
if d := math.Abs(u.FloatAt(i) - field(x, y)); d > worst {
|
||
worst = d
|
||
}
|
||
}
|
||
t.Logf("m=%2d: max nodal error %.3g", m, worst)
|
||
if previous > 0 && previous/worst < 1.4 {
|
||
t.Fatalf("m=%d: refinement ratio %.2f, want the first-order flux rate", m, previous/worst)
|
||
}
|
||
if m == 20 && worst > 5e-3 {
|
||
t.Fatalf("m=20: error %.3g too large", worst)
|
||
}
|
||
previous = worst
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DVariableKappa runs the manufactured solution
|
||
// with a spatially varying conductivity evaluated at the element
|
||
// centroids: f must carry the analytic divergence terms, and the
|
||
// P1 convergence rate must survive the varying coefficient.
|
||
func TestSolvePoissonFEM2DVariableKappa(t *testing.T) {
|
||
sin, cos := math.Pi, math.Pi
|
||
u := func(x, y float64) float64 { return math.Sin(sin*x) * math.Sin(sin*y) }
|
||
kappaF := func(x, y float64) float64 { return 1 + x*y }
|
||
ux := func(x, y float64) float64 { return cos * math.Cos(cos*x) * math.Sin(cos*y) }
|
||
uy := func(x, y float64) float64 { return cos * math.Sin(cos*x) * math.Cos(cos*y) }
|
||
lap := func(x, y float64) float64 { return -2 * math.Pi * math.Pi * u(x, y) }
|
||
source := func(x, y float64) float64 {
|
||
k := kappaF(x, y)
|
||
return -(y*ux(x, y) + x*uy(x, y) + k*lap(x, y))
|
||
}
|
||
previous := 0.0
|
||
for _, m := range []int{8, 16, 32} {
|
||
mesh, boundary := gridMesh(t, m)
|
||
values := make([]float64, len(boundary))
|
||
for p, node := range boundary {
|
||
values[p] = u(mesh.Vertices[2*node], mesh.Vertices[2*node+1])
|
||
}
|
||
uk, err := SolvePoissonFEM2D(mesh, source,
|
||
FEMPoissonOptions{KappaFunc: kappaF, DirichletNodes: boundary, DirichletValues: values})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D(m=%d): %v", m, err)
|
||
}
|
||
worst := 0.0
|
||
for i := range mesh.Vertices2() {
|
||
if d := math.Abs(uk.FloatAt(i) - u(mesh.Vertices[2*i], mesh.Vertices[2*i+1])); d > worst {
|
||
worst = d
|
||
}
|
||
}
|
||
t.Logf("m=%2d: max nodal error %.3g", m, worst)
|
||
if previous > 0 && previous/worst < 2.5 {
|
||
t.Fatalf("m=%d: refinement ratio %.2f, want the O(h²) rate", m, previous/worst)
|
||
}
|
||
previous = worst
|
||
}
|
||
}
|
||
|
||
// TestTriangleMesh2DBoundaryEdges pins the boundary-edge detection:
|
||
// the m by n grid carries exactly 2(m+n) boundary edges, every one of
|
||
// them with both endpoints on the boundary vertex ring.
|
||
func TestTriangleMesh2DBoundaryEdges(t *testing.T) {
|
||
mesh, err := GridTriangleMesh2D(0, 0, 1, 1, 5, 3)
|
||
if err != nil {
|
||
t.Fatalf("GridTriangleMesh2D: %v", err)
|
||
}
|
||
edges := mesh.BoundaryEdges()
|
||
if len(edges) != 2*2*(5+3) {
|
||
t.Fatalf("boundary edge count %d, want %d", len(edges), 2*(5+3))
|
||
}
|
||
onBoundary := func(v int) bool {
|
||
i := v % 6
|
||
j := v / 6
|
||
return i == 0 || i == 5 || j == 0 || j == 3
|
||
}
|
||
for p := 0; p < len(edges); p += 2 {
|
||
if !onBoundary(edges[p]) || !onBoundary(edges[p+1]) {
|
||
t.Fatalf("edge [%d,%d] is not on the boundary", edges[p], edges[p+1])
|
||
}
|
||
}
|
||
// The generator's vertex positions are exact.
|
||
mesh2, err := GridTriangleMesh2D(-1, 2, 2, 4, 2, 2)
|
||
if err != nil {
|
||
t.Fatalf("GridTriangleMesh2D: %v", err)
|
||
}
|
||
if mesh2.Vertices[0] != -1 || mesh2.Vertices[1] != 2 {
|
||
t.Fatalf("vertex 0 = [%g %g], want [-1 2]", mesh2.Vertices[0], mesh2.Vertices[1])
|
||
}
|
||
if mesh2.Vertices[2*(2*3+2)] != 1 || mesh2.Vertices[2*(2*3+2)+1] != 6 {
|
||
t.Fatalf("vertex (2,2) = [%g %g], want [1 6]",
|
||
mesh2.Vertices[2*(2*3+2)], mesh2.Vertices[2*(2*3+2)+1])
|
||
}
|
||
if _, err := GridTriangleMesh2D(0, 0, 1, 1, 0, 3); err == nil {
|
||
t.Fatal("a zero cell count was accepted")
|
||
}
|
||
if _, err := GridTriangleMesh2D(0, 0, -1, 1, 2, 2); err == nil {
|
||
t.Fatal("a negative extent was accepted")
|
||
}
|
||
}
|
||
|
||
// TestTriangleMesh2DBoundaryEdgesAreMeshEdges pins the pair contract of
|
||
// BoundaryEdges: every returned pair must be an edge the mesh actually
|
||
// carries, and the pairs must come out sorted, which a sort of the flat
|
||
// index list cannot deliver (it interleaves unrelated endpoints).
|
||
func TestTriangleMesh2DBoundaryEdgesAreMeshEdges(t *testing.T) {
|
||
mesh, err := GridTriangleMesh2D(0, 0, 1, 1, 1, 1)
|
||
if err != nil {
|
||
t.Fatalf("GridTriangleMesh2D: %v", err)
|
||
}
|
||
edges := mesh.BoundaryEdges()
|
||
// The square's four sides: (0,1), (0,2), (1,3), (2,3) in sorted
|
||
// pair order.
|
||
want := []int{0, 1, 0, 2, 1, 3, 2, 3}
|
||
if len(edges) != len(want) {
|
||
t.Fatalf("boundary edge count %d, want %d", len(edges)/2, len(want)/2)
|
||
}
|
||
for p := 0; p < len(edges); p += 2 {
|
||
if edges[p] > edges[p+1] {
|
||
t.Fatalf("edge [%d,%d] is not sorted as a pair", edges[p], edges[p+1])
|
||
}
|
||
}
|
||
for p := range len(want) {
|
||
if edges[p] != want[p] {
|
||
t.Fatalf("boundary edges %v, want %v", edges, want)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonFEM2DDuplicateDirichletNode pins the documented rule
|
||
// for a node listed more than once: the last value is the prescribed
|
||
// one and the node enters the assembled system exactly once, so the
|
||
// repeated listing answers what the single listing with that value
|
||
// answers. Recording it twice appends a second unit row at the same
|
||
// coordinate, which the sparse conversion merges by summing, so the
|
||
// node's diagonal doubles and the solve halves its prescribed value.
|
||
func TestSolvePoissonFEM2DDuplicateDirichletNode(t *testing.T) {
|
||
solution := func(x, y float64) float64 { return math.Sin(math.Pi*x) * math.Sin(math.Pi*y) }
|
||
source := func(x, y float64) float64 { return 2 * math.Pi * math.Pi * solution(x, y) }
|
||
mesh, boundary := gridMesh(t, 4)
|
||
values := make([]float64, len(boundary))
|
||
for p, node := range boundary {
|
||
values[p] = solution(mesh.Vertices[2*node], mesh.Vertices[2*node+1])
|
||
}
|
||
// The list names the second boundary node again at the end, with a
|
||
// different value: the last one wins and the node stays single.
|
||
const extra = 0.5
|
||
nodes := append(append([]int(nil), boundary...), boundary[1])
|
||
dupValues := append(append([]float64(nil), values...), values[1]+extra)
|
||
u, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: nodes, DirichletValues: dupValues})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D with a repeated node: %v", err)
|
||
}
|
||
single := append([]float64(nil), values...)
|
||
single[1] += extra
|
||
want, err := SolvePoissonFEM2D(mesh, source, FEMPoissonOptions{Kappa: 1, DirichletNodes: boundary, DirichletValues: single})
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonFEM2D with the node once: %v", err)
|
||
}
|
||
if got := u.FloatAt(boundary[1]); math.Abs(got-(values[1]+extra)) > 1e-12 {
|
||
t.Fatalf("the repeated node answered %g, want the last prescribed value %g", got, values[1]+extra)
|
||
}
|
||
worst := 0.0
|
||
for i := range mesh.Vertices2() {
|
||
worst = math.Max(worst, math.Abs(u.FloatAt(i)-want.FloatAt(i)))
|
||
}
|
||
if worst > 1e-12 {
|
||
t.Fatalf("the repeated listing differs from the single listing by %g, want the node recorded once", worst)
|
||
}
|
||
}
|