assemble 2D triangle mesh

src/bases/lagrange.jl

@@ -744,15 +744,15 @@ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}})
   geometry = mesh2
   cellids2, ncells2 = vertextocellmap(mesh2)
 
-  fns = Vector{Vector{Shape{T}}}(undef,num_cells1)
+  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[segment_coarse]*4) -2
+    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 segmesnts in the finer mesh
+    #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}
@@ -799,7 +799,7 @@ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}})
         shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,0.5)
     end
     # Now assign all of these shapes to the relevent segment in the coarse mesh
-    fns[segment_coarse]=shapes
+    push!(fns, shapes)
     push!(pos, cartesian(center(chart(mesh, segment_coarse))))
   end
 

src/bases/local/laglocal.jl

@@ -211,7 +211,19 @@ function (f::LagrangeRefSpace{T,10,2})(mp) where T
 end
 
 # Evaluete linear lagrange elements on a triangle
-function (f::LagrangeRefSpace{T,1,3})(t) where T
+function (f::LagrangeRefSpace{T,1,3})(t::CompScienceMeshes.MeshPointNM{T,C,D,2}) where {T,C,D}
+    u,v,w, = barycentric(t)
+
+    j = jacobian(t)
+    p = t.patch
+    SVector(
+        (value=u, curl=(p[3]-p[2])/j),
+        (value=v, curl=(p[1]-p[3])/j),
+        (value=w, curl=(p[2]-p[1])/j))
+end
+
+# Evaluete linear lagrange elements on a triangle
+function (f::LagrangeRefSpace{T,1,3})(t::CompScienceMeshes.MeshPointNM{T,C,D,3}) where {T,C,D}
     u,v,w, = barycentric(t)
 
     j = jacobian(t)

src/bases/local/ndlocal.jl

@@ -12,8 +12,27 @@ function valuetype(ref::NDRefSpace{T}, charttype::Type) where {T}
     SVector{universedimension(charttype),T}
 end
 
+function (ϕ::NDRefSpace)(nbd::CompScienceMeshes.MeshPointNM{T,C,D,2}) where {T,C,D}
 
-function (ϕ::NDRefSpace)(nbd)
+    u, v = parametric(nbd)
+    j = jacobian(nbd)
+
+    tu = tangents(nbd,1)
+    tv = tangents(nbd,2)
+
+    d = 2/j
+
+    rotate_left(t) = StaticArrays.SVector{2, Float64}(t[2], -t[1])
+
+    return SVector((
+        (value = rotate_left(tu*(u-1) + tv*v    ) / j, curl = d),
+        (value = rotate_left(tu*u     + tv*(v-1)) / j, curl = d),
+        (value = rotate_left(tu*u     + tv*v    ) / j, curl = d)
+    ))
+
+end
+
+function (ϕ::NDRefSpace)(nbd::CompScienceMeshes.MeshPointNM{T,C,D,3}) where {T,C,D}
 
     u, v = parametric(nbd)
     n = normal(nbd)

test/runtests.jl

@@ -41,6 +41,7 @@ include("test_interpolate_and_restrict.jl")
 include("test_rt3d.jl")
 include("test_gradient.jl")
 include("test_mult.jl")
+include("test_assemble_2D_surface.jl")
 
 include("test_gram.jl")
 include("test_vector_gram.jl")

test/test_assemble_2D_surface.jl

@@ -0,0 +1,110 @@
+using Test
+
+import Permutations
+import StaticArrays
+using SparseArrays
+using CompScienceMeshes
+using BEAST
+
+"""
+Test if the assemble methode for a 2D mesh with udim=2 works. 
+Test are constructing a rectangle with udim=2 and udim=3 with langrange and nedelec interpolation. 
+And assemble these two end checking that matrix from rectangle of udim=2 is equal to the matrix from the ractangle of udim=3.
+"""
+
+
+
+# overwrite the existing method of CompScienceMeshes for the positive value of the volume.
+function CompScienceMeshes._normals(tangents::StaticArrays.SVector{2,StaticArrays.SVector{2,T}}, ::Type{Val{0}}) where {T}
+
+    t = tangents[1]
+    s = tangents[2]
+    v = abs(t[1]*s[2] - t[2]*s[1])/2
+    # n[3] = tangents[1] × tangents[2]
+    # l = norm(n)
+
+    P = StaticArrays.SVector{2,T}
+    StaticArrays.SVector{0,P}(), v
+end
+
+
+function permutate_vector(X, Y)
+    tol = sqrt(eps(eltype(X[1])))
+
+    permut = Vector{Int32}()  
+    temp = collect(1:length(X))
+
+    for p in Y
+        index = findfirst(isapprox(p;atol = tol), X)
+        @assert !isnothing(index)   # vertex from Y exist in X
+        @assert temp[index] != 0    # vertex from Y unique in X
+
+        temp[index] = 0
+        append!(permut, index)
+    end
+
+    for i in temp
+        if i !=0
+            push!(permut, i)
+    end end
+
+    return permut
+end
+
+
+
+
+
+
+@testset "assemble 2D langrange" begin
+    h = 1/5
+    Ω2D = meshrectangle(1.0, 1.0, h, 2) # udim = 2
+    Ω3D = meshrectangle(1.0, 1.0, h, 3) # udim = 3
+
+    #permutate_mesh(Ω2D, Ω3D)
+
+    # lagrange spaces
+    X2D = lagrangec0d1(Ω2D, skeleton(Ω2D,0))
+    X3D = lagrangec0d1(Ω3D, skeleton(Ω3D,0))
+
+    σ =  permutate_vector(X2D.pos, map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos))
+
+    @test map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos) == X2D.pos[σ]
+
+    #@test X2D.geo.vertices == map(x-> x[1:2], X3D.geo.vertices)
+
+    I = BEAST.Identity()
+
+    A2D = assemble(I, X2D, X2D)
+    A3D = assemble(I, X3D, X3D)
+
+    @test sum(abs.(permute(sparse(A2D), σ, σ) - sparse(A3D))) ≈ 0 atol=1e-8
+
+end
+
+@testset "assemble 2D nedelec" begin
+    h = 1/5
+    Ω2D = meshrectangle(1.0, 1.0, h, 2) # udim = 2
+    Ω3D = meshrectangle(1.0, 1.0, h, 3) # udim = 3
+
+    #permutate_mesh(Ω2D, Ω3D)
+
+    # lagrange spaces
+    X2D = lagrangec0d1(Ω2D, skeleton(Ω2D,0))
+    X3D = lagrangec0d1(Ω3D, skeleton(Ω3D,0))
+
+    σ =  permutate_vector(X2D.pos, map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos))
+
+    @test map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos) == X2D.pos[σ]
+
+    @hilbertspace u
+    @hilbertspace v
+
+    I = BEAST.Identity()
+    IGG = @varform I[gradient(u), gradient(v)]
+
+    A2D = assemble(IGG, X2D, X2D)
+    A3D = assemble(IGG, X3D, X3D)
+
+    @test sum(abs.(permute(sparse(A2D), σ, σ) - sparse(A3D))) ≈ 0 atol=1e-8
+end