From 594b052f85a2fc6f2499e936a71f3cc6311c4e5c Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Fri, 4 Oct 2024 12:08:02 +1000 Subject: [PATCH 1/9] Allowing JLD2 v0.5 in Project.toml compats --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 310df8f9b..f63436fae 100644 --- a/Project.toml +++ b/Project.toml @@ -42,7 +42,7 @@ FastGaussQuadrature = "0.4.2, 1" FileIO = "1.2.2, 1.3, 1.4" FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 0.13, 1" ForwardDiff = "0.10.10" -JLD2 = "0.1.11, 0.3, 0.4" +JLD2 = "0.1.11, 0.3, 0.4, 0.5" JSON = "0.21.0" LineSearches = "7.0.1" NLsolve = "4.3.0" From 4df80c09e345f2a6de91532e3d653168f185c4fd Mon Sep 17 00:00:00 2001 From: ThePurox <28731015+ThePurox@users.noreply.github.com> Date: Tue, 8 Oct 2024 10:52:06 +0200 Subject: [PATCH 2/9] Fix broken GMSH link in Readme --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 851451681..76ee4c739 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,7 @@ pkg> add Gridap - [GridapDistributed](https://github.com/gridap/GridapDistributed.jl) Distributed-memory extension of Gridap. - [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) Embedded finite elements in Julia. -- [GridapGmsh](https://github.com/gridap/GridapGmsh.jl) Generate a FE mesh with [GMSH](www.gmsh.info) and use it in Gridap. +- [GridapGmsh](https://github.com/gridap/GridapGmsh.jl) Generate a FE mesh with [GMSH](https://www.gmsh.info) and use it in Gridap. - [GridapMakie](https://github.com/gridap/GridapMakie.jl) Makie plotting recipes for Gridap. - [GridapPardiso](https://github.com/gridap/GridapPardiso.jl) Use the [Intel Pardiso MKL direct sparse solver](https://software.intel.com/en-us/mkl-developer-reference-fortran-intel-mkl-pardiso-parallel-direct-sparse-solver-interface) in Gridap. - [GridapPETSc](https://github.com/gridap/GridapPETSc.jl) Use [PETSc](https://petsc.org/) linear and nonlinear solvers in Gridap. From f7941daa159824421226de641c41b0dac37a15df Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 11 Oct 2024 18:02:15 +1100 Subject: [PATCH 3/9] Minor --- test/AdaptivityTests/MacroFEStokesTests.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/AdaptivityTests/MacroFEStokesTests.jl b/test/AdaptivityTests/MacroFEStokesTests.jl index b961afb12..1287d64c4 100644 --- a/test/AdaptivityTests/MacroFEStokesTests.jl +++ b/test/AdaptivityTests/MacroFEStokesTests.jl @@ -11,15 +11,14 @@ function main(Dc,reftype) @assert reftype ∈ [:barycentric,:powellsabin] u_sol(x) = (Dc == 2) ? VectorValue(x[1],-x[2]) : VectorValue(x[1],-x[2],0.0) - p_sol(x) = (x[1] - 1.0/2.0) + p_sol(x) = x[1] - 1.0/2.0 domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) nc = (Dc == 2) ? (2,2) : (1,1,1) model = simplexify(CartesianDiscreteModel(domain,nc)) - min_order = (reftype == :barycentric) ? Dc : Dc-1 - order = max(2,min_order) poly = (Dc == 2) ? TRI : TET + order = (reftype == :barycentric) ? Dc : Dc-1 rrule = (reftype == :barycentric) ? Adaptivity.BarycentricRefinementRule(poly) : Adaptivity.PowellSabinRefinementRule(poly) subreffes_u = Fill(LagrangianRefFE(VectorValue{Dc,Float64},poly,order),Adaptivity.num_subcells(rrule)) @@ -62,7 +61,7 @@ function main(Dc,reftype) end end -# NOTE: Powell-Sabin split not working yet. The issue is we woudl need a global cellmap +# NOTE: Powell-Sabin split not working yet. The issue is we would need a global cellmap # directly from the sub-cells to the physical domain (due to how the split is built). # This is something I may do in the future. From eae50d24b043c96fb855e6db003d319094ac0ba3 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 12 Oct 2024 12:44:06 +1100 Subject: [PATCH 4/9] Minor --- test/AdaptivityTests/MacroFEStokesTests.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/AdaptivityTests/MacroFEStokesTests.jl b/test/AdaptivityTests/MacroFEStokesTests.jl index 1287d64c4..2dfe2c53b 100644 --- a/test/AdaptivityTests/MacroFEStokesTests.jl +++ b/test/AdaptivityTests/MacroFEStokesTests.jl @@ -27,7 +27,7 @@ function main(Dc,reftype) subreffes_p = Fill(LagrangianRefFE(Float64,poly,order-1),Adaptivity.num_subcells(rrule)) reffe_p = Adaptivity.MacroReferenceFE(rrule,subreffes_p;conformity=L2Conformity()) - qdegree = 2*order + qdegree = 2*(order-1) quad = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,qdegree) V = FESpace(model,reffe_u,dirichlet_tags=["boundary"]) @@ -53,8 +53,6 @@ function main(Dc,reftype) uh, ph = xh eh_u = uh - u_sol eh_p = ph - p_sol - println(sum(∫(eh_u⋅eh_u)dΩ)) - println(sum(∫(eh_p⋅eh_p)dΩ)) @test sum(∫(eh_u⋅eh_u)dΩ) < 1.e-10 if reftype != :powellsabin @test sum(∫(eh_p*eh_p)dΩ) < 1.e-10 From 05eb50660aba62c682dddcae57af1a08973115ff Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 12 Oct 2024 15:59:42 +1100 Subject: [PATCH 5/9] More tests --- test/AdaptivityTests/FineToCoarseFieldsTests.jl | 14 +++++++------- test/AdaptivityTests/MacroFETests.jl | 5 +++++ test/AdaptivityTests/RefinementRulesTests.jl | 12 ++++++++++++ 3 files changed, 24 insertions(+), 7 deletions(-) diff --git a/test/AdaptivityTests/FineToCoarseFieldsTests.jl b/test/AdaptivityTests/FineToCoarseFieldsTests.jl index e09b54f8d..89beb196b 100644 --- a/test/AdaptivityTests/FineToCoarseFieldsTests.jl +++ b/test/AdaptivityTests/FineToCoarseFieldsTests.jl @@ -91,13 +91,13 @@ eh2 = u_c - u_fc eh3 = u_c - u_fc2 @test sum(∫(eh3⋅eh3)*dΩ_c) < 1.e-12 -modelH=CartesianDiscreteModel((0,1,0,1),(1,1)) -modelh=refine(modelH,2) -reffe=LagrangianRefFE(Float64,QUAD,1) -XH = TestFESpace(modelH,reffe) -xH = get_fe_basis(XH) +modelH = CartesianDiscreteModel((0,1,0,1),(1,1)) +modelh = refine(modelH,2) +reffe = LagrangianRefFE(Float64,QUAD,1) +XH = TestFESpace(modelH,reffe) +xH = get_fe_basis(XH) xHh = change_domain(xH,get_triangulation(modelh),ReferenceDomain()) -evaluate(Gridap.CellData.get_data(xHh)[1], - [Point(0.0,0.0),Point(0.5,0.5)]) +evaluate(Gridap.CellData.get_data(xHh)[1],[Point(0.0,0.0),Point(0.5,0.5)]) +evaluate(Gridap.CellData.get_data(xHh)[1],Point(0.5,0.5)) end \ No newline at end of file diff --git a/test/AdaptivityTests/MacroFETests.jl b/test/AdaptivityTests/MacroFETests.jl index 09d455f6b..e7ec83f7a 100644 --- a/test/AdaptivityTests/MacroFETests.jl +++ b/test/AdaptivityTests/MacroFETests.jl @@ -13,7 +13,12 @@ function test_macro_reffe(model,fmodel,rrule,order) reffe = LagrangianRefFE(Float64,poly,order) sub_reffes = Fill(reffe,num_subcells(rrule)) macro_reffe = Adaptivity.MacroReferenceFE(rrule,sub_reffes) + ReferenceFEs.test_reference_fe(macro_reffe) + macro_quad = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,2*order) + macro_quad_bis = Adaptivity.CompositeQuadrature(Quadrature(poly,4*order),rrule) + ReferenceFEs.test_quadrature(macro_quad) + ReferenceFEs.test_quadrature(macro_quad_bis) Ω = Triangulation(model) Ωf = Triangulation(fmodel) diff --git a/test/AdaptivityTests/RefinementRulesTests.jl b/test/AdaptivityTests/RefinementRulesTests.jl index a93e903a2..3a100d322 100644 --- a/test/AdaptivityTests/RefinementRulesTests.jl +++ b/test/AdaptivityTests/RefinementRulesTests.jl @@ -24,4 +24,16 @@ for poly in polys end end +rr_bc2 = Adaptivity.BarycentricRefinementRule(TRI) +Adaptivity.test_refinement_rule(rr_bc2) + +rr_bc3 = Adaptivity.BarycentricRefinementRule(TET) +Adaptivity.test_refinement_rule(rr_bc3) + +rr_ps2 = Adaptivity.PowellSabinRefinementRule(TRI) +Adaptivity.test_refinement_rule(rr_ps2) + +rr_ps3 = Adaptivity.PowellSabinRefinementRule(TET) +Adaptivity.test_refinement_rule(rr_ps3) + end \ No newline at end of file From 60d6ed8c27ad3cd7e1fc1e30b5301a03b9436099 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 12 Oct 2024 18:10:11 +1100 Subject: [PATCH 6/9] Minor --- src/Adaptivity/CompositeQuadratures.jl | 2 -- test/AdaptivityTests/RefinementRulesTests.jl | 3 --- 2 files changed, 5 deletions(-) diff --git a/src/Adaptivity/CompositeQuadratures.jl b/src/Adaptivity/CompositeQuadratures.jl index c3cc9a793..be1c59a13 100644 --- a/src/Adaptivity/CompositeQuadratures.jl +++ b/src/Adaptivity/CompositeQuadratures.jl @@ -71,8 +71,6 @@ the quadrature `quad` into the subcells of the RefinementRule. function CompositeQuadrature( quad::Quadrature,rr::RefinementRule ) - @check ReferenceFEs.get_polytope(quad) === ReferenceFEs.get_polytope(rr) - weights = get_weights(quad) cpoints = get_coordinates(quad) diff --git a/test/AdaptivityTests/RefinementRulesTests.jl b/test/AdaptivityTests/RefinementRulesTests.jl index 3a100d322..4cc8f63a5 100644 --- a/test/AdaptivityTests/RefinementRulesTests.jl +++ b/test/AdaptivityTests/RefinementRulesTests.jl @@ -31,9 +31,6 @@ rr_bc3 = Adaptivity.BarycentricRefinementRule(TET) Adaptivity.test_refinement_rule(rr_bc3) rr_ps2 = Adaptivity.PowellSabinRefinementRule(TRI) -Adaptivity.test_refinement_rule(rr_ps2) - rr_ps3 = Adaptivity.PowellSabinRefinementRule(TET) -Adaptivity.test_refinement_rule(rr_ps3) end \ No newline at end of file From 28babb6da251224e4f709764744f86b0b88a99b0 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 12 Oct 2024 22:20:57 +1100 Subject: [PATCH 7/9] Minor --- src/Adaptivity/MacroFEs.jl | 1 + test/3d_eval.jl | 41 +++++++ test/AdaptivityTests/MacroFETests.jl | 1 - test/integration_benchmark.jl | 46 ++++++++ test/kron_vector.jl | 44 +++++++ test/mytest.jl | 92 +++++++++++++++ test/mytest2.jl | 52 +++++++++ test/raviart_thomas.jl | 57 +++++++++ test/reusable_inserters.jl | 165 +++++++++++++++++++++++++++ 9 files changed, 498 insertions(+), 1 deletion(-) create mode 100644 test/3d_eval.jl create mode 100644 test/integration_benchmark.jl create mode 100644 test/kron_vector.jl create mode 100644 test/mytest.jl create mode 100644 test/mytest2.jl create mode 100644 test/raviart_thomas.jl create mode 100644 test/reusable_inserters.jl diff --git a/src/Adaptivity/MacroFEs.jl b/src/Adaptivity/MacroFEs.jl index dc56e8089..4fd09c1a2 100644 --- a/src/Adaptivity/MacroFEs.jl +++ b/src/Adaptivity/MacroFEs.jl @@ -348,6 +348,7 @@ function MacroReferenceFE( ndofs = num_free_dofs(space) poly = get_polytope(rrule) + # This is a hack to be able to compute the orders prebasis = FineToCoarseArray(rrule,collect(map(get_prebasis,reffes))) metadata = (rrule,conn,face_own_dofs,face_own_perms) diff --git a/test/3d_eval.jl b/test/3d_eval.jl new file mode 100644 index 000000000..0640db9aa --- /dev/null +++ b/test/3d_eval.jl @@ -0,0 +1,41 @@ + +using Gridap +using Gridap.FESpaces, Gridap.CellData, Gridap.Fields, Gridap.Arrays + +using FillArrays + +# Function that goes from 2D points to 2-component vectors +u2D(x::VectorValue{2}) = VectorValue(x[1],x[2]) + +function point_3D_to_2D(x::VectorValue{3}) + return VectorValue(x[1],x[2]) +end + +function result_2D_to_3D(v::VectorValue{2}) + return VectorValue(v[1],v[2],0.0) +end + +function embed_3D(u2D::CellField,trian_3D::Triangulation) + u22_data = CellData.get_data(u2D) + ncells = length(u22_data) + + u32_data = lazy_map(Broadcasting(∘),u22_data,Fill(GenericField(point_3D_to_2D),ncells)) + u33_data = lazy_map(Broadcasting(∘),Fill(GenericField(result_2D_to_3D),ncells),u32_data) + u3D = CellData.similar_cell_field(u2D,u33_data,trian_3D,PhysicalDomain()) + return u3D +end + +model_2D = CartesianDiscreteModel((0,1,0,1),(2,2)) +model_3D = CartesianDiscreteModel((0,1,0,1,0,1),(2,2,2)) + +reffe_2D = ReferenceFE(lagrangian,VectorValue{2,Float64},1) +reffe_3D = ReferenceFE(lagrangian,VectorValue{3,Float64},1) + +V_2D = TestFESpace(model_2D,reffe_2D) +V_3D = TestFESpace(model_3D,reffe_3D) + +uh_2D = interpolate(u2D,V_2D) +uh_3D = embed_3D(uh_2D,get_triangulation(V_3D)) + +pts = get_cell_points(get_triangulation(V_3D)) +uh_3D(pts) diff --git a/test/AdaptivityTests/MacroFETests.jl b/test/AdaptivityTests/MacroFETests.jl index e7ec83f7a..b97da8a7c 100644 --- a/test/AdaptivityTests/MacroFETests.jl +++ b/test/AdaptivityTests/MacroFETests.jl @@ -13,7 +13,6 @@ function test_macro_reffe(model,fmodel,rrule,order) reffe = LagrangianRefFE(Float64,poly,order) sub_reffes = Fill(reffe,num_subcells(rrule)) macro_reffe = Adaptivity.MacroReferenceFE(rrule,sub_reffes) - ReferenceFEs.test_reference_fe(macro_reffe) macro_quad = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,2*order) macro_quad_bis = Adaptivity.CompositeQuadrature(Quadrature(poly,4*order),rrule) diff --git a/test/integration_benchmark.jl b/test/integration_benchmark.jl new file mode 100644 index 000000000..2e3932d4d --- /dev/null +++ b/test/integration_benchmark.jl @@ -0,0 +1,46 @@ + +using BenchmarkTools, ProfileView +using SparseArrays, LinearAlgebra +using Gridap, Gridap.FESpaces, Gridap.Algebra + +model = CartesianDiscreteModel((0,1,0,1),(100,100)) + +reffe_u = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1) +reffe_p = ReferenceFE(lagrangian, Float64, 1) +reffe_j = ReferenceFE(raviart_thomas, Float64, 0) +U = TestFESpace(model, reffe_u) +P = TestFESpace(model, reffe_p) +J = TestFESpace(model, reffe_j) + +X = MultiFieldFESpace([U,U]) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 2) + +Γ = Boundary(model) +dΓ = Measure(Γ, 2) + +a1(u,v) = ∫(∇(u)⊙∇(v))dΩ +A1 = assemble_matrix(a1,U,U); +@benchmark assemble_matrix!($a1,$A1,$U,$U) + +a2(u,v) = ∫(∇(u)⊙∇(v))dΩ + ∫(u⋅v)dΓ +A2 = assemble_matrix(a2,U,U); +@benchmark assemble_matrix!($a2,$A2,$U,$U) + +a3(u,v) = ∫(divergence(u)⊙divergence(v))dΩ +A3 = assemble_matrix(a3,U,U); +@benchmark assemble_matrix!($a3,$A3,$U,$U) + +a4(u,v) = ∫(divergence(u)⊙divergence(v))dΩ +A4 = assemble_matrix(a4,J,J); +@benchmark assemble_matrix!($a4,$A4,$J,$J) + +a5((u1,u2),(v1,v2)) = ∫(∇(u1)⊙∇(v1))dΩ + ∫(∇(u2)⊙∇(v2))dΩ +A5 = assemble_matrix(a5,X,X); +@benchmark assemble_matrix!($a5,$A5,$X,$X) + +ProfileView.@profview @benchmark assemble_matrix!($a1,$A1,$U,$U) +ProfileView.@profview @benchmark assemble_matrix!($a2,$A2,$U,$U) +ProfileView.@profview @benchmark assemble_matrix!($a3,$A3,$U,$U) +ProfileView.@profview @benchmark assemble_matrix!($a4,$A4,$J,$J) diff --git a/test/kron_vector.jl b/test/kron_vector.jl new file mode 100644 index 000000000..3d88c9a87 --- /dev/null +++ b/test/kron_vector.jl @@ -0,0 +1,44 @@ + +using Gridap +using Gridap.Arrays, Gridap.Geometry, Gridap.FESpaces + +model = CartesianDiscreteModel((0,1,0,1),(5,5)) + +reffe = ReferenceFE(lagrangian,Float64,1) +V = FESpace(model,reffe) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2) + +Γ = BoundaryTriangulation(model) +dΓ = Measure(Γ,2) + +a(u) = ∫(u)dΓ + +u = get_trial_fe_basis(V) +v = get_fe_basis(V) + +cu = get_array(a(u)) +cv = get_array(a(v)) + +nu = length(cu) +nv = length(cv) +cu_ids = vcat([collect(1:nu) for i in 1:nv]...) # Fast index +cv_ids = vcat([fill(i,nu) for i in 1:nv]...) # Slow index +dofs = get_cell_dof_ids(V) + +cu_exp = CompressedArray(cu,cu_ids) +cv_exp = CompressedArray(cv,cv_ids) + +uids_exp = CompressedArray(dofs,cu_ids) +vids_exp = CompressedArray(dofs,cv_ids) + +δ = fill(1.0,nv*nu) + +f(δ,v1,v2) = δ .* (v1 .- v2).^2 +arr = lazy_map(f,δ,cv_exp,cu_exp) +matdata = ([arr],[vids_exp],[uids_exp]) + +assem = SparseMatrixAssembler(V,V) +assemble_matrix(assem,matdata) + diff --git a/test/mytest.jl b/test/mytest.jl new file mode 100644 index 000000000..9ac0964a0 --- /dev/null +++ b/test/mytest.jl @@ -0,0 +1,92 @@ + +using BenchmarkTools +using SparseArrays, LinearAlgebra +using Gridap, Gridap.FESpaces, Gridap.Algebra +using Gridap.CellData, Gridap.Arrays, Gridap.Fields + +model = CartesianDiscreteModel((0,1,0,1),(2,2)) + +order = 1 +reffe = ReferenceFE(lagrangian, Float64, order) +V = TestFESpace(model, reffe) +X = MultiFieldFESpace([V,V]) + +Ω = Triangulation(model) +Γ = BoundaryTriangulation(model) +Λ = SkeletonTriangulation(model) + +degree = 2*order +dΩ = Measure(Ω,degree) +dΓ = Measure(Γ,degree) +dΛ = Measure(Λ,degree) + +n_Γ = get_normal_vector(Γ) +n_Λ = get_normal_vector(Λ) +γ = 1.0 +γ0 = 1.0 +h = 1.0 +f = 1.0 +g = 10.0 + +a((u,p),(v,q)) = + ∫( ∇(v)⊙∇(u) - ∇(q)⋅u + v⋅∇(p) )*dΩ + + ∫( (γ/h)*v⋅u - v⋅(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))⋅u + 2*(q*n_Γ)⋅u )*dΓ + + ∫( + (γ/h)*jump(v⊗n_Λ)⊙jump(u⊗n_Λ) - + jump(v⊗n_Λ)⊙mean(∇(u)) - + mean(∇(v))⊙jump(u⊗n_Λ) + + (γ0*h)*jump(q*n_Λ)⋅jump(p*n_Λ) + + jump(q*n_Λ)⋅mean(u) - + mean(v)⋅jump(p*n_Λ) + )*dΛ + +l((v,q)) = + ∫( v⋅f + q*g )*dΩ + + ∫( (γ/h)*v⋅u - (n_Γ⋅∇(v))⋅u + (q*n_Γ)⋅u )*dΓ + +op = AffineFEOperator(a,l,X,X) + + +u, p = get_trial_fe_basis(X) +v, q = get_fe_basis(X) + +u = get_trial_fe_basis(V) +v = get_fe_basis(V) + +cf1 = jump(v⊗n_Λ) +cf2 = jump(u⊗n_Λ) +cf3 = cf1 ⊙ cf2 + +cf4 = mean(∇(u)) +cf5 = cf1 ⊙ cf4 + +cf1_data = testitem(get_data(cf1)) +cf4_data = testitem(get_data(cf4)) + +_cf4 = change_domain(cf4,Λ,ReferenceDomain()) +_cf4_data = testitem(get_data(_cf4)) +return_value(Broadcasting(Operation(inner)),cf1_data,_cf4_data) + +get_data(cf5) + +args = map(get_data,cf5.args) +lazy_map(Broadcasting(cf5.op),args...) + +pts = CellData._get_cell_points(cf1,cf2) +x = testitem(get_data(pts)) +f1 = testitem(get_data(cf1)) +f2 = testitem(get_data(cf2)) + +fx1 = return_value(f1,x) +fx2 = return_value(f2,x) +r = Fields.BroadcastingFieldOpMap(⊙)(fx1,fx2) + +cf4 = v⊗n_Λ +cf4_p = testitem(get_data(cf4.plus)) +cf4_m = testitem(get_data(cf4.minus)) + +get_data(jump(n_Λ)) + +fkx = map(cfk -> evaluate(testitem(get_data(cfk)),x),cf1.args) +r = Fields.BroadcastingFieldOpMap(+)(fkx...) + diff --git a/test/mytest2.jl b/test/mytest2.jl new file mode 100644 index 000000000..f3aecda97 --- /dev/null +++ b/test/mytest2.jl @@ -0,0 +1,52 @@ + +using SparseArrays, LinearAlgebra +using Gridap, Gridap.FESpaces, Gridap.Algebra, Gridap.Adaptivity + +using GridapDistributed, PartitionedArrays + +sol(x) = x + +np = (2,2) +ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) +end + +cmodel = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4)) +fmodel = refine(cmodel) + +reffe_u = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1) +reffe_j = ReferenceFE(raviart_thomas, Float64, 0) + +Vh = TestFESpace(fmodel, reffe_u, dirichlet_tags="boundary") +Dh = TestFESpace(fmodel, reffe_j, dirichlet_tags="boundary") +Uh = TrialFESpace(Vh, sol) +Jh = TrialFESpace(Dh, sol) + +VH = TestFESpace(cmodel, reffe_u, dirichlet_tags="boundary") +DH = TestFESpace(cmodel, reffe_j, dirichlet_tags="boundary") +UH = TrialFESpace(VH, sol) +JH = TrialFESpace(DH, sol) + +f(x) = VectorValue(x[1]-x[2], x[1]+x[2]) + +Xh,XH,Yh,YH = Uh,UH,Vh,VH +#Xh,XH,Yh,YH = Jh,JH,Dh,DH + +fv_h = zero_free_values(Xh) +fv_H = zero_free_values(XH) + +dv_h = zero_dirichlet_values(Xh) +dv_H = zero_dirichlet_values(XH) + +fill!(fv_H,1.0) +uH = FEFunction(XH,fv_H,dv_H) +uh = interpolate!(uH,fv_h,Xh) + +fv_h_bis = zero_free_values(Yh) +fv_H_bis = zero_free_values(YH) + +fill!(fv_H_bis,1.0) +vH = FEFunction(XH,fv_H_bis,dv_H) +vh = interpolate!(uH,fv_h_bis,Xh) + +fv_H_bis == fv_H diff --git a/test/raviart_thomas.jl b/test/raviart_thomas.jl new file mode 100644 index 000000000..ec7906842 --- /dev/null +++ b/test/raviart_thomas.jl @@ -0,0 +1,57 @@ +using Gridap +using Gridap.ReferenceFEs, Gridap.Geometry, Gridap.FESpaces, Gridap.Arrays +using Gridap.Fields + +using BenchmarkTools + +using Profile +using ProfileView + +using FillArrays + +model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(40,40))) + + +order = 2 +reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_RT = ReferenceFE(raviart_thomas,Float64,order-1) + +V = FESpace(model,reffe) +V_RT = FESpace(model,reffe_RT) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) + +a(u,v) = ∫(u⋅v)dΩ +b(u,v) = ∫((∇⋅u)*(∇⋅v))dΩ +c(u,v) = ∫(∇(u)⊙∇(v))dΩ + +A1 = assemble_matrix(a,V,V) +A2 = assemble_matrix(a,V_RT,V_RT) + +B1 = assemble_matrix(b,V,V) +B2 = assemble_matrix(b,V_RT,V_RT) + +C1 = assemble_matrix(c,V,V) +C2 = assemble_matrix(c,V_RT,V_RT) + +@benchmark assemble_matrix!(a,A1,V,V) +@benchmark assemble_matrix!(a,A2,V_RT,V_RT) +@benchmark assemble_matrix!(b,B1,V,V) +@benchmark assemble_matrix!(b,B2,V_RT,V_RT) +@benchmark assemble_matrix!(c,C1,V,V) +@benchmark assemble_matrix!(c,C2,V_RT,V_RT) + +function p1() + for i in 1:100 + assemble_matrix!(b,B2,V_RT,V_RT) + end +end +ProfileView.@profview p1() + + +mon = V.fe_basis.cell_basis.values[1].fields +hdiv = V_RT.fe_basis.cell_basis.args[1].values[1].fields.qgrad + + + diff --git a/test/reusable_inserters.jl b/test/reusable_inserters.jl new file mode 100644 index 000000000..4922accd7 --- /dev/null +++ b/test/reusable_inserters.jl @@ -0,0 +1,165 @@ + +using SparseArrays, LinearAlgebra +using Gridap, Gridap.Algebra, Gridap.FESpaces, Gridap.Helpers + +struct MinReuse end + +struct ReusableCounterCSC{Tv,Ti} + tv::Type{Tv} + nrows::Int + ncols::Int + colnnzmax::Vector{Ti} +end + +Algebra.LoopStyle(::Type{<:ReusableCounterCSC}) = Loop() + +@inline function Algebra.add_entry!(::Function,a::ReusableCounterCSC{Tv,Ti},v,i,j) where {Tv,Ti} + a.colnnzmax[j] += Ti(1) + nothing +end + +mutable struct ReusableInserterCSC{Tv,Ti} + nrows::Int + ncols::Int + colptr::Vector{Ti} + colnnz::Vector{Ti} + rowval::Vector{Ti} + nzval::Vector{Tv} + nnz::Ti + nnzptr::Vector{Ti} +end + +Algebra.LoopStyle(::Type{<:ReusableInserterCSC}) = Loop() + +@noinline function Algebra.add_entry!(::typeof(+),a::ReusableInserterCSC{Tv,Ti},v,i,j) where {Tv,Ti} + pini = Int(a.colptr[j]) + pend = pini + Int(a.colnnz[j]) - 1 + p = searchsortedfirst(a.rowval,i,pini,pend,Base.Order.Forward) + if (p>pend) + # add new entry + a.colnnz[j] += 1 + a.rowval[p] = i + a.nzval[p] = v + elseif a.rowval[p] != i + # shift one forward from p to pend + @check pend+1 < Int(a.colptr[j+1]) + for k in pend:-1:p + o = k + 1 + a.rowval[o] = a.rowval[k] + a.nzval[o] = a.nzval[k] + end + # add new entry + a.colnnz[j] += 1 + a.rowval[p] = i + a.nzval[p] = v + else + # update existing entry + a.nzval[p] += v + end + # Save the position of the entry for reuse + a.nnz += 1 + a.nnzptr[a.nnz] = p + nothing +end + +############################################################################################ + +function Algebra.nz_counter( + ::SparseMatrixBuilder{SparseMatrixCSC{Tv,Ti},<:MinReuse}, + axes +) where {Tv,Ti} + nrows = length(axes[1]) + ncols = length(axes[2]) + colnnzmax = zeros(Ti,ncols) + ReusableCounterCSC(Tv,nrows,ncols,colnnzmax) +end + +function Algebra.nz_allocation(a::ReusableCounterCSC{Tv,Ti}) where {Tv,Ti} + colptr = Vector{Ti}(undef,a.ncols+1) + @inbounds for i in 1:a.ncols + colptr[i+1] = a.colnnzmax[i] + end + length_to_ptrs!(colptr) + ndata = colptr[end] - one(Ti) + rowval = Vector{Ti}(undef,ndata) + nzval = zeros(Tv,ndata) + nnzmax = sum(a.colnnzmax) + nnzptr = fill(zero(Ti),nnzmax) + colnnz = a.colnnzmax + fill!(colnnz,zero(Ti)) + ReusableInserterCSC(a.nrows,a.ncols,colptr,colnnz,rowval,nzval,0,nnzptr) +end + +mutable struct InserterCache{Ti} + nnz :: Ti + nnzptr :: Vector{Ti} +end + +function Algebra.create_from_nz(a::ReusableInserterCSC{Tv,Ti}) where {Tv,Ti} + k = 1 + nnzperm = zeros(Ti,a.nnz) + for j in 1:a.ncols + pini = Int(a.colptr[j]) + pend = pini + Int(a.colnnz[j]) - 1 + for p in pini:pend + a.nzval[k] = a.nzval[p] + a.rowval[k] = a.rowval[p] + nnzperm[p] = k + k += 1 + end + end + @inbounds for j in 1:a.ncols + a.colptr[j+1] = a.colnnz[j] + end + length_to_ptrs!(a.colptr) + nnz = a.colptr[end]-1 + resize!(a.rowval,nnz) + resize!(a.nzval,nnz) + nnzptr = collect(lazy_map(Reindex(nnzperm),a.nnzptr)) + SparseMatrixCSC(a.nrows,a.ncols,a.colptr,a.rowval,a.nzval), InserterCache(a.nnz,nnzptr) +end + +struct CachedInserterCSC{Tv,Ti} + mat::SparseMatrixCSC{Tv,Ti} + cache::InserterCache{Ti} +end + +Algebra.LoopStyle(::Type{<:CachedInserterCSC}) = Loop() + +@noinline function Algebra.add_entry!(::typeof(+),a::CachedInserterCSC,v,i,j) + a.cache.nnz = a.cache.nnz + 1 + k = a.cache.nnzptr[a.cache.nnz] + a.mat.nzval[k] += v + nothing +end + +function FESpaces.assemble_matrix!(mat,a::SparseMatrixAssembler,matdata,cache) + LinearAlgebra.fillstored!(mat,zero(eltype(mat))) + cache.nnz = 0 + m = CachedInserterCSC(mat,cache) + numeric_loop_matrix!(m,a,matdata) + mat +end + +############################################################################################ + +model = CartesianDiscreteModel((0,1,0,1),(10,10)) + +reffe = ReferenceFE(lagrangian, Float64, 1) +V = TestFESpace(model, reffe) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 2) +a(u,v) = ∫(∇(u)⋅∇(v))dΩ + +u = get_trial_fe_basis(V) +v = get_fe_basis(V) +matdata = collect_cell_matrix(V,V,a(u,v)) + +assem = SparseMatrixAssembler(SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinReuse()),V,V) +A1, cache = assemble_matrix(assem,matdata) + +A2 = deepcopy(A1) +assemble_matrix!(A2,assem,matdata,cache) + +A1 ≈ A2 From 91fbe0d73141a16b42f1f27744eb7ae9819adf8d Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 12 Oct 2024 22:26:03 +1100 Subject: [PATCH 8/9] Updated NEWS.md --- NEWS.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 8ea78bcd1..98d2d211f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,13 +5,20 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added + +- Added MacroFElements. These are defined as having the basis/dof-basis of a FESpace created on top of a RefinementRule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024). +- Added Barycentric refinement rule in 2D and 3D. Added Simplexify refinement rule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024). + ## [0.18.6] - 2024-08-29 ### Fixed - Improved performance of PR[#967](https://github.com/gridap/Gridap.jl/pull/967). Along the way, opened the door to Triangulations of different type in SkeletonTriangulation. Since PR[#1026](https://github.com/gridap/Gridap.jl/pull/1026). -## [0.18.5] - 2024-08-28 +## [0.18.5] - 2024-08-28 ### Changed From 309221bd097e483e42613f2b7d29a7abc614a2b8 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 14 Oct 2024 10:47:48 +1100 Subject: [PATCH 9/9] Cleanup --- test/3d_eval.jl | 41 --------- test/integration_benchmark.jl | 46 ---------- test/kron_vector.jl | 44 --------- test/mytest.jl | 92 ------------------- test/mytest2.jl | 52 ----------- test/raviart_thomas.jl | 57 ------------ test/reusable_inserters.jl | 165 ---------------------------------- 7 files changed, 497 deletions(-) delete mode 100644 test/3d_eval.jl delete mode 100644 test/integration_benchmark.jl delete mode 100644 test/kron_vector.jl delete mode 100644 test/mytest.jl delete mode 100644 test/mytest2.jl delete mode 100644 test/raviart_thomas.jl delete mode 100644 test/reusable_inserters.jl diff --git a/test/3d_eval.jl b/test/3d_eval.jl deleted file mode 100644 index 0640db9aa..000000000 --- a/test/3d_eval.jl +++ /dev/null @@ -1,41 +0,0 @@ - -using Gridap -using Gridap.FESpaces, Gridap.CellData, Gridap.Fields, Gridap.Arrays - -using FillArrays - -# Function that goes from 2D points to 2-component vectors -u2D(x::VectorValue{2}) = VectorValue(x[1],x[2]) - -function point_3D_to_2D(x::VectorValue{3}) - return VectorValue(x[1],x[2]) -end - -function result_2D_to_3D(v::VectorValue{2}) - return VectorValue(v[1],v[2],0.0) -end - -function embed_3D(u2D::CellField,trian_3D::Triangulation) - u22_data = CellData.get_data(u2D) - ncells = length(u22_data) - - u32_data = lazy_map(Broadcasting(∘),u22_data,Fill(GenericField(point_3D_to_2D),ncells)) - u33_data = lazy_map(Broadcasting(∘),Fill(GenericField(result_2D_to_3D),ncells),u32_data) - u3D = CellData.similar_cell_field(u2D,u33_data,trian_3D,PhysicalDomain()) - return u3D -end - -model_2D = CartesianDiscreteModel((0,1,0,1),(2,2)) -model_3D = CartesianDiscreteModel((0,1,0,1,0,1),(2,2,2)) - -reffe_2D = ReferenceFE(lagrangian,VectorValue{2,Float64},1) -reffe_3D = ReferenceFE(lagrangian,VectorValue{3,Float64},1) - -V_2D = TestFESpace(model_2D,reffe_2D) -V_3D = TestFESpace(model_3D,reffe_3D) - -uh_2D = interpolate(u2D,V_2D) -uh_3D = embed_3D(uh_2D,get_triangulation(V_3D)) - -pts = get_cell_points(get_triangulation(V_3D)) -uh_3D(pts) diff --git a/test/integration_benchmark.jl b/test/integration_benchmark.jl deleted file mode 100644 index 2e3932d4d..000000000 --- a/test/integration_benchmark.jl +++ /dev/null @@ -1,46 +0,0 @@ - -using BenchmarkTools, ProfileView -using SparseArrays, LinearAlgebra -using Gridap, Gridap.FESpaces, Gridap.Algebra - -model = CartesianDiscreteModel((0,1,0,1),(100,100)) - -reffe_u = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1) -reffe_p = ReferenceFE(lagrangian, Float64, 1) -reffe_j = ReferenceFE(raviart_thomas, Float64, 0) -U = TestFESpace(model, reffe_u) -P = TestFESpace(model, reffe_p) -J = TestFESpace(model, reffe_j) - -X = MultiFieldFESpace([U,U]) - -Ω = Triangulation(model) -dΩ = Measure(Ω, 2) - -Γ = Boundary(model) -dΓ = Measure(Γ, 2) - -a1(u,v) = ∫(∇(u)⊙∇(v))dΩ -A1 = assemble_matrix(a1,U,U); -@benchmark assemble_matrix!($a1,$A1,$U,$U) - -a2(u,v) = ∫(∇(u)⊙∇(v))dΩ + ∫(u⋅v)dΓ -A2 = assemble_matrix(a2,U,U); -@benchmark assemble_matrix!($a2,$A2,$U,$U) - -a3(u,v) = ∫(divergence(u)⊙divergence(v))dΩ -A3 = assemble_matrix(a3,U,U); -@benchmark assemble_matrix!($a3,$A3,$U,$U) - -a4(u,v) = ∫(divergence(u)⊙divergence(v))dΩ -A4 = assemble_matrix(a4,J,J); -@benchmark assemble_matrix!($a4,$A4,$J,$J) - -a5((u1,u2),(v1,v2)) = ∫(∇(u1)⊙∇(v1))dΩ + ∫(∇(u2)⊙∇(v2))dΩ -A5 = assemble_matrix(a5,X,X); -@benchmark assemble_matrix!($a5,$A5,$X,$X) - -ProfileView.@profview @benchmark assemble_matrix!($a1,$A1,$U,$U) -ProfileView.@profview @benchmark assemble_matrix!($a2,$A2,$U,$U) -ProfileView.@profview @benchmark assemble_matrix!($a3,$A3,$U,$U) -ProfileView.@profview @benchmark assemble_matrix!($a4,$A4,$J,$J) diff --git a/test/kron_vector.jl b/test/kron_vector.jl deleted file mode 100644 index 3d88c9a87..000000000 --- a/test/kron_vector.jl +++ /dev/null @@ -1,44 +0,0 @@ - -using Gridap -using Gridap.Arrays, Gridap.Geometry, Gridap.FESpaces - -model = CartesianDiscreteModel((0,1,0,1),(5,5)) - -reffe = ReferenceFE(lagrangian,Float64,1) -V = FESpace(model,reffe) - -Ω = Triangulation(model) -dΩ = Measure(Ω,2) - -Γ = BoundaryTriangulation(model) -dΓ = Measure(Γ,2) - -a(u) = ∫(u)dΓ - -u = get_trial_fe_basis(V) -v = get_fe_basis(V) - -cu = get_array(a(u)) -cv = get_array(a(v)) - -nu = length(cu) -nv = length(cv) -cu_ids = vcat([collect(1:nu) for i in 1:nv]...) # Fast index -cv_ids = vcat([fill(i,nu) for i in 1:nv]...) # Slow index -dofs = get_cell_dof_ids(V) - -cu_exp = CompressedArray(cu,cu_ids) -cv_exp = CompressedArray(cv,cv_ids) - -uids_exp = CompressedArray(dofs,cu_ids) -vids_exp = CompressedArray(dofs,cv_ids) - -δ = fill(1.0,nv*nu) - -f(δ,v1,v2) = δ .* (v1 .- v2).^2 -arr = lazy_map(f,δ,cv_exp,cu_exp) -matdata = ([arr],[vids_exp],[uids_exp]) - -assem = SparseMatrixAssembler(V,V) -assemble_matrix(assem,matdata) - diff --git a/test/mytest.jl b/test/mytest.jl deleted file mode 100644 index 9ac0964a0..000000000 --- a/test/mytest.jl +++ /dev/null @@ -1,92 +0,0 @@ - -using BenchmarkTools -using SparseArrays, LinearAlgebra -using Gridap, Gridap.FESpaces, Gridap.Algebra -using Gridap.CellData, Gridap.Arrays, Gridap.Fields - -model = CartesianDiscreteModel((0,1,0,1),(2,2)) - -order = 1 -reffe = ReferenceFE(lagrangian, Float64, order) -V = TestFESpace(model, reffe) -X = MultiFieldFESpace([V,V]) - -Ω = Triangulation(model) -Γ = BoundaryTriangulation(model) -Λ = SkeletonTriangulation(model) - -degree = 2*order -dΩ = Measure(Ω,degree) -dΓ = Measure(Γ,degree) -dΛ = Measure(Λ,degree) - -n_Γ = get_normal_vector(Γ) -n_Λ = get_normal_vector(Λ) -γ = 1.0 -γ0 = 1.0 -h = 1.0 -f = 1.0 -g = 10.0 - -a((u,p),(v,q)) = - ∫( ∇(v)⊙∇(u) - ∇(q)⋅u + v⋅∇(p) )*dΩ + - ∫( (γ/h)*v⋅u - v⋅(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))⋅u + 2*(q*n_Γ)⋅u )*dΓ + - ∫( - (γ/h)*jump(v⊗n_Λ)⊙jump(u⊗n_Λ) - - jump(v⊗n_Λ)⊙mean(∇(u)) - - mean(∇(v))⊙jump(u⊗n_Λ) + - (γ0*h)*jump(q*n_Λ)⋅jump(p*n_Λ) + - jump(q*n_Λ)⋅mean(u) - - mean(v)⋅jump(p*n_Λ) - )*dΛ - -l((v,q)) = - ∫( v⋅f + q*g )*dΩ + - ∫( (γ/h)*v⋅u - (n_Γ⋅∇(v))⋅u + (q*n_Γ)⋅u )*dΓ - -op = AffineFEOperator(a,l,X,X) - - -u, p = get_trial_fe_basis(X) -v, q = get_fe_basis(X) - -u = get_trial_fe_basis(V) -v = get_fe_basis(V) - -cf1 = jump(v⊗n_Λ) -cf2 = jump(u⊗n_Λ) -cf3 = cf1 ⊙ cf2 - -cf4 = mean(∇(u)) -cf5 = cf1 ⊙ cf4 - -cf1_data = testitem(get_data(cf1)) -cf4_data = testitem(get_data(cf4)) - -_cf4 = change_domain(cf4,Λ,ReferenceDomain()) -_cf4_data = testitem(get_data(_cf4)) -return_value(Broadcasting(Operation(inner)),cf1_data,_cf4_data) - -get_data(cf5) - -args = map(get_data,cf5.args) -lazy_map(Broadcasting(cf5.op),args...) - -pts = CellData._get_cell_points(cf1,cf2) -x = testitem(get_data(pts)) -f1 = testitem(get_data(cf1)) -f2 = testitem(get_data(cf2)) - -fx1 = return_value(f1,x) -fx2 = return_value(f2,x) -r = Fields.BroadcastingFieldOpMap(⊙)(fx1,fx2) - -cf4 = v⊗n_Λ -cf4_p = testitem(get_data(cf4.plus)) -cf4_m = testitem(get_data(cf4.minus)) - -get_data(jump(n_Λ)) - -fkx = map(cfk -> evaluate(testitem(get_data(cfk)),x),cf1.args) -r = Fields.BroadcastingFieldOpMap(+)(fkx...) - diff --git a/test/mytest2.jl b/test/mytest2.jl deleted file mode 100644 index f3aecda97..000000000 --- a/test/mytest2.jl +++ /dev/null @@ -1,52 +0,0 @@ - -using SparseArrays, LinearAlgebra -using Gridap, Gridap.FESpaces, Gridap.Algebra, Gridap.Adaptivity - -using GridapDistributed, PartitionedArrays - -sol(x) = x - -np = (2,2) -ranks = with_debug() do distribute - distribute(LinearIndices((prod(np),))) -end - -cmodel = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4)) -fmodel = refine(cmodel) - -reffe_u = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1) -reffe_j = ReferenceFE(raviart_thomas, Float64, 0) - -Vh = TestFESpace(fmodel, reffe_u, dirichlet_tags="boundary") -Dh = TestFESpace(fmodel, reffe_j, dirichlet_tags="boundary") -Uh = TrialFESpace(Vh, sol) -Jh = TrialFESpace(Dh, sol) - -VH = TestFESpace(cmodel, reffe_u, dirichlet_tags="boundary") -DH = TestFESpace(cmodel, reffe_j, dirichlet_tags="boundary") -UH = TrialFESpace(VH, sol) -JH = TrialFESpace(DH, sol) - -f(x) = VectorValue(x[1]-x[2], x[1]+x[2]) - -Xh,XH,Yh,YH = Uh,UH,Vh,VH -#Xh,XH,Yh,YH = Jh,JH,Dh,DH - -fv_h = zero_free_values(Xh) -fv_H = zero_free_values(XH) - -dv_h = zero_dirichlet_values(Xh) -dv_H = zero_dirichlet_values(XH) - -fill!(fv_H,1.0) -uH = FEFunction(XH,fv_H,dv_H) -uh = interpolate!(uH,fv_h,Xh) - -fv_h_bis = zero_free_values(Yh) -fv_H_bis = zero_free_values(YH) - -fill!(fv_H_bis,1.0) -vH = FEFunction(XH,fv_H_bis,dv_H) -vh = interpolate!(uH,fv_h_bis,Xh) - -fv_H_bis == fv_H diff --git a/test/raviart_thomas.jl b/test/raviart_thomas.jl deleted file mode 100644 index ec7906842..000000000 --- a/test/raviart_thomas.jl +++ /dev/null @@ -1,57 +0,0 @@ -using Gridap -using Gridap.ReferenceFEs, Gridap.Geometry, Gridap.FESpaces, Gridap.Arrays -using Gridap.Fields - -using BenchmarkTools - -using Profile -using ProfileView - -using FillArrays - -model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(40,40))) - - -order = 2 -reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) -reffe_RT = ReferenceFE(raviart_thomas,Float64,order-1) - -V = FESpace(model,reffe) -V_RT = FESpace(model,reffe_RT) - -Ω = Triangulation(model) -dΩ = Measure(Ω,2*order) - -a(u,v) = ∫(u⋅v)dΩ -b(u,v) = ∫((∇⋅u)*(∇⋅v))dΩ -c(u,v) = ∫(∇(u)⊙∇(v))dΩ - -A1 = assemble_matrix(a,V,V) -A2 = assemble_matrix(a,V_RT,V_RT) - -B1 = assemble_matrix(b,V,V) -B2 = assemble_matrix(b,V_RT,V_RT) - -C1 = assemble_matrix(c,V,V) -C2 = assemble_matrix(c,V_RT,V_RT) - -@benchmark assemble_matrix!(a,A1,V,V) -@benchmark assemble_matrix!(a,A2,V_RT,V_RT) -@benchmark assemble_matrix!(b,B1,V,V) -@benchmark assemble_matrix!(b,B2,V_RT,V_RT) -@benchmark assemble_matrix!(c,C1,V,V) -@benchmark assemble_matrix!(c,C2,V_RT,V_RT) - -function p1() - for i in 1:100 - assemble_matrix!(b,B2,V_RT,V_RT) - end -end -ProfileView.@profview p1() - - -mon = V.fe_basis.cell_basis.values[1].fields -hdiv = V_RT.fe_basis.cell_basis.args[1].values[1].fields.qgrad - - - diff --git a/test/reusable_inserters.jl b/test/reusable_inserters.jl deleted file mode 100644 index 4922accd7..000000000 --- a/test/reusable_inserters.jl +++ /dev/null @@ -1,165 +0,0 @@ - -using SparseArrays, LinearAlgebra -using Gridap, Gridap.Algebra, Gridap.FESpaces, Gridap.Helpers - -struct MinReuse end - -struct ReusableCounterCSC{Tv,Ti} - tv::Type{Tv} - nrows::Int - ncols::Int - colnnzmax::Vector{Ti} -end - -Algebra.LoopStyle(::Type{<:ReusableCounterCSC}) = Loop() - -@inline function Algebra.add_entry!(::Function,a::ReusableCounterCSC{Tv,Ti},v,i,j) where {Tv,Ti} - a.colnnzmax[j] += Ti(1) - nothing -end - -mutable struct ReusableInserterCSC{Tv,Ti} - nrows::Int - ncols::Int - colptr::Vector{Ti} - colnnz::Vector{Ti} - rowval::Vector{Ti} - nzval::Vector{Tv} - nnz::Ti - nnzptr::Vector{Ti} -end - -Algebra.LoopStyle(::Type{<:ReusableInserterCSC}) = Loop() - -@noinline function Algebra.add_entry!(::typeof(+),a::ReusableInserterCSC{Tv,Ti},v,i,j) where {Tv,Ti} - pini = Int(a.colptr[j]) - pend = pini + Int(a.colnnz[j]) - 1 - p = searchsortedfirst(a.rowval,i,pini,pend,Base.Order.Forward) - if (p>pend) - # add new entry - a.colnnz[j] += 1 - a.rowval[p] = i - a.nzval[p] = v - elseif a.rowval[p] != i - # shift one forward from p to pend - @check pend+1 < Int(a.colptr[j+1]) - for k in pend:-1:p - o = k + 1 - a.rowval[o] = a.rowval[k] - a.nzval[o] = a.nzval[k] - end - # add new entry - a.colnnz[j] += 1 - a.rowval[p] = i - a.nzval[p] = v - else - # update existing entry - a.nzval[p] += v - end - # Save the position of the entry for reuse - a.nnz += 1 - a.nnzptr[a.nnz] = p - nothing -end - -############################################################################################ - -function Algebra.nz_counter( - ::SparseMatrixBuilder{SparseMatrixCSC{Tv,Ti},<:MinReuse}, - axes -) where {Tv,Ti} - nrows = length(axes[1]) - ncols = length(axes[2]) - colnnzmax = zeros(Ti,ncols) - ReusableCounterCSC(Tv,nrows,ncols,colnnzmax) -end - -function Algebra.nz_allocation(a::ReusableCounterCSC{Tv,Ti}) where {Tv,Ti} - colptr = Vector{Ti}(undef,a.ncols+1) - @inbounds for i in 1:a.ncols - colptr[i+1] = a.colnnzmax[i] - end - length_to_ptrs!(colptr) - ndata = colptr[end] - one(Ti) - rowval = Vector{Ti}(undef,ndata) - nzval = zeros(Tv,ndata) - nnzmax = sum(a.colnnzmax) - nnzptr = fill(zero(Ti),nnzmax) - colnnz = a.colnnzmax - fill!(colnnz,zero(Ti)) - ReusableInserterCSC(a.nrows,a.ncols,colptr,colnnz,rowval,nzval,0,nnzptr) -end - -mutable struct InserterCache{Ti} - nnz :: Ti - nnzptr :: Vector{Ti} -end - -function Algebra.create_from_nz(a::ReusableInserterCSC{Tv,Ti}) where {Tv,Ti} - k = 1 - nnzperm = zeros(Ti,a.nnz) - for j in 1:a.ncols - pini = Int(a.colptr[j]) - pend = pini + Int(a.colnnz[j]) - 1 - for p in pini:pend - a.nzval[k] = a.nzval[p] - a.rowval[k] = a.rowval[p] - nnzperm[p] = k - k += 1 - end - end - @inbounds for j in 1:a.ncols - a.colptr[j+1] = a.colnnz[j] - end - length_to_ptrs!(a.colptr) - nnz = a.colptr[end]-1 - resize!(a.rowval,nnz) - resize!(a.nzval,nnz) - nnzptr = collect(lazy_map(Reindex(nnzperm),a.nnzptr)) - SparseMatrixCSC(a.nrows,a.ncols,a.colptr,a.rowval,a.nzval), InserterCache(a.nnz,nnzptr) -end - -struct CachedInserterCSC{Tv,Ti} - mat::SparseMatrixCSC{Tv,Ti} - cache::InserterCache{Ti} -end - -Algebra.LoopStyle(::Type{<:CachedInserterCSC}) = Loop() - -@noinline function Algebra.add_entry!(::typeof(+),a::CachedInserterCSC,v,i,j) - a.cache.nnz = a.cache.nnz + 1 - k = a.cache.nnzptr[a.cache.nnz] - a.mat.nzval[k] += v - nothing -end - -function FESpaces.assemble_matrix!(mat,a::SparseMatrixAssembler,matdata,cache) - LinearAlgebra.fillstored!(mat,zero(eltype(mat))) - cache.nnz = 0 - m = CachedInserterCSC(mat,cache) - numeric_loop_matrix!(m,a,matdata) - mat -end - -############################################################################################ - -model = CartesianDiscreteModel((0,1,0,1),(10,10)) - -reffe = ReferenceFE(lagrangian, Float64, 1) -V = TestFESpace(model, reffe) - -Ω = Triangulation(model) -dΩ = Measure(Ω, 2) -a(u,v) = ∫(∇(u)⋅∇(v))dΩ - -u = get_trial_fe_basis(V) -v = get_fe_basis(V) -matdata = collect_cell_matrix(V,V,a(u,v)) - -assem = SparseMatrixAssembler(SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinReuse()),V,V) -A1, cache = assemble_matrix(assem,matdata) - -A2 = deepcopy(A1) -assemble_matrix!(A2,assem,matdata,cache) - -A1 ≈ A2