methods for assembly (Id,GWP,GWP) added

Project.toml

@@ -37,7 +37,7 @@ AbstractTrees = "0.4.4"
 BlockArrays = "0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 1"
 CollisionDetection = "0.1.6"
 Combinatorics = "0.7, 1"
-CompScienceMeshes = "0.8.2"
+CompScienceMeshes = "0.8.3"
 Compat = "2, 3, 4"
 ConvolutionOperators = "0.4"
 ExtendableSparse = "1.4"

examples/mfie_gwp2.jl

@@ -1,14 +1,9 @@
 using CompScienceMeshes, BEAST
 
 # Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
-Γ = meshsphere(radius=1.0, h=0.1)
-# X, Y = raviartthomas(Γ), buffachristiansen(Γ)
-# X = brezzidouglasmarini(Γ)
-# Y = brezzidouglasmarini(Γ)
-# X = raviartthomas(Γ)
-# Y = raviartthomas(Γ)
-X = BEAST.gwpdiv(Γ; order=2)
-Y = BEAST.gwpdiv(Γ; order=2)
+Γ = meshsphere(radius=1.0, h=0.2)
+X = BEAST.gwpdiv(Γ; order=3)
+Y = BEAST.gwpdiv(Γ; order=3)
 
 ϵ, μ, ω = 1.0, 1.0, 1.0; κ = ω * √(ϵ*μ)
 # κ = 3.0

src/bases/local/gwplocal.jl

@@ -32,8 +32,8 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
     nd3 = point(T, -v, u)
 
     P = SVector{2,T}
-    vals = Vector{P}(undef,NF)
-    crls = Vector{T}(undef,NF)
+    NT = @NamedTuple{value::P, curl::T}
+    nts = Vector{NT}(undef, NF)
     idx = 1
 
     i = 0
@@ -50,10 +50,7 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
             dv = Rᵢ*dRⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ
             curl = (du*nd1[2] - dv*nd1[1]) + 2*Rᵢ*Rⱼ*Rₖ
 
-            vals[idx] = Rᵢ*Rⱼ*Rₖ*nd1
-            crls[idx] = curl
-            # push!(vals, Rᵢ*Rⱼ*Rₖ*nd1)
-            # push!(crls, curl)
+            nts[idx] = (value=Rᵢ*Rⱼ*Rₖ*nd1, curl=curl)
             idx += 1
     end
 
@@ -71,11 +68,7 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
         dv = Rᵢ*dRⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ
         curl = (du*nd2[2] - dv*nd2[1]) + 2*Rᵢ*Rⱼ*Rₖ
 
-        # push!(vals, Rᵢ*Rⱼ*Rₖ*nd2)
-        # push!(crls, curl)
-
-        vals[idx] = Rᵢ*Rⱼ*Rₖ*nd2
-        crls[idx] = curl
+        nts[idx] = (value=Rᵢ*Rⱼ*Rₖ*nd2, curl=curl)
         idx += 1
     end
 
@@ -93,12 +86,8 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
         du = dRᵢ*Rⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ
         dv = Rᵢ*dRⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ
         curl = (du*nd3[2] - dv*nd3[1]) + 2*Rᵢ*Rⱼ*Rₖ
-
-        # push!(vals, Rᵢ*Rⱼ*Rₖ*nd3)
-        # push!(crls, curl)
 
-        vals[idx] = Rᵢ*Rⱼ*Rₖ*nd3
-        crls[idx] = curl
+        nts[idx] = (value=Rᵢ*Rⱼ*Rₖ*nd3, curl=curl)
         idx += 1
     end
 
@@ -135,22 +124,14 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
             dv = Rsᵢ*dRsⱼ*Rₖ - Rsᵢ*Rsⱼ*dRₖ
             curlS3 = du*nd3[2] - dv*nd3[1] + 2*Rsᵢ*Rsⱼ*Rₖ
 
-            # push!(vals, S2-S3)
-            # push!(crls, curlS2 - curlS3)
-            vals[idx] = S2-S3
-            crls[idx] = curlS2 - curlS3
+            nts[idx] = (value=S2-S3, curl=curlS2 - curlS3)
             idx += 1
 
-            # push!(vals, S3-S1)
-            # push!(crls, curlS3 - curlS1)
-            vals[idx] = S3-S1
-            crls[idx] = curlS3 - curlS1
+            nts[idx] = (value=S3-S1, curl=curlS3 - curlS1)
             idx += 1
     end end
 
-    # NF = length(vals)
-    SVector{NF}([(value=f, curl=c) for (f,c) in zip(vals, crls)])
-    # [(value=f, curl=c) for (f,c) in zip(vals, crls)]
+    return SVector{NF}(nts)
 end
 
 

src/identityop.jl

@@ -87,6 +87,21 @@ function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg},
 end
 
 
+defaultquadstrat(::LocalOperator, ::GWPDivRefSpace{<:Real,D1},
+    ::GWPDivRefSpace{<:Real,D2}) where {D1,D2} = SingleNumQStrat(7)
+
+
+function quaddata(op::LocalOperator, g::GWPDivRefSpace, f::GWPDivRefSpace,
+    tels::Vector, bels::Vector,
+    qs::SingleNumQStrat)
+
+    u, w = trgauss(qs.quad_rule)
+    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
+    A = _alloc_workspace(qd, g, f, tels, bels)
+    return qd, A
+end
+
+
 function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs::SingleNumQStrat)
     for i in eachindex(qd)
         q = qd[i]