From 5a5acc393dfa88de7df7a90f8b60cf18715db262 Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Thu, 14 Nov 2024 09:02:01 +1100 Subject: [PATCH] Solved error in the BB space projection of the divergence of a flux FE function --- BulkGhostPenaltyAssembleMaps.jl | 41 ++++++++++++++++--------- bulk_ghost_penalty_canvas_towardsdiv.jl | 27 ++++++++-------- 2 files changed, 38 insertions(+), 30 deletions(-) diff --git a/BulkGhostPenaltyAssembleMaps.jl b/BulkGhostPenaltyAssembleMaps.jl index e18749a..9c85630 100644 --- a/BulkGhostPenaltyAssembleMaps.jl +++ b/BulkGhostPenaltyAssembleMaps.jl @@ -26,7 +26,6 @@ function Gridap.Fields.evaluate!(cache,m::BulkGhostPenaltyAssembleLhsMap,cells) end result.array end -#TODO: add {RANK,A,B} struct BulkGhostPenaltyAssembleRhsMap{A,B} <: Gridap.Fields.Map agg_cells_local_dof_ids::A agg_cells_rhs_contribs::B @@ -35,8 +34,7 @@ end function Gridap.Fields.return_cache(m::BulkGhostPenaltyAssembleRhsMap,aggregate_local_cells) cache_agg_cells_local_dof_ids=array_cache(m.agg_cells_local_dof_ids) cache_unassembled_rhs=array_cache(m.agg_cells_rhs_contribs) - rank=length(size(m.agg_cells_rhs_contribs[1])) - evaluate_result=Gridap.Arrays.CachedArray(eltype(eltype(m.agg_cells_rhs_contribs)),rank) + evaluate_result=Gridap.Arrays.CachedArray(eltype(eltype(m.agg_cells_rhs_contribs)),2) cache_agg_cells_local_dof_ids,cache_unassembled_rhs,evaluate_result end @@ -52,24 +50,37 @@ function Gridap.Fields.evaluate!(cache,m::BulkGhostPenaltyAssembleRhsMap,aggrega end end - rank=length(size(result)) - if rank==2 - Gridap.Arrays.setsize!(result,(size(contrib,1),max_local_dof_id)) - else - @assert rank==1 - Gridap.Arrays.setsize!(result,(max_local_dof_id,)) - end + Gridap.Arrays.setsize!(result,(size(contrib,1),max_local_dof_id)) + result.array .= 0.0 for (i,cell) in enumerate(aggregate_local_cells) current_cell_local_dof_ids = getindex!(cache_agg_cells_local_dof_ids,m.agg_cells_local_dof_ids,cell) contrib = getindex!(cache_unassembled_rhs,m.agg_cells_rhs_contribs,cell) for (j,local_dof) in enumerate(current_cell_local_dof_ids) - if rank==2 - result.array[:,local_dof] += contrib[:,j] - else - result.array[local_dof] += contrib[j] - end + result.array[:,local_dof] += contrib[:,j] end end result.array end + +struct BulkGhostPenaltyAssembleRhsFEFunctionMap{A} <: Gridap.Fields.Map + agg_cells_rhs_contribs::A +end + +function Gridap.Fields.return_cache(m::BulkGhostPenaltyAssembleRhsFEFunctionMap,aggregate_local_cells) + cache_unassembled_rhs=array_cache(m.agg_cells_rhs_contribs) + evaluate_result=Gridap.Arrays.CachedArray(eltype(eltype(m.agg_cells_rhs_contribs)),1) + cache_unassembled_rhs,evaluate_result +end + +function Gridap.Fields.evaluate!(cache,m::BulkGhostPenaltyAssembleRhsFEFunctionMap,aggregate_local_cells) + cache_unassembled_rhs,result=cache + contrib = getindex!(cache_unassembled_rhs,m.agg_cells_rhs_contribs,1) + Gridap.Arrays.setsize!(result,(size(contrib,1),)) + result.array .= 0.0 + for (i,cell) in enumerate(aggregate_local_cells) + contrib = getindex!(cache_unassembled_rhs,m.agg_cells_rhs_contribs,cell) + result.array .+= contrib + end + result.array +end \ No newline at end of file diff --git a/bulk_ghost_penalty_canvas_towardsdiv.jl b/bulk_ghost_penalty_canvas_towardsdiv.jl index d6cc400..05a395c 100644 --- a/bulk_ghost_penalty_canvas_towardsdiv.jl +++ b/bulk_ghost_penalty_canvas_towardsdiv.jl @@ -489,9 +489,9 @@ function divergence_in_pressure_space_fe_function(uh::Gridap.FESpaces.SingleFiel end testxh=FEFunction(X, rand(num_free_dofs(X))) testuh,testph = testxh -testFEfunction = divergence_in_pressure_space_fe_function(testuh,P, dq) -eh = testFEfunction-(∇⋅testuh) -norm_testFEfunction = sum(∫(eh⋅eh)dΩhact) +testdivuh_pressure = divergence_in_pressure_space_fe_function(testuh,P, dq) +eh = testdivuh_pressure-(∇⋅testuh) +@assert sum(∫(eh⋅eh)dΩhact) < 1.0e-12 # TODO: the following functions should be tested and can then replace the lines below these commented lines. # function divergence_in_pressure_space_basis(uh::FEBasis,P,dq) @@ -530,28 +530,25 @@ dqbb_Ωagg_cells =change_domain_bb_to_agg_cells(qbb,ref_agg_cell_to_ref_bb_map, # agg_cells_rhs_contribs=get_array(∫(dqbb_Ωagg_cells⋅div_du_in_pressure_space_single_field)dΩagg_cells) # ass_rhs_map=BulkGhostPenaltyAssembleRhsMap(P_agg_cells_local_dof_ids,agg_cells_rhs_contribs) ## TODO: prep to assemble vectors # rhs=lazy_map(ass_rhs_map,aggregate_to_local_cells) -test_agg_cells_rhs_contribs=get_array(∫(dqbb_Ωagg_cells⋅testFEfunction)dΩagg_cells) -test_ass_rhs_map=BulkGhostPenaltyAssembleRhsMap(P_agg_cells_local_dof_ids,test_agg_cells_rhs_contribs) +test_agg_cells_rhs_contribs=get_array(∫(dqbb_Ωagg_cells⋅testdivuh_pressure)dΩagg_cells) +test_ass_rhs_map=BulkGhostPenaltyAssembleRhsFEFunctionMap(test_agg_cells_rhs_contribs) test_rhs=lazy_map(test_ass_rhs_map,aggregate_to_local_cells) test_div_dv_l2_proj_bb_dofs=lazy_map(\,p_lhs,test_rhs) -# TODO: can be removed if we fix the BulkGhostPenaltyAssembleRhsMap. -# transposed_test_div_dv_l2_proj_bb_dofs= lazy_map(x->collect(reshape(x,(size(x,1)))),lazy_map(transpose,test_div_dv_l2_proj_bb_dofs)) -# test_div_dv_l2_proj_bb_array=lazy_map(Gridap.Fields.linear_combination, -# transposed_test_div_dv_l2_proj_bb_dofs, -# Gridap.CellData.get_data(qbb)) - test_div_dv_l2_proj_bb_array=lazy_map(Gridap.Fields.linear_combination, test_div_dv_l2_proj_bb_dofs, Gridap.CellData.get_data(qbb)) test_div_dv_l2_proj_bb_array_agg_cells=lazy_map(Broadcasting(∘), lazy_map(Reindex(test_div_dv_l2_proj_bb_array),agg_cells_to_aggregate), ref_agg_cell_to_ref_bb_map) -div_uh_proj_bb=Gridap.CellData.GenericCellField(test_div_dv_l2_proj_bb_array_agg_cells,dΩagg_cells.quad.trian,ReferenceDomain()) +div_uh_proj_bb=Gridap.CellData.GenericCellField(test_div_dv_l2_proj_bb_array_agg_cells, + dΩagg_cells.quad.trian, + ReferenceDomain()) x=get_cell_points(dΩagg_cells.quad.trian) -div_uh_proj_bb(x)[1] +div_uh_proj_bb(x)[8] +testdivuh_pressure(x)[7] -res = testFEfunction - div_uh_proj_bb -norm = sum(∫((res)*(res))*dΩagg_cells) \ No newline at end of file +res = testdivuh_pressure - div_uh_proj_bb +sum(∫(res*res)*dΩagg_cells) \ No newline at end of file