diff --git a/src/FESpaces.jl b/src/FESpaces.jl index 053d376..8b451c6 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -49,19 +49,28 @@ function _generate_hanging_faces_owner_face_dofs(num_hanging_faces, ptrs[1] = 1 for fid_hanging = 1:num_hanging_faces glue = hanging_faces_glue[fid_hanging] - ocell_lface = glue[2] - ptrs[fid_hanging+1] = ptrs[fid_hanging] + length(face_dofs[ocell_lface]) + if (glue[1]!=-1) + ocell_lface = glue[2] + ptrs[fid_hanging+1] = ptrs[fid_hanging] + length(face_dofs[ocell_lface]) + else + ptrs[fid_hanging+1] = ptrs[fid_hanging] + 1 + end end data_owner_face_dofs = Vector{Int}(undef, ptrs[num_hanging_faces+1] - 1) for fid_hanging = 1:num_hanging_faces glue = hanging_faces_glue[fid_hanging] ocell = glue[1] - ocell_lface = glue[2] - s = ptrs[fid_hanging] - e = ptrs[fid_hanging+1] - 1 - current_cell_dof_ids = getindex!(cache, cell_dof_ids, ocell) - for (j, ldof) in enumerate(face_dofs[ocell_lface]) + if (ocell!=-1) + ocell_lface = glue[2] + s = ptrs[fid_hanging] + e = ptrs[fid_hanging+1] - 1 + current_cell_dof_ids = getindex!(cache, cell_dof_ids, ocell) + for (j, ldof) in enumerate(face_dofs[ocell_lface]) data_owner_face_dofs[s+j-1] = current_cell_dof_ids[ldof] + end + else + s = ptrs[fid_hanging] + data_owner_face_dofs[s] = 1 end end Gridap.Arrays.Table(data_owner_face_dofs, ptrs) @@ -141,61 +150,72 @@ function _generate_constraints!(Df, current_cell_dof_ids = getindex!(cache_dof_ids, cell_dof_ids, cell) lface = hanging_faces_to_lface[fid_hanging] ocell, ocell_lface, subface = hanging_faces_glue[fid_hanging] - ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) - oface_dim = face_dim(Val{Dc}, ocell_lface) - - if (Df == 0) # Am I a vertex? - hanging_lvertex_within_first_subface = 2^oface_dim - cur_subface_own_dofs = subface_own_dofs[oface_dim][hanging_lvertex_within_first_subface] - elseif (Df == 1 && Dc == 3) # Am I an edge? - if (subface < 0) # Edge hanging in the interior of a face - @assert subface == -1 || subface == -2 || subface == -3 || subface == -4 - @assert oface_dim == Dc - 1 - abs_subface = abs(subface) - if (abs_subface == 1) - subface = 1 - edge = 4 + 4 # num_vertices+edge_id - elseif (abs_subface == 2) - subface = 3 - edge = 4 + 4 # num_vertices+edge_id - elseif (abs_subface == 3) - subface = 3 - edge = 4 + 1 # num_vertices+edge_id - elseif (abs_subface == 4) - subface = 4 - edge = 4 + 1 # num_vertices+edge_id - end - cur_subface_own_dofs = subface_own_dofs[oface_dim][edge] - else - @assert subface == 1 || subface == 2 - @assert oface_dim == 1 + if (ocell!=-1) + ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) + oface_dim = face_dim(Val{Dc}, ocell_lface) + + if (Df == 0) # Am I a vertex? hanging_lvertex_within_first_subface = 2^oface_dim + cur_subface_own_dofs = subface_own_dofs[oface_dim][hanging_lvertex_within_first_subface] + elseif (Df == 1 && Dc == 3) # Am I an edge? + if (subface < 0) # Edge hanging in the interior of a face + @assert subface == -1 || subface == -2 || subface == -3 || subface == -4 + @assert oface_dim == Dc - 1 + abs_subface = abs(subface) + if (abs_subface == 1) + subface = 1 + edge = 4 + 4 # num_vertices+edge_id + elseif (abs_subface == 2) + subface = 3 + edge = 4 + 4 # num_vertices+edge_id + elseif (abs_subface == 3) + subface = 3 + edge = 4 + 1 # num_vertices+edge_id + elseif (abs_subface == 4) + subface = 4 + edge = 4 + 1 # num_vertices+edge_id + end + cur_subface_own_dofs = subface_own_dofs[oface_dim][edge] + else + @assert subface == 1 || subface == 2 + @assert oface_dim == 1 + hanging_lvertex_within_first_subface = 2^oface_dim + cur_subface_own_dofs = subface_own_dofs[oface_dim][end] + end + elseif (Df == Dc - 1) # Am I a face? + @assert oface_dim == Dc - 1 cur_subface_own_dofs = subface_own_dofs[oface_dim][end] end - elseif (Df == Dc - 1) # Am I a face? - @assert oface_dim == Dc - 1 - cur_subface_own_dofs = subface_own_dofs[oface_dim][end] - end - oface = getindex!(cache_cell_faces[oface_dim+1], cell_faces[oface_dim+1], ocell)[ocell_lface_within_dim] - oface_lid, _ = owner_faces_lids[oface_dim][oface] - pindex = owner_faces_pindex[oface_dim][oface_lid] - for ((ldof, dof), ldof_subface) in zip(enumerate(face_own_dofs[offset+lface]), cur_subface_own_dofs) - push!(sDOF_to_dof, current_cell_dof_ids[dof]) - push!(sDOF_to_dofs, hanging_faces_owner_face_dofs[fid_hanging]) - coeffs = Vector{Float64}(undef, length(hanging_faces_owner_face_dofs[fid_hanging])) - # Go over dofs of ocell_lface - for (ifdof, icdof) in enumerate(face_dofs[ocell_lface]) - pifdof = node_permutations[oface_dim][pindex][ifdof] - println("XXXX: $(ifdof) $(pifdof)") - ldof_coarse = face_dofs[ocell_lface][pifdof] - coeffs[ifdof] = - ref_constraints[face_subface_ldof_to_cell_ldof[oface_dim][ocell_lface_within_dim][subface][ldof_subface], ldof_coarse] + oface = getindex!(cache_cell_faces[oface_dim+1], cell_faces[oface_dim+1], ocell)[ocell_lface_within_dim] + oface_lid, _ = owner_faces_lids[oface_dim][oface] + pindex = owner_faces_pindex[oface_dim][oface_lid] + for ((ldof, dof), ldof_subface) in zip(enumerate(face_own_dofs[offset+lface]), cur_subface_own_dofs) + push!(sDOF_to_dof, current_cell_dof_ids[dof]) + push!(sDOF_to_dofs, hanging_faces_owner_face_dofs[fid_hanging]) + coeffs = Vector{Float64}(undef, length(hanging_faces_owner_face_dofs[fid_hanging])) + # Go over dofs of ocell_lface + for (ifdof, icdof) in enumerate(face_dofs[ocell_lface]) + pifdof = node_permutations[oface_dim][pindex][ifdof] + println("XXXX: $(ifdof) $(pifdof)") + ldof_coarse = face_dofs[ocell_lface][pifdof] + coeffs[ifdof] = + ref_constraints[face_subface_ldof_to_cell_ldof[oface_dim][ocell_lface_within_dim][subface][ldof_subface], ldof_coarse] + end + push!(sDOF_to_coeffs, coeffs) end - push!(sDOF_to_coeffs, coeffs) - end + else + for (ldof, dof) in enumerate(face_own_dofs[offset+lface]) + push!(sDOF_to_dof, current_cell_dof_ids[dof]) + push!(sDOF_to_dofs, hanging_faces_owner_face_dofs[fid_hanging]) + push!(sDOF_to_coeffs, [1.0]) + end + end end + println("sDOF_to_dof [$(Df)]= $(sDOF_to_dof)") + println("sDOF_to_dofs [$(Df)]= $(sDOF_to_dofs)") + println("sDOF_to_coeffs [$(Df)]= $(sDOF_to_coeffs)") end # count how many different owner faces @@ -217,13 +237,15 @@ function _compute_owner_faces_pindex_and_lids(Df, owner_faces_lids = Dict{Int,Tuple{Int,Int,Int}}() for fid_hanging = 1:num_hanging_faces ocell, ocell_lface, _ = hanging_faces_glue[fid_hanging] - ocell_dim = face_dim(Val{Dc}, ocell_lface) - if (ocell_dim == Df) - ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) - owner_face = cell_faces[ocell][ocell_lface_within_dim] - if !(haskey(owner_faces_lids, owner_face)) - num_owner_faces += 1 - owner_faces_lids[owner_face] = (num_owner_faces, ocell, ocell_lface) + if (ocell!=-1) + ocell_dim = face_dim(Val{Dc}, ocell_lface) + if (ocell_dim == Df) + ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) + owner_face = cell_faces[ocell][ocell_lface_within_dim] + if !(haskey(owner_faces_lids, owner_face)) + num_owner_faces += 1 + owner_faces_lids[owner_face] = (num_owner_faces, ocell, ocell_lface) + end end end end @@ -236,16 +258,18 @@ function _compute_owner_faces_pindex_and_lids(Df, for fid_hanging = 1:num_hanging_faces ocell, ocell_lface, subface = hanging_faces_glue[fid_hanging] - ocell_dim = face_dim(Val{Dc}, ocell_lface) - if (ocell_dim == Df) - cell = hanging_faces_to_cell[fid_hanging] - lface = hanging_faces_to_lface[fid_hanging] - cvertex = lface_to_cvertices[lface][subface] - vertex = cell_vertices[cell][cvertex] - ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) - owner_face = cell_faces[ocell][ocell_lface_within_dim] - owner_face_lid, _ = owner_faces_lids[owner_face] - owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+subface] = vertex + if (ocell!=-1) + ocell_dim = face_dim(Val{Dc}, ocell_lface) + if (ocell_dim == Df) + cell = hanging_faces_to_cell[fid_hanging] + lface = hanging_faces_to_lface[fid_hanging] + cvertex = lface_to_cvertices[lface][subface] + vertex = cell_vertices[cell][cvertex] + ocell_lface_within_dim = face_lid_within_dim(Val{Dc}, ocell_lface) + owner_face = cell_faces[ocell][ocell_lface_within_dim] + owner_face_lid, _ = owner_faces_lids[owner_face] + owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+subface] = vertex + end end end @@ -314,11 +338,11 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, hanging_faces_glue, dmodel.models, spaces_wo_constraints) do gridap_cell_faces, - num_regular_faces, - num_hanging_faces, - hanging_faces_glue, - model, - V + num_regular_faces, + num_hanging_faces, + hanging_faces_glue, + model, + V hanging_faces_to_cell = Vector{Vector{Int}}(undef, Dc) hanging_faces_to_lface = Vector{Vector{Int}}(undef, Dc) @@ -502,7 +526,7 @@ end # An auxiliary function which we use in order to generate a version of # get_cell_dof_ids() for FE spaces with linear constraints which is suitable # for the algorithm which generates the global DoFs identifiers -function fe_space_with_linear_constraints_cell_dof_ids(Uc::FESpaceWithLinearConstraints) +function fe_space_with_linear_constraints_cell_dof_ids(Uc::FESpaceWithLinearConstraints,sDOF_to_dof) U_cell_dof_ids = Gridap.Arrays.Table(get_cell_dof_ids(Uc.space)) ndata = U_cell_dof_ids.ptrs[end] - 1 Uc_cell_dof_ids_data = zeros(eltype(U_cell_dof_ids.data), ndata) @@ -515,6 +539,7 @@ function fe_space_with_linear_constraints_cell_dof_ids(Uc::FESpaceWithLinearCons n_cells = length(U_cell_dof_ids) n_fdofs = num_free_dofs(Uc.space) n_fmdofs = Uc.n_fmdofs + dof_to_sDOF = Dict(val=>key for (key,val) in enumerate(sDOF_to_dof)) for cell in 1:n_cells pini = U_cell_dof_ids.ptrs[cell] pend = U_cell_dof_ids.ptrs[cell+1] - 1 @@ -522,13 +547,13 @@ function fe_space_with_linear_constraints_cell_dof_ids(Uc::FESpaceWithLinearCons dof = U_cell_dof_ids.data[p] DOF = Gridap.FESpaces._dof_to_DOF(dof, n_fdofs) qini = Uc.DOF_to_mDOFs.ptrs[DOF] - qend = Uc.DOF_to_mDOFs.ptrs[DOF+1] - 1 - if (qend - qini == 0) # master DOF + qend = Uc.DOF_to_mDOFs.ptrs[DOF+1]-1 + if (!haskey(dof_to_sDOF,dof)) # master DoF + @assert qend-qini==0 mDOF = Uc.DOF_to_mDOFs.data[qini] mdof = Gridap.FESpaces._DOF_to_dof(mDOF, n_fmdofs) Uc_cell_dof_ids_data[p] = mdof else # slave DoF - @assert qend - qini > 0 Uc_cell_dof_ids_data[p] = max_negative_minus_one end end @@ -563,14 +588,18 @@ function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc}, reff Vc = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V) end - local_cell_dof_ids = map(spaces_w_constraints) do Vc - result = fe_space_with_linear_constraints_cell_dof_ids(Vc) + local_cell_dof_ids = map(spaces_w_constraints,sDOF_to_dof) do Vc,sDOF_to_dof + result = fe_space_with_linear_constraints_cell_dof_ids(Vc,sDOF_to_dof) println("result=", result) result end nldofs = map(num_free_dofs,spaces_w_constraints) cell_gids = get_cell_gids(model) gids=GridapDistributed.generate_gids(cell_gids,local_cell_dof_ids,nldofs) + map(partition(gids)) do indices + println("[$(part_id(indices))]: l2g=$(local_to_global(indices))") + println("[$(part_id(indices))]: l2o=$(local_to_owner(indices))") + end vector_type = GridapDistributed._find_vector_type(spaces_w_constraints,gids) GridapDistributed.DistributedSingleFieldFESpace(spaces_w_constraints,gids,vector_type) end diff --git a/src/OctreeDistributedDiscreteModels.jl b/src/OctreeDistributedDiscreteModels.jl index 64ab486..4c0d759 100644 --- a/src/OctreeDistributedDiscreteModels.jl +++ b/src/OctreeDistributedDiscreteModels.jl @@ -246,12 +246,33 @@ end function pXest_partition!(::Type{Val{Dc}}, ptr_pXest) where Dc if (Dc==2) + # The 1 here is required to avoid that the children of the + # same parent are assigned to different partitions p4est_partition(ptr_pXest, 0, C_NULL) else @assert false end end +function pXest_balance!(::Type{Val{Dc}}, ptr_pXest; k_2_1_balance=0) where Dc + if (Dc==2) + if (k_2_1_balance==0) + p4est_balance(ptr_pXest, P4est_wrapper.P4EST_CONNECT_FULL, C_NULL) + else + p4est_balance(ptr_pXest, P4est_wrapper.P4EST_CONNECT_FACE, C_NULL) + end + else + if (k_2_1_balance==0) + p8est_balance(ptr_pXest, P4est_wrapper.P8EST_CONNECT_FULL, C_NULL) + elseif (k_2_1_balance==1) + p8est_balance(ptr_pXest, P4est_wrapper.P8EST_CONNECT_EDGE, C_NULL) + else + @assert k_2_1_balance==2 + p8est_balance(ptr_pXest, P4est_wrapper.P8EST_CONNECT_FACE, C_NULL) + end + end +end + function pXest_partition_given!(::Type{Val{Dc}}, ptr_pXest, new_num_cells_per_part) where Dc if (Dc==2) p4est_partition_given(ptr_pXest, new_num_cells_per_part) @@ -340,6 +361,163 @@ const rrule_f_to_c_dim_2D=Vector{Vector{Vector{UInt8}}}( [[1,2,1,2], [1,2,2,1], [2,1,1,2], [2,1,2,1]] # e ]) + function p4est_update_flags!(ptr_p4est_old, + ptr_p4est_new) + p4est_old = ptr_p4est_old[] + p4est_new = ptr_p4est_new[] + flags=unsafe_wrap(Array, + Ptr{Cint}(p4est_old.user_pointer), + p4est_old.local_num_quadrants) + + num_trees = Cint(p4est_old.connectivity[].num_trees) + @assert num_trees == Cint(p4est_new.connectivity[].num_trees) + + for itree = 0:num_trees-1 + tree_old = _p4est_tree_array_index(Val{2},p4est_old.trees,itree)[] + tree_new = _p4est_tree_array_index(Val{2},p4est_new.trees,itree)[] + num_quads_old = Cint(tree_old.quadrants.elem_count) + iquad_new = 0 + iquad_old = 0 + while iquad_old < num_quads_old + q_old = _p4est_quadrant_array_index(Val{2},tree_old.quadrants,iquad_old) + q_new = _p4est_quadrant_array_index(Val{2},tree_new.quadrants,iquad_new) + if (p4est_quadrant_compare(q_old,q_new) == 0) # q_old was not refined nor coarsened + flags[iquad_old+1] = nothing_flag + iquad_new += 1 + iquad_old += 1 + elseif (p4est_quadrant_is_parent(q_old,q_new)!=0) # q_old was refined + flags[iquad_old+1] = refine_flag + iquad_new += 4 + iquad_old += 1 + elseif (p4est_quadrant_is_parent(q_new,q_old)!=0) # q_old and its siblings were coarsened + for i=0:3 + flags[iquad_old+i+1] = coarsen_flag + end + iquad_old += 4 + iquad_new += 1 + else + @assert false + end + end + end + end + + # void F90_p4est_update_refinement_and_coarsening_flags(p4est_t * p4est_old, p4est_t * p4est_new) + # { + # p4est_tree_t *tree_old; + # p4est_quadrant_t *q_old; + # sc_array_t *quadrants_old; + # int old_quadrant_index; + + # p4est_tree_t *tree_new; + # p4est_quadrant_t *q_new; + # sc_array_t *quadrants_new; + # int i, new_quadrant_index; + + # int * user_pointer; + + # P4EST_ASSERT(p4est_old->user_pointer == p4est_new->user_pointer); + + # user_pointer = (int *) p4est_old->user_pointer; + + # // Extract references to the first (and uniquely allowed) trees + # tree_old = p4est_tree_array_index (p4est_old->trees,0); + # tree_new = p4est_tree_array_index (p4est_new->trees,0); + + # quadrants_old = &(tree_old->quadrants); + # quadrants_new = &(tree_new->quadrants); + + # new_quadrant_index = 0; + # for (old_quadrant_index=0; old_quadrant_index < quadrants_old->elem_count;) + # { + # q_old = p4est_quadrant_array_index(quadrants_old, old_quadrant_index); + # q_new = p4est_quadrant_array_index(quadrants_new, new_quadrant_index); + # if ( p4est_quadrant_compare(q_old,q_new) == 0 ) //q_old was not refined nor coarsened + # { + # user_pointer[old_quadrant_index] = FEMPAR_do_nothing_flag; + # old_quadrant_index++; + # new_quadrant_index++; + # } + # else if ( p4est_quadrant_is_parent(q_old,q_new) ) //q_old was refined + # { + # user_pointer[old_quadrant_index] = FEMPAR_refinement_flag; + # old_quadrant_index++; + # new_quadrant_index = new_quadrant_index + P4EST_CHILDREN; + # } + # else if ( p4est_quadrant_is_parent(q_new,q_old) ) //q_old and its siblings were coarsened + # { + # for (i=0; i < P4EST_CHILDREN; i++) + # { + # user_pointer[old_quadrant_index] = FEMPAR_coarsening_flag; + # old_quadrant_index++; + # } + # new_quadrant_index++; + # } + # else + # { + # P4EST_ASSERT(0); + # } + # } + # } + + # void F90_p8est_update_refinement_and_coarsening_flags(p8est_t * p8est_old, p8est_t * p8est_new) + # { + # p8est_tree_t *tree_old; + # p8est_quadrant_t *q_old; + # sc_array_t *quadrants_old; + # int old_quadrant_index; + + # p8est_tree_t *tree_new; + # p8est_quadrant_t *q_new; + # sc_array_t *quadrants_new; + # int i, new_quadrant_index; + + # int * user_pointer; + + # P4EST_ASSERT(p8est_old->user_pointer == p8est_new->user_pointer); + + # user_pointer = (int *) p8est_old->user_pointer; + + # // Extract references to the first (and uniquely allowed) trees + # tree_old = p8est_tree_array_index (p8est_old->trees,0); + # tree_new = p8est_tree_array_index (p8est_new->trees,0); + + # quadrants_old = &(tree_old->quadrants); + # quadrants_new = &(tree_new->quadrants); + + # new_quadrant_index = 0; + # for (old_quadrant_index=0; old_quadrant_index < quadrants_old->elem_count;) + # { + # q_old = p8est_quadrant_array_index(quadrants_old, old_quadrant_index); + # q_new = p8est_quadrant_array_index(quadrants_new, new_quadrant_index); + # if ( p8est_quadrant_compare(q_old,q_new) == 0 ) //q_old was not refined nor coarsened + # { + # user_pointer[old_quadrant_index] = FEMPAR_do_nothing_flag; + # old_quadrant_index++; + # new_quadrant_index++; + # } + # else if ( p8est_quadrant_is_parent(q_old,q_new) ) //q_old was refined + # { + # user_pointer[old_quadrant_index] = FEMPAR_refinement_flag; + # old_quadrant_index++; + # new_quadrant_index = new_quadrant_index + P8EST_CHILDREN; + # } + # else if ( p8est_quadrant_is_parent(q_new,q_old) ) //q_old and its siblings were coarsened + # { + # for (i=0; i < P8EST_CHILDREN; i++) + # { + # user_pointer[old_quadrant_index] = FEMPAR_coarsening_flag; + # old_quadrant_index++; + # } + # new_quadrant_index++; + # } + # else + # { + # P4EST_ASSERT(0); + # } + # } + + # } function _compute_fine_to_coarse_model_glue( cparts, @@ -812,9 +990,11 @@ function Gridap.Adaptivity.refine(model::OctreeDistributedDiscreteModel{Dc,Dp}, refine_callback_2d_c = @cfunction($refine_callback_2d, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t})) - # Copy and refine input p4est + # Copy input p4est, refine and balance ptr_new_pXest = pXest_copy(Val{Dc}, model.ptr_pXest) pXest_refine!(Val{Dc}, ptr_new_pXest, refine_callback_2d_c) + pXest_balance!(Val{Dc}, ptr_new_pXest) + p4est_update_flags!(model.ptr_pXest,ptr_new_pXest) # Extract ghost and lnodes ptr_pXest_ghost = setup_pXest_ghost(Val{Dc}, ptr_new_pXest) @@ -1966,6 +2146,7 @@ function generate_cell_faces_and_non_conforming_glue(::Type{Val{Dc}}, gridap_cell_edges[PXEST_2_GRIDAP_EDGE[p4est_ledge+1]] = hanging_edge_from_face_code end end + end end @@ -2034,6 +2215,7 @@ function generate_cell_faces_and_non_conforming_glue(::Type{Val{Dc}}, end end end + end end end @@ -2041,8 +2223,29 @@ function generate_cell_faces_and_non_conforming_glue(::Type{Val{Dc}}, # Go over all touched hanging faces and start # assigning IDs from the last num_regular_faces ID # For each hanging face, keep track of (owner_cell,lface) - println("regular_faces_p4est_to_gridap: $(regular_faces_p4est_to_gridap)") - println("hanging_faces_pairs_to_owner_face: $(hanging_faces_pairs_to_owner_face)") + println("regular_faces_p4est_to_gridap [$(part_id(partition(cell_prange).item))]: $(regular_faces_p4est_to_gridap)") + println("hanging_faces_pairs_to_owner_face[$(part_id(partition(cell_prange).item))]: $(hanging_faces_pairs_to_owner_face)") + + # Go over all hanging faces + # Detect if the owner face is in a ghost cell. + # If not in a ghost cell or touched + # Else + # The current face becomes a regular face + # end + function is_ghost(cell) + cell>own_length(indices) + end + # for key in keys(hanging_faces_pairs_to_owner_face) + # (cell, lface) = key + # (owner_p4est_gface, half) = hanging_faces_pairs_to_owner_face[key] + # if (is_ghost(cell) && !(haskey(regular_faces_p4est_to_gridap,owner_p4est_gface))) + # # Current face becomes a regular face (locally) + # num_regular_faces[Dc]+=1 + # start_gridap_faces = (cell - 1) * n_cell_faces + # gridap_cell_faces_data[start_gridap_faces+lface] = num_regular_faces[Dc] + # end + # end + hanging_faces_owner_cell_and_lface = Vector{Tuple{Int,Int,Int}}(undef, length(keys(hanging_faces_pairs_to_owner_face))) @@ -2050,38 +2253,62 @@ function generate_cell_faces_and_non_conforming_glue(::Type{Val{Dc}}, for key in keys(hanging_faces_pairs_to_owner_face) (cell, lface) = key (owner_p4est_gface, half) = hanging_faces_pairs_to_owner_face[key] - owner_gridap_gface = regular_faces_p4est_to_gridap[owner_p4est_gface] num_hanging_faces[Dc] += 1 start_gridap_faces = (cell - 1) * n_cell_faces gridap_cell_faces_data[start_gridap_faces+lface] = num_regular_faces[Dc] + num_hanging_faces[Dc] - (owner_cell, p4est_lface) = p4est_gface_to_gcell_p4est_lface[owner_p4est_gface] - hanging_faces_owner_cell_and_lface[num_hanging_faces[Dc]] = - (owner_cell, n_cell_vertices+n_cell_edges+PXEST_2_GRIDAP_FACE[p4est_lface], half) + if (!(is_ghost(cell)) || haskey(regular_faces_p4est_to_gridap,owner_p4est_gface)) + (owner_cell, p4est_lface) = p4est_gface_to_gcell_p4est_lface[owner_p4est_gface] + hanging_faces_owner_cell_and_lface[num_hanging_faces[Dc]] = + (owner_cell, n_cell_vertices+n_cell_edges+PXEST_2_GRIDAP_FACE[p4est_lface], half) + else + # Glue info cannot be computed for this hanging face + hanging_faces_owner_cell_and_lface[num_hanging_faces[Dc]] = (-1,-1,-1) + end end println("AAA: $(gridap_cell_faces_data)") + # owner_p4est_gface_to_regular_vertex = Dict{Int,Int}() + # for key in keys(hanging_vertices_pairs_to_owner_face) + # (cell, lvertex) = key + # owner_p4est_gface = hanging_vertices_pairs_to_owner_face[key] + # if (is_ghost(cell) && !(haskey(regular_faces_p4est_to_gridap,owner_p4est_gface))) + # if !(haskey(owner_p4est_gface_to_regular_vertex, owner_p4est_gface)) + # num_regular_faces[1] += 1 + # owner_p4est_gface_to_regular_vertex[owner_p4est_gface] = num_regular_faces[1] + # end + # start_gridap_vertices = (cell - 1) * n_cell_vertices + # gridap_cell_vertices_data[start_gridap_vertices+lvertex] = + # owner_p4est_gface_to_regular_vertex[owner_p4est_gface] + # end + # end + + # Go over all touched hanging vertices and start # assigning IDs from the last num_regular_vertices ID # For each hanging vertex, keep track of (owner_cell,lface) num_hanging_faces[1] = 0 hanging_vertices_owner_cell_and_lface = Tuple{Int,Int,Int}[] - owner_gridap_gface_to_hanging_vertex = Dict{Int,Int}() half=1 + owner_p4est_gface_to_hanging_vertex = Dict{Int,Int}() for key in keys(hanging_vertices_pairs_to_owner_face) (cell, lvertex) = key owner_p4est_gface = hanging_vertices_pairs_to_owner_face[key] - owner_gridap_gface = regular_faces_p4est_to_gridap[owner_p4est_gface] - if !(haskey(owner_gridap_gface_to_hanging_vertex, owner_gridap_gface)) + if !(haskey(owner_p4est_gface_to_hanging_vertex, owner_p4est_gface)) num_hanging_faces[1] += 1 - owner_gridap_gface_to_hanging_vertex[owner_gridap_gface] = num_hanging_faces[1] - (owner_cell, p4est_lface) = p4est_gface_to_gcell_p4est_lface[owner_p4est_gface] - push!(hanging_vertices_owner_cell_and_lface, - (owner_cell, n_cell_vertices+n_cell_edges+PXEST_2_GRIDAP_FACE[p4est_lface],half)) - end + owner_p4est_gface_to_hanging_vertex[owner_p4est_gface] = num_hanging_faces[1] + if (!is_ghost(cell) || (haskey(regular_faces_p4est_to_gridap,owner_p4est_gface))) + (owner_cell, p4est_lface) = p4est_gface_to_gcell_p4est_lface[owner_p4est_gface] + owner_gridap_gface = regular_faces_p4est_to_gridap[owner_p4est_gface] + push!(hanging_vertices_owner_cell_and_lface, + (owner_cell, n_cell_vertices+n_cell_edges+PXEST_2_GRIDAP_FACE[p4est_lface],half)) + else + push!(hanging_vertices_owner_cell_and_lface,(-1, -1,-1)) + end + end start_gridap_vertices = (cell - 1) * n_cell_vertices gridap_cell_vertices_data[start_gridap_vertices+lvertex] = num_regular_faces[1] + - owner_gridap_gface_to_hanging_vertex[owner_gridap_gface] + owner_p4est_gface_to_hanging_vertex[owner_p4est_gface] end if (Dc==3) diff --git a/src/RedistributeTools.jl b/src/RedistributeTools.jl index 3d7e121..32a6a94 100644 --- a/src/RedistributeTools.jl +++ b/src/RedistributeTools.jl @@ -114,6 +114,7 @@ function _allocate_cell_wise_dofs(cell_to_ldofs) end end end + function get_glue_components(glue::GridapDistributed.RedistributeGlue,reverse::Val{false}) return glue.lids_rcv, glue.lids_snd, glue.parts_rcv, glue.parts_snd, glue.new2old @@ -170,14 +171,32 @@ function _allocate_cell_wise_dofs(cell_to_ldofs) snd_data, rcv_data, cell_dof_values_new = caches lids_rcv, lids_snd, parts_rcv, parts_snd, new2old = get_glue_components(glue,Val(reverse)) + + map(model_new.parts,lids_rcv,lids_snd,parts_rcv,parts_snd,new2old) do rank, + lids_rcv, + lids_snd, + parts_rcv, + parts_snd, + new2old + print("$(rank): [$(lids_rcv)]: $(parts_rcv)"); print("\n") + print("$(rank): [$(lids_snd)]: $(parts_snd)"); print("\n") + print("$(rank): [$(new2old)]"); print("\n") + end cell_dof_values_old = change_parts(cell_dof_values_old,get_parts(glue);default=[]) cell_dof_ids_new = change_parts(cell_dof_ids_new,get_parts(glue);default=[[]]) _pack_snd_data!(snd_data,cell_dof_values_old,lids_snd) + map(model_new.parts,snd_data,cell_dof_values_old,rcv_data) do rank,snd_data,cell_dof_values_old, rcv_data + print("$(rank): [$(snd_data)]"); print("\n") + print("$(rank): [$(cell_dof_values_old)]"); print("\n") + print("$(rank) XXX: [$(rcv_data)]"); print("\n") + end + graph=ExchangeGraph(parts_snd,parts_rcv) t=exchange!(rcv_data,snd_data,graph) + wait(t) # We have to build the owned part of "cell_dof_values_new" out of # 1. cell_dof_values_old (for those cells s.t. new2old[:]!=0) @@ -186,21 +205,34 @@ function _allocate_cell_wise_dofs(cell_to_ldofs) cell_dof_values_old, new2old) - wait(t) - _unpack_rcv_data!(cell_dof_values_new,rcv_data,lids_rcv) # Now that every part knows it's new owned dofs, exchange ghosts new_parts = model_new.parts cell_dof_values_new = change_parts(cell_dof_values_new,new_parts) if i_am_in(new_parts) + map(partition(get_cell_gids(model_new))) do indices + println("!!!$(part_id(indices))!!!: $(local_to_global(indices))") + println("!!!$(part_id(indices))!!!: $(local_to_owner(indices))") + end cache = fetch_vector_ghost_values_cache(cell_dof_values_new, partition(get_cell_gids(model_new))) + map(partition(get_cell_gids(model_new)),cache) do indices,cache + println("!!!$(part_id(indices))!!!: !!!$(cache.cache.local_indices_snd)!!!: $(cache.cache.local_indices_rcv)") + end fetch_vector_ghost_values!(cell_dof_values_new,cache) |> wait end return cell_dof_values_new end + function _get_cell_dof_ids_inner_space(s::FESpace) + get_cell_dof_ids(s) + end + + function _get_cell_dof_ids_inner_space(s::FESpaceWithLinearConstraints) + get_cell_dof_ids(s.space) + end + function get_redistribute_free_values_cache(fv_new::Union{PVector,Nothing}, Uh_new::Union{GridapDistributed.DistributedSingleFieldFESpace,VoidDistributedFESpace}, fv_old::Union{PVector,Nothing}, @@ -210,7 +242,7 @@ function _allocate_cell_wise_dofs(cell_to_ldofs) glue::GridapDistributed.RedistributeGlue; reverse=false) cell_dof_values_old = !isa(fv_old,Nothing) ? map(scatter_free_and_dirichlet_values,local_views(Uh_old),local_views(fv_old),dv_old) : nothing - cell_dof_ids_new = !isa(fv_new,Nothing) ? map(get_cell_dof_ids, local_views(Uh_new)) : nothing + cell_dof_ids_new = !isa(fv_new,Nothing) ? map(_get_cell_dof_ids_inner_space, local_views(Uh_new)) : nothing caches = get_redistribute_cell_dofs_cache(cell_dof_values_old,cell_dof_ids_new,model_new,glue;reverse=reverse) return caches @@ -240,7 +272,7 @@ function _allocate_cell_wise_dofs(cell_to_ldofs) reverse=false) cell_dof_values_old = !isa(fv_old,Nothing) ? map(scatter_free_and_dirichlet_values,local_views(Uh_old),local_views(fv_old),dv_old) : nothing - cell_dof_ids_new = !isa(fv_new,Nothing) ? map(get_cell_dof_ids, local_views(Uh_new)) : nothing + cell_dof_ids_new = !isa(fv_new,Nothing) ? map(_get_cell_dof_ids_inner_space, local_views(Uh_new)) : nothing cell_dof_values_new = redistribute_cell_dofs!(caches,cell_dof_values_old,cell_dof_ids_new,model_new,glue;reverse=reverse) # Gather the new free dofs @@ -257,7 +289,7 @@ function _allocate_cell_wise_dofs(cell_to_ldofs) reverse=false) cell_dof_values_old = !isa(uh_old,Nothing) ? map(get_cell_dof_values,local_views(uh_old)) : nothing - cell_dof_ids_new = !isa(Uh_new,VoidDistributedFESpace) ? map(get_cell_dof_ids,local_views(Uh_new)) : nothing + cell_dof_ids_new = !isa(Uh_new,VoidDistributedFESpace) ? map(_get_cell_dof_ids_inner_space,local_views(Uh_new)) : nothing map(cell_dof_values_old,partition(Uh_new.gids)) do cell_dof_values_old, indices print("[$(part_id(indices))]: $(cell_dof_values_old)"); print("\n") diff --git a/test.jl b/test.jl index 825e267..38fb671 100644 --- a/test.jl +++ b/test.jl @@ -6,143 +6,164 @@ using MPI # Define integration mesh and quadrature -order=1 +order=2 # Define manufactured functions u(x) = x[1]+x[2]^order f(x) = -Δ(u)(x) degree = 2*order+1 MPI.Init() -ranks=distribute_with_mpi(LinearIndices((2,))) +ranks=distribute_with_mpi(LinearIndices((4,))) coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,1) -map(ranks,partition(dmodel.dmodel.face_gids[end])) do rank, indices - print("$(rank): $(local_to_owner(indices))"); print("\n") -end -reffe=ReferenceFE(lagrangian,Float64,order) -VH=FESpace(dmodel,reffe,conformity=:H1;dirichlet_tags="boundary") -UH=TrialFESpace(VH,u) -ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices - flags=zeros(Cint,length(indices)) - flags.=-1 - if (rank==1) - flags[1:2].=[refine_flag,nothing_flag] - elseif (rank==2) - flags[1:2].=[nothing_flag,nothing_flag] +function test(ranks,dmodel) + map(ranks,partition(dmodel.dmodel.face_gids[end])) do rank, indices + print("$(rank): $(local_to_owner(indices))"); print("\n") + end + reffe=ReferenceFE(lagrangian,Float64,order) + VH=FESpace(dmodel,reffe,conformity=:H1;dirichlet_tags="boundary") + UH=TrialFESpace(VH,u) + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices + flags=zeros(Cint,length(indices)) + flags.=nothing_flag + + flags[1]=refine_flag + flags[own_length(indices)]=refine_flag + + # To create some unbalance + if (rank%2==0 && own_length(indices)>1) + flags[div(own_length(indices),2)]=refine_flag + end + + print("rank: $(rank) flags: $(flags)"); print("\n") + flags + end + fmodel,glue=refine(dmodel,ref_coarse_flags); + map(ranks,glue) do rank, glue + if rank==2 + print(glue.n2o_faces_map[end]); print("\n") + print(glue.refinement_rules); print("\n") + end end - print("rank: $(rank) flags: $(flags)"); print("\n") - flags -end -fmodel,glue=refine(dmodel,ref_coarse_flags); -map(ranks,glue) do rank, glue - if rank==2 - print(glue.n2o_faces_map[end]); print("\n") - print(glue.refinement_rules); print("\n") - end -end - -Vh=FESpace(fmodel,reffe,conformity=:H1;dirichlet_tags="boundary") -map(ranks,partition(Vh.gids)) do rank, indices - print("$(rank): $(local_to_owner(indices))"); print("\n") - print("$(rank): $(local_to_global(indices))"); print("\n") -end -Uh=TrialFESpace(Vh,u) - -ΩH = Triangulation(dmodel) -dΩH = Measure(ΩH,degree) - -aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH -bH(v) = ∫(v*f)*dΩH - -op = AffineFEOperator(aH,bH,UH,VH) -uH = solve(op) -e = u - uH - -# # Compute errors -el2 = sqrt(sum( ∫( e*e )*dΩH )) -eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩH )) - -tol=1e-8 -@assert el2 < tol -@assert eh1 < tol - - -Ωh = Triangulation(fmodel) -dΩh = Measure(Ωh,degree) - -ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh -bh(v) = ∫(v*f)*dΩh - -op = AffineFEOperator(ah,bh,Uh,Vh) -uh = solve(op) -e = u - uh - -# # Compute errors -el2 = sqrt(sum( ∫( e*e )*dΩh )) -eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩh )) - -tol=1e-8 -@assert el2 < tol -@assert eh1 < tol - -# prolongation via interpolation -uHh=interpolate(uH,Uh) -e = uh - uHh -el2 = sqrt(sum( ∫( e*e )*dΩh )) -tol=1e-8 -@assert el2 < tol - -# prolongation via L2-projection -# Coarse FEFunction -> Fine FEFunction, by projection -ah(u,v) = ∫(v⋅u)*dΩh -lh(v) = ∫(v⋅uH)*dΩh -oph = AffineFEOperator(ah,lh,Uh,Vh) -uHh = solve(oph) -e = uh - uHh -el2 = sqrt(sum( ∫( e*e )*dΩh )) -tol=1e-8 -@assert el2 < tol - -# restriction via interpolation -uhH=interpolate(uh,UH) -e = uH - uhH -el2 = sqrt(sum( ∫( e*e )*dΩh )) -tol=1e-8 -@assert el2 < tol - -# restriction via L2-projection -dΩhH = Measure(ΩH,Ωh,2*order) -aH(u,v) = ∫(v⋅u)*dΩH -lH(v) = ∫(v⋅uh)*dΩhH -oph = AffineFEOperator(aH,lH,UH,VH) -uhH = solve(oph) -e = uH - uhH -el2 = sqrt(sum( ∫( e*e )*dΩH )) - -fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); -Vhred=FESpace(fmodel_red,reffe,conformity=:H1;dirichlet_tags="boundary") -Uhred=TrialFESpace(Vhred,u) - -Ωhred = Triangulation(fmodel_red) -dΩhred = Measure(Ωhred,degree) - -ahred(u,v) = ∫( ∇(v)⊙∇(u) )*dΩhred -bhred(v) = ∫(v*f)*dΩhred - -op = AffineFEOperator(ahred,bhred,Uhred,Vhred) -uhred = solve(op) -e = u - uhred -el2 = sqrt(sum( ∫( e*e )*dΩhred )) -@assert el2 < tol - -uhred2 = GridapP4est.redistribute_fe_function(uh,Vhred,fmodel_red,red_glue) -e = u - uhred2 -el2 = sqrt(sum( ∫( e*e )*dΩhred )) -tol=1e-8 -println(el2) -@assert el2 < tol + Vh=FESpace(fmodel,reffe,conformity=:H1;dirichlet_tags="boundary") + map(ranks,partition(Vh.gids)) do rank, indices + print("$(rank): $(local_to_owner(indices))"); print("\n") + print("$(rank): $(local_to_global(indices))"); print("\n") + end + Uh=TrialFESpace(Vh,u) + + ΩH = Triangulation(dmodel) + dΩH = Measure(ΩH,degree) + + aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH + bH(v) = ∫(v*f)*dΩH + + op = AffineFEOperator(aH,bH,UH,VH) + uH = solve(op) + e = u - uH + + # # Compute errors + el2 = sqrt(sum( ∫( e*e )*dΩH )) + eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩH )) + + tol=1e-8 + @assert el2 < tol + @assert eh1 < tol + + + Ωh = Triangulation(fmodel) + dΩh = Measure(Ωh,degree) + + ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh + bh(v) = ∫(v*f)*dΩh + + op = AffineFEOperator(ah,bh,Uh,Vh) + uh = solve(op) + e = u - uh + + # # Compute errors + el2 = sqrt(sum( ∫( e*e )*dΩh )) + eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩh )) + + tol=1e-8 + @assert el2 < tol + @assert eh1 < tol + + # prolongation via interpolation + uHh=interpolate(uH,Uh) + e = uh - uHh + el2 = sqrt(sum( ∫( e*e )*dΩh )) + tol=1e-8 + @assert el2 < tol + + # prolongation via L2-projection + # Coarse FEFunction -> Fine FEFunction, by projection + ahp(u,v) = ∫(v⋅u)*dΩh + lhp(v) = ∫(v⋅uH)*dΩh + oph = AffineFEOperator(ahp,lhp,Uh,Vh) + uHh = solve(oph) + e = uh - uHh + el2 = sqrt(sum( ∫( e*e )*dΩh )) + tol=1e-8 + @assert el2 < tol + + # restriction via interpolation + uhH=interpolate(uh,UH) + e = uH - uhH + el2 = sqrt(sum( ∫( e*e )*dΩh )) + tol=1e-8 + @assert el2 < tol + + # restriction via L2-projection + dΩhH = Measure(ΩH,Ωh,2*order) + aHp(u,v) = ∫(v⋅u)*dΩH + lHp(v) = ∫(v⋅uh)*dΩhH + oph = AffineFEOperator(aHp,lHp,UH,VH) + uhH = solve(oph) + e = uH - uhH + el2 = sqrt(sum( ∫( e*e )*dΩH )) + + fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); + Vhred=FESpace(fmodel_red,reffe,conformity=:H1;dirichlet_tags="boundary") + Uhred=TrialFESpace(Vhred,u) + + Ωhred = Triangulation(fmodel_red) + dΩhred = Measure(Ωhred,degree) + + ahred(u,v) = ∫( ∇(v)⊙∇(u) )*dΩhred + bhred(v) = ∫(v*f)*dΩhred + + op = AffineFEOperator(ahred,bhred,Uhred,Vhred) + uhred = solve(op) + e = u - uhred + el2 = sqrt(sum( ∫( e*e )*dΩhred )) + @assert el2 < tol + + + map(Vhred.spaces) do space + if (MPI.Comm_rank(MPI.COMM_WORLD)==1) + println("XXXX: $(get_cell_dof_ids(space.space))") + end + end + + uhred2 = GridapP4est.redistribute_fe_function(uh,Vhred,fmodel_red,red_glue) + e = u - uhred2 + el2 = sqrt(sum( ∫( e*e )*dΩhred )) + tol=1e-8 + println(el2) + @assert el2 < tol + + fmodel_red +end +function f() + rdmodel=dmodel + for i=1:5 + rdmodel=test(ranks,rdmodel) + end +end +f()