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

fix topology getneighborhood call in 2D #941

Closed
koehlerson opened this issue May 23, 2024 · 18 comments · Fixed by #942
Closed

fix topology getneighborhood call in 2D #941

koehlerson opened this issue May 23, 2024 · 18 comments · Fixed by #942

Comments

@koehlerson
Copy link
Member

The renaming of the entities leads to a wrong code branch for getneighborhood(top,grid,facetidx) if the grid is in 2D. This is due to the fact that the facetidx is converted to an edgeidx and then an edge neighborhood reconstruction is built which was only designed for the three dimensional case.

cellid, local_edgeidx = edgeidx[1], edgeidx[2]
cell_edges = edges(getcells(grid,cellid))
nonlocal_edgeid = cell_edges[local_edgeidx]
cell_neighbors = getneighborhood(top,grid,CellIndex(cellid))
self_reference_local = EdgeIndex[]
for cellid in cell_neighbors
local_neighbor_edgeid = findfirst(x->issubset(x,nonlocal_edgeid),edges(getcells(grid,cellid)))
local_neighbor_edgeid === nothing && continue
local_edge = EdgeIndex(cellid,local_neighbor_edgeid)
push!(self_reference_local, local_edge)
end

What should be done is this:

if include_self
return [top.face_face_neighbor[faceidx[1],faceidx[2]].neighbor_info; faceidx]
else
return top.face_face_neighbor[faceidx[1],faceidx[2]].neighbor_info
end

but with the edge_edge table in 2D. Some dispatch distinction for 2D and 3D for an EdgeIndex should be enough. Some regression testing should be included. It could be that the first, wrong code branch returns the same as the intended one, problem for me is that there are assumptions about the grid type that I don't fulfill in the p4est branch, and that's why I noticed this. It should be changed anyway for performance reasons.

@koehlerson
Copy link
Member Author

koehlerson commented May 23, 2024

Context for the first code branch is that some edge neighborhood is reconstructed, since the table only stores the exclusively connected element via the edge. However, there are more elements connected by this edge (but also by a face) in 3D

@KnutAM
Copy link
Member

KnutAM commented May 23, 2024

I'm not sure if I understand when this will fail. Which assumptions are you talking about?
Do you have a specific case showing the that we could test for?

Regarding the performance, that is clear

grid = generate_grid(Quadrilateral, (2,2));
top = ExclusiveTopology(grid);

Latest release (btiming: 55.127 ns (2 allocations: 128 bytes))

julia> getneighborhood(top, grid, FaceIndex(1, 2), true)
2-element Vector{FaceIndex}:
 FaceIndex((2, 4))
 FaceIndex((1, 2))

Master (btiming: 363.636 ns (19 allocations: 1.50 KiB))

julia> getneighborhood(top, grid, FacetIndex(1, 2), true)
2-element Vector{EdgeIndex}:
 EdgeIndex((2, 4))
 EdgeIndex((1, 2))

@koehlerson
Copy link
Member Author

I'm not sure if I understand when this will fail. Which assumptions are you talking about? Do you have a specific case showing the that we could test for?

It fails if you cannot fulfill what is expected of a cell, because edges is called on a cell. The adaptive cells don't know what their global edge ids are. This doesn't matter too much, since the computations done in the code branch don't make sense anyway. The problem is how to test it

@bplcn
Copy link

bplcn commented Jun 12, 2024

Hi! Recently I have tryed the adaptive mesh. However I find that there is something wrong with the hanging nodes. The balance tree failures at some edges. Is it something to do with the getneighborhood function?

屏幕截图 2024-06-12 143119

@koehlerson
Copy link
Member Author

Have you tried on the latest mk/p4est commit ?

@koehlerson
Copy link
Member Author

Hi! Recently I have tryed the adaptive mesh. However I find that there is something wrong with the hanging nodes. The balance tree failures at some edges. Is it something to do with the getneighborhood function?
屏幕截图 2024-06-12 143119

