Add support for adaptive mesh refinement with hanging-nodes to Gridap.jl (serial runs) #21
Comments
|
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:
|
|
Hi @principejavier ! Which is the status of this development? |
|
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. |
|
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. |
|
cc @JordiManyer ... highly related to his current work. |
|
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. |
|
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 |
|
For the records, and to help understanding the
|
|
The same as in #21 (comment) 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.
|
|
Note that the Therefore replacing the function |
Yes, I also noted that. Unfortunately, we cannot run the current algorithm for conforming meshes verbatim to generate the model produced by In 2D, the logic behind the definition of
[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. |
|
Just documenting in this post (to not forget) an scenario we will have to deal with correctly when working with 3D non-conforming meshes:
|
|
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.,
|
|
@principejavier ... please note:
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. |
|
After e21df63 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 On a separate note, I do believe that the function |
|
@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 #21 (comment). Thanks to @JordiManyer that helped us via this Gridap PR to provide some missing functionality needed to build the constraints (gridap/Gridap.jl#886). Immediate lines of action:
|
|
for the records ... I have already completed the following in 2D:
I am going to start next with the 3D case. To this end, I will need at some point an extension of
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. |
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! |
|
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 #35 for a more up-to-date list of pending tasks. |
Potential tasks (to be filled as we go):
CONCERNS
get_cell_permutations(...)may not provide the correct value for inter-octree owner (coarser) faces for non-oriented meshes.TO-DO
DONE
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]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
FESpaceWithLinearConstraintsencima 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 anUnstructuredDiscreteModel.[Step 1] Familiarize yourself with the so-called
p4est_lnodes_tstruct. The subroutine that builds this struct in p4est providesdegree=-2option. 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 fedFESpaceWithLinearConstraints. Note that in FEMPAR we usedp4est_mesh_tinstead ofp4est_lnodes_t struct. Note also that inGridapP4est.jlI am already usingp4est_lnodes_tto 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 thisP4est_wrapper.jlscript 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
UnconstrainedFESpaceand 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.The text was updated successfully, but these errors were encountered: