Add support for adaptive mesh refinement with hanging-nodes to Gridap.jl (serial runs)

amartinhuertas commented 2 years ago

Potential tasks (to be filled as we go):

CONCERNS

[x] A remainder/concern to take into account while developing the code. Working with broken meshes like the ones in Step 0 below may lead to inconsistencies/violations of invariants inherent to Gridap's Geometry module.
- One that immediately comes into my mind is that hanging faces are solely in one cell. IF the code uses this as a criterion to determine faces sitting on the boundary, we may end having hanging faces in the boundary of the domain. This is obviously not possible/inconsistent.
- get_cell_permutations(...) may not provide the correct value for inter-octree owner (coarser) faces for non-oriented meshes.
- do you anticipate any other inconsistency that may arise along these lines?
[x] To think about how to implement facet integration on the boundary of cells with different levels of refinement and Neumann boundary faces. Does the approach based on broken meshes + glue versus meshes aware of non-conformity precludes us from implementing this? THINK NOW. [long-term] dev

TO-DO DONE

[x] @principejavier ... could you please help me out and extend the current P4est_wrapper.jl script so that we can start familiarizing ourselves with the logic behind p4est_lnodes_t in 3D? I think it might be a good idea to eairly stress/think deeply about the approach that we pursue to handle non-conforming meshes also in 3D, even if we defer the devs for 3D case for a later stage. Let me know what you think. [UPDATE by @amartinhuertas]: I added the P4est_wrapper.jl script that generates and shows the internals of p8est_lnodes_t. See: https://github.com/gridap/P4est_wrapper.jl/blob/adding_lnodes_example/test/example_lnodes_3d.jl and https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1399232101
[x] Immediate lines of action (2D code): (a) test with finer meshes/non-oriented meshes. (b) support for order>1 elements (currently only works for order=1 elements) (c) clean/abstract software. 26e34e130595f69f40a93cd9520be29ef811ac3f 0a290f3371e79f3202fc6a676d9368ea9cd92781 c5fe5a43c8658fade53bfe74fd8a027cdc6c8cbb 6d5fc79bee4a2718837220aa16a186411bda107e a099a54420671c8ac6a8cac3f6450b106cb09762 5debb9d91ff7a0a4b29fe5babe96a147522021cb e75910bca5d31f576c6e87779df9d51004236add
[x] To extend the program in https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384635611 for order>1 FEs. [short-term]
[x] To think how to extend to the nonconforming case the algorithms that generate the geometric entity IDs of vertices, edges, faces on the boundary of the adapted mesh from the ones of the coarsest mesh. See here for the case of conforming meshes. [short-term]
[x] To transform this TMP code such that it returns a fully consistent OctreeDistributedDiscreteModel. To this end, I expect we will have to to adapt/rewrite this function to the non-conforming case. In any case, it should be possible to re-use some of the building blocks, such as, generate_node_coordinates, but not all. [short-term]
[x] To extend this TMP code to 3D. [mid-term]
[x] [ONLY IN 3d] To take into account https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1386402744 while developing the code. Related to previous task. [mid-term] 8c6c278b80159a21e1114c6be699d9c4731b8ccc
[x] On a separate note, I do believe that the function generate_cell_vertices_and_faces should return dimension-indexed arrays (i.e., 1 for corners, 2 for edges, etc.), instead of separated arrays for vertices, and faces, etc. This would avoid code replication in the generation of symbolic constraint dependencies. cec3510e6349ef97b964e6e46b5b8fa036b67fc1

cc @santiagobadia @fverdugo

I create this issue as requested by @principejavier. The issue develops a bit further the project abstract available at

https://github.com/gridap/GSoC/blob/main/2022/ideas-list.md#add-support-for-adaptive-mesh-refinement-with-hanging-nodes-to-gridapjl

with additional remarks/code/details.