To answer the question: I don't think it is related to the function, since I explicitly call the neighborhood tables now. It's more likely that this is a bug from an old commit, that is hopefully fixed on a more recent commit. There were some leftover cases in the balancing code which didn't take care of some hanging nodes. Maybe it's a bug that hasn't been fixed. Is that the adaptivity example from the docs?

@bplcn
Copy link

bplcn commented Jun 12, 2024

Hi! Recently I have tryed the adaptive mesh. However I find that there is something wrong with the hanging nodes. The balance tree failures at some edges. Is it something to do with the getneighborhood function?
屏幕截图 2024-06-12 143119

To answer the question: I don't think it is related to the function, since I explicitly call the neighborhood tables now. It's more likely that this is a bug from an old commit, that is hopefully fixed on a more recent commit. There were some leftover cases in the balancing code which didn't take care of some hanging nodes. Maybe it's a bug that hasn't been fixed. Is that the adaptivity example from the docs?

No, it not from the doc examples. The examples works all right.

@koehlerson
Copy link
Member Author

Interesting, do you have a MWE?

@bplcn
Copy link

bplcn commented Jun 12, 2024

Interesting, do you have a MWE?

No. What is the MWE ? :)

@koehlerson
Copy link
Member Author

ah sorry! minimum working example such that I can try to reproduce the non-balanced case :)

@bplcn
Copy link

bplcn commented Jun 12, 2024

ah sorry! minimum working example such that I can try to reproduce the non-balanced case :)

OK, I will prepare for it as quick as possible.

@bplcn
Copy link

bplcn commented Jun 12, 2024

Following is the script for the AMR process

using Ferrite, FerriteMeshParser
using DelimitedFiles

function iniatial_refine!(adaptive_grid, region_func; lmin=0.01)

    ip = Lagrange{RefQuadrilateral, 1}()
    qr = QuadratureRule{RefQuadrilateral}(2)
    cv = CellValues(qr, ip)

    ifVec = [false,false,false,false]

    for k = 1:100

        @info k

        marked_cells = Vector{Int64}()

        grid = Ferrite.creategrid(adaptive_grid) 

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

        cc = CellCache(dh)

        for cellid in 1:getncells(grid)
            reinit!(cc, cellid)
            reinit!(cv, cc)

            coords = getcoordinates(cc)

            for (i, x) in enumerate(coords)
                ifVec[i] = region_func(x)
            end

            if_close = true in ifVec

            for q_point in 1:getnquadpoints(cv)
                dΩ = getdetJdV(cv, q_point) #wᵢ * Jξ
                le = sqrt(dΩ)*2 
                ifVec[q_point] = le > lmin
            end

            if (true in ifVec) && if_close
                push!(marked_cells, cellid)
            end
        end

        if_refine = ~isempty(marked_cells)

        if ~if_refine
            break
        end

        # update the grid and buffer infomation
        Ferrite.refine!(adaptive_grid, marked_cells)
        Ferrite.balanceforest!(adaptive_grid)

    end

end

