Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

p4est datastructures #780

Open
wants to merge 159 commits into
base: master
Choose a base branch
from
Open

p4est datastructures #780

wants to merge 159 commits into from

Conversation

koehlerson
Copy link
Member

@koehlerson koehlerson commented Jul 31, 2023

to be discussed at FerriteCon 23, register now: https://ferrite-fem.github.io/FerriteCon/

TODO

  • p4est types
    • OctantBWG morton index encoded octant in 2 and 3D
    • OctreeBWG linear tree with homogeneous OctantBWG leaves
    • ForestBWG forest of homogeneous OctreeBWG
  • Octant Operations
  • Octree operations
    • balance intertree
    • balance intratree
    • coarsen
    • refine
  • materializing a Ferrite.Grid for unstructured meshes
  • hanging nodes
  • error estimator interface
    • ZZ type error estimator slow version (probably working in docs/srcs/literate-tutorials/adaptivity.jl
  • Docs
    • Docs to the p4est types and how they interplay
    • Docs to the yet to be defined interface
    • quick and dirty implementation of elasticity docs/src/literate-tutorials/adaptivity.jl with (almost?maybe?) working but slow ZZ error estimator
    • manufactured solution example for heat equation at docs/src/literate-tutorials/heat-adaptivity.jl
    • adjust quick and dirty docs/src/literate-tutorials/adaptivity.jl to a proper tutorial
  • More debug statements to track down issues
  • multiple neighbors attached to edge and vertex
  • rotation test in 3D

Followup PRs

@termi-official
Copy link
Member

1. Why doesn't `Ferrite.refine!` include `Ferrite.balanceforest!(grid)`? (Of course, internally it makes sense to have these separated, but I would have assumed that the API function `refine!` would do both.

There is a performance trap hiding here. Users might refine single elements one by one and if this triggers rebalancing the refinement can become quickly very expensive. Hence the users is responsible to call the balance function once he is done.

4. [..] I guess at least some old to new cells mapping is required to do this efficiently?

Yes, that is missing and necessary for any time-dependent problem. We add more features later (see follow-up PRs). This PR concentrates on the core implementation with documentation.

Error Estimation
Note that this is up the the user to implement? (If it is, I just assume here).

Here we run into the issue that Ferrite is "too low-level". The majority of error estimators need direct access to an evaluation of the energy or flux. So for now the only option is that the packages on top of Ferrite provide these.

Comment on lines +107 to +119
# Maps from entity to dofs
# `vertexdict` keeps track of the visited vertices. The first dof added to vertex v is
# stored in vertexdict[v].
vertexdicts::Vector{Vector{Int}}
# `edgedict` keeps track of the visited edges, this will only be used for a 3D problem.
# An edge is uniquely determined by two global vertices, with global direction going
# from low to high vertex number.
edgedicts::Vector{Dict{Tuple{Int,Int}, Int}}
# `facedict` keeps track of the visited faces. We only need to store the first dof we
# add to the face since currently more dofs per face isn't supported. In
# 2D a face (i.e. a line) is uniquely determined by 2 vertices, and in 3D a face (i.e. a
# surface) is uniquely determined by 3 vertices.
facedicts::Vector{Dict{NTuple{3,Int}, Int}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reminder to me: Try to remove this again.

@koehlerson
Copy link
Member Author

koehlerson commented Jun 11, 2024

I assumed it doesn't follow the AbstractGrid interface, otherwise you could have used that directly in the DofHandler I assumed. If it is, would be nice to document what it does include and what it doesn't from the interface.

It fulfills the AbstractGrid interface, but the problem is that the AbstractCell interface is too strict for adaptive/distributed grids (cells respectively). Thus, the DofHandler fails in distributing dofs, because of calls like

for (vi, vertex) in pairs(vertices(cell))
which assume that a cell knows about its global vertex ids, which isn't the case for the OctantBWG cells. Even the OctreeBWG won't know about its global ids, instead given o::OctantBWG in some tree::OctreeBWG part of forest::ForestBWG you can determine o's global vertex ids. I think it's the same with distributed but just one layer removed and instead you can imagine to track the shared dofs somehow

Yeah for the algorithm or interface in general I need to discuss with @termi-official a bit.

@koehlerson
Copy link
Member Author

@fredrikekre what do you think in terms of minimal abstraction for gauss point data which eases state variable transfer and error estimation vs giving it up to the user to do this?

@KnutAM
Copy link
Member

KnutAM commented Jun 11, 2024

There is a performance trap hiding here. Users might refine single elements one by one and if this triggers rebalancing the refinement can become quickly very expensive. Hence the users is responsible to call the balance function once he is done.

OK, if the intention is that you should be able to call refine! for one cell at a time, and then do balance, then let's leave that to later than. I thought the interface refine!(grid, cell_numbers::Vector) was intended as a single step refine all cells.

I guess we can at some point introduce a new function remesh that does refine + balance for P4est AMR and only refine for e.g. "triangle-splitting" AMR. But to not take the namespace (going into 1.0), wouldn't it make sense to call this refine_tree or smth since the split is specific to P4est?

@KnutAM
Copy link
Member

KnutAM commented Jun 11, 2024

Side-note that I discovered when trying to understand the above:
I think users could be surprised that refine!(grid, cellids) is not equivalent to

for cellid in cellids
    refine!(grid, cellid)
end

Perhaps the refine!(grid, ::Int) should be an internal function to encourage users to always use a collection of cells to refine?

@fredrikekre
Copy link
Member

Because the first call generates new cell numbers or?

@koehlerson
Copy link
Member Author

koehlerson commented Jun 11, 2024

Because the first call generates new cell numbers or?

Yes, as soon as the integer dispatch is called 4 (2D) or 8 (3D) new cells are inserted therefore shifting the cellids of the other cells

@koehlerson
Copy link
Member Author

I guess we can at some point introduce a new function remesh that does refine + balance for P4est AMR and only refine for e.g. "triangle-splitting" AMR. But to not take the namespace (going into 1.0), wouldn't it make sense to call this refine_tree or smth since the split is specific to P4est?

I'd vote for structs that act as a flag whether or not after refinement should be balanced. That should be okay for all approaches. Balancing is not something super tightly connected to p4est. It's useful in cases where you don't want nested non-conformity constraints. I would rather go for refine since more or less every adaptive approach can be seen as a tree approach (also triangle splitting)

@KnutAM
Copy link
Member

KnutAM commented Jun 12, 2024

I'd vote for structs that act as a flag whether or not after refinement should be balanced.

That sounds good to me!

@fredrikekre
Copy link
Member

@fredrikekre what do you think in terms of minimal abstraction for gauss point data which eases state variable transfer and error estimation vs giving it up to the user to do this?

No idea, but I guess there have to be some way to transfer things between the meshes, also the solution, right? A good start would probably be to have a cellmapping and then you can interpolate as you wish for quadrature data perhaps?

@koehlerson
Copy link
Member Author

koehlerson commented Jun 13, 2024

@fredrikekre what do you think in terms of minimal abstraction for gauss point data which eases state variable transfer and error estimation vs giving it up to the user to do this?

No idea, but I guess there have to be some way to transfer things between the meshes, also the solution, right? A good start would probably be to have a cellmapping and then you can interpolate as you wish for quadrature data perhaps?

Yes the solution transfer is not so problematic since we have the knowledge about the fields, their underlying interpolation and the previous grid they lived on. However, what about the Gauss point data, e.g. internal variables. We need to transfer those as well in the long run. Further, we also need some interface for stress/flux computation since most error estimator do something with the stress/flux. If its just linear elastic, then it's simple you can just have a function somehow, but what about the case with internal variables? I think dealii also doesn't have a solution to it, since they only provide the kelly estimator for the poisson problem AFAIU https://www.dealii.org/current/doxygen/deal.II/classKellyErrorEstimator.html

@fredrikekre
Copy link
Member

Yea at least for internal variables I don't think there is a good default solution so you would just have to provide the information for users to do it in a way that fits the problem at hand.

@termi-official
Copy link
Member

Yea at least for internal variables I don't think there is a good default solution so you would just have to provide the information for users to do it in a way that fits the problem at hand.

xref Nonlinear Finite Elements (Ch 8.6.2) by Wriggers

@bplcn
Copy link

bplcn commented Jun 14, 2024

I tried to construct a mapping dictionary of the cell ids from a old grid/forest to the new one for 2D grid.

function cell_id_map(forest_old::f, forest_new::f) where f
    
    cell_id_old = 0
    cell_id_new = 0

    cellid_old2new_Dict = Dict()

    for (ktree, tree) in enumerate(forest_old.cells)
        kleaf_new = 1

        for kleaf in 1:length(tree.leaves)
            
            cell_id_old += 1

            if forest_old.cells[ktree].leaves[kleaf].l == forest_new.cells[ktree].leaves[kleaf_new].l
                # the element is not refine  
                cell_id_new += 1
                cellid_old2new_Dict[cell_id_old] = [cell_id_new]
                kleaf_new += 1
            else 
                # refined
                cellid_old2new_Dict[cell_id_old] = [cell_id_new+1, cell_id_new+2, cell_id_new+3, cell_id_new+4]
                cell_id_new += 4
                kleaf_new += 4
            end
        end
    end
    return cellid_old2new_Dict
end

For each time of refinement, the old element may split or not change. Thus the values in the dictionary is a vector contain one or four elements. The implementation compares the old and new forest. Thus both of them should be stored.

@koehlerson
Copy link
Member Author

I tried to construct a mapping dictionary of the cell ids from a old grid/forest to the new one for 2D grid.

function cell_id_map(forest_old::f, forest_new::f) where f
    
    cell_id_old = 0
    cell_id_new = 0

    cellid_old2new_Dict = Dict()

    for (ktree, tree) in enumerate(forest_old.cells)
        kleaf_new = 1

        for kleaf in 1:length(tree.leaves)
            
            cell_id_old += 1

            if forest_old.cells[ktree].leaves[kleaf].l == forest_new.cells[ktree].leaves[kleaf_new].l
                # the element is not refine  
                cell_id_new += 1
                cellid_old2new_Dict[cell_id_old] = [cell_id_new]
                kleaf_new += 1
            else 
                # refined
                cellid_old2new_Dict[cell_id_old] = [cell_id_new+1, cell_id_new+2, cell_id_new+3, cell_id_new+4]
                cell_id_new += 4
                kleaf_new += 4
            end
        end
    end
    return cellid_old2new_Dict
end

For each time of refinement, the old element may split or not change. Thus the values in the dictionary is a vector contain one or four elements. The implementation compares the old and new forest. Thus both of them should be stored.

Cool thanks! That's a good starting point. You should consider making PRs 🚀

@koehlerson
Copy link
Member Author

I tried to construct a mapping dictionary of the cell ids from a old grid/forest to the new one for 2D grid.

function cell_id_map(forest_old::f, forest_new::f) where f
    
    cell_id_old = 0
    cell_id_new = 0

    cellid_old2new_Dict = Dict()

    for (ktree, tree) in enumerate(forest_old.cells)
        kleaf_new = 1

        for kleaf in 1:length(tree.leaves)
            
            cell_id_old += 1

            if forest_old.cells[ktree].leaves[kleaf].l == forest_new.cells[ktree].leaves[kleaf_new].l
                # the element is not refine  
                cell_id_new += 1
                cellid_old2new_Dict[cell_id_old] = [cell_id_new]
                kleaf_new += 1
            else 
                # refined
                cellid_old2new_Dict[cell_id_old] = [cell_id_new+1, cell_id_new+2, cell_id_new+3, cell_id_new+4]
                cell_id_new += 4
                kleaf_new += 4
            end
        end
    end
    return cellid_old2new_Dict
end

For each time of refinement, the old element may split or not change. Thus the values in the dictionary is a vector contain one or four elements. The implementation compares the old and new forest. Thus both of them should be stored.

We have to think about how we can construct the mapping without having two instantiated ForestBWGs so either we have refine! functions that alter some mapping struct as e.g. your dict and then we forward it to refine! and balance! or we have to return the mappings from the functions

@termi-official
Copy link
Member

Thanks @bplcn !

I tried to construct a mapping dictionary of the cell ids from a old grid/forest to the new one for 2D grid.

function cell_id_map(forest_old::f, forest_new::f) where f
    
    cell_id_old = 0
    cell_id_new = 0

    cellid_old2new_Dict = Dict()

    for (ktree, tree) in enumerate(forest_old.cells)
        kleaf_new = 1

        for kleaf in 1:length(tree.leaves)
            
            cell_id_old += 1

            if forest_old.cells[ktree].leaves[kleaf].l == forest_new.cells[ktree].leaves[kleaf_new].l
                # the element is not refine  
                cell_id_new += 1
                cellid_old2new_Dict[cell_id_old] = [cell_id_new]
                kleaf_new += 1
            else 
                # refined
                cellid_old2new_Dict[cell_id_old] = [cell_id_new+1, cell_id_new+2, cell_id_new+3, cell_id_new+4]
                cell_id_new += 4
                kleaf_new += 4
            end
        end
    end
    return cellid_old2new_Dict
end

For each time of refinement, the old element may split or not change. Thus the values in the dictionary is a vector contain one or four elements. The implementation compares the old and new forest. Thus both of them should be stored.

We have to think about how we can construct the mapping without having two instantiated ForestBWGs so either we have refine! functions that alter some mapping struct as e.g. your dict and then we forward it to refine! and balance! or we have to return the mappings from the functions

This is similar to what I had in mind.

amr_map = initialize_amr_map(forest, instanced_grid)
refine!(forest, elements_to_refine, amr_map)
coarsen!(forest, elements_to_coarsen, amr_map)
balance!(forest, amp_map)

where we can now query all refinement and coarsening information through amr_map.

@fredrikekre
Copy link
Member

refine!(forest, elements_to_refine, amr_map)
coarsen!(forest, elements_to_coarsen, amr_map)

This would have the same problem as #780 (comment)? Need to refine and coarsen in a single call instead?

@termi-official
Copy link
Member

refine!(forest, elements_to_refine, amr_map)
coarsen!(forest, elements_to_coarsen, amr_map)

This would have the same problem as #780 (comment)? Need to refine and coarsen in a single call instead?

Indeed. This was just a sketch assuming that we have all indices properly determined a-priori. In general the loop body is a bit more involved

  1. build leaf mesh
  2. setup problem on leaf mesh
  3. solve problem
  4. estimate error
  5. mark cells
  6. refine marked cells
  7. balance
  8. build leaf mesh
  9. setup problem on leaf mesh
  10. solve problem
  11. estimate error again
  12. mark cells with the error estimate
  13. coarsen marked cells
  14. balance again
  15. if not converged goto 1

And now you can have different variations of

We could also provide a function

enum WatDo
   Coarsen
   Refine
   Keep
end
elements_to_adapt = Vector{WatDo}(...)
adapt_mesh!(forest, elements_to_adapt, amr_map)

@termi-official
Copy link
Member

Found 2 more issues in 3D with Maxi. Last commit has test coverage for them.

  1. Not all hanging nodes are detected.
  2. For rotated elements the unique node detection fails.

@termi-official
Copy link
Member

Failing analytical example:

using Ferrite, FerriteGmsh, SparseArrays

const ref_geometry = RefHexahedron
const geometry_type = Hexahedron
const sdim = 3

grid = generate_grid(geometry_type, ntuple(_->4, sdim));

function random_deformation_field(x::Vec{dim}) where dim
    if any(x .≈ -1.0) || any(x .≈ 1.0)
        return x
    else
        Vec{dim}(x .+ (rand(dim).-0.5)*0.1)
    end
end
transform_coordinates!(grid, random_deformation_field)
grid  = ForestBWG(grid,10)

analytical_solution(x) = atan(2*(norm(x)-0.5)/0.02)
analytical_rhs(x) = -laplace(analytical_solution,x)

function assemble_cell!(ke, fe, cellvalues, ue, coords)
    fill!(ke, 0.0)
    fill!(fe, 0.0)

    n_basefuncs = getnbasefunctions(cellvalues)
    for q_point in 1:getnquadpoints(cellvalues)
        x = spatial_coordinate(cellvalues, q_point, coords)
        dΩ = getdetJdV(cellvalues, q_point)
        for i in 1:n_basefuncs
            Nᵢ = shape_value(cellvalues, q_point, i)
            ∇Nᵢ = shape_gradient(cellvalues, q_point, i)
            fe[i] += analytical_rhs(x) * Nᵢ * dΩ
            for j in 1:n_basefuncs
                ∇Nⱼ = shape_gradient(cellvalues, q_point, j)
                ke[i, j] += ∇Nⱼ ⋅ ∇Nᵢ * dΩ
            end
        end
    end
end

function assemble_global!(K, f, a, dh, cellvalues)
    ## Allocate the element stiffness matrix and element force vector
    n_basefuncs = getnbasefunctions(cellvalues)
    ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    ## Create an assembler
    assembler = start_assemble(K, f)
    ## Loop over all cells
    for cell in CellIterator(dh)
        reinit!(cellvalues, cell)
        @views ue = a[celldofs(cell)]
        ## Compute element contribution
        coords = getcoordinates(cell)
        assemble_cell!(ke, fe, cellvalues, ue, coords)
        ## Assemble ke and fe into K and f
        assemble!(assembler, celldofs(cell), ke, fe)
    end
    return K, f
end

function solve(grid)
    dim = 3
    order = 1
    ip = Lagrange{ref_geometry, order}()
    qr = QuadratureRule{ref_geometry}(2)
    cellvalues = CellValues(qr, ip);

    dh = DofHandler(grid)
    add!(dh, :u, ip)
    close!(dh);

    ch = ConstraintHandler(dh)
    add!(ch, ConformityConstraint(:u))
    add!(ch, Dirichlet(:u, getfacetset(grid, "top"), (x, t) -> 0.0))
    add!(ch, Dirichlet(:u, getfacetset(grid, "right"), (x, t) -> 0.0))
    add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0))
    add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0))
    close!(ch);

    K = create_sparsity_pattern(dh,ch)
    f = zeros(ndofs(dh))
    a = zeros(ndofs(dh))
    assemble_global!(K, f, a, dh, cellvalues);
    apply!(K, f, ch)
    u = K \ f;
    apply!(u,ch)
    return u,dh,ch,cellvalues