[Step 0]. A mano (no hace falta p4est), construir una malla con dos celdas, una refinada, y la otra no (un hanging node), construir un FE space non-conforming encima (el hanging node recibira un id de dof positivo), y luego construir el FESpaceWithLinearConstraints encima de este ultimo, y probar que todo ok para resolver un Problema de Poisson con solucion manufacturada con la maquinaria que tenemos ahora mismo en Gridap. Para construir la malla con las cinco celdas podemos utilizar, e.g., UnstructuredDiscreteModel. See, e.g., here for an example of manually building an UnstructuredDiscreteModel.
[Step 1] Familiarize yourself with the so-called p4est_lnodes_t struct. The subroutine that builds this struct in p4est provides degree=-2 option. I am positive this topological data structure is going to provide all the information that we need in order to glue vertices, edges, faces, etc. on the boundary of the cells in a parallel distributed context, find hanging vertices, edges, and faces, their coarser cells around, etc. In other words, to be able to build the linear multipoint constraints which are required to fed FESpaceWithLinearConstraints. Note that in FEMPAR we used p4est_mesh_t instead of p4est_lnodes_t struct. Note also that in GridapP4est.jl I am already using p4est_lnodes_t to glue together vertices, edges, faces on the boundary of the cells in a distributed context. (See here). In order to help you out understanding this struct, I have written this P4est_wrapper.jl script that decodes the contents of the struct. Read the comments, and execute it to see the output on the terminal!!!
[Step 2] Combine the DoF numbering from UnconstrainedFESpace and the info extracted out of [Step 1] to automatically build the symbolic CSR part of the linear multipoint constraints. For the coefficients, you can use a similar strategy that the one in FEMPAR, i.e., basis functions evaluated in reference space at the vertices of a grid which represents the uniform refinement rule 1:2^d of the reference n-cube. To this end this piece of code can be useful as a starting point.

principejavier commented 2 years ago

Thanks @amartinhuertas! I finally had time to start with this. So far I only completed Step0 and I got a fair first understanding of Gridap.Geometry. Everything seems to work, see the pictures without and with constrains:

non_conforming

constrained

amartinhuertas commented 1 year ago

Hi @principejavier ! Which is the status of this development?

principejavier commented 1 year ago

Hi @amartinhuertas I was at step 1 when shifted focus to other issues. I expect to have some time for this at the end of the month.

amartinhuertas commented 1 year ago

Ok. Thanks for the update. I may have some time the following two weeks to contribute to this dev. I am just asking so that we do not overlap work. My goal would be ideally to come up with the linear constraints built for the serial case + 2D, and be able to solve a 2D Poisson problem with a non-conforming mesh top-bottom.

amartinhuertas commented 1 year ago

cc @JordiManyer ... highly related to his current work.

JordiManyer commented 1 year ago

I have time to put into this PR, although I don't really have experience with the low-level P4est api. But I'm glad to take any work off your hands if you need me to.

principejavier commented 1 year ago

Here is the code for Step 0

module poisson

using FillArrays
using Gridap
using Gridap.Arrays
using Gridap.FESpaces
import Gridap: ∇

# Define manufactured functions
u(x) = x[1] # + x[2]
∇u(x) = VectorValue(1,1)
∇(::typeof(u)) = ∇u
f(x) = 0.0

function setup_model()
#
#  7---8---9------11
#  | 3 | 4 |       |
#  4---5---6   5   |
#  | 1 | 2 |       |
#  1---2---3------10
#
    ptr  = [ 1, 5, 9, 13, 17, 21 ]
    data = [ 1,2,4,5, 2,3,5,6, 4,5,7,8, 5,6,8,9, 3,10,9,11 ]
    cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
    node_coordinates = Vector{Point{2,Float64}}(undef,11)
    for j in 1:3
        for i in 1:3
            node_coordinates[3*(j-1)+i]=Point{2,Float64}((i-1)/2.0,(j-1)/2.0)
        end
    end
    node_coordinates[10]=Point{2,Float64}(2.0,0.0)
    node_coordinates[11]=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())
    Gridap.Geometry.UnstructuredDiscreteModel(grid)
end

function run()

