test interpolation for RTRefSpace

examples/efie.jl

@@ -1,9 +1,9 @@
 using CompScienceMeshes
 using BEAST
 
-# Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
+Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
 # Γ = meshsphere(radius=1.0, h=0.1)
-Γ = CompScienceMeshes.meshmobius(h=0.035)
+# Γ = CompScienceMeshes.meshmobius(h=0.035)
 X = raviartthomas(Γ)
 
 κ, η = 1.0, 1.0

src/bases/local/rtlocal.jl

@@ -99,7 +99,52 @@ const _dof_perms_rt = [
 ]
 
 function dof_permutation(::RTRefSpace, vert_permutation)
-    i = findfirst(==(tuple(vert_permutation...)), _vert_perms_rt)
-    @assert i != nothing
+    i = something(findfirst(==(tuple(vert_permutation...)), _vert_perms_rt),0)
+    @assert i != 0
     return _dof_perms_rt[i]
+end
+
+
+"""
+    interpolate(interpolant::RefSpace, chart1, interpolee::RefSpace, chart2)
+
+Computes by interpolation approximations of the local shape functions for
+`interpolee` on `chart2` in terms of the local shape functions for `interpolant`
+on `chart1`. The returned value is a matrix `Q` such that
+
+```math
+\\phi_i \\approx \\sum_j Q_{ij} \\psi_j
+```
+
+with ``\\phi_i`` the i-th local shape function for `interpolee` and ``\\psi_j`` the
+j-th local shape function for `interpolant`.
+"""
+function interpolate(interpolant::RTRefSpace, chart1, interpolee::RefSpace, chart2)
+    function fields(p)
+        x = cartesian(p)
+        v = carttobary(chart2, x)
+        r = neighborhood(chart2, v)
+        fieldvals = [f.value for f in interpolee(r)]
+    end
+
+    interpolate(fields, interpolant, chart1)
+end
+
+function interpolate(fields, interpolant::RTRefSpace, chart)
+    Q = map(faces(chart)) do face
+        p = center(face)
+        x = cartesian(p)
+        u = carttobary(chart, x)
+        q = neighborhood(chart, u)
+        n = normal(q)
+
+        # minus because in CSM the tangent points towards vertex[1]
+        t = -tangents(p,1)
+        m = cross(t,n)
+
+        fieldvals = fields(q)
+        q = [dot(fv,m) for fv in fieldvals]
+    end
+
+    return hcat(Q...)
 end

test/runtests.jl

@@ -32,6 +32,7 @@ include("test_ndjunction.jl")
 include("test_ndspace.jl")
 include("test_restrict.jl")
 include("test_ndlcd_restrict.jl")
+include("test_interpolate_and_restrict.jl")
 include("test_rt3d.jl")
 include("test_gradient.jl")
 include("test_mult.jl")

test/test_interpolate_and_restrict.jl

@@ -0,0 +1,29 @@
+using Test
+
+using BEAST
+using CompScienceMeshes
+
+chart1 = simplex(
+    point(1,0,0),
+    point(0,1,0),
+    point(0,0,0))
+
+chart2 = simplex(
+    point(1/2,0,0),
+    point(0,1/2,0),
+    point(1,1,0))
+
+X = BEAST.RTRefSpace{Float64}()
+@time Q1 = BEAST.restrict(X, chart1, chart2)
+@time Q2 = BEAST.interpolate(X, chart2, X, chart1)
+@test Q1 ≈ Q2
+
+constant_vector_field = point(1,2,0)
+Q3 = BEAST.interpolate(X, chart2) do p
+    return [constant_vector_field]
+end
+
+ctr = center(chart2)
+vals = [f.value for f in X(ctr)]
+itpol = sum(w*val for (w,val) in zip(Q3,vals))
+@test itpol ≈ constant_vector_field