From b7737ab745534a771e42072a10c2be72a6ac5100 Mon Sep 17 00:00:00 2001 From: bboeyken Date: Thu, 5 Feb 2026 10:31:14 +0100 Subject: [PATCH 1/5] assemble 2D triangle mesh With this assemble works for a triangulated mesh with udim = 2 for lagrange and nedelec interpolation. --- src/bases/local/laglocal.jl | 14 +++- src/bases/local/ndlocal.jl | 21 +++++- test/runtests.jl | 1 + test/test_assemble_2D_surface.jl | 110 +++++++++++++++++++++++++++++++ 4 files changed, 144 insertions(+), 2 deletions(-) create mode 100644 test/test_assemble_2D_surface.jl diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 2b140c88..7a6f4e5a 100644 --- a/src/bases/local/laglocal.jl +++ b/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) diff --git a/src/bases/local/ndlocal.jl b/src/bases/local/ndlocal.jl index a7f0f977..617d3e9b 100644 --- a/src/bases/local/ndlocal.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index 198fa68d..1d8b600a 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_assemble_2D_surface.jl b/test/test_assemble_2D_surface.jl new file mode 100644 index 00000000..94be3bc4 --- /dev/null +++ b/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 \ No newline at end of file From 6e21c050440db2fdee8a420518a74c8d3437aeea Mon Sep 17 00:00:00 2001 From: bboeyken Date: Mon, 16 Mar 2026 21:35:06 +0100 Subject: [PATCH 2/5] duallagrangec0d1 fix First attempt to fix the duallagrangec0d1. --- src/bases/lagrange.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index ef26ac56..f285e873 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -744,12 +744,13 @@ 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 < 1 && continue # skip detached vertices 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 @@ -799,7 +800,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 From cfb83263ee7b2f08d6daf9674495d86a3dfda93a Mon Sep 17 00:00:00 2001 From: bboeyken Date: Mon, 23 Mar 2026 13:21:17 +0100 Subject: [PATCH 3/5] fix duallagrange c0d1 --- src/bases/lagrange.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index f285e873..729a8ab1 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -749,11 +749,10 @@ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}}) # 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 < 1 && continue # skip detached vertices + 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} From 6d5f38ceb10004289095415198c5113dbbae9a4b Mon Sep 17 00:00:00 2001 From: bboeyken Date: Mon, 1 Jun 2026 09:19:53 +0200 Subject: [PATCH 4/5] add gradien for Direct productspaces --- src/bases/divergence.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/bases/divergence.jl b/src/bases/divergence.jl index a4b32add..f66d4b43 100644 --- a/src/bases/divergence.jl +++ b/src/bases/divergence.jl @@ -48,6 +48,11 @@ function curl(x::Space) curl(x, geo, crl) end +function gradient(X::DirectProductSpace{T}) where T + x = Space{T}[gradient(s) for s in X.factors] + DirectProductSpace(x) +end + function gradient(x::Space) ref = refspace(x) From ad80baff7bf935291611909c30014c27a6c56829 Mon Sep 17 00:00:00 2001 From: bboeyken Date: Mon, 1 Jun 2026 09:20:58 +0200 Subject: [PATCH 5/5] Allow handling bases with no functions --- src/excitation.jl | 4 +++- src/localop.jl | 7 +++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/excitation.jl b/src/excitation.jl index 3c7d2c52..f3fc8982 100644 --- a/src/excitation.jl +++ b/src/excitation.jl @@ -39,7 +39,9 @@ function assemble!(field::Functional, tfs::Space, store; quadstrat=defaultquadstrat(field, tfs)) qs = quadstrat(field, tfs) - tels, tad = assemblydata(tfs) + + tr = assemblydata(tfs); tr == nothing && return + tels, tad = tr trefs = refspace(tfs) qd = quaddata(field, trefs, tels, qs) diff --git a/src/localop.jl b/src/localop.jl index 8b3ad69a..ada78f1a 100644 --- a/src/localop.jl +++ b/src/localop.jl @@ -300,8 +300,11 @@ function assemble_local_mixed!(biop::LocalOperator, tfs::Space{T}, bfs::Space{T} num_trefs = numfunctions(trefs, tdom) num_brefs = numfunctions(brefs, bdom) - tels, tad = assemblydata(tfs) - bels, bad = assemblydata(bfs) + tr = assemblydata(tfs); tr == nothing && return + br = assemblydata(bfs); br == nothing && return + + tels, tad = tr + bels, bad = br qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)