# Read the mesh file
grid_ini = get_ferrite_grid("snep.inp");
adaptive_grid = ForestBWG(grid_ini,#=maxlevel=#15);

if_in_reg_fun(x) =  abs(x[2]) < 0.01 && abs(x[1]) < 0.01

iniatial_refine!(adaptive_grid, if_in_reg_fun; lmin=0.0025)

grid_ini_fine = Ferrite.creategrid(adaptive_grid) 
ip = Lagrange{RefQuadrilateral, 1}()
dh = DofHandler(grid_ini_fine)
    add!(dh, :p, ip)
close!(dh)

VTKFile("ini_mesh", dh) do vtk
end

data = readdlm("marked_cell.txt",',')
n_refine = size(data, 1)
n_max_elem = size(data, 2)
marked_cell_Vec = Vector{Vector{Int64}}()
for k_refine in 1:n_refine
    mark_cells_now = Vector{Int64}()
    for kelem in 1:n_max_elem
        if ~isempty(data[k_refine, kelem])
            push!(mark_cells_now, data[k_refine, kelem])
        end
    end
    push!(marked_cell_Vec, mark_cells_now)
end

pvd = VTKFileCollection("ada.pvd", grid_ini_fine);
for k_refine in 1:n_refine
    @info k_refine
    Ferrite.refine!(adaptive_grid, marked_cell_Vec[k_refine])
    Ferrite.balanceforest!(adaptive_grid)
    grid = Ferrite.creategrid(adaptive_grid) 
    dh = DofHandler(grid)
    add!(dh, :p, ip)
    close!(dh)
    
    vtk = VTKFile("test/ada-$(k_refine)", grid) 
    addstep!(pvd,vtk,k_refine)
end
close(pvd)

awe.zip

The whole package is attached, which includes the mesh for each iteration. The edge belongs to cell 207(+1). I think that there may be some problem with the balance_tree! function.

@koehlerson
Copy link
Member Author

Interesting! I can reproduce it. There should be some issue with the balanceforest! code, because balance_tree only balances within a single tree and if I understand the example correctly, then the problem is between two trees:

image
image

I think the usage of the functions are correct but I need to take a closer look. Anyhow, thanks for reporting :)

@koehlerson
Copy link
Member Author

@termi-official and me found the problem (at least we think so). The last picture of my latest post is your initial grid, is that correct? If so, then you already start from a non-conforming grid, i.e. included hanging nodes. You need to start from a grid without hanging nodes. The topology can only handle conforming meshes and the adaptive approach relies quite heavily on the topology object

@koehlerson
Copy link
Member Author

@KnutAM maybe worth to include in #974 when working on the ExclusiveTopology docs a note or warning that the stuff only works for conforming grids

@bplcn
Copy link

bplcn commented Jun 12, 2024

@termi-official and me found the problem (at least we think so). The last picture of my latest post is your initial grid, is that correct? If so, then you already start from a non-conforming grid, i.e. included hanging nodes. You need to start from a grid without hanging nodes. The topology can only handle conforming meshes and the adaptive approach relies quite heavily on the topology object

image

The initial grid here is the .inp file snep.inp, which is a conforming one. Then the mesh is read and the mesh of center region (|x| < 0.01 and |y| <0.01 ) is refiened at first. After that the refined area will expand. So I think the initial mesh may not lead to the unbalance. The unbalanced mesh only appears at the last frame. I noticed the problem because it will lead to the error of nested affine constrain.

@koehlerson
Copy link
Member Author

koehlerson commented Jun 13, 2024

sorry for the confusion @bplcn ! I took a closer look and I think I resolved the issue on the latest commit of mk/p4est. It was some remaining logical error when the balancing happens. Specifically, the problem was that in case of refinement due to balancing another non-conformity is introduced, the check was wrong that something is "still happening" in the balancing. Could you please try the latest commit again with your code? The example works now for me:

pre balance
image

Now a vertex balance is triggered where the arrow tip points to, which leads on the previous version to
image

and now on the latest version you should get
image

@bplcn
Copy link

bplcn commented Jun 13, 2024

sorry for the confusion @bplcn ! I took a closer look and I think I resolved the issue on the latest commit of mk/p4est. It was some remaining logical error when the balancing happens. Specifically, the problem was that in case of refinement due to balancing another non-conformity is introduced, the check was wrong that something is "still happening" in the balancing. Could you please try the latest commit again with your code? The example works now for me:

pre balance image

Now a vertex balance is triggered where the arrow tip points to, which leads on the previous version to image

and now on the latest version you should get image

I have tried the last commit. It works all right now and no more errors now. Many thanks!

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

Successfully merging a pull request may close this issue.

3 participants