# Construct the discrete model
model=setup_model()
labels = get_face_labeling(model)
labels.d_to_dface_to_entity[1]=[1,0,0,1,0,0,1,0,0,1,1]
labels.d_to_dface_to_entity[2]=[0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,1] # this is not needed...
add_tag!(labels,"dirichlet",[1])

# Construct the FEspace
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="dirichlet")
U = TrialFESpace(V0,u)

sDOF_to_dof = [4]
sDOF_to_dofs = Table([[2,6]])
sDOF_to_coeffs = Table([[0.5,0.5]])
Vc = FESpaceWithLinearConstraints(
  sDOF_to_dof,
  sDOF_to_dofs,
  sDOF_to_coeffs,
  V0)
Uc = TrialFESpace(Vc,u)

# Define integration mesh and quadrature
degree = 2
Ω = 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

# Write the numerical solution, the manufactured solution, and the error
# in a vtu file
#writevtk(trian,"results",nref=1,cellfields=["uh"=>uh,"u"=>u,"e"=>e])
writevtk(Ω,"results",cellfields=["uh"=>uh,"u"=>u,"e"=>e])

end

end

amartinhuertas commented 1 year ago

For the records, and to help understanding the p4est_lnodes_t data structure, let me post below, assuming that we run the sample snippet program (in its current status: https://github.com/gridap/P4est_wrapper.jl/blob/a0b27a292f4684ccd485661b7ac05d9bccac2e32/test/example_lnodes.jl) in a single MPI task:

A picture of the element_nodes member variable, along with the values of other relevant member variables.
The output that the program generates on screen.

cell=1 faces=Int32[0, 1, 2, 3]
lface=0 gface=0
lface=1 gface=1
lface=2 gface=2
lface=3 gface=3
cell=1 vertices=Int32[4, 5, 6, 7]
cell=1 does not have hanging faces
cell=2 faces=Int32[1, 8, 9, 10]
lface=1 gface=1
lface=8 gface=8
lface=9 gface=9
lface=10 gface=10
cell=2 vertices=Int32[5, 11, 7, 16]
cell=2 has hanging faces
cell=2 lface=1 is NOT hanging
cell=2 lface=2 is hanging from half=0
cell=2 lface=3 is NOT hanging
cell=2 lface=4 is NOT hanging
cell=3 faces=Int32[12, 13, 3, 14]
lface=12 gface=12
lface=13 gface=13
lface=3 gface=3
lface=14 gface=14
cell=3 vertices=Int32[6, 7, 15, 16]
cell=3 has hanging faces
cell=3 lface=1 is NOT hanging
cell=3 lface=2 is NOT hanging
cell=3 lface=3 is NOT hanging
cell=3 lface=4 is hanging from half=0
cell=4 faces=Int32[13, 8, 10, 14]
lface=13 gface=13
lface=8 gface=8
lface=10 gface=10
lface=14 gface=14
cell=4 vertices=Int32[7, 11, 15, 16]
cell=4 has hanging faces
cell=4 lface=1 is NOT hanging
cell=4 lface=2 is hanging from half=1
cell=4 lface=3 is NOT hanging
cell=4 lface=4 is hanging from half=1
cell=5 faces=Int32[8, 17, 18, 19]
lface=8 gface=8
lface=17 gface=17
lface=18 gface=18
lface=19 gface=19
cell=5 vertices=Int32[11, 20, 16, 21]
cell=5 does not have hanging faces
cell=6 faces=Int32[22, 23, 14, 24]
lface=22 gface=22
lface=23 gface=23
lface=14 gface=14
lface=24 gface=24
cell=6 vertices=Int32[15, 16, 25, 26]
cell=6 does not have hanging faces
cell=7 faces=Int32[23, 27, 19, 28]
lface=23 gface=23
lface=27 gface=27
lface=19 gface=19
lface=28 gface=28
cell=7 vertices=Int32[16, 21, 26, 29]
cell=7 does not have hanging faces

amartinhuertas commented 1 year ago

The same as in https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384819749 but with an updated version of the test driver that just prints additional output on screen (https://github.com/gridap/P4est_wrapper.jl/blob/568fc013a991bde9f619886474cc35e5ae2f26ea/test/example_lnodes.jl) run on two MPI tasks: EDIT: (1) the upper left corner in the 4th cell in P0 should be labeled as 15 not 16; (2) the upper right corner of the 4th cell in P0 is labeled as 16.

[1,0]<stdout>:num_local_nodes=17
[1,0]<stdout>:owned_count=17
[1,0]<stdout>:num_nonlocal_nodes=0
[1,0]<stdout>:nonlocal_nodes=Int64[]
[1,0]<stdout>:global_owned_count=Int32[17, 13]
[1,0]<stdout>:global_offset=0
[1,0]<stdout>:num_nonlocal_nodes=0
[1,0]<stdout>:nonlocal_nodes=Int64[]
[1,0]<stdout>:cell=1 faces=Int32[0, 1, 2, 3]
[1,0]<stdout>:lface=0 gface=0
[1,0]<stdout>:lface=1 gface=1
[1,0]<stdout>:lface=2 gface=2
[1,0]<stdout>:lface=3 gface=3
[1,0]<stdout>:cell=1 vertices=Int32[4, 5, 6, 7]
[1,0]<stdout>:cell=1 does not have hanging faces
[1,0]<stdout>:cell=2 faces=Int32[1, 8, 9, 10]
[1,0]<stdout>:lface=1 gface=1
[1,0]<stdout>:lface=8 gface=8
[1,0]<stdout>:lface=9 gface=9
[1,0]<stdout>:lface=10 gface=10
[1,0]<stdout>:cell=2 vertices=Int32[5, 11, 7, 16]
[1,0]<stdout>:cell=2 has hanging faces
[1,0]<stdout>:cell=2 lface=1 is NOT hanging
[1,0]<stdout>:cell=2 lface=2 is hanging from half=0
[1,0]<stdout>:cell=2 lface=3 is NOT hanging
[1,0]<stdout>:cell=2 lface=4 is NOT hanging
[1,0]<stdout>:cell=3 faces=Int32[12, 13, 3, 14]
[1,0]<stdout>:lface=12 gface=12
[1,0]<stdout>:lface=13 gface=13
[1,0]<stdout>:lface=3 gface=3
[1,0]<stdout>:lface=14 gface=14
[1,0]<stdout>:cell=3 vertices=Int32[6, 7, 15, 16]
[1,0]<stdout>:cell=3 has hanging faces
[1,0]<stdout>:cell=3 lface=1 is NOT hanging
[1,0]<stdout>:cell=3 lface=2 is NOT hanging
[1,0]<stdout>:cell=3 lface=3 is NOT hanging
[1,0]<stdout>:cell=3 lface=4 is hanging from half=0
[1,0]<stdout>:cell=4 faces=Int32[13, 8, 10, 14]
[1,0]<stdout>:lface=13 gface=13
[1,0]<stdout>:lface=8 gface=8
[1,0]<stdout>:lface=10 gface=10
[1,0]<stdout>:lface=14 gface=14
[1,0]<stdout>:cell=4 vertices=Int32[7, 11, 15, 16]
[1,0]<stdout>:cell=4 has hanging faces
[1,0]<stdout>:cell=4 lface=1 is NOT hanging
[1,0]<stdout>:cell=4 lface=2 is hanging from half=1
[1,0]<stdout>:cell=4 lface=3 is NOT hanging
[1,0]<stdout>:cell=4 lface=4 is hanging from half=1

[1,1]<stdout>:num_local_nodes=18
[1,1]<stdout>:owned_count=13
[1,1]<stdout>:num_nonlocal_nodes=5
[1,1]<stdout>:nonlocal_nodes=[8, 11, 14, 15, 16]
[1,1]<stdout>:global_owned_count=Int32[17, 13]
[1,1]<stdout>:global_offset=17
[1,1]<stdout>:num_nonlocal_nodes=5
[1,1]<stdout>:nonlocal_nodes=[8, 11, 14, 15, 16]
[1,1]<stdout>:cell=1 faces=Int32[13, 0, 1, 2]
[1,1]<stdout>:lface=13 gface=30
[1,1]<stdout>:lface=0 gface=17
[1,1]<stdout>:lface=1 gface=18
[1,1]<stdout>:lface=2 gface=19
[1,1]<stdout>:cell=1 vertices=Int32[14, 3, 17, 4]
[1,1]<stdout>:cell=1 does not have hanging faces
[1,1]<stdout>:cell=2 faces=Int32[5, 6, 15, 7]
[1,1]<stdout>:lface=5 gface=22
[1,1]<stdout>:lface=6 gface=23
[1,1]<stdout>:lface=15 gface=32
[1,1]<stdout>:lface=7 gface=24
[1,1]<stdout>:cell=2 vertices=Int32[16, 17, 8, 9]
[1,1]<stdout>:cell=2 does not have hanging faces
[1,1]<stdout>:cell=3 faces=Int32[6, 10, 2, 11]
[1,1]<stdout>:lface=6 gface=23
[1,1]<stdout>:lface=10 gface=27
[1,1]<stdout>:lface=2 gface=19
[1,1]<stdout>:lface=11 gface=28
[1,1]<stdout>:cell=3 vertices=Int32[17, 4, 9, 12]
[1,1]<stdout>:cell=3 does not have hanging faces

principejavier commented 1 year ago

Note that the p4est does not generate an ID for hanging nodes but it fills the cell with the constraining one (that is not already in the cell), e.g upper right corner of cell 3 in https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384819749 is printed as 16.

Therefore replacing the function setup_model in the code of Step 0 above will not be immediate, we need to generate a new numbering (adding IDs to hanging nodes) to create the space V0. One option would be to just using a counter starting at num_local_nodes. In the parallel case it would require a reduction to generate global IDs.

amartinhuertas commented 1 year ago

Note that the p4est does not generate an ID for hanging nodes but it fills the cell with the constraining one (that is not already in the cell), e.g upper right corner of cell 3 in https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384819749 is printed as 16. Therefore replacing the function setup_model in the code of Step 0 above will not be immediate, we need to generate a new numbering (adding IDs to hanging nodes) to create the space V0.

Yes, I also noted that. Unfortunately, we cannot run the current algorithm for conforming meshes verbatim to generate the model produced by setup_model() in step 0 (e.g., for faces, this would imply gluing together the hanging faces and their owner face, which i think does not make much sense). As you say, it seems we need an ad-hoc numbering of vertices, edges, and faces algorithm for non-conforming meshes.

In 2D, the logic behind the definition of p4est_lnodes_t seems to be:

For hanging faces, it provides the "owner" face.
For hanging vertices, it provides the "opposite" (regular assuming 2:1 balance) vertex at the boundary of the owner face. As a result, p4est_lnodes_t does not "glue" together hanging vertices with the same global ID in all cells where the vertex belongs.

[UPDATE]: the logic behind p4est_lnodes is to provide, for a given cell, all the dofs ids which are ultimately touched when assembling the cell + the constraints on its hanging dofs.

amartinhuertas commented 1 year ago

Just documenting in this post (to not forget) an scenario we will have to deal with correctly when working with 3D non-conforming meshes:

amartinhuertas commented 1 year ago

For the records, below you can find the output of https://github.com/gridap/P4est_wrapper.jl/blob/51530b6d8a88b2a6b14228eeffcd910b652a19a6/test/example_lnodes_3d.jl (i.e., p8est_lnodes_t for degree=-3) along with an illustration that covers all possible cases:

cell=1 faces=Int32[0, 1, 2, 3, 4, 5]
lface=0 gface=0
lface=1 gface=1
lface=2 gface=2
lface=3 gface=3
lface=4 gface=4
lface=5 gface=5
cell=1 edges=Int32[6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
cell=1 corners=Int32[18, 19, 20, 21, 22, 23, 24, 25]
cell=1 does not have hanging faces
cell=2 faces=Int32[1, 26, 27, 28, 29, 30]
lface=1 gface=1
lface=26 gface=26
lface=27 gface=27
lface=28 gface=28
lface=29 gface=29
lface=30 gface=30
cell=2 edges=Int32[31, 32, 33, 34, 11, 35, 13, 36, 15, 37, 17, 38]
cell=2 corners=Int32[19, 39, 21, 40, 23, 41, 25, 42]
cell=2 has hanging faces
hanging_face=Int32[0, -1, -1, -1, -1, -1]
hanging_edge=Int32[-1, -1, -1, -1, 0, -1, 4, -1, 0, -1, 4, -1]
cell=3 faces=Int32[26, 43, 44, 45, 46, 47]
lface=26 gface=26
lface=43 gface=43
lface=44 gface=44
lface=45 gface=45
lface=46 gface=46
lface=47 gface=47
cell=3 edges=Int32[48, 49, 50, 51, 35, 52, 36, 53, 37, 54, 38, 55]
cell=3 corners=Int32[39, 56, 40, 57, 41, 58, 42, 59]
cell=3 does not have hanging faces
cell=4 faces=Int32[1, 60, 28, 61, 62, 63]
lface=1 gface=1
lface=60 gface=60
lface=28 gface=28
lface=61 gface=61
lface=62 gface=62
lface=63 gface=63
cell=4 edges=Int32[32, 64, 34, 89, 11, 65, 13, 66, 15, 38, 17, 72]
cell=4 corners=Int32[19, 40, 21, 73, 23, 42, 25, 91]
cell=4 has hanging faces
hanging_face=Int32[1, -1, -1, 0, -1, -1]
hanging_edge=Int32[-1, 0, -1, 4, 3, -1, 4, -1, 4, -1, 0, 4]
cell=5 faces=Int32[60, 67, 45, 61, 68, 69]
lface=60 gface=60
lface=67 gface=67
lface=45 gface=45
lface=61 gface=61
lface=68 gface=68
lface=69 gface=69
cell=5 edges=Int32[49, 64, 51, 89, 65, 70, 66, 71, 38, 55, 17, 72]
cell=5 corners=Int32[40, 57, 21, 73, 42, 59, 25, 91]
cell=5 has hanging faces
hanging_face=Int32[-1, -1, -1, 1, -1, -1]
hanging_edge=Int32[-1, 1, -1, 4, -1, -1, -1, -1, -1, -1, 4, 2]
cell=6 faces=Int32[1, 74, 75, 76, 30, 77]
lface=1 gface=1
lface=74 gface=74
lface=75 gface=75
lface=76 gface=76
lface=30 gface=30
lface=77 gface=77
cell=6 edges=Int32[33, 34, 78, 89, 11, 36, 13, 84, 15, 79, 17, 80]
cell=6 corners=Int32[19, 41, 21, 42, 23, 87, 25, 91]
cell=6 has hanging faces
hanging_face=Int32[2, -1, -1, -1, -1, 0]
hanging_edge=Int32[-1, -1, 0, 4, 4, -1, 0, 4, 3, -1, 4, -1]
cell=7 faces=Int32[74, 81, 82, 83, 47, 77]
lface=74 gface=74
lface=81 gface=81
lface=82 gface=82
lface=83 gface=83
lface=47 gface=47
lface=77 gface=77
cell=7 edges=Int32[50, 51, 78, 89, 36, 53, 13, 84, 79, 85, 80, 86]
cell=7 corners=Int32[41, 58, 42, 59, 23, 87, 25, 91]
cell=7 has hanging faces
hanging_face=Int32[-1, -1, -1, -1, -1, 1]
hanging_edge=Int32[-1, -1, 1, 4, -1, -1, 4, 2, -1, -1, -1, -1]
cell=8 faces=Int32[1, 88, 76, 61, 63, 77]
lface=1 gface=1
lface=88 gface=88
lface=76 gface=76
lface=61 gface=61
lface=63 gface=63
lface=77 gface=77
cell=8 edges=Int32[34, 64, 78, 89, 11, 66, 13, 84, 15, 80, 17, 72]
cell=8 corners=Int32[19, 42, 21, 73, 23, 87, 25, 91]
cell=8 has hanging faces
hanging_face=Int32[3, -1, -1, 2, -1, 2]
hanging_edge=Int32[-1, 4, 4, 2, 4, -1, 1, 4, 4, -1, 1, 4]
cell=9 faces=Int32[88, 90, 83, 61, 69, 77]
lface=88 gface=88
lface=90 gface=90
lface=83 gface=83
lface=61 gface=61
lface=69 gface=69
lface=77 gface=77
cell=9 edges=Int32[51, 64, 78, 89, 66, 71, 13, 84, 80, 86, 17, 72]
cell=9 corners=Int32[42, 59, 21, 73, 23, 87, 25, 91]
cell=9 has hanging faces
hanging_face=Int32[-1, -1, -1, 3, -1, 3]
hanging_edge=Int32[-1, 4, 4, 3, -1, -1, 4, 3, -1, -1, 4, 3]
cell=10 faces=Int32[92, 93, 3, 94, 95, 96]
lface=92 gface=92
lface=93 gface=93
lface=3 gface=3
lface=94 gface=94
lface=95 gface=95
lface=96 gface=96
cell=10 edges=Int32[7, 97, 9, 98, 99, 100, 101, 102, 16, 17, 103, 104]
cell=10 corners=Int32[20, 21, 105, 106, 24, 25, 107, 108]
cell=10 does not have hanging faces
cell=11 faces=Int32[93, 109, 61, 110, 111, 112]
lface=93 gface=93
lface=109 gface=109
lface=61 gface=61
lface=110 gface=110
lface=111 gface=111
lface=112 gface=112
cell=11 edges=Int32[64, 113, 89, 114, 100, 115, 102, 116, 17, 72, 104, 117]
cell=11 corners=Int32[21, 73, 106, 118, 25, 91, 108, 119]
cell=11 does not have hanging faces
cell=12 faces=Int32[120, 121, 122, 123, 5, 124]
lface=120 gface=120
lface=121 gface=121
lface=122 gface=122
lface=123 gface=123
lface=5 gface=5
lface=124 gface=124
cell=12 edges=Int32[8, 9, 125, 126, 12, 13, 127, 128, 129, 130, 131, 132]
cell=12 corners=Int32[22, 23, 24, 25, 133, 134, 135, 136]
cell=12 does not have hanging faces
cell=13 faces=Int32[121, 137, 138, 139, 77, 140]
lface=121 gface=121
lface=137 gface=137
lface=138 gface=138
lface=139 gface=139
lface=77 gface=77
lface=140 gface=140
cell=13 edges=Int32[78, 89, 141, 142, 13, 84, 128, 143, 130, 144, 132, 145]
cell=13 corners=Int32[23, 87, 25, 91, 134, 146, 136, 147]
cell=13 does not have hanging faces
cell=14 faces=Int32[148, 149, 123, 150, 96, 151]
lface=148 gface=148
lface=149 gface=149
lface=123 gface=123
lface=150 gface=150
lface=96 gface=96
lface=151 gface=151
cell=14 edges=Int32[9, 98, 126, 152, 101, 102, 153, 154, 131, 132, 155, 156]
cell=14 corners=Int32[24, 25, 107, 108, 135, 136, 157, 158]
cell=14 does not have hanging faces
cell=15 faces=Int32[149, 159, 139, 160, 112, 161]
lface=149 gface=149
lface=159 gface=159
lface=139 gface=139
lface=160 gface=160
lface=112 gface=112
lface=161 gface=161
cell=15 edges=Int32[89, 114, 142, 162, 102, 116, 154, 163, 132, 145, 156, 164]
cell=15 corners=Int32[25, 91, 108, 119, 136, 147, 158, 165]
cell=15 does not have hanging faces

amartinhuertas commented 1 year ago

@principejavier ... please note:

The following commit on GridapP4est.jl, p4est-migration-exploring-non-conformity branch. https://github.com/gridap/GridapP4est.jl/commit/7d7fc38eca3549fe46466a992a7e8499a8499734
I have updated the issue description (https://github.com/gridap/GridapP4est.jl/issues/21#issue-1333404245) such that that I have added some additional tasks that I believe have to be completed after the work in 1.

I dont have much time the following weeks to continue with this dev (teaching starts). Can you please take the lead? We can have a meeting if you like, I can also help you out, but not too much time to write code.

principejavier commented 1 year ago

After https://github.com/gridap/GridapP4est.jl/commit/e21df63f15b18e65d95141941161929ba23bb457 the code is working in 2D producing a DistributedDiscreteModel as in the uniform case, including geometric IDs. I need to change the focus (to produce some results with GridapMHD.jl) but I will continue with this development as much as I can.

amartinhuertas commented 1 year ago

After https://github.com/gridap/GridapP4est.jl/commit/e21df63f15b18e65d95141941161929ba23bb457 the code is working in 2D producing a DistributedDiscreteModel as in the uniform case, including geometric IDs. I need to change the focus (to produce some results with GridapMHD.jl) but I will continue with this development as much as I can.

FYI ... I could work on this the past couple of days. I have now to leave it here for this week. I moved part of the code to src and implemented the generation of symbolic constraints dependencies (the coefficient values are still missing).

On a separate note, I do believe that the function generate_cell_vertices_and_faces should return dimension-indexed arrays (i.e., 1 for corners, 2 for edges, etc.), instead of separated arrays for vertices, and faces, etc. This would avoid code replication in the generation of symbolic constraint dependencies.

amartinhuertas commented 1 year ago

@principejavier ... FYI, this week, I could work a day on this. I just pushed into the branch. The current code (https://github.com/gridap/GridapP4est.jl/blob/61b08ac74e8e560a762532ef8a1d896a02ca3d70/test/non_conforming_octrees_wip.jl) already automates the generation of constraints and reproduces https://github.com/gridap/GridapP4est.jl/issues/21#issuecomment-1384635611.

Thanks to @JordiManyer that helped us via this Gridap PR to provide some missing functionality needed to build the constraints (https://github.com/gridap/Gridap.jl/pull/886).

Immediate lines of action:

test with finer meshes/non-oriented meshes.
support for order>1 elements (currently only works for order=1 elements)
clean/abstract software.

amartinhuertas commented 1 year ago

for the records ... I have already completed the following in 2D:

Immediate lines of action: test with finer meshes/non-oriented meshes. support for order>1 elements (currently only works for order=1 elements) clean/abstract software.

I am going to start next with the 3D case. To this end, I will need at some point an extension of

Thanks to @JordiManyer that helped us via this Gridap PR to provide some missing functionality needed to build the constraints (https://github.com/gridap/Gridap.jl/pull/886).

to the 3D case. No hurries. There is still a lot of work in the geometric/topological part before starting with the construction of the hanging node constraints.

amartinhuertas commented 1 year ago

I am going to start next with the 3D case.

For the records, I have now a preliminary implementation of 3D AMR with forest-of-octrees that works for scalar Lagrangian FEs and arbitrary order. See https://github.com/gridap/GridapP4est.jl/blob/a65eb67e165e1489cef922dafa797bd29d196277/test/NonConformingOctreeDistributedDiscreteModelsTests.jl! Thanks to @JordiManyer that helped us to extending the functionality required to compute the constraints in 3D!

amartinhuertas commented 1 year ago

Closing this issue. At present, the code works in 2D/3D and in parallel with non-conforming meshes. all tasks except facet integration have been completed. See https://github.com/gridap/GridapP4est.jl/issues/35 for a more up-to-date list of pending tasks.

gridap / GridapP4est.jl

Add support for adaptive mesh refinement with hanging-nodes to Gridap.jl (serial runs) #21