end

function solve_adaptive(initial_grid::ForestBWG{sdim}) where sdim
    ip = Lagrange{ref_geometry, 1}()
    qr_sc = QuadratureRule{ref_geometry}(2)
    cellvalues = CellValues(qr_sc, ip);
    finished = false
    i = 1
    grid = deepcopy(initial_grid)
    pvd = VTKFileCollection("heat_amr.pvd",grid);
    while !finished && i<=10
        @show i
        transfered_grid = Ferrite.creategrid(grid)
        u,dh,ch,cv = solve(transfered_grid)
        # σ_gp, σ_gp_sc = compute_fluxes(u,dh)
        # projector = L2Projector(Lagrange{ref_geometry, 1}(), transfered_grid)
        # σ_dof = project(projector, σ_gp, QuadratureRule{ref_geometry}(2))
        cells_to_refine = Int[]
        error_arr = Float64[]
        for (cellid,cell) in enumerate(CellIterator(dh))
            reinit!(cellvalues, cell)
            @views ue = u[celldofs(cell)]
            error = 0.0
            for q_point in 1:getnquadpoints(cellvalues)
                x = spatial_coordinate(cv, q_point, getcoordinates(cell))
                dΩ = getdetJdV(cellvalues,q_point)
                error += abs(analytical_solution(x) - function_value(cellvalues, q_point, ue)) * dΩ
            end
            if error > 0.001
                push!(cells_to_refine,cellid)
            end
            push!(error_arr,error)
        end

        addstep!(pvd, i, dh) do vtk
            write_solution(vtk, dh, u)
            # write_projection(vtk, projector, σ_dof, "flux")
            # write_cell_data(vtk, getindex.(collect(Iterators.flatten(σ_gp_sc)),1), "flux sc x")
            # write_cell_data(vtk, getindex.(collect(Iterators.flatten(σ_gp_sc)),2), "flux sc y")
            write_cell_data(vtk, error_arr, "error")
        end

        Ferrite.refine!(grid, cells_to_refine)
        Ferrite.balanceforest!(grid)

        i += 1
        if isempty(cells_to_refine)
            finished = true
        end
    end
    close(pvd);
    transfered_grid = Ferrite.creategrid(grid)
    u,dh,ch,cv = solve(transfered_grid)
end

u,dh,ch,cv = solve_adaptive(grid)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants