quadrature over triangles

Project.toml

@@ -1,7 +1,7 @@
 name = "DomainIntegrals"
 uuid = "cc6bae93-f070-4015-88fd-838f9505a86c"
 authors = ["Daan Huybrechs <[email protected]>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"

src/functional/measure.jl

@@ -155,6 +155,23 @@ function unsafe_weight(μ::DiscreteWeight, i)
 end
 
 
+####################
+## Generic weights
+####################
+
+"A generic continuous weight function."
+struct GenericWeight{T,F,S} <: Weight{T}
+    fun     ::  F
+    support ::  S
+end
+
+GenericWeight(fun, support = ChebyshevInterval()) =
+    GenericWeight{eltype(support)}(fun, support)
+GenericWeight{T}(fun::F, support::S) where {T,F,S} =
+    GenericWeight{T,F,S}(fun, support)
+
+support(μ::GenericWeight) = μ.support
+unsafe_weightfun(μ::GenericWeight, x) = μ.fun(x)
 
 "A generic discrete weight that stores points and weights."
 struct GenericDiscreteWeight{T,P,W} <: DiscreteWeight{T}

src/processing/domain.jl

@@ -63,7 +63,8 @@ function process_domain(qs, integrand, domain::Domain, measure::Lebesgue, proper
     if has_integration_domain(domain)
         paramdomain = integration_domain(domain)
         fmap = mapfrom_integration_domain(domain)
-        (DomainSets.diffvolume(fmap) * (integrand ∘ fmap), paramdomain, Lebesgue{eltype(paramdomain)}(), properties...)
+        process_domain(qs, diffvolume(fmap) * (integrand ∘ fmap), paramdomain, Lebesgue{eltype(paramdomain)}(), properties...)
+        # (diffvolume(fmap) * (integrand ∘ fmap), paramdomain, Lebesgue{eltype(paramdomain)}(), properties...)
     else
         (integrand, domain, measure, properties...)
     end
@@ -82,6 +83,12 @@ integrate_domain(qs, integrand, domain::Point, measure, properties...) =
 #     process_domain(qs, (t->jacdet(m,t)) * (integrand ∘ m), superdomain(domain), measure, properties...)
 # end
 
+# TODO: implement triangles in terms of UnitSimplex, rather than the other way around
+function process_domain(qs, integrand, domain::EuclideanUnitSimplex{2,T}, measure::Lebesgue, properties...) where {T}
+    m = AffineMap(SA[-1 zero(T); 0 1], SA[one(T); 0])
+    process_domain(qs, diffvolume(m) * (integrand ∘ m), LowerRightTriangle(zero(T),one(T)), measure, properties...)
+end
+
 function process_domain(qs, integrand, domain::LowerRightTriangle{T}, measure::Lebesgue, properties...) where {T}
     a = domain.a
     b = domain.b