Bugfix: Dual Lagrange spaces

src/bases/lagrange.jl

@@ -601,7 +601,7 @@ function unitfunctioncxd0(mesh)
     # in case of a global function. Might be replaced by something
     # more useful.
     # For now, we fill it with the average position of the shape functions
-    p = vertextype(mesh)(0.0, 0.0, 0.0)
+    p = zero(vertextype(mesh))
     for cell in mesh
         p += cartesian(center(chart(mesh, cell)))
     end
@@ -814,11 +814,20 @@ the area so that overall the integral over the dual function is one.
 When `interpolatory=true` is used, the function value is one on the support, and thus,
 it gives rise to a partition of unity.
 """
-function duallagrangecxd0(mesh, jct=CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1); interpolatory=false)
-    vertexlist = interior_and_junction_vertices(mesh, jct)
-    duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory)
+function duallagrangecxd0(mesh; interpolatory=false)
+    jct = CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1)
+    duallagrangecxd0(mesh, jct; interpolatory=interpolatory)
 end
 
+function duallagrangecxd0(mesh, jct::CompScienceMeshes.AbstractMesh; interpolatory=false)
+    if dimension(jct) == dimension(mesh) - 1
+        vertexlist = interior_and_junction_vertices(mesh, jct)
+    else
+        @assert dimension(jct) == 0
+        vertexlist = Int[CompScienceMeshes.indices(jct, v)[1] for v in jct]
+    end
+    duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory)
+end
 
 function duallagrangecxd0(mesh, vertexlist::Vector{Int}; interpolatory=false)
 
@@ -841,14 +850,6 @@ function duallagrangecxd0(mesh, vertexlist::Vector{Int}; interpolatory=false)
     LagrangeBasis{0,-1,NF}(fine, fns, pos)
 end
 
-
-function duallagrangecxd0(mesh, vertices::CompScienceMeshes.AbstractMesh{U,1}; interpolatory=false) where {U}
-    # vertexlist = Int[v[1] for v in vertices]
-    vertexlist =Int[CompScienceMeshes.indices(vertices, v)[1] for v in vertices]
-    return duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory)
-end
-
-
 """
     singleduallagd0(fine, F, v; interpolatory=false)
 
@@ -1192,6 +1193,7 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})
     uv_ctr = ones(dimension(mesh))/(dimension(mesh)+1)
 
     vtoc, vton = vertextocellmap(refined)
+    fine_cells = collect(refined)
     for (i,p) in enumerate(mesh)
         coarse_idcs = CompScienceMeshes.indices(mesh, p)
         coarse_chart = chart(mesh,p)
@@ -1208,6 +1210,7 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})
         I = CollisionDetection.find(fine_vertices, centroid)
         @assert length(I) == 1
         centroid_id = I[1]
+        children = vtoc[centroid_id, 1:vton[centroid_id]]
 
         # get the indx in fine.vertices of the centroid of the faces of coarse_cell
         face_center_ids = Vector{Int}(undef,num_faces)
@@ -1224,39 +1227,16 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})
             face_center_ids[f] = I[1]
         end
 
-        n = vton[centroid_id]
-        for c in vtoc[centroid_id,1:n]
-            fine_idcs = cells(refined)[c].indices
-            local_id = something(findfirst(isequal(centroid_id), fine_idcs), 0)
-            @assert local_id != 0
-            shape = Shape(c, local_id, 1.0)
-            push!(fns[i], shape)
-        end
-
-        for f in 1:num_faces
-            v = face_center_ids[f]
-            jct_pred(vertices(refined)[v]) && continue
-            n = vton[v]
-            for c in vtoc[v,1:n]
-                fine_idcs = cells(refined)[c].indices
-                local_id = something(findfirst(isequal(v), fine_idcs),0)
-                @assert local_id != 0
-                #shape = Shape(c, local_id, 1/n/2)
-                shape = Shape(c, local_id, 1/2)
-
-                push!(fns[i], shape)
-            end
-        end
-
-        for f in 1:length(coarse_idcs)
-            v = coarse_idcs[f]
-            jct_pred(vertices(refined)[v]) && continue
+        for v in (centroid_id, face_center_ids..., coarse_idcs...)
+            v != centroid_id && jct_pred(vertices(refined)[v]) && continue
             n = vton[v]
