-
Notifications
You must be signed in to change notification settings - Fork 85
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
Comments
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 |
I'm not sure if I understand when this will fail. Which assumptions are you talking about? Regarding the performance, that is clear grid = generate_grid(Quadrilateral, (2,2));
top = ExclusiveTopology(grid); Latest release (btiming: julia> getneighborhood(top, grid, FaceIndex(1, 2), true)
2-element Vector{FaceIndex}:
FaceIndex((2, 4))
FaceIndex((1, 2)) Master (btiming: julia> getneighborhood(top, grid, FacetIndex(1, 2), true)
2-element Vector{EdgeIndex}:
EdgeIndex((2, 4))
EdgeIndex((1, 2)) |
It fails if you cannot fulfill what is expected of a |
Have you tried on the latest |
Interesting, do you have a MWE? |
No. What is the MWE ? :) |
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. |
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)
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. |
Interesting! I can reproduce it. There should be some issue with the I think the usage of the functions are correct but I need to take a closer look. Anyhow, thanks for reporting :) |
@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 |
The initial grid here is the .inp file |
sorry for the confusion @bplcn ! I took a closer look and I think I resolved the issue on the latest commit of Now a vertex balance is triggered where the arrow tip points to, which leads on the previous version to |
I have tried the last commit. It works all right now and no more errors now. Many thanks! |
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 thefacetidx
is converted to anedgeidx
and then an edge neighborhood reconstruction is built which was only designed for the three dimensional case.Ferrite.jl/src/Grid/topology.jl
Lines 282 to 292 in 68e1ab1
What should be done is this:
Ferrite.jl/src/Grid/topology.jl
Lines 259 to 263 in 68e1ab1
but with the
edge_edge
table in 2D. Some dispatch distinction for 2D and 3D for anEdgeIndex
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.The text was updated successfully, but these errors were encountered: