From af0f6bcda2fd384e5ddb643ed0b50aa63ba6cbc0 Mon Sep 17 00:00:00 2001
From: Thomas Kemmer <thomas@binaryden.de>
Date: Tue, 23 Aug 2022 14:11:37 +0200
Subject: [PATCH] [Base] Extend quadraturepoints to compute points on triangles

---
 docs/src/lib/quadrature.md |  7 +++++
 src/NESSie.jl              |  2 +-
 src/base/quadrature.jl     | 56 +++++++++++++++++++++++++++++++++++++-
 test/base/quadrature.jl    | 30 ++++++++++++++++++++
 4 files changed, 93 insertions(+), 2 deletions(-)

diff --git a/docs/src/lib/quadrature.md b/docs/src/lib/quadrature.md
index 8d5d7c03..3701d643 100644
--- a/docs/src/lib/quadrature.md
+++ b/docs/src/lib/quadrature.md
@@ -7,9 +7,16 @@
     QuadPts3D
 ```
 
+## Quadrature points on elements
+```@docs
+    ElementQuad
+    TriangleQuad
+```
+
 ## Generators
 ```@docs
     quadraturepoints
+    quadraturepoints(::Vector{Triangle})
 ```
 
 ## Laplace potential
diff --git a/src/NESSie.jl b/src/NESSie.jl
index 3d0dc809..0950e3d0 100644
--- a/src/NESSie.jl
+++ b/src/NESSie.jl
@@ -15,7 +15,7 @@ include("base/model.jl")
 export Element, SurfaceElement, VolumeElement, Triangle, Tetrahedron, Charge, Model
 
 include("base/quadrature.jl")
-export QuadraturePoints, QuadPts2D, QuadPts3D, quadraturepoints
+export QuadraturePoints, QuadPts2D, QuadPts3D, ElementQuad, TriangleQuad, quadraturepoints
 
 include("base/util.jl")
 export meshunion, obspoints_line, obspoints_plane
diff --git a/src/base/quadrature.jl b/src/base/quadrature.jl
index a43b2ca0..2d29c2f8 100644
--- a/src/base/quadrature.jl
+++ b/src/base/quadrature.jl
@@ -67,7 +67,7 @@ Generator function for quadrature points:
  * *Triangles*: 7 points per element [[Rad48]](@ref Bibliography)
  * *Tetrahedra*: 5 points per element [[Kea86]](@ref Bibliography)
 
-## Return type
+# Return type
 `QuadPts2D` or `QuadPts3D`
 """ quadraturepoints
 for T in [:Float64, :Float32]
@@ -97,3 +97,57 @@ for T in [:Float64, :Float32]
         @inline quadraturepoints(::Type{Tetrahedron{$(T)}}) = $(varname)
     end
 end
+
+
+# =========================================================================================
+"""
+    abstract type ElementQuad{T <: AbstractFloat} end
+
+Abstract base type for quadrature points on specific elements
+"""
+abstract type ElementQuad{T <: AbstractFloat} end
+
+
+# =========================================================================================
+"""
+    struct TriangleQuad{T} <: ElementQuad{T}
+        elem   ::Triangle{T}   # given element
+        qpts   ::Matrix{T}     # quadrature points on the element
+        weights::Vector{T}     # weights for the quadrature points
+    end
+
+Quadrature points on a specific surface triangle, including weights.
+"""
+struct TriangleQuad{T} <: ElementQuad{T}
+    """Given element"""
+    elem   ::Triangle{T}
+    """Quadrature points on the given element as 3xN matrix"""
+    qpts   ::Matrix{T}
+    """Weights for the quadrature points"""
+    weights::Vector{T}
+end
+
+
+# =========================================================================================
+"""
+    quadraturepoints{T}(elements::Vector{Triangle{T}})
+
+Computes and returns quadrature points on all given elements, including weights.
+
+# Return type
+`Vector{TriangleQuad{T}}`
+"""
+function quadraturepoints(elements::Vector{Triangle{T}}) where T
+    qpts = quadraturepoints(Triangle{T})
+    ret  = Vector{TriangleQuad{T}}(undef, length(elements))
+
+    for i in eachindex(elements)
+        elem = elements[i]
+        mat  = Matrix{T}(undef, 3, qpts.num)
+        for j in 1:qpts.num
+            mat[:, j] .= qpts.x[j] .* (elem.v2 .- elem.v1) .+ qpts.y[j] .* (elem.v3 .- elem.v1) .+ elem.v1
+        end
+        ret[i] = TriangleQuad{T}(elem, mat, qpts.weight)
+    end
+    ret
+end
diff --git a/test/base/quadrature.jl b/test/base/quadrature.jl
index 55493631..c10b198f 100644
--- a/test/base/quadrature.jl
+++ b/test/base/quadrature.jl
@@ -1,6 +1,36 @@
 @testset "quadraturepoints" begin
     for T in testtypes
+        # simple generators
         @test typeof(quadraturepoints(Triangle{T})) == QuadPts2D{T}
         @test typeof(quadraturepoints(Tetrahedron{T})) == QuadPts3D{T}
+
+        # on-triangle generators
+        elements = Triangle{T}[]
+        ret = quadraturepoints(elements)
+        @test typeof(ret) == Vector{TriangleQuad{T}}
+        @test isempty(ret)
+
+        push!(elements, Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1]))
+        qpts = quadraturepoints(Triangle{T})
+        ret  = quadraturepoints(elements)
+        @test typeof(ret) == Vector{TriangleQuad{T}}
+        @test length(ret) == 1
+        @test ret[1].elem === elements[1]
+        @test ret[1].weights === qpts.weight
+        @test size(ret[1].qpts) == (3, qpts.num)
+        for i in 1:qpts.num
+            @test ret[1].qpts[:, i] ≈ T[0, qpts.x[i], qpts.y[i]]
+        end
+
+        push!(elements, Triangle(T[1, 1, 1], T[1, 1, 3], T[1, 3, 1]))
+        ret  = quadraturepoints(elements)
+        @test typeof(ret) == Vector{TriangleQuad{T}}
+        @test length(ret) == 2
+        @test ret[2].elem === elements[2]
+        @test ret[2].weights === qpts.weight
+        @test size(ret[2].qpts) == (3, qpts.num)
+        for i in 1:qpts.num
+            @test ret[2].qpts[:, i] ≈ T[1, 1 + 2 * qpts.y[i], 1 + 2 * qpts.x[i]]
+        end
     end
 end