-            for c in vtoc[v,1:n]
-                fine_idcs = cells(refined)[c].indices
-                local_id = something(findfirst(isequal(v), fine_idcs),0)
+            cells_at_v = vtoc[v, 1:n]
+            coeff = T(count(in(children), cells_at_v)) / n
+            for c in cells_at_v
+                fine_idcs = CompScienceMeshes.indices(refined, fine_cells[c])
+                local_id = something(findfirst(isequal(v), fine_idcs), 0)
                 @assert local_id != 0
-                shape = Shape(c, local_id, 1/(n/2))
+                shape = Shape(c, local_id, coeff)
                 push!(fns[i], shape)
             end
         end
@@ -1279,79 +1259,50 @@ This basis function creats the dual Lagrange basis function and return an object
 It also return a gemoetry containing the refined mesh.
 """
 function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}})
-  T = coordtype(mesh)
-  U = universedimension(mesh)
-  # get the information about number of vertices, number of faces , and the maping between vertices and faces for the original mesh
-  numverts1 = numvertices(mesh)
-  num_cells1 = numcells(mesh)
-  cellids1, ncells1=vertextocellmap(mesh)
-  # get the information about number of vertices, number of faces , and the maping between vertices and faces for the refined mesh
-  num_cells2 = numcells(mesh2)
-  numverts2 = numvertices(mesh2)
-  geometry = mesh2
-  cellids2, ncells2 = vertextocellmap(mesh2)
-
-  fns = Vector{Vector{Shape{T}}}()
-  pos = Vector{vertextype(mesh)}()
-  # We will iterate over the coarse mesh segments to assign all the functions to it.
-  for segment_coarse in 1 : num_cells1
-    # For the dual Lagrange there is a 6 shapes per segment
-    numshapes = (ncells1[mesh.faces[segment_coarse].indices[1]]*4) -2
-    shapes = Vector{Shape{T}}(undef,numshapes)
-    # Now we will get all the smaller faces within the coarse segment
-    #i.e The coose segment will have two points, and these tow points are connected to two segments in the finer mesh
-    # This will give us a 4 smaller faces per Dual lagrange basis, we store them first in all_faces
-    # all_faces= Array{SVector{2,Int}}(undef,4)                      # faces in the original segment (4)
-    CT = CompScienceMeshes.SimplexGraph{2}
-    all_faces = Vector{CT}(undef, 4)
-    # the follwoing code get the verteciec for each coarse segment
-    # then it looks for the two faces connected to each point in the finer mesh
-    # if for example the segment would connect to more than two faces we will have
-    # mesh.faces[segment_coarse][1],n] and iterate over how many segments [n] are connected
-    all_faces[1]=mesh2.faces[cellids2[mesh.faces[segment_coarse][1],1]]
-    all_faces[2]=mesh2.faces[cellids2[mesh.faces[segment_coarse][1],2]]
-    all_faces[3]=mesh2.faces[cellids2[mesh.faces[segment_coarse][2],1]]
-    all_faces[4]=mesh2.faces[cellids2[mesh.faces[segment_coarse][2],2]]
-
-    # now we now the first point of the corse segment will have the left hand side basis
-    # and we know that that point is connected to faces 1,2 in the array all_faces
-    # but now we want to which one of them is inner segment and which one is at the edge
-    # as both of these two faces are in left hand side
-    # The inner face will have the left hand side corase point in its first place for example
-    # we have original point 2,3 , and the coarse segment (2,3)
-    # after refinment we will have for example  points 20,2,21,3,22 with four faces (20,2),(2,21),(21,3),(3,22)
-    # so the faces to the left (20,2),(2,21) we can decide the inner one if the original point(2) is at first index
-    # same for right hand side faces (21,3),(3,22) if the original point(3) is in the second index
-    # then the inner faces are (2,21) and (21,3)
-    # the following code does the same checking process and assign the shapes for the dual lagragne right away
-
-    # For the Left hand side faces
-    if(all_faces[1][1] == mesh.faces[segment_coarse][1])
-        shapes[1]= Shape(cellids2[mesh.faces[segment_coarse][1],2],1,0.5)
-        shapes[2]= Shape(cellids2[mesh.faces[segment_coarse][1],1],2,0.5)
-        shapes[3]= Shape(cellids2[mesh.faces[segment_coarse][1],1],1,1.0)
-    elseif(all_faces[2][1] == mesh.faces[segment_coarse][1])
-        shapes[1]= Shape(cellids2[mesh.faces[segment_coarse][1],2],1,1.0)
-        shapes[2]= Shape(cellids2[mesh.faces[segment_coarse][1],2],2,0.5)
-        shapes[3]= Shape(cellids2[mesh.faces[segment_coarse][1],1],1,0.5)
-    end
-    # For the Right hand side faces
-    if(all_faces[3][2] == mesh.faces[segment_coarse][2])
-        shapes[4]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,1.0)
-        shapes[5]= Shape(cellids2[mesh.faces[segment_coarse][2],1],1,0.5)
-        shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],2],2,0.5)
-    elseif(all_faces[4][2] == mesh.faces[segment_coarse][2])
-        shapes[4]= Shape(cellids2[mesh.faces[segment_coarse][2],2],2,1.0)
-        shapes[5]= Shape(cellids2[mesh.faces[segment_coarse][2],2],1,0.5)
-        shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,0.5)
+    T = coordtype(mesh)
+
+    fns = Vector{Vector{Shape{T}}}(undef, numcells(mesh))
+    pos = Vector{vertextype(mesh)}()
+
+    # store the fine mesh's vertices in an octree for fast retrieval
+    fine_vertices = Octree(vertices(mesh2))
+    vtoc, vton = vertextocellmap(mesh2)
+    fine_cells = collect(mesh2)
+
+    for (i, p) in enumerate(mesh)
+        coarse_idcs = CompScienceMeshes.indices(mesh, p)
+        centroid = cartesian(center(chart(mesh, p)))
+
+        fns[i] = Vector{Shape{T}}()
+        push!(pos, centroid)
+
+        # midpoint of a coarse segment is a vertex of the refinement, and the
+        # fine cells meeting there are the two children of that segment
+        I = CollisionDetection.find(fine_vertices, centroid)
+        @assert length(I) == 1
+        centroid_id = I[1]
+        children = vtoc[centroid_id, 1:vton[centroid_id]]
+
+        # The coefficient at a fine vertex is the fraction of the fine cells meeting
+        # there that belong to this coarse segment: one at the midpoint, one half at
+        # an interior coarse vertex, one at the end of an open curve.
+        for v in (centroid_id, coarse_idcs...)
+            v != centroid_id && pred(vertices(mesh2)[v]) && continue
+            n = vton[v]
+            cells_at_v = vtoc[v, 1:n]
+            coeff = T(count(in(children), cells_at_v)) / n
+            for c in cells_at_v
+                fine_idcs = CompScienceMeshes.indices(mesh2, fine_cells[c])
+                local_id = something(findfirst(isequal(v), fine_idcs), 0)
+                @assert local_id != 0
+                push!(fns[i], Shape(c, local_id, coeff))
+            end
+        end
     end
-    # Now assign all of these shapes to the relevent segment in the coarse mesh
-    push!(fns, shapes)
-    push!(pos, cartesian(center(chart(mesh, segment_coarse))))
-  end
+    @assert length(fns) == length(pos)
 
-  NF = 2
-  LagrangeBasis{1,0,NF}(geometry, fns, pos)
+    NF = 2
+    LagrangeBasis{1,0,NF}(mesh2, fns, pos)
 end
 
 # """
