diff --git a/src/GridapAPIExtensions.jl b/src/GridapAPIExtensions.jl index f1ba7ca..039d473 100644 --- a/src/GridapAPIExtensions.jl +++ b/src/GridapAPIExtensions.jl @@ -933,7 +933,7 @@ function Gridap.Arrays.evaluate!(cache, Gridap.Fields.TransposeFieldIndices(gx) end -# + # When we do local operation like uK_u∂K = Reconstruction_operator(trial::MultiFieldCellField) where, # uK_u∂K is associated with a different object id for example, a view of the triangulation or # distributed setup, then cell_field operations like vK*ur∂K break, where vK, v∂K = test::MultiFieldCellField. diff --git a/src/HybridAffineFEOperators.jl b/src/HybridAffineFEOperators.jl index 3019397..f88596c 100644 --- a/src/HybridAffineFEOperators.jl +++ b/src/HybridAffineFEOperators.jl @@ -350,8 +350,8 @@ function _find_skeleton(dc::DomainContribution) end function _find_bulk(D, dc::DomainContribution) - [trian for trian in keys(dc.dict) - if isa(trian, Triangulation{D,D}) && !(isa(trian, SkeletonTriangulation))] + [trian for trian in keys(dc.dict) + if ( isa(trian, Triangulation{D}) && ( !(isa(trian, SkeletonTriangulation)) || !(isa(trian, SkeletonView)) ) )] end function _find_boundary(dc::DomainContribution) @@ -363,7 +363,7 @@ function _add_static_condensation(matvec, bulk_fields, skeleton_fields) Gridap.Helpers.@check length(keys(matvec.dict)) == 1 _matvec = DomainContribution() for (trian, t) in matvec.dict - Gridap.Helpers.@check isa(trian, SkeletonTriangulation) + Gridap.Helpers.@check isa(trian, SkeletonTriangulation) || isa(trian, GridapHybrid.SkeletonView) _matvec.dict[trian] = lazy_map(StaticCondensationMap(bulk_fields, skeleton_fields), t) end _matvec diff --git a/src/ProjectionFEOperators.jl b/src/ProjectionFEOperators.jl index 57aef08..ee68018 100644 --- a/src/ProjectionFEOperators.jl +++ b/src/ProjectionFEOperators.jl @@ -28,7 +28,9 @@ # end # end -# Adding the flexibility to pass an operator to the constructor. +# Adding the flexibility to pass an operator to the constructor. Without this change, explicitly computng +# (in the distributed setting) xK_x∂K = P(xh), where xh is the basis associaed with the hybrid space +# and P is the projection operator will not work. struct ProjectionFEOperator{T1,T2,T3,T4,T5} <: GridapType LHS_form :: T1 RHS_form :: T2 diff --git a/test/PoissonHHODistributedTestsSerial.jl b/test/PoissonHHODistributedTestsSerial.jl new file mode 100644 index 0000000..3ecf310 --- /dev/null +++ b/test/PoissonHHODistributedTestsSerial.jl @@ -0,0 +1,294 @@ +module PoissonHHODistributedTests + using Gridap + using Gridap.Geometry, Gridap.ReferenceFEs, Gridap.Arrays, Gridap.CellData, Gridap.FESpaces + using GridapHybrid + using GridapDistributed + using PartitionedArrays + using Plots + + + function setup_reconstruction_operator(model, order, dΩ, d∂K) + nK = get_cell_normal_vector(d∂K.quad.trian) + refferecᵤ = ReferenceFE(orthogonal_basis, Float64, order+1) + reffe_c = ReferenceFE(monomial_basis , Float64, order+1; subspace=:OnlyConstant) + + Ω = dΩ.quad.trian + + VKR = TestFESpace(Ω , refferecᵤ; conformity=:L2) + UKR = TrialFESpace(VKR) + VKR_C = TestFESpace(Ω, reffe_c ; conformity=:L2, vector_type=Vector{Float64}) + UKR_C = TrialFESpace(VKR_C) + + V = MultiFieldFESpace([VKR,VKR_C]) + U = MultiFieldFESpace([UKR,UKR_C]) + + m( (u,u_c), (v,v_c) ) = ∫(∇(v)⋅∇(u))dΩ + ∫(v_c*u)dΩ + ∫(v*u_c)dΩ + n( (uK,u∂K), (v,v_c) ) = ∫(∇(v)⋅∇(uK))dΩ + ∫(v_c*uK)dΩ + ∫((∇(v)⋅nK)*u∂K)d∂K - ∫((∇(v)⋅nK)*uK)d∂K + + ReconstructionFEOperator((m,n), U, V) + end + + + function setup_projection_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K) + m( (uK,u∂K), (vK,v∂K) ) = ∫(vK*uK)dΩ + ∫(v∂K*u∂K)d∂K + function n(uK_u∂K, (vK,v∂K)) + urK_ur∂K = R(uK_u∂K) + urK,ur∂K = urK_ur∂K + uK,u∂K = uK_u∂K + # ∫(vK*(urK-uK))dΩ+∫(vK*ur∂K)dΩ - # bulk projection terms + # ∫(v∂K*urK)d∂K+∫(v∂K*u∂K)d∂K-∫(v∂K*ur∂K)d∂K # skeleton projection terms + ∫(vK*(urK-uK))dΩ+∫(vK*ur∂K)dΩ + # bulk projection terms + ∫(v∂K*urK)d∂K-∫(v∂K*u∂K)d∂K+∫(v∂K*ur∂K)d∂K # skeleton projection terms + end + function n(uK_u∂K::Gridap.MultiField.MultiFieldFEFunction, (vK,v∂K)) + uK,u∂K = uK_u∂K + urK = R(uK_u∂K) + ∫(vK*(urK-uK))dΩ + # bulk projection terms + ∫(v∂K*(urK-u∂K))d∂K # skeleton projection terms + end + GridapHybrid.ProjectionFEOperator((m,n),UK_U∂K,VK_V∂K,R) + end + + + + # Compute length of diagonals in the reference domain + function foo(m) + p0 = evaluate(m,Point(0.0,0.0)) + p1 = evaluate(m,Point(1.0,1.0)) + norm(p1-p0) + # To Do: incorporate the other diagonal + # p2 = evaluate(m,Point(1.0,0.0)) + # p3 = evaluate(m,Point(0.0,1.0)) + # maximum( norm(p1-p0), norm(p2-p3) ) + end + + + function setup_reconstruction_operator(model::GridapDistributed.GenericDistributedDiscreteModel, + order, dΩ::GridapDistributed.DistributedMeasure, d∂K::GridapDistributed.DistributedMeasure) + R = map(local_views(model), local_views(dΩ), local_views(d∂K)) do m, dΩi, d∂Ki + setup_reconstruction_operator(m, order, dΩi, d∂Ki) + end + return R + end + + function setup_projection_operator(trial::GridapDistributed.DistributedMultiFieldFESpace, + test::GridapDistributed.DistributedMultiFieldFESpace, R::AbstractArray{<:ReconstructionFEOperator}, + dΩ::GridapDistributed.DistributedMeasure, d∂K::GridapDistributed.DistributedMeasure) + P = map(local_views(trial), local_views(test), local_views(R), local_views(dΩ), local_views(d∂K)) do triali, testi, Ri, dΩi, d∂Ki + setup_projection_operator(triali, testi, Ri, dΩi, d∂Ki) + end + return P + end + + function r(u,v,R,dΩ) + uK_u∂K=R(u) + vK_v∂K=R(v) + uK,u∂K = uK_u∂K + vK,v∂K = vK_v∂K + ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + + ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ + end + + function s(u,v,P,d∂K) + # Currently, we cannot use this CellField in the integrand of the skeleton integrals below. + # we think we need to develop a specific version for _transform_face_to_cell_lface_expanded_array + # we will be using h_T_1 in the meantime + + # Ω = get_triangulation(u) + # cmaps = collect1d(get_cell_map(Ω)) + # h_array = lazy_map(foo,cmaps) + # # h_T = CellField(h_array,Ω,PhysicalDomain()) + # h_T = CellField(h_array,Ω) + # h_T_1 = 1.0/h_T + + lΩ = get_triangulation(u) + dlΩ = Measure(lΩ, 3) + h_T=CellField(get_array(∫(1)dlΩ),lΩ) + h_T_1=1.0/h_T + + uK_u∂K_ΠK,uK_u∂K_Π∂K=P(u) + vK_v∂K_ΠK,vK_v∂K_Π∂K=P(v) + + uK_ΠK , u∂K_ΠK = uK_u∂K_ΠK + uK_Π∂K , u∂K_Π∂K = uK_u∂K_Π∂K + + vK_ΠK , v∂K_ΠK = vK_v∂K_ΠK + vK_Π∂K , v∂K_Π∂K = vK_v∂K_Π∂K + + vK_Π∂K_vK_ΠK=vK_Π∂K-vK_ΠK + v∂K_Π∂K_v∂K_ΠK=v∂K_Π∂K-v∂K_ΠK + uK_Π∂K_uK_ΠK=uK_Π∂K-uK_ΠK + u∂K_Π∂K_u∂K_ΠK=u∂K_Π∂K-u∂K_ΠK + + ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(uK_Π∂K_uK_ΠK))d∂K + + ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + + ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + + ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(uK_Π∂K_uK_ΠK))d∂K + end + + function r(u::GridapDistributed.DistributedMultiFieldCellField, + v::GridapDistributed.DistributedMultiFieldCellField, + R::AbstractArray{<:ReconstructionFEOperator}, + dΩ::GridapDistributed.DistributedMeasure) + a_r = map(local_views(u), local_views(v), local_views(R), local_views(dΩ)) do ui, vi, Ri, dΩi + a_r = r(ui,vi,Ri,dΩi) + return a_r + end + end + + function s(u::GridapDistributed.DistributedMultiFieldCellField, + v::GridapDistributed.DistributedMultiFieldCellField, + P::AbstractArray{<:ProjectionFEOperator}, + d∂K::GridapDistributed.DistributedMeasure) + a_s = map(local_views(u), local_views(v), local_views(P), local_views(d∂K)) do ui, vi, Pi, d∂Ki + a_s = s(ui,vi,Pi,d∂Ki) + return a_s + end + end + + # ranks = with_debug() do distribute + # distribute(LinearIndices((prod(np),))) + # end + + function solve_hho(domain, mesh_partition, np, ranks, order, u) + model = CartesianDiscreteModel(ranks, np, domain, mesh_partition) + + D = num_cell_dims(model) + Ω = Triangulation(ReferenceFE{D},model) + Γ = Triangulation(ReferenceFE{D-1},model) + ∂K= GridapHybrid.Skeleton(model) + + # Local Triangulations to be used for local problems + ∂K_loc = GridapDistributed.DistributedTriangulation( + map(t -> t.parent,local_views(∂K)), ∂K.model) + Ω_loc = GridapDistributed.add_ghost_cells(Ω) + + reffe = ReferenceFE(lagrangian, Float64, order; space=:P) + VK = TestFESpace(Ω , reffe; conformity=:L2) + V∂K = TestFESpace(Γ , reffe; conformity=:L2, dirichlet_tags=collect(5:8)) + UK = TrialFESpace(VK) + U∂K= TrialFESpace(V∂K,u) + V = MultiFieldFESpace([VK,V∂K]) + U = MultiFieldFESpace([UK,U∂K]) + + degree = 2*(order+1) + dΩ = Measure(Ω, degree) + d∂K= Measure(∂K, degree) + d∂K_loc = Measure(∂K_loc, degree) + dΩ_loc = Measure(Ω_loc, degree) + # Ω_loc = get_triangulation(VK) # alternative way to define Ω_loc which ensures that the object id is same as that of VK. + # ndofs= num_free_dofs(U) + + R = setup_reconstruction_operator(model, order, dΩ_loc, d∂K_loc) + P = setup_projection_operator(U, V, R, dΩ_loc, d∂K_loc) + + a(u,v,dΩ,d∂K) = map(+, local_views(r(u,v,R,dΩ)), local_views(s(u,v,P,d∂K))) + l((vK,),dΩ) = ∫(vK*f)dΩ + + weakform = (u,v)->(a(u,v,dΩ,d∂K),l(v,dΩ)) + op=HybridAffineFEOperator(weakform, U, V, [1], [2]) + + xh = solve(op); + uhK,uh∂K = xh + e = uhK-u + + return sqrt(sum(∫(e⋅e)dΩ)) + end + + function convg_test(domain,ns,np,order,u) + + ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) + end + + el2 = Float64[] + hs = Float64[] + for n in ns + l2 = solve_hho(domain,(n,n),np,ranks,order,u) + println(l2) + push!(el2,l2) + h = 1/n + push!(hs,h) + end + println(el2) + println(hs) + el2, hs + end + + function slope(hs,errors) + x = log10.(hs) + y = log10.(errors) + linreg = hcat(fill!(similar(x), 1), x) \ y + linreg[2] + end + + + u(x) = sin(2π*x[1])*sin(2π*x[2]) + f(x) = -Δ(u)(x) + + domain = (0,1,0,1) + + np = (2,2) + ns = [np[1]*2^i for i in 2:5] + # mesh_partition = (2,2) + + order = 1 + el, hs = convg_test(domain,ns,np,order,u) + println("Slope L2-norm u: $(slope(hs,el))") + slopekp1=[Float64(ni)^(-(order+1)) for ni in ns] + display(plot(hs,[el slopekp1], + xaxis=:log, yaxis=:log, + label=["L2u (measured)" "slope k+1"], + shape=:auto, + xlabel="h",ylabel="L2 error",legend=:bottomright)) + +end + +# xh = get_trial_fe_basis(U); +# yh = get_fe_basis(V); + +# xh1 = local_views(xh).items[1] +# yh1 = local_views(yh).items[1] + +# Ω1 = get_triangulation(xh1) + +# cmaps = collect1d(get_cell_map(Ω1)) +# h_arr = lazy_map(foo,cmaps) +# # h_T = CellField(h_array,Ω,PhysicalDomain()) +# H_T = CellField(h_arr,Ω1) +# H_T_1 = 1.0/H_T + +# Ω1_loc = local_views(Ω_loc).items[1] +# dΩ1_loc = local_views(dΩ_loc).items[1] + +# h_T=CellField(get_array(∫(1)dΩ1_loc),Ω1_loc) +# h_T_1=1.0/h_T + + +# P1 = local_views(P).items[1] + +# uK_u∂K_ΠK,uK_u∂K_Π∂K=P1(xh1) +# vK_v∂K_ΠK,vK_v∂K_Π∂K=P1(yh1) + +# uK_ΠK , u∂K_ΠK = uK_u∂K_ΠK +# uK_Π∂K , u∂K_Π∂K = uK_u∂K_Π∂K + +# vK_ΠK , v∂K_ΠK = vK_v∂K_ΠK +# vK_Π∂K , v∂K_Π∂K = vK_v∂K_Π∂K + +# vK_Π∂K_vK_ΠK=vK_Π∂K-vK_ΠK +# v∂K_Π∂K_v∂K_ΠK=v∂K_Π∂K-v∂K_ΠK +# uK_Π∂K_uK_ΠK=uK_Π∂K-uK_ΠK +# u∂K_Π∂K_u∂K_ΠK=u∂K_Π∂K-u∂K_ΠK + +# d∂K_loc1 = local_views(d∂K_loc).items[1] + +# (h_T_1*(vK_Π∂K_vK_ΠK)*(uK_Π∂K_uK_ΠK)) + +# ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(uK_Π∂K_uK_ΠK))d∂K_loc1 + +# ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K_loc1 + +# ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K_loc1 + +# ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(uK_Π∂K_uK_ΠK))d∂K_loc1 + + + diff --git a/test/_dev/Distributed_testsC.jl b/test/_dev/Distributed_testsC.jl deleted file mode 100644 index 97ce437..0000000 --- a/test/_dev/Distributed_testsC.jl +++ /dev/null @@ -1,213 +0,0 @@ -using Gridap -using Gridap.Geometry, Gridap.ReferenceFEs, Gridap.Arrays, Gridap.CellData, Gridap.FESpaces -using GridapHybrid -using GridapDistributed -using PartitionedArrays - -function setup_reconstruction_operator(model, order, dΩ, d∂K) - nK = get_cell_normal_vector(d∂K.quad.trian) - refferecᵤ = ReferenceFE(orthogonal_basis, Float64, order+1) - reffe_c = ReferenceFE(monomial_basis , Float64, order+1; subspace=:OnlyConstant) - - Ω = dΩ.quad.trian - - VKR = TestFESpace(Ω , refferecᵤ; conformity=:L2) - UKR = TrialFESpace(VKR) - VKR_C = TestFESpace(Ω, reffe_c ; conformity=:L2, vector_type=Vector{Float64}) - UKR_C = TrialFESpace(VKR_C) - - V = MultiFieldFESpace([VKR,VKR_C]) - U = MultiFieldFESpace([UKR,UKR_C]) - - m( (u,u_c), (v,v_c) ) = ∫(∇(v)⋅∇(u))dΩ + ∫(v_c*u)dΩ + ∫(v*u_c)dΩ - n( (uK,u∂K), (v,v_c) ) = ∫(∇(v)⋅∇(uK))dΩ + ∫(v_c*uK)dΩ + ∫((∇(v)⋅nK)*u∂K)d∂K - ∫((∇(v)⋅nK)*uK)d∂K - - ReconstructionFEOperator((m,n), U, V) -end - - -function setup_projection_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K) - m( (uK,u∂K), (vK,v∂K) ) = ∫(vK*uK)dΩ + ∫(v∂K*u∂K)d∂K - function n(uK_u∂K, (vK,v∂K)) - urK_ur∂K = R(uK_u∂K) - urK,ur∂K = urK_ur∂K - uK,u∂K = uK_u∂K - # ∫(vK*(urK-uK))dΩ+∫(vK*ur∂K)dΩ - # bulk projection terms - # ∫(v∂K*urK)d∂K+∫(v∂K*u∂K)d∂K-∫(v∂K*ur∂K)d∂K # skeleton projection terms - ∫(vK*(urK-uK))dΩ+∫(vK*ur∂K)dΩ + # bulk projection terms - ∫(v∂K*urK)d∂K-∫(v∂K*u∂K)d∂K+∫(v∂K*ur∂K)d∂K # skeleton projection terms - end - function n(uK_u∂K::Gridap.MultiField.MultiFieldFEFunction, (vK,v∂K)) - uK,u∂K = uK_u∂K - urK = R(uK_u∂K) - ∫(vK*(urK-uK))dΩ + # bulk projection terms - ∫(v∂K*(urK-u∂K))d∂K # skeleton projection terms - end - GridapHybrid.ProjectionFEOperator((m,n),UK_U∂K,VK_V∂K,R) -end - - - -# Compute length of diagonals in the reference domain -function foo(m) -p0 = evaluate(m,Point(0.0,0.0)) -p1 = evaluate(m,Point(1.0,1.0)) -norm(p1-p0) -# To Do: incorporate the other diagonal -# p2 = evaluate(m,Point(1.0,0.0)) -# p3 = evaluate(m,Point(0.0,1.0)) -# maximum( norm(p1-p0), norm(p2-p3) ) -end - -u(x) = x[1]^p+x[2]^p -f(x) = -Δ(u)(x) - - np = (2,2) - ranks = with_debug() do distribute - distribute(LinearIndices((prod(np),))) - end - - domain = (0,1,0,1) - mesh_partition = (2,2) - model = CartesianDiscreteModel(ranks,np,domain, mesh_partition) - - D = num_cell_dims(model) - Ω = Triangulation(ReferenceFE{D},model) - Γ = Triangulation(ReferenceFE{D-1},model) - ∂K= GridapHybrid.Skeleton(model) - - # Local Triangulations to be used for local problems - ∂K_loc = GridapDistributed.DistributedTriangulation( - map(t -> t.parent,local_views(∂K)), ∂K.model - ) - Ω_loc = GridapDistributed.add_ghost_cells(Ω) - - order = 0 - reffe = ReferenceFE(lagrangian, Float64, order; space=:P) - VK = TestFESpace(Ω , reffe; conformity=:L2) - V∂K = TestFESpace(Γ , reffe; conformity=:L2) - UK = TrialFESpace(VK) - U∂K= TrialFESpace(V∂K) - - # Ω_loc = get_triangulation(VK) # alternative way to define Ω_loc which ensures that the object id is same as that of VK. - - degree = 2*(order+1) - dΩ = Measure(Ω, degree) - d∂K= Measure(∂K, degree) - d∂K_loc = Measure(∂K_loc, degree) - dΩ_loc = Measure(Ω_loc, degree) - - V = MultiFieldFESpace([VK,V∂K]) - U = MultiFieldFESpace([UK,U∂K]) - - function setup_reconstruction_operator(model::GridapDistributed.GenericDistributedDiscreteModel, - order, dΩ::GridapDistributed.DistributedMeasure, d∂K::GridapDistributed.DistributedMeasure) - R = map(local_views(model), local_views(dΩ), local_views(d∂K)) do m, dΩi, d∂Ki - setup_reconstruction_operator(m, order, dΩi, d∂Ki) - end - return R - end - - function setup_projection_operator(trial::GridapDistributed.DistributedMultiFieldFESpace, - test::GridapDistributed.DistributedMultiFieldFESpace, R::AbstractArray{<:ReconstructionFEOperator}, - dΩ::GridapDistributed.DistributedMeasure, d∂K::GridapDistributed.DistributedMeasure) - P = map(local_views(trial), local_views(test), local_views(R), local_views(dΩ), local_views(d∂K)) do triali, testi, Ri, dΩi, d∂Ki - setup_projection_operator(triali, testi, Ri, dΩi, d∂Ki) - end - return P - end - - R = setup_reconstruction_operator(model, order, dΩ_loc, d∂K_loc) - P = setup_projection_operator(U, V, R, dΩ_loc, d∂K_loc) - - function r(u,v,R,dΩ) - uK_u∂K=R(u) - vK_v∂K=R(v) - uK,u∂K = uK_u∂K - vK,v∂K = vK_v∂K - ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + - ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ - end - - function r(u::GridapDistributed.DistributedMultiFieldCellField, - v::GridapDistributed.DistributedMultiFieldCellField, - R::AbstractArray{<:ReconstructionFEOperator}, - dΩ::GridapDistributed.DistributedMeasure) - consistency = map(local_views(u), local_views(v), local_views(R), local_views(dΩ)) do ui, vi, Ri, dΩi - consistency = r(ui,vi,Ri,dΩi) - return consistency - end - end - - - - - xh = get_trial_fe_basis(U); - yh = get_fe_basis(V); - - a_consistency = r(xh,yh,R,dΩ) - - - xK_x∂K_ΠK,xK_x∂K_Π∂K=P(xh) - xK_x∂K_ΠK_xK_x∂K_Π∂K=P(yh) - - P1 = local_views(P).items[1]; - xh1 = local_views(xh).items[1]; - - basis_style = GridapHybrid._get_basis_style(xh1) - LHSf, RHSf = GridapHybrid._evaluate_forms(P1, xh1) - cell_dofs= GridapHybrid._compute_cell_dofs(LHSf,RHSf) - O = P1.test_space - GridapHybrid._generate_image_space_span(P1,O,xh1,cell_dofs,basis_style) - - # function s(u,v,d∂K) - # Ω = get_triangulation(u) - # cmaps = collect1d(get_cell_map(Ω)) - # h_array = lazy_map(foo,cmaps) - # # h_T = CellField(h_array,Ω,PhysicalDomain()) - # h_T = CellField(h_array,Ω) - # h_T_1 = 1.0/h_T - - - # uK_u∂K_ΠK,uK_u∂K_Π∂K=P(u) - # vK_v∂K_ΠK,vK_v∂K_Π∂K=P(v) - - # uK_ΠK , u∂K_ΠK = uK_u∂K_ΠK - # uK_Π∂K , u∂K_Π∂K = uK_u∂K_Π∂K - - # vK_ΠK , v∂K_ΠK = vK_v∂K_ΠK - # vK_Π∂K , v∂K_Π∂K = vK_v∂K_Π∂K - - # vK_Π∂K_vK_ΠK=vK_Π∂K-vK_ΠK - # v∂K_Π∂K_v∂K_ΠK=v∂K_Π∂K-v∂K_ΠK - # uK_Π∂K_uK_ΠK=uK_Π∂K-uK_ΠK - # u∂K_Π∂K_u∂K_ΠK=u∂K_Π∂K-u∂K_ΠK - - # ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(uK_Π∂K_uK_ΠK))d∂K + - # ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + - # ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + - # ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(uK_Π∂K_uK_ΠK))d∂K - # end - - # a(u,v,dΩ,d∂K) = r(u,v,dΩ) + s(u,v,d∂K) - # l((vK,),dΩ) = ∫(vK*f)dΩ - - # weakform = (u,v)->(a(u,v,dΩ,d∂K),l(v,dΩ)) - # @enter op=HybridAffineFEOperator(weakform, U, V, [1], [2]) - - - -# return A, h, ndofsΩ -# end - -# function domain_operation(operation, l::DomainContribution, r::DomainContribution) -# @check get_domains(r) == get_domains(l) -# o = DomainContribution() -# for trian in get_domains(l) -# cl = get_contribution(l, trian) -# cr = get_contribution(r, trian) -# c = lazy_map(Broadcasting(operation), cl, cr) -# add_contribution!(o, trian, c) -# end -# return o -# end \ No newline at end of file diff --git a/test/_dev/Serial.jl b/test/_dev/Serial.jl index 717c765..4ce4210 100644 --- a/test/_dev/Serial.jl +++ b/test/_dev/Serial.jl @@ -4,6 +4,16 @@ using GridapHybrid using GridapDistributed using PartitionedArrays +function foo(m) + p0 = evaluate(m,Point(0.0,0.0)) + p1 = evaluate(m,Point(1.0,1.0)) + norm(p1-p0) + # To Do: incorporate the other diagonal + # p2 = evaluate(m,Point(1.0,0.0)) + # p3 = evaluate(m,Point(0.0,1.0)) + # maximum( norm(p1-p0), norm(p2-p3) ) +end + function setup_reconstruction_operator(model, order, dΩ, d∂K) nK = get_cell_normal_vector(d∂K.quad.trian) refferecᵤ = ReferenceFE(orthogonal_basis, Float64, order+1) @@ -60,6 +70,9 @@ norm(p1-p0) end +p = 3 +u(x) = x[1]^p+x[2]^p +f(x) = -Δ(u)(x) np = (2,2) # ranks = with_debug() do distribute @@ -79,9 +92,9 @@ end order = 0 reffe = ReferenceFE(lagrangian, Float64, order; space=:P) VK = TestFESpace(Ω , reffe; conformity=:L2) - V∂K = TestFESpace(Γ , reffe; conformity=:L2) + V∂K = TestFESpace(Γ , reffe; conformity=:L2, dirichlet_tags=collect(5:8)) UK = TrialFESpace(VK) - U∂K= TrialFESpace(V∂K) + U∂K= TrialFESpace(V∂K,u) degree = 2*(order+1) @@ -99,136 +112,108 @@ end # h_array = lazy_map(foo,cmaps) R = setup_reconstruction_operator(model, order, dΩ, d∂K) - P = setup_projection_operator(U,V,R,dΩ,d∂K) + projection_op = setup_projection_operator(U,V,R,dΩ,d∂K) xh = get_trial_fe_basis(U) yh = get_fe_basis(V) - - basis_style = GridapHybrid._get_basis_style(xh) - - LHSf, RHSf = GridapHybrid._evaluate_forms(P, xh) - - trial = GridapHybrid._to_trial_basis(xh) - test = GridapHybrid._get_test_fe_basis(P) - - urK_ur∂K = R(trial) - urK,ur∂K = urK_ur∂K - uK,u∂K = trial - vK,v∂K = test - # ∫(vK*(urK-uK))dΩ+∫(vK*ur∂K)dΩ - # bulk projection terms - # ∫(v∂K*urK)d∂K+∫(v∂K*u∂K)d∂K-∫(v∂K*ur∂K)d∂K # skeleton projection terms - ∫(vK*(urK-uK))dΩ - ∫(vK*ur∂K)dΩ # bulk projection terms - ∫(v∂K*urK)d∂K - ∫(v∂K*u∂K)d∂K - ∫(v∂K*ur∂K)d∂K - - - - - - - - - xK_x∂K = R(xh) - - @enter xK_x∂K_ΠK,xK_x∂K_Π∂K=P(xh) - - -trial = GridapHybrid._to_trial_basis(xh) -test = GridapHybrid._get_test_fe_basis(P) -@enter RHS_contribs = P.RHS_form(trial, test) - - @enter LHSf,RHSf = GridapHybrid._evaluate_forms(P,xh) - cell_dofs = GridapHybrid._compute_cell_dofs(LHSf, RHSf) - O = R.test_space - GridapHybrid._generate_image_space_span(R,O,xh,cell_dofs,basis_style) - + # Definition of r bilinear form whenever u is a TrialBasis + function r(u,v) + uK_u∂K=R(u) + vK_v∂K=R(v) + uK,u∂K = uK_u∂K + vK,v∂K = vK_v∂K + ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + + ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ + end + + function s(u,v) + Ω = get_triangulation(u) + cmaps = collect1d(get_cell_map(Ω)) + h_array = lazy_map(foo,cmaps) + # h_T = CellField(h_array,Ω,PhysicalDomain()) + h_T = CellField(h_array,Ω) + h_T_1 = 1.0/h_T + # h_T=CellField(get_array(∫(1)dΩ),Ω) + # h_T_1=1.0/h_T + # h_T_2=1.0/(h_T*h_T) + + # # Currently, I cannot use this CellField in the integrand of the skeleton integrals below. + # # I think we need to develop a specific version for _transform_face_to_cell_lface_expanded_array + # # I will be using h_T_1 in the meantime + # h_F=CellField(get_array(∫(1)dΓ),Γ) + # h_F=1.0/h_F + + uK_u∂K_ΠK,uK_u∂K_Π∂K=projection_op(u) + vK_v∂K_ΠK,vK_v∂K_Π∂K=projection_op(v) + + uK_ΠK , u∂K_ΠK = uK_u∂K_ΠK + uK_Π∂K , u∂K_Π∂K = uK_u∂K_Π∂K + + vK_ΠK , v∂K_ΠK = vK_v∂K_ΠK + vK_Π∂K , v∂K_Π∂K = vK_v∂K_Π∂K + + vK_Π∂K_vK_ΠK=vK_Π∂K-vK_ΠK + v∂K_Π∂K_v∂K_ΠK=v∂K_Π∂K-v∂K_ΠK + uK_Π∂K_uK_ΠK=uK_Π∂K-uK_ΠK + u∂K_Π∂K_u∂K_ΠK=u∂K_Π∂K-u∂K_ΠK + + ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(uK_Π∂K_uK_ΠK))d∂K + + ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + + ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + + ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(uK_Π∂K_uK_ΠK))d∂K + end + aa(u,v) = r(u,v) + s(u,v) + l((vK,)) = ∫(vK*f)dΩ + weakform = (u,v)->(aa(u,v),l(v)) + op=HybridAffineFEOperator((u,v)->(aa(u,v),l(v,)), U, V, [1], [2]) + lh = solve(op.skeleton_op) + + A = op.skeleton_op.op.matrix + xh = get_trial_fe_basis(U) + yh = get_fe_basis(V) +xK_x∂K = R(xh) +xK, x∂K = xK_x∂K +quad_cell_data = get_data(d∂K.quad) +cp_∂K = get_cell_points(d∂K.quad) +x∂K_at_cp_∂K = evaluate(x∂K, cp_∂K) +xK_at_cp_∂K = evaluate(xK, cp_∂K) +cp_Ω = get_cell_points(dΩ.quad) +x∂K_at_cp_Ω = evaluate(x∂K, cp_Ω) +xK_at_cp_Ω = evaluate(xK, cp_Ω) +xK_x∂K_ΠK_xK_x∂K_Π∂K = projection_op(xh) +xK_x∂K_ΠK, xK_x∂K_Π∂K = xK_x∂K_ΠK_xK_x∂K_Π∂K +xK_ΠK , x∂K_ΠK = xK_x∂K_ΠK +xK_ΠK_at_cp_Ω = evaluate(xK_ΠK, cp_Ω) +xK_ΠK_at_cp_∂K = evaluate(xK_ΠK, cp_∂K) - function r(u,v) - uK_u∂K=R(u) - vK_v∂K=R(v) - uK,u∂K = uK_u∂K - vK,v∂K = vK_v∂K - ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + - ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ - end - - function s(u,v) - Ω = get_triangulation(u) - cmaps = collect1d(get_cell_map(Ω)) - h_array = lazy_map(foo,cmaps) - # h_T = CellField(h_array,Ω,PhysicalDomain()) - h_T = CellField(h_array,Ω) - h_T_1 = 1.0/h_T - - - uK_u∂K_ΠK,uK_u∂K_Π∂K=P(u) - vK_v∂K_ΠK,vK_v∂K_Π∂K=P(v) - - uK_ΠK , u∂K_ΠK = uK_u∂K_ΠK - uK_Π∂K , u∂K_Π∂K = uK_u∂K_Π∂K - - vK_ΠK , v∂K_ΠK = vK_v∂K_ΠK - vK_Π∂K , v∂K_Π∂K = vK_v∂K_Π∂K - - vK_Π∂K_vK_ΠK=vK_Π∂K-vK_ΠK - v∂K_Π∂K_v∂K_ΠK=v∂K_Π∂K-v∂K_ΠK - uK_Π∂K_uK_ΠK=uK_Π∂K-uK_ΠK - u∂K_Π∂K_u∂K_ΠK=u∂K_Π∂K-u∂K_ΠK - - ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(uK_Π∂K_uK_ΠK))d∂K + - ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + - ∫(h_T_1*(vK_Π∂K_vK_ΠK)*(u∂K_Π∂K_u∂K_ΠK))d∂K + - ∫(h_T_1*(v∂K_Π∂K_v∂K_ΠK)*(uK_Π∂K_uK_ΠK))d∂K - end - - u(x) = x[1]^p+x[2]^p - f(x) = -Δ(u)(x) - - a(u,v)=r(u,v)+s(u,v) - l((vK,))=∫(vK*f)dΩ +x∂K_ΠK_at_cp_Ω = evaluate(x∂K_ΠK, cp_Ω) +x∂K_ΠK_at_cp_∂K = evaluate(x∂K_ΠK, cp_∂K) - uh = get_trial_fe_basis(U) - @enter Ruh = R(uh) +xK_Π∂K, x∂K_Π∂K = xK_x∂K_Π∂K - @enter op=HybridAffineFEOperator((u,v)->(a(u,v),l(v,)), U, V, [1], [2]) - -# return A, h, ndofsΩ -# end - -function domain_operation(operation, l::DomainContribution, r::DomainContribution) - @check get_domains(r) == get_domains(l) - o = DomainContribution() - for trian in get_domains(l) - cl = get_contribution(l, trian) - cr = get_contribution(r, trian) - c = lazy_map(Broadcasting(operation), cl, cr) - add_contribution!(o, trian, c) - end - return o -end \ No newline at end of file +xK_Π∂K_at_cp_∂K = evaluate(xK_Π∂K, cp_∂K) +x∂K_Π∂K_at_cp_∂K = evaluate(x∂K_Π∂K, cp_∂K) \ No newline at end of file diff --git a/test/_dev/Distributed/DarcyHDGTests.jl b/test/_dev/_deprecated/Distributed/DarcyHDGTests.jl similarity index 100% rename from test/_dev/Distributed/DarcyHDGTests.jl rename to test/_dev/_deprecated/Distributed/DarcyHDGTests.jl diff --git a/test/_dev/_deprecated/Distributed/HybridAffineFEOperators.jl b/test/_dev/_deprecated/Distributed/HybridAffineFEOperators.jl deleted file mode 100644 index 4a7282f..0000000 --- a/test/_dev/_deprecated/Distributed/HybridAffineFEOperators.jl +++ /dev/null @@ -1,85 +0,0 @@ -function _merge_bulk_and_skeleton_contributions( - matcontribs::GridapDistributed.DistributedDomainContribution, - veccontribs::GridapDistributed.DistributedDomainContribution) - dmat,dvec=map_parts(matcontribs.contribs,veccontribs.contribs) do matseq, vecseq - _merge_bulk_and_skeleton_contributions(matseq,vecseq) - end - GridapDistributed.DistributedDomainContribution(dmat), - GridapDistributed.DistributedDomainContribution(dvec) -end - -function Gridap.FESpaces._pair_contribution_when_possible( - obiform::GridapDistributed.DistributedDomainContribution, - oliform::GridapDistributed.DistributedDomainContribution) - dmv,dm,dv=map_parts(obiform.contribs,oliform.contribs) do obiform, oliform - Gridap.FESpaces._pair_contribution_when_possible(obiform,oliform) - end - GridapDistributed.DistributedDomainContribution(dmv), - GridapDistributed.DistributedDomainContribution(dm), - GridapDistributed.DistributedDomainContribution(dv) -end - -function Gridap.FESpaces._collect_cell_matrix_and_vector( - trial :: GridapDistributed.DistributedFESpace, - test :: GridapDistributed.DistributedFESpace, - matvec :: GridapDistributed.DistributedDomainContribution, - mat :: GridapDistributed.DistributedDomainContribution, - vec :: GridapDistributed.DistributedDomainContribution) - map_parts(GridapDistributed.local_views(trial), - GridapDistributed.local_views(test), - matvec.contribs,mat.contribs,vec.contribs) do trial, test, matvec, mat, vec - Gridap.FESpaces._collect_cell_matrix_and_vector(trial,test,matvec,mat,vec) - end -end - -function _add_static_condensation( - matvec::GridapDistributed.DistributedDomainContribution, - bulk_fields, skeleton_fields) - dmv=map_parts(matvec.contribs) do matvec - _add_static_condensation(matvec,bulk_fields,skeleton_fields) - end - GridapDistributed.DistributedDomainContribution(dmv) -end - -function Gridap.FESpaces._attach_dirichlet( - matvec:: GridapDistributed.DistributedDomainContribution, - mat :: GridapDistributed.DistributedDomainContribution, - uhd) - dmv,dm=map_parts(matvec.contribs,mat.contribs,uhd.fields) do matvec, mat, uhd - m,v=Gridap.FESpaces._attach_dirichlet(matvec,mat, uhd) - m,v - end - GridapDistributed.DistributedDomainContribution(dmv), - GridapDistributed.DistributedDomainContribution(dm) -end - - -function _compute_hybridizable_from_skeleton_free_dof_values( - skeleton_fe_function :: GridapDistributed.DistributedCellField, - trial_hybridizable :: GridapDistributed.DistributedFESpace, - test_hybridizable :: GridapDistributed.DistributedFESpace, - trial_skeleton :: GridapDistributed.DistributedFESpace, - matvec :: GridapDistributed.DistributedDomainContribution, - bulk_fields, skeleton_fields) - - function f(skel_fe, trial_hyb, test_hyb, trial_skel, matvec) - _compute_hybridizable_from_skeleton_free_dof_values(skel_fe, - trial_hyb, - test_hyb, - trial_skel, - matvec, - bulk_fields, - skeleton_fields) - end - # f(skeleton_fe_function.fields.parts[1], - # trial_hybridizable.part_fe_space.parts[1], - # test_hybridizable.part_fe_space.parts[1], - # trial_skeleton.spaces.parts[1], - # matvec.contribs.parts[1]) - values=map_parts(f,GridapDistributed.local_views(skeleton_fe_function), - GridapDistributed.local_views(trial_hybridizable), - GridapDistributed.local_views(test_hybridizable), - GridapDistributed.local_views(trial_skeleton), - GridapDistributed.local_views(matvec)) - PVector(values,trial_hybridizable.gids) -end diff --git a/test/_dev/_deprecated/Distributed/Skeleton.jl b/test/_dev/_deprecated/Distributed/Skeleton.jl deleted file mode 100644 index 76c92a4..0000000 --- a/test/_dev/_deprecated/Distributed/Skeleton.jl +++ /dev/null @@ -1,75 +0,0 @@ - -function SkeletonTriangulation(model::GridapDistributed.DistributedDiscreteModel) - SkeletonTriangulation(no_ghost,model::GridapDistributed.DistributedDiscreteModel) -end - - -# function _sign_flips(model::GridapDistributed.DistributedDiscreteModel) -# cell_reffes = map_parts(model.models) do m -# basis,reffe_args,reffe_kwargs = ReferenceFE(raviart_thomas,Float64,0) -# reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...) -# end -# GridapDistributed._generate_sign_flips(model,cell_reffes) -# end - -function _sign_flips(model::GridapDistributed.DistributedDiscreteModel) - cell_reffes = map(local_views(model)) do m - basis,reffe_args,reffe_kwargs = ReferenceFE(raviart_thomas,Float64,0) - ReferenceFE(m,basis,reffe_args...;reffe_kwargs...) - end - GridapDistributed._generate_sign_flips(model,cell_reffes) -end - -# function SkeletonTriangulation(portion::GridapDistributed.WithGhost, -# model::GridapDistributed.DistributedDiscreteModel) -# dtrians=map_parts(model.models,_sign_flips(model)) do model, sign_flip -# Skeleton(model,sign_flip) -# end -# GridapDistributed.DistributedTriangulation(dtrians,model) -# end - -function SkeletonTriangulation(portion::GridapDistributed.WithGhost, - model::GridapDistributed.DistributedDiscreteModel) - dtrians=map(local_views(model),_sign_flips(model)) do model, sign_flip - Skeleton(model,sign_flip) - end - GridapDistributed.DistributedTriangulation(dtrians,model) -end - -# function SkeletonTriangulation(portion::GridapDistributed.NoGhost, -# model::GridapDistributed.DistributedDiscreteModel) -# cell_gids=get_cell_gids(model) -# dtrians=map_parts(model.models, -# _sign_flips(model), -# cell_gids.partition) do model, sign_flip, partition -# cell_to_parent_cell= -# findall([partition.lid_to_part[cell]==partition.part -# for cell=1:length(partition.lid_to_part)]) -# Skeleton(model,cell_to_parent_cell,sign_flip) -# end -# GridapDistributed.DistributedTriangulation(dtrians,model) -# end - -function SkeletonTriangulation(portion::GridapDistributed.NoGhost, - model::GridapDistributed.DistributedDiscreteModel) - cell_gids=get_cell_gids(model) - dtrians=map(local_views(model),_sign_flips(model),partition(cell_gids)) do model, sign_flip, partition - cell_to_parent_cell=own_to_local(partition.indices) - GridapHybrid.Skeleton(model,cell_to_parent_cell,sign_flip) - end - GridapDistributed.DistributedTriangulation(dtrians,model) -end - -function get_cell_normal_vector(dskeleton::GridapDistributed.DistributedTriangulation) - fields=map_parts(dskeleton.trians) do skeleton - get_cell_normal_vector(skeleton) - end - GridapDistributed.DistributedCellField(fields) -end - -function get_cell_owner_normal_vector(dskeleton::GridapDistributed.DistributedTriangulation) - fields=map_parts(dskeleton.trians) do skeleton - get_cell_owner_normal_vector(skeleton) - end - GridapDistributed.DistributedCellField(fields) -end diff --git a/src/HybridDistributed.jl b/test/_dev/_deprecated/HybridDistributed.jl similarity index 100% rename from src/HybridDistributed.jl rename to test/_dev/_deprecated/HybridDistributed.jl