diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index d90abca9..81e9e4ec 100644 --- a/src/bases/lagrange.jl +++ b/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 \ No newline at end of file +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 diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 61eacffa..eca8e224 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -15,7 +15,7 @@ valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T # Evaluate constant lagrange elements on anything (ϕ::LagrangeRefSpace{T,0,2})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) - +(ϕ:: LagrangeRefSpace{T,0,4})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) # Evaluate linear Lagrange elements on a segment # The derivative denotes the tangential derivative function (f::LagrangeRefSpace{T,1,2})(mp) where T @@ -596,53 +596,6 @@ function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T return Q end - -function restrict(f::LagrangeRefSpace{T,2}, dom1, dom2) where T - - D = numfunctions(f) - Q = zeros(T, D, D) - - # for each point of the new domain - for i in 1:3 - - #vertices - v = dom2.vertices[i] - - # find the barycentric coordinates in dom1 - uvn = carttobary(dom1, v) - - # evaluate the shape functions in this point - x = neighborhood(dom1, uvn) - fx = f(x) - - for j in 1:D - Q[j,i] = fx[j][1] - end - - - #edges - # find the center of edge i of dom2 - a = dom2.vertices[mod1(i+1,3)] - b = dom2.vertices[mod1(i+2,3)] - v = (a + b) / 2 - - # find the barycentric coordinates in dom1 - uvn = carttobary(dom1, v) - - # evaluate the shape functions in this point - x = neighborhood(dom1, uvn) - fx = f(x) - - for j in 4:D - Q[j,i+3] = fx[j][1] - end - end - - return Q -end - - - # Quadratic Lagrange element on a triangle # function (f::LagrangeRefSpace{T,2,3})(t) where T # u,v,w, = barycentric(t) @@ -905,4 +858,59 @@ function interpolate!(out, fields, interpolant::LagrangeRefSpace{T,0,3}, chart) vals = fields(p) for (g, val) in zip(axes(out, 1), vals) out[g,1] = val -end end \ No newline at end of file +end end + +# segment methods +function interpolate(fields, interpolant::LagrangeRefSpace{T, Degree, 2},chart) where {T, Degree} + vals = Vector{Vector{T}}() + if Degree > 0 + push!(vals, fields(neighborhood(chart, (T(1),)))) + push!(vals, fields(neighborhood(chart, (T(0),)))) + for k in 1:Degree-1 + push!(vals, fields(neighborhood(chart, (T(Degree-k)/Degree,)))) + end + else + push!(vals, fields(center(chart))) + end + Q = Matrix{T}(undef, length(vals[1]), length(vals)) + for i in eachindex(vals) + Q[:,i] .= vals[i] + end + return Q +end + +function interpolate!(out, fields, interpolant::LagrangeRefSpace{T, Degree, 2}, chart) where {T, Degree} + us = Degree > 0 ? (T(1), T(0), (T(Degree-k)/Degree for k in 1:Degree-1)...) : (T(1)/2,) + for (idx, u) in enumerate(us) + vals = fields(neighborhood(chart, (u,))) + for (g, val) in zip(axes(out, 1), vals) + out[g, idx] = val + end + end +end + +# tetrahedron methods +function interpolate(fields, interpolant::LagrangeRefSpace{T, 1, 4}) where {T} + vals = [fields(neighborhood(chart, uvw)) for uvw in ((T(1),T(0),T(0)),(T(0),T(1),T(0)),(T(0),T(0),T(1)),(T(0),T(0),T(0)))] + Q = Matrix{T}(undef, length(vals[1]), length(vals)) + for i in eachindex(vals) + Q[:,i] .= vals[i] + end + return Q +end + +function interpolate!(out, fields, interpolant::LagrangeRefSpace{T, 1, 4}, chart) where {T} + for (idx, uvw) in enumerate(((T(1),T(0),T(0)),(T(0),T(1),T(0)),(T(0),T(0),T(1)),(T(0),T(0),T(0)))) + vals = fields(neighborhood(chart, uvw)) + for (g, val) in zip(axes(out, 1), vals) + out[g, idx] = val + end + end +end + +function interpolate!(out, fields, interpolant::LagrangeRefSpace{T,0,4}, chart) where {T} + vals = fields(center(chart)) + for (g, val) in zip(axes(out, 1), vals) + out[g, 1] = val + end +end diff --git a/src/localop.jl b/src/localop.jl index ada78f1a..e2bca943 100644 --- a/src/localop.jl +++ b/src/localop.jl @@ -188,6 +188,7 @@ function assemble_local_refines!(biop::LocalOperator, tfs::Space, bfs::Space, st P = ta2g[p] Q = CompScienceMeshes.parent(tgeo, P) q = bg2a[Q] + q == 0 && continue bcell = bels[q] @assert overlap(tcell, bcell) diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl index 0da2f922..35d36856 100644 --- a/src/quadrature/rules/testrefinestrialqrule.jl +++ b/src/quadrature/rules/testrefinestrialqrule.jl @@ -37,10 +37,7 @@ function integrate!(out, op, trial_charts = [chart(test_mesh, p) for p in CompScienceMeshes.children(parent_mesh, trial_cell)] trial_overlaps = map(trial_charts) do chart - simplex( - carttobary(trial_chart, chart.vertices[1]), - carttobary(trial_chart, chart.vertices[2]), - carttobary(trial_chart, chart.vertices[3])) + simplex(map(v -> carttobary(trial_chart, v), chart.vertices)) end quadstrat = qr.conforming_qstrat @@ -64,3 +61,35 @@ function integrate!(out, op, for k in 1:size(Q, 2) out[i,j] += zlocal[i,k] * Q[j,k] end end end end end + +@testitem "TestRefinesTrialQRule on a curve" begin + using BEAST, Test + using CompScienceMeshes + using LinearAlgebra + + Γ = meshcircle(1.0, 2π/20) + X = duallagrangec0d1(Γ) + Y = lagrangecxd0(Γ) + @test CompScienceMeshes.refines(BEAST.geometry(X), Γ) + + 𝒱 = Helmholtz2D.singlelayer(wavenumber=1.0) + A = Matrix(assemble(𝒱, X, Y)) + @test all(isfinite, A) + + # reference: lift the coarse constants onto the fine mesh, so both spaces live on the + # same geometry and the conforming rules apply + fine = BEAST.geometry(X) + fns = [[BEAST.Shape(2*(i-1)+1, 1, 1.0), BEAST.Shape(2*(i-1)+2, 1, 1.0)] for i in 1:numcells(Γ)] + pos = [cartesian(CompScienceMeshes.center(chart(Γ,p))) for p in Γ] + Yf = BEAST.LagrangeBasis{0,-1,1}(fine, fns, pos) + @test Matrix(assemble(BEAST.Identity(), X, Yf)) ≈ Matrix(assemble(BEAST.Identity(), X, Y)) + + B = Matrix(assemble(𝒱, X, Yf)) + @test norm(A-B)/norm(B) < 1e-12 + + # opposite orientation goes through TrialRefinesTestQRule + A2 = Matrix(assemble(𝒱, Y, X)) + B2 = Matrix(assemble(𝒱, Yf, X)) + @test all(isfinite,A2) + @test norm(A2 - B2)/norm(B2) < 1e-12 +end diff --git a/test/test_interpolate_and_restrict.jl b/test/test_interpolate_and_restrict.jl index 9671c0e6..bfa53bde 100644 --- a/test/test_interpolate_and_restrict.jl +++ b/test/test_interpolate_and_restrict.jl @@ -58,4 +58,59 @@ itpol = sum(w*val for (w,val) in zip(Q3,vals)) end end +@testitem "restrict Lagrange on a refined cell" begin + using CompScienceMeshes + using LinearAlgebra + + coarse = simplex(point(0,0,0), point(1,0,0), point(0,1,0)) + # one child of the barycentric refinement of coarse: + child = simplex(point(1,0,0),point(0.5,0.5,0), point(1/3,1/3,0)) + + tobary(chart, ch) = simplex(map(v -> carttobary(chart,v), ch.vertices)) + + # the three routes to a restriction matrix must agree on a genuine + # restriciton; self_restriction alone passes under any node ordering + for (degree, n) in ((1,3),(2,6)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 3, n}() + Q1 = BEAST.restrict(rs, coarse, child) + Q2 = BEAST.restrict(rs, coarse, child, tobary(coarse, child)) + Q3 = zeros(n, n) + BEAST.restrict!(Q3, rs, coarse, child, tobary(coarse, child)) + @test Q1 ≈ Q2 + @test Q1 ≈ Q3 + end + for (degree, n) in ((0,1),(1,3),(2,6),(3,10)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 3, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, coarse, coarse, tobary(coarse, coarse)) + @test Q ≈ I + end + + # segments: nodes are the two ends followed by interior points descending + seg = simplex(point(1,0,0), point(0,0,0)) + for (degree, n) in ((0,1),(1,2),(2,3),(3,4),(4,5)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 2, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, seg, seg, tobary(seg, seg)) + @test Q ≈ I + end + + # genuine restriction on a segment: child is the left half + half = simplex(point(1,0,0), point(0.5,0,0)) + for ( degree, n) in ((1,2),(2,3),(3,4)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 2, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, seg, half, tobary(seg, half)) + @test Q ≈ BEAST.restrict(rs, seg, half) + end + + # tetrahedra: only degrees zero and one have an evaluator + tet = simplex(point(0,0,0), point(1,0,0),point(0,1,0), point(0,0,1)) + for (degree, n) in ((0,1),(1,4)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 4, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, tet, tet, tobary(tet, tet)) + @test Q ≈ I + end +end diff --git a/test/test_local_assembly.jl b/test/test_local_assembly.jl index 49946454..866ee6e9 100644 --- a/test/test_local_assembly.jl +++ b/test/test_local_assembly.jl @@ -1,6 +1,7 @@ @info "Executing test_local_assembly.jl" using CompScienceMeshes using BEAST +using LinearAlgebra using Test @@ -46,6 +47,22 @@ BEAST.assemble_local_mixed!(Id, BC, RT, st2) Q2 = fr2() @test isapprox(Q1, Q2, atol=1e-8) +# Dual functions reach into cells of the coarse mesh that carry no trial function. +# Refines algorithm has to skip those cells. +Xd = duallagrangec0d1(m) +Yp = BEAST.subset(BEAST.lagrangecxd0(m), collect(1:2:length(m))) + +fr1, st1 = BEAST.allocatestorage(Id, Xd, Yp, Val{:bandedstorage}, BEAST.LongDelays{:ignore}) +BEAST.assemble_local_refines!(Id, Xd, Yp, st1) +Q1 = fr1() + +fr2, st2 = BEAST.allocatestorage(Id, Xd, Yp, Val{:bandedstorage}, BEAST.LongDelays{:ignore}) +BEAST.assemble_local_mixed!(Id, Xd, Yp, st2) +Q2 = fr2() + +@test isapprox(Q1, Q2, atol=1e-8) +@test norm(Q1) > 0 +@test isapprox(Matrix(assemble(Id, Xd, Yp)), Q2, atol=1e-8) @testitem "localop assembly returns sparse" begin using CompScienceMeshes @@ -71,4 +88,4 @@ Q2 = fr2() # T = typeof(Z1.lmap.A.lmap) @test Z1.lmap.A.lmap isa SparseMatrixCSC -end \ No newline at end of file +end