@@ -1958,8 +1909,103 @@ end
     @test BEAST._surface(duallagrangec0d1(Γ)) ≈ 6.0
 end
 
+@testitem "duallagrangec0d1 on open surfaces" begin
+    using CompScienceMeshes
+    # dual P1 must reproduce constants on an open surface too: a boundary edge
+    # belongs to one coarse cell, so its coefficient is one, not one half
+    Γ = meshrectangle(1.0, 1.0, 0.25, 3)
+    @test BEAST._surface(duallagrangec0d1(Γ)) ≈ 1.0
+    # closed surfaces are unchanged
+    @test BEAST._surface(duallagrangec0d1(meshcuboid(1.0,1.0,1.0,0.5))) ≈ 6.0
+end
+
 @testitem "duallagrangecxd0" begin
     using CompScienceMeshes
     Γ = meshcuboid(1.0,1.0,1.0,0.5)
     @test BEAST._surface(duallagrangecxd0(Γ)) ≈ 6.0
-end
+end
+
+@testitem "duallagrangecxd0 on curves" begin
+    using CompScienceMeshes
+
+    # closed curve, vertex numbering inherited from the parent surface
+    Γ = boundary(meshrectangle(1.0,1.0,0.5,3))
+    X = duallagrangecxd0(Γ)
+    @test numfunctions(X) == numcells(Γ)
+    @test BEAST._surface(X) ≈ 4.0
+
+    # open curve: the two end vertices carry no dual function
+    Λ = meshsegment(1.0, 1/4, 3)
+    Y = duallagrangecxd0(Λ)
+    @test numfunctions(Y) == numcells(Λ) - 1
+    @test BEAST._surface(Y) ≈ 0.75
+
+    # junction supplied, half-pyramids on end vertices are kept
+    Z = duallagrangecxd0(Λ, boundary(Λ))
+    @test numfunctions(Z) == numcells(Λ) + 1
+    @test BEAST._surface(Z) ≈ 1.0
+end
+
+@testitem "duallagrangec0d1 on curves" begin
+    using CompScienceMeshes
+    using StaticArrays
+
+    function maxjump(X) #function to find maximum jump of basis X
+        fine = BEAST.geometry(X)
+        fcells = collect(fine)
+        incident = Dict{Int, Vector{Int}}()
+        for c in 1:numcells(fine), v in CompScienceMeshes.indices(fine, fcells[c])
+            push!(get!(incident, v, Int[]), c)
+        end
+        jump = 0.0
+        for fn in X.fns
+            val = Dict{Tuple{Int,Int},Float64}()
+            for s in fn
+                v = CompScienceMeshes.indices(fine, fcells[s.cellid])[s.refid]
+                val[(s.cellid, v)] = get(val, (s.cellid, v), 0.0) + s.coeff
+            end
+            for (c, v) in keys(val)
+                xs = [get(val, (c2, v), 0.0) for c2 in incident[v]]
+                jump = max(jump, maximum(xs) - minimum(xs))
+            end
+        end
+        jump
+    end
+
+    len(m) = sum(volume(chart(m, p)) for p in m)
+
+    #close curve
+    n = 8
+    verts = [SVector(cos(2π*k/n), sin(2π*k/n), 0.0) for k in 0:n-1]
+    faces = [CompScienceMeshes.SimplexGraph{2}(SVector(k, mod1(k+1, n))) for k in 1:n]
+    Γ = Mesh(verts, faces)
+    X = duallagrangec0d1(Γ)
+    @test numfunctions(X) == numcells(Γ)
+    @test all(==(6), length.(X.fns))
+    @test maxjump(X) < 1e-12 # check if space is C0
+    @test BEAST._surface(X) ≈ len(Γ)
+
+    # open curve: end segments carry one shape less
+    Λ = meshsegment(1.0,1/4, 3)
+    Y = duallagrangec0d1(Λ)
+    @test numfunctions(Y) == numcells(Λ)
+    @test sort(unique(length.(Y.fns))) == [5, 6]
+    @test maxjump(Y) < 1e-12
+    @test BEAST._surface(Y) ≈ 1.0
+
+    # vertex numbering inherited from a parent surface (github issue #195)
+    Γb = boundary(meshrectangle(1.0,1.0,1/4, 2))
+    Z = duallagrangec0d1(Γb)
+    @test numfunctions(Z) == numcells(Γb)
+    @test maxjump(Z) < 1e-12
+end
+
+@testitem "unitfunctioncxd0 in a two dimensional universe" begin
+    using CompScienceMeshes
+
+    Γ = meshcircle(1.0, 2π/20) # a curve in the plane
+    len = sum(volume(chart(Γ, p)) for p in Γ)
+    @test numfunctions(BEAST.unitfunctioncxd0(Γ)) == 1
+    @test BEAST._surface(duallagrangec0d1(Γ)) ≈ len
+    @test BEAST._surface(duallagrangecxd0(Γ)) ≈ len
+end