From 26e34e130595f69f40a93cd9520be29ef811ac3f Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Tue, 11 Apr 2023 23:05:10 +1000 Subject: [PATCH] Adding non-conforming octree tests --- Manifest.toml | 26 +- ...ingOctreeDistributedDiscreteModelsTests.jl | 521 ++++++++++++++++++ test/OctreeDistributedDiscreteModelsTests.jl | 2 +- test/non_conforming_octrees_wip.jl | 518 ----------------- test/runtests.jl | 3 + 5 files changed, 538 insertions(+), 532 deletions(-) create mode 100644 test/NonConformingOctreeDistributedDiscreteModelsTests.jl delete mode 100644 test/non_conforming_octrees_wip.jl diff --git a/Manifest.toml b/Manifest.toml index 284457c..0bc8822 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -38,9 +38,9 @@ version = "1.1.1" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "a89acc90c551067cd84119ff018619a1a76c6277" +git-tree-sha1 = "38911c7737e123b28182d89027f4216cfc8a9da7" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.2.1" +version = "7.4.3" [[deps.ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra", "SparseArrays"] @@ -57,9 +57,9 @@ uuid = "15f4f7f2-30c1-5605-9d31-71845cf9641f" version = "0.2.0" [[deps.BSON]] -git-tree-sha1 = "86e9781ac28f4e80e9b98f7f96eae21891332ac2" +git-tree-sha1 = "2208958832d6e1b59e49f53697483a84ca8d664e" uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" -version = "0.3.6" +version = "0.3.7" [[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" @@ -215,7 +215,7 @@ uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" [[deps.Gridap]] deps = ["AbstractTrees", "BSON", "BlockArrays", "Combinatorics", "DataStructures", "DocStringExtensions", "FastGaussQuadrature", "FileIO", "FillArrays", "ForwardDiff", "JLD2", "JSON", "LineSearches", "LinearAlgebra", "NLsolve", "NearestNeighbors", "PolynomialBases", "QuadGK", "Random", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "Test", "WriteVTK"] -git-tree-sha1 = "414c21416961116fe48f7dccd65419b0495e83a8" +git-tree-sha1 = "c06fafb41527457d8be668a82e7985a8bc31bf0e" repo-rev = "refinement-rule-subface-glue" repo-url = "git@github.com:gridap/Gridap.jl.git" uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" @@ -376,9 +376,9 @@ version = "2.28.0+0" [[deps.MicrosoftMPI_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "a16aa086d335ed7e0170c5265247db29172af2f9" +git-tree-sha1 = "a8027af3d1743b3bfae34e54872359fdebb31422" uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf" -version = "10.1.3+2" +version = "10.1.3+4" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -438,9 +438,9 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" [[deps.OrderedCollections]] -git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" +git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.4.1" +version = "1.6.0" [[deps.P4est_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Pkg", "TOML", "Zlib_jll"] @@ -584,9 +584,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[deps.StatsAPI]] deps = ["LinearAlgebra"] -git-tree-sha1 = "f9af7f195fb13589dd2e2d57fdb401717d2eb1f6" +git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7" uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.5.0" +version = "1.6.0" [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] @@ -613,9 +613,9 @@ version = "1.0.1" [[deps.TranscodingStreams]] deps = ["Random", "Test"] -git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f" +git-tree-sha1 = "0b829474fed270a4b0ab07117dce9b9a2fa7581a" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.9.11" +version = "0.9.12" [[deps.UUIDs]] deps = ["Random", "SHA"] diff --git a/test/NonConformingOctreeDistributedDiscreteModelsTests.jl b/test/NonConformingOctreeDistributedDiscreteModelsTests.jl new file mode 100644 index 0000000..b4bab67 --- /dev/null +++ b/test/NonConformingOctreeDistributedDiscreteModelsTests.jl @@ -0,0 +1,521 @@ +module NonConformingOctreeDistributedDiscreteModelsTests + using P4est_wrapper + using GridapP4est + using Gridap + using PartitionedArrays + using GridapDistributed + using MPI + using Gridap.FESpaces + using FillArrays + # Generate a local numbering of vertices that includes hanging vertices + # Generate a local numbering of faces out of the one generated by vertices (automatic? to confirm) + + # Establish the correspondence among local numbering of vertices and p4est global numbering + # Establish the correspondence among local numbering of faces and p4est global numbering + + # Generate a global numbering of (regular,hanging) vertices? + # Generate a global numbering of (regular,hanging) faces? + + function setup_model(perm) + @assert perm ∈ (1,2,3,4) + # + # 3-------4-------6 + # | | | + # | | | + # | | | + # 1-------2-------5 + # + ptr = [ 1, 5, 9 ] + if (perm==1) + data = [ 1,2,3,4, 2,5,4,6 ] + elseif (perm==2) + data = [ 1,2,3,4, 6,4,5,2 ] + elseif (perm==3) + data = [ 4,3,2,1, 2,5,4,6 ] + elseif (perm==4) + data = [ 4,3,2,1, 6,4,5,2 ] + end + cell_vertex_lids = Gridap.Arrays.Table(data,ptr) + node_coordinates = Vector{Point{2,Float64}}(undef,6) + node_coordinates[1]=Point{2,Float64}(0.0,0.0) + node_coordinates[2]=Point{2,Float64}(1.0,0.0) + node_coordinates[3]=Point{2,Float64}(0.0,1.0) + node_coordinates[4]=Point{2,Float64}(1.0,1.0) + node_coordinates[5]=Point{2,Float64}(2.0,0.0) + node_coordinates[6]=Point{2,Float64}(2.0,1.0) + + polytope=QUAD + scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1) + cell_types=collect(Fill(1,length(cell_vertex_lids))) + cell_reffes=[scalar_reffe] + grid = Gridap.Geometry.UnstructuredGrid(node_coordinates, + cell_vertex_lids, + cell_reffes, + cell_types, + Gridap.Geometry.NonOriented()) + m=Gridap.Geometry.UnstructuredDiscreteModel(grid) + labels = get_face_labeling(m) + labels.d_to_dface_to_entity[1]=[7,7,7,7,7,7] + if (perm==1 || perm==2) + labels.d_to_dface_to_entity[2]=[7,7,7,0,7,7,7] + elseif (perm==3 || perm==4) + labels.d_to_dface_to_entity[2]=[7,7,0,7,7,7,7] + end + add_tag!(labels,"boundary",[7]) + m + end + + + ## Better to use a C-enum. But I did not use it in order to keep the Julia + ## version of this C example as simple as possible + const nothing_flag = Cint(0) + const refine_flag = Cint(1) + + ## Refine those cells with even identifier (0,2,4,6,8,...) + ## Leave untouched cells with odd identifier (1,3,5,7,9,...) + function allocate_and_set_refinement_and_coarsening_flags(forest_ptr::Ptr{p4est_t}) + forest = forest_ptr[] + tree = p4est_tree_array_index(forest.trees, 0)[] + return [i != 1 ? nothing_flag : refine_flag for i = 1:tree.quadrants.elem_count] + end + + MPI.Init() + parts = get_part_ids(MPIBackend(), 1) + # run(parts, (1, 1)) + # MPI.Finalize() + + function test(perm,order) + # This is for debuging + Dc=2 + domain = (0, 1, 0, 1) + subdomains = (1, 1) + coarse_model = setup_model(perm) + model = OctreeDistributedDiscreteModel(parts, coarse_model, 0) + + ref_coarse_flags=map_parts(parts) do _ + [refine_flag,nothing_flag] + #allocate_and_set_refinement_and_coarsening_flags(model.ptr_pXest) + end + dmodel,non_conforming_glue=refine(model,ref_coarse_flags) + + p4est_vtk_write_file(dmodel.ptr_pXest, C_NULL, string("adapted_forest")) + + # FE Spaces + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(dmodel,reffe,dirichlet_tags="boundary") + U = TrialFESpace(V) + + function _h_refined_reffe(reffe::Tuple{<:Lagrangian,Any,Any}) + (reffe[1],(reffe[2][1],2*reffe[2][2]),reffe[3]) + end + + basis, reffe_args,reffe_kwargs = reffe + cell_reffe = ReferenceFE(QUAD,basis,reffe_args...;reffe_kwargs...) + reffe_cell = cell_reffe + + h_refined_reffe = _h_refined_reffe(reffe) + basis, reffe_args,reffe_kwargs = h_refined_reffe + cell_reffe_h_refined = ReferenceFE(QUAD,basis,reffe_args...;reffe_kwargs...) + reffe_cell_h_refined = cell_reffe_h_refined + + dof_basis_h_refined = Gridap.CellData.get_dof_basis(reffe_cell_h_refined) + + coarse_shape_funs=Gridap.ReferenceFEs.get_shapefuns(reffe_cell) + ref_constraints=evaluate(dof_basis_h_refined,coarse_shape_funs) + + + rr = Gridap.Adaptivity.RedRefinementRule(QUAD) + face_subface_ldof_to_cell_ldof = Gridap.Adaptivity.coarse_nodes_above_fine_nodes(rr,(order,order),1) + + + # To-think: might this info go to the glue? + # If it is required in different scenarios, I would say it may make sense + function _generate_hanging_faces_to_cell_and_lface(num_regular_faces, + num_hanging_faces, + gridap_cell_faces) + # Locate for each hanging vertex a cell to which it belongs + # and local position within that cell + hanging_faces_to_cell = Vector{Int}(undef, num_hanging_faces) + hanging_faces_to_lface = Vector{Int}(undef, num_hanging_faces) + for cell=1:length(gridap_cell_faces) + s=gridap_cell_faces.ptrs[cell] + e=gridap_cell_faces.ptrs[cell+1] + l=e-s + for j=1:l + fid=gridap_cell_faces.data[s+j-1] + if fid>num_regular_faces + fid_hanging=fid-num_regular_faces + hanging_faces_to_cell[fid_hanging]=cell + hanging_faces_to_lface[fid_hanging]=j + end + end + end + hanging_faces_to_cell, hanging_faces_to_lface + end + + + function _generate_hanging_faces_owner_face_dofs(num_hanging_faces, + face_dofs, + hanging_faces_glue, + cell_dof_ids) + + cache = array_cache(cell_dof_ids) + ptrs=Vector{Int}(undef, num_hanging_faces+1) + 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]) + 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]) + data_owner_face_dofs[s+j-1]=current_cell_dof_ids[ldof] + end + end + Gridap.Arrays.Table(data_owner_face_dofs, ptrs) + end + + function _generate_constraints!(Df, + Dc, + cell_faces, + num_hanging_faces, + hanging_faces_to_cell, + hanging_faces_to_lface, + hanging_faces_owner_face_dofs, + hanging_faces_glue, + face_subface_ldof_to_cell_ldof, + face_dofs, + face_own_dofs, + subface_own_dofs, + cell_dof_ids, + node_permutations, + owner_faces_pindex, + owner_faces_lids, + sDOF_to_dof, + sDOF_to_dofs, + sDOF_to_coeffs) + + @assert Dc==2 || Dc==3 + @assert 0 ≤ Df < Dc + num_vertices = 2^Dc + num_edges = (Dc==2 ? 0 : 12) + num_faces = (2*Dc) + + offset=0 + if (Df ≥ 1) + offset+=num_vertices + if (Df == 2) + offset+=num_edges + end + end + + cache_dof_ids=array_cache(cell_dof_ids) + cache_faces=array_cache(cell_faces) + for fid_hanging=1:num_hanging_faces + cell=hanging_faces_to_cell[fid_hanging] + 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 = ocell_lface - num_vertices + oface=getindex!(cache_faces,cell_faces,ocell)[ocell_lface] + oface_lid,_=owner_faces_lids[oface] + pindex=owner_faces_pindex[oface_lid] + for ((ldof,dof),ldof_subface) in zip(enumerate(face_own_dofs[offset+lface]),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+num_vertices]) + pifdof=node_permutations[pindex][ifdof] + println("XXXX: $(ifdof) $(pifdof)") + ldof_coarse=face_dofs[ocell_lface+num_vertices][pifdof] + coeffs[ifdof]=ref_constraints[face_subface_ldof_to_cell_ldof[ocell_lface][subface][ldof_subface],ldof_coarse] + end + push!(sDOF_to_coeffs,coeffs) + end + end + end + + + # count how many different owner faces + # for each owner face + # track the global IDs of its face vertices from the perspective of the subfaces + # for each owner face + # compute permutation id + function _compute_owner_faces_pindex_and_lids(Dc, + num_hanging_faces, + hanging_faces_glue, + hanging_faces_to_cell, + hanging_faces_to_lface, + cell_vertices, + cell_faces, + lface_to_cvertices, + pindex_to_cfvertex_to_fvertex) + num_vertices=2^Dc + num_owner_faces=0 + owner_faces_lids=Dict{Int,Tuple{Int,Int,Int}}() + for fid_hanging=1:num_hanging_faces + ocell,ocell_lface,_=hanging_faces_glue[fid_hanging] + owner_face=cell_faces[ocell][ocell_lface-num_vertices] + if !(haskey(owner_faces_lids,owner_face)) + num_owner_faces+=1 + owner_faces_lids[owner_face]=(num_owner_faces,ocell,ocell_lface) + end + end + + num_face_vertices=length(first(lface_to_cvertices)) + owner_face_vertex_ids=Vector{Int}(undef,num_face_vertices*num_owner_faces) + + for fid_hanging=1:num_hanging_faces + ocell,ocell_lface,subface=hanging_faces_glue[fid_hanging] + 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] + owner_face=cell_faces[ocell][ocell_lface-num_vertices] + owner_face_lid,_=owner_faces_lids[owner_face] + owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+subface]=vertex + end + + owner_faces_pindex=Vector{Int}(undef,num_owner_faces) + for owner_face in keys(owner_faces_lids) + (owner_face_lid,ocell,ocell_lface)=owner_faces_lids[owner_face] + # Compute permutation id by comparing + # 1. cell_vertices[ocell][ocell_lface] + # 2. owner_face_vertex_ids + pindexfound = false + cfvertex_to_cvertex = lface_to_cvertices[ocell_lface-num_vertices] + for (pindex, cfvertex_to_fvertex) in enumerate(pindex_to_cfvertex_to_fvertex) + found = true + for (cfvertex,fvertex) in enumerate(cfvertex_to_fvertex) + vertex1 = owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fvertex] + cvertex = cfvertex_to_cvertex[cfvertex] + vertex2 = cell_vertices[ocell][cvertex] + if vertex1 != vertex2 + found = false + break + end + end + if found + owner_faces_pindex[owner_face_lid] = pindex + pindexfound = true + break + end + end + @assert pindexfound "Valid pindex not found" + end + owner_faces_pindex, owner_faces_lids + end + + function generate_constraints(dmodel, + V, + reffe, + non_conforming_glue, + ref_constraints, + face_subface_ldof_to_cell_ldof) + num_regular_faces, + num_hanging_faces, + gridap_cell_faces, + hanging_faces_glue = non_conforming_glue + + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = map_parts(gridap_cell_faces[1], + num_regular_faces[1], + num_hanging_faces[1], + hanging_faces_glue[1], + gridap_cell_faces[2], + num_regular_faces[2], + num_hanging_faces[2], + hanging_faces_glue[2], + dmodel.dmodel.models, V.spaces) do gridap_cell_vertices, + num_regular_vertices, num_hanging_vertices, + hanging_vertices_owner_cell_and_lface, + gridap_cell_faces, + num_regular_faces, num_hanging_faces, + hanging_faces_owner_cell_and_lface, + model, V + + # Locate for each hanging vertex a cell to which it belongs + # and local position within that cell + hanging_vertices_to_cell, + hanging_vertices_to_lvertex = _generate_hanging_faces_to_cell_and_lface(num_regular_vertices, + num_hanging_vertices, + gridap_cell_vertices) + + + # Locate for each hanging facet a cell to which it belongs + # and local position within that cell + hanging_faces_to_cell, + hanging_faces_to_lface = _generate_hanging_faces_to_cell_and_lface(num_regular_faces, + num_hanging_faces, + gridap_cell_faces) + + basis, reffe_args,reffe_kwargs = reffe + cell_reffe = ReferenceFE(QUAD,basis,reffe_args...;reffe_kwargs...) + reffe_cell = cell_reffe + + cell_dof_ids = get_cell_dof_ids(V) + face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(reffe_cell) + face_dofs = Gridap.ReferenceFEs.get_face_dofs(reffe_cell) + + hanging_vertices_owner_face_dofs = _generate_hanging_faces_owner_face_dofs(num_hanging_vertices, + face_dofs, + hanging_vertices_owner_cell_and_lface, + cell_dof_ids) + + hanging_faces_owner_face_dofs = _generate_hanging_faces_owner_face_dofs(num_hanging_faces, + face_dofs, + hanging_faces_owner_cell_and_lface, + cell_dof_ids) + + sDOF_to_dof = Int[] + sDOF_to_dofs = Vector{Int}[] + sDOF_to_coeffs = Vector{Float64}[] + + Dc = 2 + facet_polytope = Dc==2 ? SEGMENT : QUAD + basis, reffe_args,reffe_kwargs = reffe + face_reffe = ReferenceFE(facet_polytope,basis,reffe_args...;reffe_kwargs...) + pindex_to_cfvertex_to_fvertex = Gridap.ReferenceFEs.get_vertex_permutations(facet_polytope) + + lface_to_cvertices = Gridap.ReferenceFEs.get_faces(QUAD,Dc-1,0) + owner_faces_pindex, owner_faces_lids=_compute_owner_faces_pindex_and_lids(Dc, + num_hanging_faces, + hanging_faces_owner_cell_and_lface, + hanging_faces_to_cell, + hanging_faces_to_lface, + gridap_cell_vertices, + gridap_cell_faces, + lface_to_cvertices, + pindex_to_cfvertex_to_fvertex) + + nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope,[reffe_args[2]]) + node_permutations=Gridap.ReferenceFEs._compute_node_permutations(facet_polytope,nodes) + + hanging_lvertex_within_first_subface = 2^(Dc-1) + subface_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(face_reffe) + vertex_subface_own_dofs = subface_own_dofs[hanging_lvertex_within_first_subface] + face_dofs = Gridap.ReferenceFEs.get_face_dofs(cell_reffe) + + _generate_constraints!(0, + Dc, + gridap_cell_faces, + num_hanging_vertices, + hanging_vertices_to_cell, + hanging_vertices_to_lvertex, + hanging_vertices_owner_face_dofs, + hanging_vertices_owner_cell_and_lface, + face_subface_ldof_to_cell_ldof, + face_dofs, + face_own_dofs, + vertex_subface_own_dofs, + cell_dof_ids, + node_permutations, + owner_faces_pindex, + owner_faces_lids, + sDOF_to_dof, + sDOF_to_dofs, + sDOF_to_coeffs) + + face_subface_own_dofs = subface_own_dofs[end] + _generate_constraints!(1, + Dc, + gridap_cell_faces, + num_hanging_faces, + hanging_faces_to_cell, + hanging_faces_to_lface, + hanging_faces_owner_face_dofs, + hanging_faces_owner_cell_and_lface, + face_subface_ldof_to_cell_ldof, + face_dofs, + face_own_dofs, + face_subface_own_dofs, + cell_dof_ids, + node_permutations, + owner_faces_pindex, + owner_faces_lids, + sDOF_to_dof, + sDOF_to_dofs, + sDOF_to_coeffs) + sDOF_to_dof, Gridap.Arrays.Table(sDOF_to_dofs), Gridap.Arrays.Table(sDOF_to_coeffs) + end + end + + num_cells(dmodel.dmodel.models.part) + + sDOF_to_dof, sDOF_to_dofs,sDOF_to_coeffs= + generate_constraints(dmodel, V, reffe, + non_conforming_glue, ref_constraints, face_subface_ldof_to_cell_ldof) + println(sDOF_to_dof) + println(sDOF_to_dofs) + println(sDOF_to_coeffs) + + # Define manufactured functions + u(x) = x[1]+x[2]^order + f(x) = -Δ(u)(x) + + map_parts(dmodel.dmodel.models,V.spaces,U.spaces,sDOF_to_dof,sDOF_to_dofs,sDOF_to_coeffs) do model,V,U,sDOF_to_dof,sDOF_to_dofs,sDOF_to_coeffs + + println(get_cell_dof_ids(V)) + fl=get_face_labeling(model) + t=GridapDistributed.get_grid_topology(model) + println(fl.d_to_dface_to_entity[1]) + println(fl.d_to_dface_to_entity[2]) + println(Gridap.Geometry.get_faces(t,2,0)) + + Vc = FESpaceWithLinearConstraints( + sDOF_to_dof, + sDOF_to_dofs, + sDOF_to_coeffs, + V) + Uc = TrialFESpace(Vc,u) + + # Define integration mesh and quadrature + degree = 2*order+1 + Ω = Triangulation(model) + dΩ = Measure(Ω,degree) + + a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ + b(v) = ∫(v*f)*dΩ + + # op = AffineFEOperator(a,b,U,V0) + op = AffineFEOperator(a,b,Uc,Vc) + uh = solve(op) + + # Define exact solution and error + e = u - uh + + # Compute errors + el2 = sqrt(sum( ∫( e*e )*dΩ )) + eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩ )) + + tol=1e-8 + @assert el2 < tol + @assert eh1 < tol + end + + end + + test(1,1) + test(1,2) + test(1,3) + + test(2,1) + test(2,2) + test(2,3) + + test(3,1) + test(3,2) + test(3,3) + + test(4,1) + test(4,2) + test(4,3) + +end \ No newline at end of file diff --git a/test/OctreeDistributedDiscreteModelsTests.jl b/test/OctreeDistributedDiscreteModelsTests.jl index 241b0a9..61da5b1 100644 --- a/test/OctreeDistributedDiscreteModelsTests.jl +++ b/test/OctreeDistributedDiscreteModelsTests.jl @@ -94,7 +94,7 @@ module OctreeDistributedDiscreteModelsTests map_parts(model.dmodel.models,model_back.dmodel.models) do m1, m2 Ωh1 = Triangulation(m1) dΩh1 = Measure(Ωh1,2) - Ωh2 = Triangulation(m2) + Ωh2OctreeDistributedDiscreteModelsTests.jl = Triangulation(m2) dΩh2 = Measure(Ωh2,2) sum(∫(1)dΩh1) ≈ sum(∫(1)dΩh2) end diff --git a/test/non_conforming_octrees_wip.jl b/test/non_conforming_octrees_wip.jl deleted file mode 100644 index 1859e9d..0000000 --- a/test/non_conforming_octrees_wip.jl +++ /dev/null @@ -1,518 +0,0 @@ -using P4est_wrapper -using GridapP4est -using Gridap -using PartitionedArrays -using GridapDistributed -using MPI -using Gridap.FESpaces -using FillArrays -# Generate a local numbering of vertices that includes hanging vertices -# Generate a local numbering of faces out of the one generated by vertices (automatic? to confirm) - -# Establish the correspondence among local numbering of vertices and p4est global numbering -# Establish the correspondence among local numbering of faces and p4est global numbering - -# Generate a global numbering of (regular,hanging) vertices? -# Generate a global numbering of (regular,hanging) faces? - -function setup_model(perm) - @assert perm ∈ (1,2,3,4) - # - # 3-------4-------6 - # | | | - # | | | - # | | | - # 1-------2-------5 - # - ptr = [ 1, 5, 9 ] - if (perm==1) - data = [ 1,2,3,4, 2,5,4,6 ] - elseif (perm==2) - data = [ 1,2,3,4, 6,4,5,2 ] - elseif (perm==3) - data = [ 4,3,2,1, 2,5,4,6 ] - elseif (perm==4) - data = [ 4,3,2,1, 6,4,5,2 ] - end - cell_vertex_lids = Gridap.Arrays.Table(data,ptr) - node_coordinates = Vector{Point{2,Float64}}(undef,6) - node_coordinates[1]=Point{2,Float64}(0.0,0.0) - node_coordinates[2]=Point{2,Float64}(1.0,0.0) - node_coordinates[3]=Point{2,Float64}(0.0,1.0) - node_coordinates[4]=Point{2,Float64}(1.0,1.0) - node_coordinates[5]=Point{2,Float64}(2.0,0.0) - node_coordinates[6]=Point{2,Float64}(2.0,1.0) - - polytope=QUAD - scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1) - cell_types=collect(Fill(1,length(cell_vertex_lids))) - cell_reffes=[scalar_reffe] - grid = Gridap.Geometry.UnstructuredGrid(node_coordinates, - cell_vertex_lids, - cell_reffes, - cell_types, - Gridap.Geometry.NonOriented()) - m=Gridap.Geometry.UnstructuredDiscreteModel(grid) - labels = get_face_labeling(m) - labels.d_to_dface_to_entity[1]=[7,7,7,7,7,7] - if (perm==1 || perm==2) - labels.d_to_dface_to_entity[2]=[7,7,7,0,7,7,7] - elseif (perm==3 || perm==4) - labels.d_to_dface_to_entity[2]=[7,7,0,7,7,7,7] - end - add_tag!(labels,"boundary",[7]) - m -end - - -## Better to use a C-enum. But I did not use it in order to keep the Julia -## version of this C example as simple as possible -const nothing_flag = Cint(0) -const refine_flag = Cint(1) - -## Refine those cells with even identifier (0,2,4,6,8,...) -## Leave untouched cells with odd identifier (1,3,5,7,9,...) -function allocate_and_set_refinement_and_coarsening_flags(forest_ptr::Ptr{p4est_t}) - forest = forest_ptr[] - tree = p4est_tree_array_index(forest.trees, 0)[] - return [i != 1 ? nothing_flag : refine_flag for i = 1:tree.quadrants.elem_count] -end - -MPI.Init() -parts = get_part_ids(MPIBackend(), 1) -# run(parts, (1, 1)) -# MPI.Finalize() - -function test(perm,order) - # This is for debuging - Dc=2 - domain = (0, 1, 0, 1) - subdomains = (1, 1) - coarse_model = setup_model(perm) - model = OctreeDistributedDiscreteModel(parts, coarse_model, 0) - - ref_coarse_flags=map_parts(parts) do _ - [refine_flag,nothing_flag] - #allocate_and_set_refinement_and_coarsening_flags(model.ptr_pXest) - end - dmodel,non_conforming_glue=refine(model,ref_coarse_flags) - - p4est_vtk_write_file(dmodel.ptr_pXest, C_NULL, string("adapted_forest")) - - # FE Spaces - reffe = ReferenceFE(lagrangian,Float64,order) - V = TestFESpace(dmodel,reffe,dirichlet_tags="boundary") - U = TrialFESpace(V) - - function _h_refined_reffe(reffe::Tuple{<:Lagrangian,Any,Any}) - (reffe[1],(reffe[2][1],2*reffe[2][2]),reffe[3]) - end - - basis, reffe_args,reffe_kwargs = reffe - cell_reffe = ReferenceFE(QUAD,basis,reffe_args...;reffe_kwargs...) - reffe_cell = cell_reffe - - h_refined_reffe = _h_refined_reffe(reffe) - basis, reffe_args,reffe_kwargs = h_refined_reffe - cell_reffe_h_refined = ReferenceFE(QUAD,basis,reffe_args...;reffe_kwargs...) - reffe_cell_h_refined = cell_reffe_h_refined - - dof_basis_h_refined = Gridap.CellData.get_dof_basis(reffe_cell_h_refined) - - coarse_shape_funs=Gridap.ReferenceFEs.get_shapefuns(reffe_cell) - ref_constraints=evaluate(dof_basis_h_refined,coarse_shape_funs) - - - rr = Gridap.Adaptivity.RedRefinementRule(QUAD) - face_subface_ldof_to_cell_ldof = Gridap.Adaptivity.coarse_nodes_above_fine_nodes(rr,(order,order),1) - - - # To-think: might this info go to the glue? - # If it is required in different scenarios, I would say it may make sense - function _generate_hanging_faces_to_cell_and_lface(num_regular_faces, - num_hanging_faces, - gridap_cell_faces) - # Locate for each hanging vertex a cell to which it belongs - # and local position within that cell - hanging_faces_to_cell = Vector{Int}(undef, num_hanging_faces) - hanging_faces_to_lface = Vector{Int}(undef, num_hanging_faces) - for cell=1:length(gridap_cell_faces) - s=gridap_cell_faces.ptrs[cell] - e=gridap_cell_faces.ptrs[cell+1] - l=e-s - for j=1:l - fid=gridap_cell_faces.data[s+j-1] - if fid>num_regular_faces - fid_hanging=fid-num_regular_faces - hanging_faces_to_cell[fid_hanging]=cell - hanging_faces_to_lface[fid_hanging]=j - end - end - end - hanging_faces_to_cell, hanging_faces_to_lface - end - - - function _generate_hanging_faces_owner_face_dofs(num_hanging_faces, - face_dofs, - hanging_faces_glue, - cell_dof_ids) - - cache = array_cache(cell_dof_ids) - ptrs=Vector{Int}(undef, num_hanging_faces+1) - 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]) - 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]) - data_owner_face_dofs[s+j-1]=current_cell_dof_ids[ldof] - end - end - Gridap.Arrays.Table(data_owner_face_dofs, ptrs) - end - - function _generate_constraints!(Df, - Dc, - cell_faces, - num_hanging_faces, - hanging_faces_to_cell, - hanging_faces_to_lface, - hanging_faces_owner_face_dofs, - hanging_faces_glue, - face_subface_ldof_to_cell_ldof, - face_dofs, - face_own_dofs, - subface_own_dofs, - cell_dof_ids, - node_permutations, - owner_faces_pindex, - owner_faces_lids, - sDOF_to_dof, - sDOF_to_dofs, - sDOF_to_coeffs) - - @assert Dc==2 || Dc==3 - @assert 0 ≤ Df < Dc - num_vertices = 2^Dc - num_edges = (Dc==2 ? 0 : 12) - num_faces = (2*Dc) - - offset=0 - if (Df ≥ 1) - offset+=num_vertices - if (Df == 2) - offset+=num_edges - end - end - - cache_dof_ids=array_cache(cell_dof_ids) - cache_faces=array_cache(cell_faces) - for fid_hanging=1:num_hanging_faces - cell=hanging_faces_to_cell[fid_hanging] - 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 = ocell_lface - num_vertices - oface=getindex!(cache_faces,cell_faces,ocell)[ocell_lface] - oface_lid,_=owner_faces_lids[oface] - pindex=owner_faces_pindex[oface_lid] - for ((ldof,dof),ldof_subface) in zip(enumerate(face_own_dofs[offset+lface]),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+num_vertices]) - pifdof=node_permutations[pindex][ifdof] - println("XXXX: $(ifdof) $(pifdof)") - ldof_coarse=face_dofs[ocell_lface+num_vertices][pifdof] - coeffs[ifdof]=ref_constraints[face_subface_ldof_to_cell_ldof[ocell_lface][subface][ldof_subface],ldof_coarse] - end - push!(sDOF_to_coeffs,coeffs) - end - end - end - - - # count how many different owner faces - # for each owner face - # track the global IDs of its face vertices from the perspective of the subfaces - # for each owner face - # compute permutation id - function _compute_owner_faces_pindex_and_lids(Dc, - num_hanging_faces, - hanging_faces_glue, - hanging_faces_to_cell, - hanging_faces_to_lface, - cell_vertices, - cell_faces, - lface_to_cvertices, - pindex_to_cfvertex_to_fvertex) - num_vertices=2^Dc - num_owner_faces=0 - owner_faces_lids=Dict{Int,Tuple{Int,Int,Int}}() - for fid_hanging=1:num_hanging_faces - ocell,ocell_lface,_=hanging_faces_glue[fid_hanging] - owner_face=cell_faces[ocell][ocell_lface-num_vertices] - if !(haskey(owner_faces_lids,owner_face)) - num_owner_faces+=1 - owner_faces_lids[owner_face]=(num_owner_faces,ocell,ocell_lface) - end - end - - num_face_vertices=length(first(lface_to_cvertices)) - owner_face_vertex_ids=Vector{Int}(undef,num_face_vertices*num_owner_faces) - - for fid_hanging=1:num_hanging_faces - ocell,ocell_lface,subface=hanging_faces_glue[fid_hanging] - 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] - owner_face=cell_faces[ocell][ocell_lface-num_vertices] - owner_face_lid,_=owner_faces_lids[owner_face] - owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+subface]=vertex - end - - owner_faces_pindex=Vector{Int}(undef,num_owner_faces) - for owner_face in keys(owner_faces_lids) - (owner_face_lid,ocell,ocell_lface)=owner_faces_lids[owner_face] - # Compute permutation id by comparing - # 1. cell_vertices[ocell][ocell_lface] - # 2. owner_face_vertex_ids - pindexfound = false - cfvertex_to_cvertex = lface_to_cvertices[ocell_lface-num_vertices] - for (pindex, cfvertex_to_fvertex) in enumerate(pindex_to_cfvertex_to_fvertex) - found = true - for (cfvertex,fvertex) in enumerate(cfvertex_to_fvertex) - vertex1 = owner_face_vertex_ids[(owner_face_lid-1)*num_face_vertices+fvertex] - cvertex = cfvertex_to_cvertex[cfvertex] - vertex2 = cell_vertices[ocell][cvertex] - if vertex1 != vertex2 - found = false - break - end - end - if found - owner_faces_pindex[owner_face_lid] = pindex - pindexfound = true - break - end - end - @assert pindexfound "Valid pindex not found" - end - owner_faces_pindex, owner_faces_lids - end - - function generate_constraints(dmodel, - V, - reffe, - non_conforming_glue, - ref_constraints, - face_subface_ldof_to_cell_ldof) - num_regular_faces, - num_hanging_faces, - gridap_cell_faces, - hanging_faces_glue = non_conforming_glue - - sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = map_parts(gridap_cell_faces[1], - num_regular_faces[1], - num_hanging_faces[1], - hanging_faces_glue[1], - gridap_cell_faces[2], - num_regular_faces[2], - num_hanging_faces[2], - hanging_faces_glue[2], - dmodel.dmodel.models, V.spaces) do gridap_cell_vertices, - num_regular_vertices, num_hanging_vertices, - hanging_vertices_owner_cell_and_lface, - gridap_cell_faces, - num_regular_faces, num_hanging_faces, - hanging_faces_owner_cell_and_lface, - model, V - - # Locate for each hanging vertex a cell to which it belongs - # and local position within that cell - hanging_vertices_to_cell, - hanging_vertices_to_lvertex = _generate_hanging_faces_to_cell_and_lface(num_regular_vertices, - num_hanging_vertices, - gridap_cell_vertices) - - - # Locate for each hanging facet a cell to which it belongs - # and local position within that cell - hanging_faces_to_cell, - hanging_faces_to_lface = _generate_hanging_faces_to_cell_and_lface(num_regular_faces, - num_hanging_faces, - gridap_cell_faces) - - basis, reffe_args,reffe_kwargs = reffe - cell_reffe = ReferenceFE(QUAD,basis,reffe_args...;reffe_kwargs...) - reffe_cell = cell_reffe - - cell_dof_ids = get_cell_dof_ids(V) - face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(reffe_cell) - face_dofs = Gridap.ReferenceFEs.get_face_dofs(reffe_cell) - - hanging_vertices_owner_face_dofs = _generate_hanging_faces_owner_face_dofs(num_hanging_vertices, - face_dofs, - hanging_vertices_owner_cell_and_lface, - cell_dof_ids) - - hanging_faces_owner_face_dofs = _generate_hanging_faces_owner_face_dofs(num_hanging_faces, - face_dofs, - hanging_faces_owner_cell_and_lface, - cell_dof_ids) - - sDOF_to_dof = Int[] - sDOF_to_dofs = Vector{Int}[] - sDOF_to_coeffs = Vector{Float64}[] - - Dc = 2 - facet_polytope = Dc==2 ? SEGMENT : QUAD - basis, reffe_args,reffe_kwargs = reffe - face_reffe = ReferenceFE(facet_polytope,basis,reffe_args...;reffe_kwargs...) - pindex_to_cfvertex_to_fvertex = Gridap.ReferenceFEs.get_vertex_permutations(facet_polytope) - - lface_to_cvertices = Gridap.ReferenceFEs.get_faces(QUAD,Dc-1,0) - owner_faces_pindex, owner_faces_lids=_compute_owner_faces_pindex_and_lids(Dc, - num_hanging_faces, - hanging_faces_owner_cell_and_lface, - hanging_faces_to_cell, - hanging_faces_to_lface, - gridap_cell_vertices, - gridap_cell_faces, - lface_to_cvertices, - pindex_to_cfvertex_to_fvertex) - - nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope,[reffe_args[2]]) - node_permutations=Gridap.ReferenceFEs._compute_node_permutations(facet_polytope,nodes) - - hanging_lvertex_within_first_subface = 2^(Dc-1) - subface_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(face_reffe) - vertex_subface_own_dofs = subface_own_dofs[hanging_lvertex_within_first_subface] - face_dofs = Gridap.ReferenceFEs.get_face_dofs(cell_reffe) - - _generate_constraints!(0, - Dc, - gridap_cell_faces, - num_hanging_vertices, - hanging_vertices_to_cell, - hanging_vertices_to_lvertex, - hanging_vertices_owner_face_dofs, - hanging_vertices_owner_cell_and_lface, - face_subface_ldof_to_cell_ldof, - face_dofs, - face_own_dofs, - vertex_subface_own_dofs, - cell_dof_ids, - node_permutations, - owner_faces_pindex, - owner_faces_lids, - sDOF_to_dof, - sDOF_to_dofs, - sDOF_to_coeffs) - - face_subface_own_dofs = subface_own_dofs[end] - _generate_constraints!(1, - Dc, - gridap_cell_faces, - num_hanging_faces, - hanging_faces_to_cell, - hanging_faces_to_lface, - hanging_faces_owner_face_dofs, - hanging_faces_owner_cell_and_lface, - face_subface_ldof_to_cell_ldof, - face_dofs, - face_own_dofs, - face_subface_own_dofs, - cell_dof_ids, - node_permutations, - owner_faces_pindex, - owner_faces_lids, - sDOF_to_dof, - sDOF_to_dofs, - sDOF_to_coeffs) - sDOF_to_dof, Gridap.Arrays.Table(sDOF_to_dofs), Gridap.Arrays.Table(sDOF_to_coeffs) - end - end - - num_cells(dmodel.dmodel.models.part) - - sDOF_to_dof, sDOF_to_dofs,sDOF_to_coeffs= - generate_constraints(dmodel, V, reffe, - non_conforming_glue, ref_constraints, face_subface_ldof_to_cell_ldof) - println(sDOF_to_dof) - println(sDOF_to_dofs) - println(sDOF_to_coeffs) - - # Define manufactured functions - u(x) = x[1]+x[2]^order - f(x) = -Δ(u)(x) - - map_parts(dmodel.dmodel.models,V.spaces,U.spaces,sDOF_to_dof,sDOF_to_dofs,sDOF_to_coeffs) do model,V,U,sDOF_to_dof,sDOF_to_dofs,sDOF_to_coeffs - - println(get_cell_dof_ids(V)) - fl=get_face_labeling(model) - t=GridapDistributed.get_grid_topology(model) - println(fl.d_to_dface_to_entity[1]) - println(fl.d_to_dface_to_entity[2]) - println(Gridap.Geometry.get_faces(t,2,0)) - - Vc = FESpaceWithLinearConstraints( - sDOF_to_dof, - sDOF_to_dofs, - sDOF_to_coeffs, - V) - Uc = TrialFESpace(Vc,u) - - # Define integration mesh and quadrature - degree = 2*order+1 - Ω = Triangulation(model) - dΩ = Measure(Ω,degree) - - a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ - b(v) = ∫(v*f)*dΩ - - # op = AffineFEOperator(a,b,U,V0) - op = AffineFEOperator(a,b,Uc,Vc) - uh = solve(op) - - # Define exact solution and error - e = u - uh - - # Compute errors - el2 = sqrt(sum( ∫( e*e )*dΩ )) - eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩ )) - - tol=1e-8 - @assert el2 < tol - @assert eh1 < tol - end - -end - -test(1,1) -test(1,2) -test(1,3) - -test(2,1) -test(2,2) -test(2,3) - -test(3,1) -test(3,2) -test(3,3) - -test(4,1) -test(4,2) -test(4,3) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index df3fad3..896a03e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,6 +38,9 @@ function run_tests(testdir) "OctreeDistributedDiscreteModelsNoEnvTests.jl"] np = 6 extra_args = "" + elseif f in ["NonConformingOctreeDistributedDiscreteModelsTests.jl"] + np = 1 + extra_args = "" else np = nprocs extra_args = ""