diff --git a/docs/make.jl b/docs/make.jl index f3a68d1689..d6375d3dd5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -68,7 +68,8 @@ bibtex_plugin = CitationBibliography( "topics/boundary_conditions.md", "topics/constraints.md", "topics/grid.md", - "topics/export.md" + "topics/export.md", + "topics/amr.md" ], "Reference" => [ "Reference overview" => "reference/index.md", diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index a8333e9665..c69a191315 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -123,3 +123,33 @@ @article{Mu:2014:IP keywords = {Discontinuous Galerkin, Finite element, Interior penalty, Second-order elliptic equations, Hybrid mesh}, abstract = {This paper provides a theoretical foundation for interior penalty discontinuous Galerkin methods for second-order elliptic equations on very general polygonal or polyhedral meshes. The mesh can be composed of any polygons or polyhedra that satisfy certain shape regularity conditions characterized in a recent paper by two of the authors, Wang and Ye (2012) [11]. The usual H1-conforming finite element methods on such meshes are either very complicated or impossible to implement in practical computation. The interior penalty discontinuous Galerkin method provides a simple and effective alternative approach which is efficient and robust. Results with such general meshes have important application in computational sciences.} } +@article{BWG2011, + title={p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees}, + author={Burstedde, Carsten and Wilcox, Lucas C and Ghattas, Omar}, + journal={SIAM Journal on Scientific Computing}, + volume={33}, + number={3}, + pages={1103--1133}, + year={2011}, + publisher={SIAM} +} +@article{IBWG2015, + title={Recursive algorithms for distributed forests of octrees}, + author={Isaac, Tobin and Burstedde, Carsten and Wilcox, Lucas C and Ghattas, Omar}, + journal={SIAM Journal on Scientific Computing}, + volume={37}, + number={5}, + pages={C497--C531}, + year={2015}, + publisher={SIAM} +} +@article{SSB2008, + title={Bottom-up construction and 2: 1 balance refinement of linear octrees in parallel}, + author={Sundar, Hari and Sampath, Rahul S and Biros, George}, + journal={SIAM Journal on Scientific Computing}, + volume={30}, + number={5}, + pages={2675--2708}, + year={2008}, + publisher={SIAM} +} diff --git a/docs/src/devdocs/AMR.md b/docs/src/devdocs/AMR.md new file mode 100644 index 0000000000..049052cc9e --- /dev/null +++ b/docs/src/devdocs/AMR.md @@ -0,0 +1,217 @@ +# Adaptive Mesh Refinement (AMR) + +## P4est + +Ferrite's P4est implementation is based on these papers: + +- [BWG2011](@citet) +- [IBWG2015](@citet) + +where almost everything is implemented in a serial way from the first paper. +Only certain specific algorithms of the second paper are implemented and there is a lot of open work to include the iterators of the second paper. +Look into the issues of Ferrite.jl and search for the AMR tag. + +### Important Concepts + +One of the most important concepts, which everything is based on, are space filling curves (SFC). +In particular, [Z-order (also named Morton order, Morton space-filling curves)](https://en.wikipedia.org/wiki/Z-order_curve) are used in p4est. +The basic idea is that each Octant (in 3D) or quadrant (in 2D) can be encoded by 2 quantities + +- the level `l` +- the lower left (front) coordinates `xyz` + +Based on them a unique identifier, the morton index, can be computed. +The mapping from (`l`, `xyz`) -> `mortonidx(l,xyz)` is bijective, meaning we can flip the approach +and can construct each octant/quadrant solely by the `mortonidx` and a given level `l`. + +The current implementation of an octant looks currently like this: +```julia +struct OctantBWG{dim, N, T} <: AbstractCell{RefHypercube{dim}} + #Refinement level + l::T + #x,y,z \in {0,...,2^b} where (0 ≤ l ≤ b)} + xyz::NTuple{dim,T} +end +``` +whenever coordinates are considered we follow the z order logic, meaning x before y before z. +Note that the acronym BWG stands for the initials of the surname of the authors of the p4est paper. +The coordinates of an octant are described in the *octree coordinate system* which goes from $[0,2^b]^{dim}$. +The parameter $b$ describes the maximum level of refinement and is set a priori. +Another important aspect of the octree coordinate system is, that it is a discrete integer coordinate system. +The size of an octant at the lowest possible level `b` is always 1, sometimes these octants are called atoms. + +The octree is implemented as: +```julia +struct OctreeBWG{dim,N,T} <: AbstractAdaptiveCell{RefHypercube{dim}} + leaves::Vector{OctantBWG{dim,N,T}} + #maximum refinement level + b::T + nodes::NTuple{N,Int} +end +``` + +So, only the leaves of the tree are stored and not any intermediate refinement level. +The field `b` is the maximum refinement level and is crucial. This parameter determines the size of the octree coordinate system. +The octree coordinate system is the coordinate system in which the coordinates `xyz` of any `octant::OctantBWG` are described. + +### Examples + +Let's say the maximum octree level is $b=3$, then the coordinate system is in 2D $[0,2^3]^2 = [0, 8]^2$. +So, our root is on level 0 of size 8 and has the lower left coordinates `(0,0)` + +```julia +# different constructors available, first one OctantBWG(dim,level,mortonid,maximumlevel) +# other possibility by giving directly level and a tuple of coordinates OctantBWG(level,(x,y)) +julia> dim = 2; level = 0; maximumlevel=3 +julia> oct = OctantBWG(dim,level,1,maximumlevel) +OctantBWG{2,4,4} + l = 0 + xy = 0,0 +``` +The size of octants at a specific level can be computed by a simple operation +```julia +julia> Ferrite.AMR._compute_size(#=b=#3,#=l=#0) +8 +``` +This computation is based on the relation $\text{size}=2^{b-l}$. +Now, to fully understand the octree coordinate system we go a level down, i.e. we cut the space in $x$ and $y$ in half. +This means, that the octants are now of size $2^{3-1}=4$. +Construct all level 1 octants based on mortonid: +```julia +# note the arguments are dim,level,mortonid,maximumlevel +julia> dim = 2; level = 1; maximumlevel = 3 +julia> oct = Ferrite.AMR.OctantBWG(dim, level, 1, maximumlevel) +OctantBWG{2,4,4} + l = 1 + xy = 0,0 + +julia> oct = Ferrite.AMR.OctantBWG(dim, level, 2, maximumlevel) +OctantBWG{2,4,4} + l = 1 + xy = 4,0 + +julia> oct = Ferrite.AMR.OctantBWG(dim, level, 3, maximumlevel) +OctantBWG{2,4,4} + l = 1 + xy = 0,4 + +julia> oct = Ferrite.AMR.OctantBWG(dim, level, 4, maximumlevel) +OctantBWG{2,4,4} + l = 1 + xy = 4,4 +``` + +So, the morton index is on **one** specific level just a x before y before z "cell" or "element" identifier +``` +x-----------x-----------x +| | | +| | | +| 3 | 4 | +| | | +| | | +x-----------x-----------x +| | | +| | | +| 1 | 2 | +| | | +| | | +x-----------x-----------x +``` + +The operation to compute octants/quadrants is cheap, since it is just bitshifting. +An important aspect of the morton index is that it's only consecutive on **one** level in this specific implementation. +Note that other implementation exists that incorporate the level integer within the morton identifier and by that have a unique identifier across levels. +If you have a tree like this below: + +``` +x-----------x-----------x +| | | +| | | +| 9 | 10 | +| | | +| | | +x-----x--x--x-----------x +| |6 |7 | | +| 3 x--x--x | +| |4 |5 | | +x-----x--x--x 8 | +| | | | +| 1 | 2 | | +x-----x-----x-----------x +``` + +you would maybe think this is the morton index, but strictly speaking it is not. +What we see above is just the `leafindex`, i.e. the index where you find this leaf in the `leaves` array of `OctreeBWG`. +Let's try to construct the lower right based on the morton index on level 1 + +```julia +julia> o = Ferrite.OctantBWG(2,1,8,3) +ERROR: AssertionError: m ≤ (one(T) + one(T)) ^ (dim * l) # 8 > 4 +Stacktrace: + [1] OctantBWG(dim::Int64, l::Int32, m::Int32, b::Int32) + @ Ferrite ~/repos/Ferrite.jl/src/Adaptivity/AdaptiveCells.jl:23 + [2] OctantBWG(dim::Int64, l::Int64, m::Int64, b::Int64) + @ Ferrite ~/repos/Ferrite.jl/src/Adaptivity/AdaptiveCells.jl:43 + [3] top-level scope + @ REPL[93]:1 +``` + +The assertion expresses that it is not possible to construct a morton index 8 octant, since the upper bound of the morton index is 4 on level 1. +The morton index of the lower right cell is 2 on level 1. + +```julia +julia> o = Ferrite.AMR.OctantBWG(2,1,2,3) +OctantBWG{2,4,4} + l = 1 + xy = 4,0 +``` + +### Octant operation + +There are multiple useful functions to compute information about an octant e.g. parent, childs, etc. + +```@docs +Ferrite.AMR.isancestor +Ferrite.AMR.morton +Ferrite.AMR.children +Ferrite.AMR.vertices +Ferrite.AMR.edges +Ferrite.AMR.faces +Ferrite.AMR.transform_pointBWG +``` + +### Intraoctree operation + +Intraoctree operation stay within one octree and compute octants that are attached in some way to a pivot octant `o`. +These operations are useful to collect unique entities within a single octree or to compute possible neighbors of `o`. +[BWG2011](@citet) Algorithm 5, 6, and 7 describe the following intraoctree operations: + +```@docs +Ferrite.AMR.corner_neighbor +Ferrite.AMR.edge_neighbor +Ferrite.AMR.facet_neighbor +Ferrite.AMR.possibleneighbors +``` + +### Interoctree operation + +Interoctree operation are in contrast to intraoctree operation by computing octant transformations across different octrees. +Thereby, one needs to account for topological connections between the octrees as well as possible rotations of the octrees. +[BWG2011](@citet) Algorithm 8, 10, and 12 explain the algorithms that are implemented in the following functions: + +```@docs +Ferrite.AMR.transform_corner +Ferrite.AMR.transform_edge +Ferrite.AMR.transform_facet +``` + +Note that we flipped the input and to expected output logic a bit to the proposed algorithms of the paper. +However, the original proposed versions are implemented as well in: + +```@docs +Ferrite.AMR.transform_corner_remote +Ferrite.AMR.transform_edge_remote +Ferrite.AMR.transform_facet_remote +``` + +despite being never used in the code base so far. diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 9c16b83d7c..aafd41f00f 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -5,5 +5,5 @@ developing the library. ```@contents Depth = 1 -Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md"] +Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md", "AMR.md"] ``` diff --git a/docs/src/literate-tutorials/adaptivity.jl b/docs/src/literate-tutorials/adaptivity.jl new file mode 100644 index 0000000000..fae90aeb5e --- /dev/null +++ b/docs/src/literate-tutorials/adaptivity.jl @@ -0,0 +1,196 @@ +using Ferrite, FerriteGmsh, SparseArrays +#cells = [Hexahedron((1, 2, 5, 4, 10, 11, 14, 13)),Hexahedron((4, 5, 8, 7, 13, 14, 17, 16)), Hexahedron((3, 4, 7, 6, 12, 13, 16, 15)), Hexahedron((12, 13, 16, 15, 21, 22, 25, 24)), Hexahedron((13, 14, 17, 16, 22, 23, 26, 25)), Hexahedron((10, 11, 14, 13, 19, 20, 23, 22)), Hexahedron((9, 10, 13, 12, 18, 19, 22, 21))] +### beefy L cells = [Hexahedron((4, 5, 8, 7, 13, 14, 17, 16)), Hexahedron((3, 4, 7, 6, 12, 13, 16, 15)), Hexahedron((12, 13, 16, 15, 21, 22, 25, 24)), Hexahedron((13, 14, 17, 16, 22, 23, 26, 25)), Hexahedron((10, 11, 14, 13, 19, 20, 23, 22)), Hexahedron((9, 10, 13, 12, 18, 19, 22, 21))] +### cells = [Hexahedron((3, 4, 7, 6, 12, 13, 16, 15)), Hexahedron((12, 13, 16, 15, 21, 22, 25, 24)), Hexahedron((9, 10, 13, 12, 18, 19, 22, 21))] +#nodes = Node{3, Float64}[Node{3, Float64}(Vec{3}((0.0, -1.0, -1.0))), Node{3, Float64}(Vec{3}((1.0, -1.0, -1.0))), Node{3, Float64}(Vec{3}((-1.0, 0.0, -1.0))), Node{3, Float64}(Vec{3}((0.0, 0.0, -1.0))), Node{3, Float64}(Vec{3}((1.0, 0.0, -1.0))), Node{3, Float64}(Vec{3}((-1.0, 1.0, -1.0))), Node{3, Float64}(Vec{3}((0.0, 1.0, -1.0))), Node{3, Float64}(Vec{3}((1.0, 1.0, -1.0))), Node{3, Float64}(Vec{3}((-1.0, -1.0, 0.0))), Node{3, Float64}(Vec{3}((0.0, -1.0, 0.0))), Node{3, Float64}(Vec{3}((1.0, -1.0, 0.0))), Node{3, Float64}(Vec{3}((-1.0, 0.0, 0.0))), Node{3, Float64}(Vec{3}((0.0, 0.0, 0.0))), Node{3, Float64}(Vec{3}((1.0, 0.0, 0.0))), Node{3, Float64}(Vec{3}((-1.0, 1.0, 0.0))), Node{3, Float64}(Vec{3}((0.0, 1.0, 0.0))), Node{3, Float64}(Vec{3}((1.0, 1.0, 0.0))), Node{3, Float64}(Vec{3}((-1.0, -1.0, 1.0))), Node{3, Float64}(Vec{3}((0.0, -1.0, 1.0))), Node{3, Float64}(Vec{3}((1.0, -1.0, 1.0))), Node{3, Float64}(Vec{3}((-1.0, 0.0, 1.0))), Node{3, Float64}(Vec{3}((0.0, 0.0, 1.0))), Node{3, Float64}(Vec{3}((1.0, 0.0, 1.0))), Node{3, Float64}(Vec{3}((-1.0, 1.0, 1.0))), Node{3, Float64}(Vec{3}((0.0, 1.0, 1.0))), Node{3, Float64}(Vec{3}((1.0, 1.0, 1.0)))] +#grid = Grid(cells,nodes) +#addfacetset!(grid,"front",x->x[2]≈-1) +#addfacetset!(grid,"back",x->x[2]≈1) +#addfacetset!(grid,"left",x->x[3]≈-1) +#addfacetset!(grid,"right",x->x[3]≈1) +#addfacetset!(grid,"top",x->x[1]≈-1) +#addfacetset!(grid,"bottom",x->x[1]≈1) +#addfacetset!(grid,"pull",x->x[1]≈0 && x[2] <= 0.5 && x[3] <= 0.5) +grid = generate_grid(Hexahedron,(2,1,1)) +grid = ForestBWG(grid,20) +#Ferrite.refine_all!(grid,1) +Ferrite.refine!(grid,[1,2]) +Ferrite.balanceforest!(grid) + +struct Elasticity + G::Float64 + K::Float64 +end + +function material_routine(material::Elasticity, ε::SymmetricTensor{dim}) where dim + (; G, K) = material + stress(ε) = 2G * dev(ε) + K * tr(ε) * one(ε) + ∂σ∂ε, σ = gradient(stress, ε, :all) + return σ, ∂σ∂ε +end + +E = 200e3 # Young's modulus [MPa] +ν = 0.2 # Poisson's ratio [-] +material = Elasticity(E/2(1+ν), E/3(1-2ν)); + +function assemble_cell!(ke, fe, cellvalues, material, ue) + fill!(ke, 0.0) + fill!(fe, 0.0) + + n_basefuncs = getnbasefunctions(cellvalues) + for q_point in 1:getnquadpoints(cellvalues) + ## For each integration point, compute strain, stress and material stiffness + ε = function_symmetric_gradient(cellvalues, q_point, ue) + σ, ∂σ∂ε = material_routine(material, ε) + + dΩ = getdetJdV(cellvalues, q_point) + for i in 1:n_basefuncs + ∇Nᵢ = shape_gradient(cellvalues, q_point, i) + fe[i] += σ ⊡ ∇Nᵢ * dΩ + for j in 1:n_basefuncs + ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j) + ke[i, j] += (∂σ∂ε ⊡ ∇ˢʸᵐNⱼ) ⊡ ∇Nᵢ * dΩ + end + end + end +end + +function assemble_global!(K, f, a, dh, cellvalues, material) + ## 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 + assemble_cell!(ke, fe, cellvalues, material, ue) + ## 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{RefHexahedron, order}()^dim + qr = QuadratureRule{RefHexahedron}(2) + cellvalues = CellValues(qr, ip); + + dh = DofHandler(grid) + add!(dh, :u, ip) + close!(dh); + + ch = ConstraintHandler(dh) + add!(ch, Ferrite.ConformityConstraint(:u)) + add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> Vec{3}((0.0,0.0,0.0)), [1,2,3])) + add!(ch, Dirichlet(:u, getfacetset(grid, "back"), (x, t) -> Vec{3}((0.0,0.0,0.0)), [1,2,3])) + add!(ch, Dirichlet(:u, getfacetset(grid, "front"), (x, t) -> Vec{3}((0.0,0.0,0.0)), [1,2,3])) + add!(ch, Dirichlet(:u, getfacetset(grid, "top"), (x, t) -> Vec{3}((0.0,0.0,0.0)), [1,2,3])) + add!(ch, Dirichlet(:u, getfacetset(grid, "right"), (x, t) -> Vec{3}((0.0,0.0,0.0)), [1,2,3])) + add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> Vec{3}((-0.1,0.0,0.0)), [1,2,3])) + close!(ch); + + K = create_sparsity_pattern(dh,ch) + f = zeros(ndofs(dh)) + a = zeros(ndofs(dh)) + assemble_global!(K, f, a, dh, cellvalues, material); + apply!(K, f, ch) + u = K \ f; + apply!(u,ch) + return u,dh,ch,cellvalues +end + +function compute_fluxes(u,dh::DofHandler{dim}) where dim + ip = Lagrange{RefHexahedron, 1}()^dim + qr = QuadratureRule{RefHexahedron}(2) + cellvalues_sc = CellValues(qr, ip); + ## "Normal" quadrature points for the fluxes + qr = QuadratureRule{RefHexahedron}(2) + cellvalues = CellValues(qr, ip); + ## Buffers + σ_gp_sc = Vector{Vector{SymmetricTensor{2,dim,Float64,dim*2}}}() + σ_gp_sc_loc = Vector{SymmetricTensor{2,dim,Float64,dim*2}}() + σ_gp = Vector{Vector{SymmetricTensor{2,dim,Float64,dim*2}}}() + σ_gp_loc = Vector{SymmetricTensor{2,dim,Float64,dim*2}}() + for (cellid,cell) in enumerate(CellIterator(dh)) + reinit!(cellvalues, cell) + reinit!(cellvalues_sc, cell) + @views ue = u[celldofs(cell)] + for q_point in 1:getnquadpoints(cellvalues) + ε = function_symmetric_gradient(cellvalues, q_point, ue) + σ, _ = material_routine(material, ε) + push!(σ_gp_loc, σ) + end + for q_point in 1:getnquadpoints(cellvalues_sc) + ε = function_symmetric_gradient(cellvalues_sc, q_point, ue) + σ, _ = material_routine(material, ε) + push!(σ_gp_sc_loc, σ) + end + push!(σ_gp,copy(σ_gp_loc)) + push!(σ_gp_sc,copy(σ_gp_sc_loc)) + ## Reset buffer for local points + empty!(σ_gp_loc) + empty!(σ_gp_sc_loc) + end + return σ_gp, σ_gp_sc +end + +function solve_adaptive(initial_grid) + ip = Lagrange{RefHexahedron, 1}() + qr = QuadratureRule{RefHexahedron}(2) + cellvalues_tensorial = CellValues(qr, ip); + finished = false + i = 1 + grid = initial_grid + transfered_grid = Ferrite.creategrid(grid) + pvd = VTKFileCollection("linear-elasticity.pvd",transfered_grid) + while !finished + transfered_grid = Ferrite.creategrid(grid) + u,dh,ch,cv = solve(transfered_grid) + σ_gp, σ_gp_sc = compute_fluxes(u,dh) + projector = L2Projector(Lagrange{RefHexahedron, 1}()^3, transfered_grid) + σ_dof = project(projector, σ_gp, QuadratureRule{RefHexahedron}(2)) + cells_to_refine = Int[] + error_arr = Float64[] + for (cellid,cell) in enumerate(CellIterator(projector.dh)) + reinit!(cellvalues_tensorial, cell) + @views σe = σ_dof[celldofs(cell)] + error = 0.0 + for q_point in 1:getnquadpoints(cellvalues_tensorial) + σ_dof_at_sc = function_value(cellvalues_tensorial, q_point, σe) + error += norm((σ_gp_sc[cellid][q_point] - σ_dof_at_sc )) * getdetJdV(cellvalues_tensorial,q_point) + end + push!(error_arr,error) + end + η = maximum(error_arr) + θ = 0.5 + for (cellid,cell_err) in enumerate(error_arr) + if cell_err > θ*η + push!(cells_to_refine,cellid) + end + end + addstep!(pvd,i,transfered_grid) do vtk + write_solution(vtk, dh, u) + write_projection(vtk, projector, σ_dof, "stress") + write_cell_data(vtk, getindex.(collect(Iterators.flatten(σ_gp_sc)),1), "stress sc") + write_cell_data(vtk, error_arr, "error") + end + + Ferrite.refine!(grid,cells_to_refine) + Ferrite.balanceforest!(grid) + + i += 1 + if isempty(cells_to_refine) || maximum(error_arr) < 0.05 + finished = true + end + break + end + close(pvd) +end + +solve_adaptive(grid) diff --git a/docs/src/literate-tutorials/heat_adaptivity.jl b/docs/src/literate-tutorials/heat_adaptivity.jl new file mode 100644 index 0000000000..955c1c2ee6 --- /dev/null +++ b/docs/src/literate-tutorials/heat_adaptivity.jl @@ -0,0 +1,171 @@ +using Ferrite, FerriteGmsh, SparseArrays +grid = generate_grid(Hexahedron, (4,4,4)); +function random_deformation_field(x) + if any(x .≈ -1.0) || any(x .≈ 1.0) + return x + else + Vec{3}(x .+ (rand(3).-0.5)*0.15) + 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{RefHexahedron, order}() + qr = QuadratureRule{RefHexahedron}(2) + cellvalues = CellValues(qr, ip); + + dh = DofHandler(grid) + add!(dh, :u, ip) + close!(dh); + + ch = ConstraintHandler(dh) + 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)) + add!(ch, ConformityConstraint(:u)) + 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 compute_fluxes(u,dh) + ip = Lagrange{RefHexahedron, 1}() + ## Normal quadrature points + qr = QuadratureRule{RefHexahedron}(2) + cellvalues = CellValues(qr, ip); + ## Superconvergent point + qr_sc = QuadratureRule{RefHexahedron}(1) + cellvalues_sc = CellValues(qr_sc, ip); + ## Buffers + σ_gp = Vector{Vector{Vec{3,Float64}}}() + σ_gp_loc = Vector{Vec{3,Float64}}() + σ_gp_sc = Vector{Vector{Vec{3,Float64}}}() + σ_gp_sc_loc = Vector{Vec{3,Float64}}() + for (cellid,cell) in enumerate(CellIterator(dh)) + @views ue = u[celldofs(cell)] + + reinit!(cellvalues, cell) + for q_point in 1:getnquadpoints(cellvalues) + gradu = function_gradient(cellvalues, q_point, ue) + push!(σ_gp_loc, gradu) + end + push!(σ_gp,copy(σ_gp_loc)) + empty!(σ_gp_loc) + + reinit!(cellvalues_sc, cell) + for q_point in 1:getnquadpoints(cellvalues_sc) + gradu = function_gradient(cellvalues_sc, q_point, ue) + push!(σ_gp_sc_loc, gradu) + end + push!(σ_gp_sc,copy(σ_gp_sc_loc)) + empty!(σ_gp_sc_loc) + end + return σ_gp, σ_gp_sc +end + +function solve_adaptive(initial_grid) + ip = Lagrange{RefHexahedron, 1}()^3 + qr_sc = QuadratureRule{RefHexahedron}(1) + cellvalues_flux = 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{RefHexahedron, 1}(), transfered_grid) + σ_dof = project(projector, σ_gp, QuadratureRule{RefHexahedron}(2)) + cells_to_refine = Int[] + error_arr = Float64[] + for (cellid,cell) in enumerate(CellIterator(projector.dh)) + reinit!(cellvalues_flux, cell) + @views σe = σ_dof[celldofs(cell)] + error = 0.0 + for q_point in 1:getnquadpoints(cellvalues_flux) + σ_dof_at_sc = function_value(cellvalues_flux, q_point, σe) + error += norm((σ_gp_sc[cellid][q_point] - σ_dof_at_sc )) + error *= getdetJdV(cellvalues_flux,q_point) + 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); +end + +solve_adaptive(grid) diff --git a/docs/src/literate-tutorials/l.geo b/docs/src/literate-tutorials/l.geo new file mode 100644 index 0000000000..b0d2bbbca1 --- /dev/null +++ b/docs/src/literate-tutorials/l.geo @@ -0,0 +1,50 @@ +//+ +Point(1) = {0, 0, 0, 1.0}; +//+ +Point(2) = {1, 0, 0, 1.0}; +//+ +Point(3) = {1, 0.5, 0, 1.0}; +//+ +Point(4) = {0.5, 0.5, 0, 1.0}; +//+ +Point(5) = {0.5, 1, 0, 1.0}; +//+ +Point(6) = {0, 1, 0, 1.0}; +//+ +Line(1) = {1, 2}; +//+ +Line(2) = {2, 3}; +//+ +Line(3) = {3, 4}; +//+ +Line(4) = {4, 5}; +//+ +Line(5) = {5, 6}; +//+ +Line(6) = {6, 1}; +//+ +Line(7) = {1, 4}; +//+ +Curve Loop(1) = {1, 2, 3, 4, 5, 6}; +//+ +Curve Loop(2) = {3, -7, 1, 2}; +//+ +Plane Surface(1) = {2}; +//+ +Curve Loop(3) = {4, 5, 6, 7}; +//+ +Plane Surface(2) = {3}; +//+ +Transfinite Surface {1}; +//+ +Transfinite Surface {2}; +//+ +Physical Curve("top", 8) = {5}; +//+ +Physical Curve("left", 9) = {6}; +//+ +Physical Curve("bottom", 10) = {1}; +//+ +Physical Curve("right", 11) = {2}; +//+ +Physical Surface("domain", 12) = {2, 1}; diff --git a/docs/src/literate-tutorials/l.msh b/docs/src/literate-tutorials/l.msh new file mode 100644 index 0000000000..f37a7b98ea --- /dev/null +++ b/docs/src/literate-tutorials/l.msh @@ -0,0 +1,71 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 8 "top" +1 9 "left" +1 10 "bottom" +1 11 "right" +2 12 "domain" +$EndPhysicalNames +$Entities +6 7 2 0 +1 0 0 0 0 +2 1 0 0 0 +3 1 0.5 0 0 +4 0.5 0.5 0 0 +5 0.5 1 0 0 +6 0 1 0 0 +1 0 0 0 1 0 0 1 10 2 1 -2 +2 1 0 0 1 0.5 0 1 11 2 2 -3 +3 0.5 0.5 0 1 0.5 0 0 2 3 -4 +4 0.5 0.5 0 0.5 1 0 0 2 4 -5 +5 0 1 0 0.5 1 0 1 8 2 5 -6 +6 0 0 0 0 1 0 1 9 2 6 -1 +7 0 0 0 0.5 0.5 0 0 2 1 -4 +1 0 0 0 1 0.5 0 1 12 4 3 -7 1 2 +2 0 0 0 0.5 1 0 1 12 4 4 5 6 7 +$EndEntities +$Nodes +12 6 1 6 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +1 0.5 0 +0 4 0 1 +4 +0.5 0.5 0 +0 5 0 1 +5 +0.5 1 0 +0 6 0 1 +6 +0 1 0 +1 1 0 0 +1 2 0 0 +1 5 0 0 +1 6 0 0 +2 1 0 0 +2 2 0 0 +$EndNodes +$Elements +6 6 1 6 +1 1 1 1 +1 1 2 +1 2 1 1 +2 2 3 +1 5 1 1 +3 5 6 +1 6 1 1 +4 6 1 +2 1 3 1 +5 3 4 1 2 +2 2 3 1 +6 4 5 6 1 +$EndElements diff --git a/docs/src/topics/amr.md b/docs/src/topics/amr.md new file mode 100644 index 0000000000..b2353c0927 --- /dev/null +++ b/docs/src/topics/amr.md @@ -0,0 +1,141 @@ +# Adaptive Mesh Refinement + +Adaptive mesh refinement (AMR) is a computational technique used in finite element analysis to enhance the accuracy and efficiency of simulations. +It involves dynamically adjusting the mesh resolution based on some criteria. +By refining the mesh in regions where the solution exhibits features of interest, AMR ensures that computational resources are concentrated where they are most needed, leading to more accurate results without a proportional increase in computational cost. +This refinement can be achieved in different fashions by either adjusting the mesh size (h-adaptivity), the polynomial order of the Ansatz functions (p-adaptivity) or the nodal positions (r-adaptivity). + +In Ferrite.jl, adaptivity is achieved through a p4est type of implementation which covers h-adaptivity. +This approach is designed to handle unstructured hexahedral (in 3D) and quadrilateral (in 2D) meshes. +A further restriction of the p4est type of implementation is isotropic refinement, meaning that elements are subdivided into smaller elements of the same shape. + +In AMR different phenomena and vocabulary emerge which we group into the following aspects + +- hanging nodes +- balancing +- error estimation + +## Hanging Nodes +### What Are Hanging Nodes? + +Hanging nodes occur during the process of mesh refinement. +When a mesh is refined, some elements may be subdivided while their neighboring elements are not. +This results in nodes that lie on the edges or faces of coarser elements but do not coincide with the existing nodes of those elements. +These intermediate nodes are referred to as hanging nodes. + +Consider the following situation: +``` +x-----------x-----------x x-----------x-----------x +|7 |8 |9 |7 |8 |9 +| | | | | | +| | | | | | +| | | | | | +| | | | | | +| | | | | | +| | | refinement | | | +x-----------x-----------x --------> x-----x-----x-----------x +|4 |5 |6 |4 |14 |5 |6 +| | | | | | | +| | | | | | | +| | | x-----x-----x | +| | | |11 |12 |13 | +| | | | | | | +| | | | | | | +x-----------x-----------x x-----x-----x-----------x + 1 2 3 1 10 2 3 +``` +The new introduced nodes 10, 11 and 12 are shared and therefore are not hanging nodes. +However, the nodes 14 and 13 are not shared with the neighboring coarser element and therefore hanging. + +### Implications of Hanging Nodes + +The presence of hanging nodes poses the challenge of non-conformity: +A mesh with hanging nodes is non-conforming because the finite element mesh no longer adheres to the requirement that elements meet only at their nodes or along entire edges or faces. +This lack of conformity can lead to difficulties in maintaining the continuity of the solution across the mesh. +However we can recover the continuity of the solution by constraining the hanging nodes. + + +### How to Treat Hanging Nodes + +To address the issues introduced by hanging nodes, specific strategies and constraints are employed. +The degrees of freedom (DoFs) associated with hanging nodes are constrained based on the surrounding coarser mesh elements. +For example, in a linear finite element method, the value at a hanging node can be constrained to be the average of the values at the adjacent vertices of the coarser element. +As for the example above node 13 could be constrained to $\boldsymbol{u}[13]=0.5\boldsymbol{u}[5]+0.5\boldsymbol{u}[2]$. +As soon as higher polynomial degrees are involved, things become more involved. +In Ferrite, a conformity constraint can be constructed with the ConstraintHandler when using a DofHandler which has been constructed with a grid passed from `Ferrite.creategrid(adaptive_grid::ForestBWG)`. +This conformity constraint ensures that each hanging node is constrained appropriately. + +```julia +ch = ConstraintHandler(dh) +add!(ch, Ferrite.ConformityConstraint(:u)) +``` + +## Balancing + +Hanging nodes can depend in a nested fashion on other hanging nodes. +This case complicates the `ConformityConstraint` massively and is therefore often prohibited by the act of balancing. +Mesh balancing refers to the process of ensuring a smooth and gradual transition in element sizes across the computational mesh. +This is especially important in adaptive mesh refinement (AMR), where different regions of the mesh may undergo varying levels of refinement based on the solution's characteristics. +The goal of mesh balancing is to prevent isolated regions of excessive refinement, which can lead to inefficiencies in terms of constructing conformity constraints and numerical instability. +A famous balancing approach is the 2:1 balance, where it is ensured that a hanging node only depends on non-hanging nodes. +Therefore, exactly one level of non-conformity is allowed. +An example of this process is visualized below. + +``` +x-----------x-----------x x-----------x-----------x +| | | | | | | +| | | | | | | +| | | | | | | +| | | |-----x-----| | +| | | | | | | +| | | | | | | +| | | balancing | | | | +x-----x--x--x-----------x --------> x-----x--x--x-----------x +| | | | | | | | | | | +| x--x--x | | x--x--x | | +| | | | | | | | | | | +x-----x--x--x | x-----x--x--x-----x-----| +| | | | | | | | | +| | | | | | | | | +| | | | | | | | | +x-----x-----x-----------x x-----x-----x-----------x +``` + +Note that in the example above, the top right element hasn't been refined. +However, in some cases it is advantageous to do so in order to have a smoother transition in element size. +Therefore, by default, the adaptive mesh is also refined over this vertex, leading to the following result: + +``` +x-----------x-----------x x-----------x-----------x +| | | | | | | | +| | | | | | | | +| | | | | | | | +| | | |-----x-----|-----x-----| +| | | | | | | | +| | | | | | | | +| | | balancing | | | | | +x-----x--x--x-----------x --------> x-----x--x--x-----------x +| | | | | | | | | | | +| x--x--x | | x--x--x | | +| | | | | | | | | | | +x-----x--x--x | x-----x--x--x-----x-----| +| | | | | | | | | +| | | | | | | | | +| | | | | | | | | +x-----x-----x-----------x x-----x-----x-----------x +``` + +In Ferrite's p4est implementation, one can call `balanceforest!` to balance the adaptive grid. + +```julia +Ferrite.balanceforest!(adaptive_grid) +``` + +## Error Estimation +Error estimation is a critical component of adaptive mesh refinement (AMR) in finite element analysis. +The primary objective of error estimation is to identify regions of the computational domain where the numerical solution is less accurate and requires further refinement. +Accurate error estimation guides the adaptive refinement process, ensuring that computational resources are concentrated in areas where they will have the most significant impact on improving solution accuracy. + +In practice, error estimators evaluate the local error in the finite element solution by comparing it to an approximation of the true solution. +Common techniques include residual-based error estimation, where the residuals of the finite element equations are used to estimate the local error, and recovery-based methods, which involve constructing a higher-order approximation of the solution and assessing the difference between this approximation and the finite element solution. +By identifying elements with high estimated errors, these methods provide a targeted approach to mesh refinement, enhancing the overall efficiency and accuracy of the simulation. diff --git a/docs/src/topics/index.md b/docs/src/topics/index.md index 6e8a1d7642..e5ae4c1ded 100644 --- a/docs/src/topics/index.md +++ b/docs/src/topics/index.md @@ -11,6 +11,7 @@ Pages = [ "boundary_conditions.md", "constraints.md", "grid.md", - "export.md" + "export.md", + "amr.md" ] ``` diff --git a/src/Adaptivity/AMR.jl b/src/Adaptivity/AMR.jl new file mode 100644 index 0000000000..ab5eaee2f4 --- /dev/null +++ b/src/Adaptivity/AMR.jl @@ -0,0 +1,21 @@ +module AMR + +using .. Ferrite +import Ferrite: @debug +using SparseArrays +using OrderedCollections + +include("BWG.jl") +include("ncgrid.jl") +include("constraints.jl") + +export ForestBWG, + refine!, + refine_all!, + coarsen!, + balanceforest!, + creategrid, + ConformityConstraint, + NonConformingGrid + +end diff --git a/src/Adaptivity/BWG.jl b/src/Adaptivity/BWG.jl new file mode 100644 index 0000000000..5bea4c4c0e --- /dev/null +++ b/src/Adaptivity/BWG.jl @@ -0,0 +1,2009 @@ +# TODO we should remove the mixture of indices. Maybe with these: +# - struct FacetIndexBWG ... end +# - struct QuadrilateralBWG ... end +# - struct HexahedronBWG ... end + +abstract type AbstractAdaptiveGrid{dim} <: Ferrite.AbstractGrid{dim} end +abstract type AbstractAdaptiveCell{refshape <: Ferrite.AbstractRefShape} <: Ferrite.AbstractCell{refshape} end + +const ncorners_face3D = 4 +const ncorners_face2D = 2 +const ncorners_edge = ncorners_face2D + +_maxlevel = [30,19] + +function set_maxlevel(dim::Integer,maxlevel::Integer) + _maxlevel[dim-1] = maxlevel +end + +struct OctantBWG{dim, N, T} <: Ferrite.AbstractCell{Ferrite.RefHypercube{dim}} + #Refinement level + l::T + #x,y,z \in {0,...,2^b} where (0 ≤ l ≤ b)} + xyz::NTuple{dim,T} +end + +""" + OctantBWG(dim::Integer, l::Integer, b::Integer, m::Integer) +Construct an `octant` based on dimension `dim`, level `l`, amount of levels `b` and morton index `m` +""" +function OctantBWG(dim::Integer, l::T1, m::T2, b::T1=_maxlevel[dim-1]) where {T1 <: Integer, T2 <: Integer} + @assert l ≤ b #maximum refinement level exceeded + @assert m ≤ (one(T1)+one(T1))^(dim*l) + x,y,z = (zero(T1),zero(T1),zero(T1)) + h = Int32(_compute_size(b,l)) + _zero = zero(T1) + _one = one(T1) + _two = _one + _one + for i in _zero:l-_one + x = x | (h*((m-_one) & _two^(dim*i))÷_two^((dim-_one)*i)) + y = y | (h*((m-_one) & _two^(dim*i+_one))÷_two^((dim-_one)*i+_one)) + z = z | (h*((m-_one) & _two^(dim*i+_two))÷_two^((dim-_one)*i+_two)) + end + if dim < 3 + OctantBWG{2,4,T1}(l,(x,y)) + else + OctantBWG{3,8,T1}(l,(x,y,z)) + end +end + +#OctantBWG(dim::Int,l::Int,m::Int,b::Int=_maxlevel[dim-1]) = OctantBWG(dim,l,m,b) +#OctantBWG(dim::Int,l::Int,m::Int,b::Int32) = OctantBWG(dim,l,m,b) +#OctantBWG(dim::Int,l::Int32,m::Int,b::Int32) = OctantBWG(dim,l,Int32(m),b) +function OctantBWG(level::T,coords::NTuple) where T <: Integer + dim = length(coords) + nnodes = 2^dim + OctantBWG{dim,nnodes,eltype(coords)}(level,coords) +end +#OctantBWG(level::Int32,coords::NTuple) = OctantBWG(level,Int32.(coords)) +#OctantBWG(level::Int32, coords::NTuple{dim,Int32}) where dim = OctantBWG{dim,2^dim,2*dim,Int32}(level,coords) + +""" +From [BWG2011](@citet); +> The octant coordinates are stored as integers of a fixed number b of bits, +> where the highest (leftmost) bit represents the first vertical level of the +> octree (counting the root as level zero), the second highest bit the second level of the octree, and so on. +A morton index can thus be constructed by interleaving the integer bits (2D): +\$m(\\text{Oct}) := (y_b,x_b,y_{b-1},x_{b-1},...y_0,x_0)_2\$ +further we assume the following +> Due to the two-complement representation of integers in practically all current hardware, +> where the highest digit denotes the negated appropriate power of two, bitwise operations as used, +> for example, in Algorithm 1 yield the correct result even for negative coordinates. +also from [BWG2011](@citet) + +TODO: use LUT method from https://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/ +""" +function morton(octant::OctantBWG{dim,N,T},l::T,b::T) where {dim,N,T<:Integer} + o = one(T) + z = zero(T) + id = zero(widen(eltype(octant.xyz))) + loop_length = (sizeof(typeof(id))*T(8)) ÷ dim - o + for i in z:loop_length + for d in z:dim-o + # first shift extract i-th bit and second shift inserts it at interleaved index + id = id | ((octant.xyz[d+o] & (o << i)) << ((dim-o)*i+d)) + end + end + # discard the bit information about deeper levels + return (id >> ((b-l)*dim))+o +end +morton(octant::OctantBWG{dim,N,T1},l::T2,b::T3) where {dim,N,T1<:Integer,T2<:Integer,T3<:Integer} = morton(octant,T1(l),T1(b)) + +Base.zero(::Type{OctantBWG{3, 8}}) = OctantBWG(3, 0, 1) +Base.zero(::Type{OctantBWG{2, 4}}) = OctantBWG(2, 0, 1) +root(dim::T) where T<:Integer = zero(OctantBWG{dim,2^dim}) +Base.eltype(::Type{OctantBWG{dim,N,T}}) where {dim,N,T} = T + +ncorners(::Type{OctantBWG{dim,N,T}}) where {dim,N,T} = N # TODO change to how many corners +ncorners(o::OctantBWG) = ncorners(typeof(o)) +nnodes(::Type{OctantBWG{dim,N,T}}) where {dim,N,T} = N +nnodes(o::OctantBWG) = ncorners(typeof(o)) +nchilds(::Type{OctantBWG{dim,N,T}}) where {dim,N,T} = N +nchilds(o::OctantBWG) = nchilds(typeof(o))# Follow z order, x before y before z for faces, edges and corners + +Base.isequal(o1::OctantBWG, o2::OctantBWG) = (o1.l == o2.l) && (o1.xyz == o2.xyz) +""" + o1::OctantBWG < o2::OctantBWG +Implements Algorithm 2.1 of [IBWG2015](@citet). +Checks first if mortonid is smaller and later if level is smaller. +Thus, ancestors precede descendants (preordering). +""" +function Base.isless(o1::OctantBWG, o2::OctantBWG) + if o1.xyz != o2.xyz + #TODO verify b=o1.l/b=o2.l as argument potential bug otherwise + return morton(o1,o1.l,o1.l) < morton(o2,o2.l,o2.l) + else + return o1.l < o2.l + end +end + +""" + children(octant::OctantBWG{dim,N,T}, b::Integer) -> NTuple{M,OctantBWG} +Computes all childern of `octant` +""" +function children(octant::OctantBWG{dim,N,T}, b::Integer) where {dim,N,T} + o = one(T) + _nchilds = nchilds(octant) + startid = morton(octant,octant.l+o,b) + endid = startid + _nchilds + o + return ntuple(i->OctantBWG(dim,octant.l+o,(startid:endid)[i],b),_nchilds) +end + +abstract type OctantIndex{T<:Integer} end +Base.isequal(i1::T,i2::T) where T<:OctantIndex = i1.idx == i2.idx #same type +Base.isequal(i1::T1,i2::T2) where {T1<:OctantIndex,T2<:OctantIndex} = false #different type + +struct OctantCornerIndex{T} <: OctantIndex{T} + idx::T +end +Base.hash(idx::OctantCornerIndex) = Base.hash((0,idx.idx)) +Base.show(io::IO, ::MIME"text/plain", c::OctantCornerIndex) = print(io, "O-Corner $(c.idx)") +Base.show(io::IO, c::OctantCornerIndex) = print(io, "O-Corner $(c.idx)") + +struct OctantEdgeIndex{T} <: OctantIndex{T} + idx::T +end +Base.hash(idx::OctantEdgeIndex) = Base.hash((1,idx.idx)) +Base.show(io::IO, ::MIME"text/plain", e::OctantEdgeIndex) = print(io, "O-Edge $(e.idx)") +Base.show(io::IO, e::OctantEdgeIndex) = print(io, "O-Edge $(e.idx)") + +struct OctantFaceIndex{T} <: OctantIndex{T} + idx::T +end +Base.hash(idx::OctantFaceIndex) = Base.hash((2,idx.idx)) +Base.show(io::IO, ::MIME"text/plain", f::OctantFaceIndex) = print(io, "O-Face $(f.idx)") +Base.show(io::IO, f::OctantFaceIndex) = print(io, "O-Face $(f.idx)") + +vertex(octant::OctantBWG, c::OctantCornerIndex, b::Integer) = vertex(octant,c.idx,b) +function vertex(octant::OctantBWG{dim,N,T}, c::Integer, b::Integer) where {dim,N,T} + h = T(_compute_size(b,octant.l)) + return ntuple(d->((c-1) & (2^(d-1))) == 0 ? octant.xyz[d] : octant.xyz[d] + h ,dim) +end + +""" + vertices(octant::OctantBWG{dim}, b::Integer) + +Computes all vertices of a given `octant`. Each vertex is encoded within the octree coordinates i.e. by integers. +""" +function vertices(octant::OctantBWG{dim},b::Integer) where {dim} + _nvertices = 2^dim + return ntuple(i->vertex(octant,i,b),_nvertices) +end + +face(octant::OctantBWG, f::OctantFaceIndex, b::Integer) = face(octant,f.idx,b) +function face(octant::OctantBWG{2}, f::Integer, b::Integer) + cornerid = view(𝒱₂,f,:) + return ntuple(i->vertex(octant, cornerid[i], b),2) +end + +function face(octant::OctantBWG{3}, f::Integer, b::Integer) + cornerid = view(𝒱₃,f,:) + return ntuple(i->vertex(octant, cornerid[i], b),4) +end + +""" + faces(octant::OctantBWG{dim}, b::Integer) + +Computes all faces of a given `octant`. Each face is encoded within the octree coordinates i.e. by integers. +Further, each face consists of either two two-dimensional integer coordinates or four three-dimensional integer coordinates. +""" +function faces(octant::OctantBWG{dim}, b::Integer) where dim + _nfaces = 2*dim + return ntuple(i->face(octant,i,b),_nfaces) +end + +edge(octant::OctantBWG, e::OctantEdgeIndex, b::Integer) = edge(octant,e.idx,b) +function edge(octant::OctantBWG{3}, e::Integer, b::Integer) + cornerid = view(𝒰,e,:) + return ntuple(i->vertex(octant,cornerid[i], b),2) +end + +""" + edges(octant::OctantBWG{dim}, b::Integer) + +Computes all edges of a given `octant`. Each edge is encoded within the octree coordinates i.e. by integers. +Further, each edge consists of two three-dimensional integer coordinates. +""" +edges(octant::OctantBWG{3}, b::Integer) = ntuple(i->edge(octant,i,b),12) + +""" + boundaryset(o::OctantBWG{2}, i::Integer, b::Integer +implements two dimensional boundaryset table from Fig.4.1 [IBWG2015](@citet) +TODO: could be done little bit less ugly +""" +function boundaryset(o::OctantBWG{2,N,T}, i::Integer, b::Integer) where {N,T} + if i==1 + return Set((OctantCornerIndex(1),OctantFaceIndex(1),OctantFaceIndex(3))) + elseif i==2 + return Set((OctantCornerIndex(2),OctantFaceIndex(2),OctantFaceIndex(3))) + elseif i==3 + return Set((OctantCornerIndex(3),OctantFaceIndex(1),OctantFaceIndex(4))) + elseif i==4 + return Set((OctantCornerIndex(4),OctantFaceIndex(2),OctantFaceIndex(4))) + else + throw("no boundary") + end +end + +""" + boundaryset(o::OctantBWG{3}, i::Integer, b::Integer +implements three dimensional boundaryset table from Fig.4.1 [IBWG2015](@citet) +TODO: could be done little bit less ugly +""" +function boundaryset(o::OctantBWG{3,N,T}, i::Integer, b::Integer) where {N,T} + if i==1 + return Set((OctantCornerIndex(1),OctantEdgeIndex(1),OctantEdgeIndex(5),OctantEdgeIndex(9), OctantFaceIndex(1),OctantFaceIndex(3),OctantFaceIndex(5))) + elseif i==2 + return Set((OctantCornerIndex(2),OctantEdgeIndex(1),OctantEdgeIndex(6),OctantEdgeIndex(10),OctantFaceIndex(2),OctantFaceIndex(3),OctantFaceIndex(5))) + elseif i==3 + return Set((OctantCornerIndex(3),OctantEdgeIndex(2),OctantEdgeIndex(5),OctantEdgeIndex(11),OctantFaceIndex(1),OctantFaceIndex(4),OctantFaceIndex(5))) + elseif i==4 + return Set((OctantCornerIndex(4),OctantEdgeIndex(2),OctantEdgeIndex(6),OctantEdgeIndex(12),OctantFaceIndex(2),OctantFaceIndex(4),OctantFaceIndex(5))) + elseif i==5 + return Set((OctantCornerIndex(5),OctantEdgeIndex(3),OctantEdgeIndex(7),OctantEdgeIndex(9), OctantFaceIndex(1),OctantFaceIndex(3),OctantFaceIndex(6))) + elseif i==6 + return Set((OctantCornerIndex(6),OctantEdgeIndex(3),OctantEdgeIndex(8),OctantEdgeIndex(10),OctantFaceIndex(2),OctantFaceIndex(3),OctantFaceIndex(6))) + elseif i==7 + return Set((OctantCornerIndex(7),OctantEdgeIndex(4),OctantEdgeIndex(7),OctantEdgeIndex(11),OctantFaceIndex(1),OctantFaceIndex(4),OctantFaceIndex(6))) + elseif i==8 + return Set((OctantCornerIndex(8),OctantEdgeIndex(4),OctantEdgeIndex(8),OctantEdgeIndex(12),OctantFaceIndex(2),OctantFaceIndex(4),OctantFaceIndex(6))) + else + throw("no boundary") + end +end + +""" + find_range_boundaries(f::OctantBWG{dim,N,T}, l::OctantBWG{dim,N,T}, s::OctantBWG{dim,N,T}, idxset, b) + find_range_boundaries(s::OctantBWG{dim,N,T}, idxset, b) +Algorithm 4.2 of [IBWG2015](@citet) +TODO: write tests +""" +function find_range_boundaries(f::OctantBWG{dim,N,T1}, l::OctantBWG{dim,N,T1}, s::OctantBWG{dim,N,T1}, idxset::Set{OctantIndex{T2}}, b) where {dim,N,T1,T2} + o = one(T1) + if isempty(idxset) || s.l == b + return idxset + end + j = ancestor_id(f,s.l+o,b); k = ancestor_id(l,s.l+o,b) + boundary_j = boundaryset(s,j,b) + kidz = children(s,b) + if j==k + return find_range_boundaries(f,l,kidz[j],idxset ∩ boundary_j,b) + end + idxset_match = Set{OctantIndex{T2}}() + for i in (j+o):(k-o) + union!(idxset_match,idxset ∩ boundaryset(s,i,b)) + end + boundary_k = boundaryset(s,k,b) + idxset_match_j = setdiff((idxset ∩ boundary_j),idxset_match) + fj, lj = descendants(kidz[j],b) + if fj != f + idxset_match_j = find_range_boundaries(f,lj,kidz[j],idxset_match_j,b) + end + idxset_match_k = setdiff(setdiff((idxset ∩ boundary_k),idxset_match),idxset_match_j) + fk, lk = descendants(kidz[k],b) + if lk != l + idxset_match_k = find_range_boundaries(fk,l,kidz[k],idxset_match_k,b) + end + return idxset_match ∪ idxset_match_j ∪ idxset_match_k +end + +#for convenience, should probably changed to parent(s) until parent(s)==root and then descendants(root) +function find_range_boundaries(s::OctantBWG, idxset, b) + f,l = descendants(s,b) + return find_range_boundaries(f,l,s,idxset,b) +end + +function isrelevant(xyz::NTuple{dim,T},leafsuppₚ::Set{<:OctantBWG}) where {dim,T} + ###### only relevant for distributed + #for all s in leafsuppₚ + # if s in 𝒪ₚ + # return true + # else + # check stuff Algorithm 5.1 line 4-5 + # end + #end + return true +end + +struct OctreeBWG{dim,N,T} <: AbstractAdaptiveCell{Ferrite.RefHypercube{dim}} + leaves::Vector{OctantBWG{dim,N,T}} + #maximum refinement level + b::T + nodes::NTuple{N,Int} +end + +function refine!(octree::OctreeBWG{dim,N,T}, pivot_octant::OctantBWG{dim,N,T}) where {dim,N,T<:Integer} + if !(pivot_octant.l + 1 <= octree.b) + return + end + o = one(T) + # TODO replace this with recursive search function + leave_idx = findfirst(x->x==pivot_octant,octree.leaves) + old_octant = popat!(octree.leaves,leave_idx) + _children = children(pivot_octant,octree.b) + for child in _children + insert!(octree.leaves,leave_idx,child) + leave_idx += 1 + end +end + +function coarsen!(octree::OctreeBWG{dim,N,T}, o::OctantBWG{dim,N,T}) where {dim,N,T<:Integer} + _two = T(2) + leave_idx = findfirst(x->x==o,octree.leaves) + shift = child_id(o,octree.b) - one(T) + if shift != zero(T) + old_morton = morton(o,o.l,octree.b) + o = OctantBWG(dim,o.l,old_morton,octree.b) + end + window_start = leave_idx - shift + window_length = _two^dim - one(T) + new_octant = parent(o, octree.b) + octree.leaves[leave_idx - shift] = new_octant + deleteat!(octree.leaves,leave_idx-shift+one(T):leave_idx-shift+window_length) +end + +OctreeBWG{3,8}(nodes::NTuple,b=_maxlevel[2]) = OctreeBWG{3,8,Int64}([zero(OctantBWG{3,8})],Int64(b),nodes) +OctreeBWG{2,4}(nodes::NTuple,b=_maxlevel[1]) = OctreeBWG{2,4,Int64}([zero(OctantBWG{2,4})],Int64(b),nodes) +OctreeBWG(cell::Quadrilateral,b=_maxlevel[2]) = OctreeBWG{2,4}(cell.nodes,b) +OctreeBWG(cell::Hexahedron,b=_maxlevel[1]) = OctreeBWG{3,8}(cell.nodes,b) + +Base.length(tree::OctreeBWG) = length(tree.leaves) +Base.eltype(::Type{OctreeBWG{dim,N,T}}) where {dim,N,T} = T + +function inside(oct::OctantBWG{dim},b) where dim + maxsize = _maximum_size(b) + outside = any(xyz -> xyz >= maxsize, oct.xyz) || any(xyz -> xyz < 0, oct.xyz) + return !outside +end + +inside(tree::OctreeBWG{dim},oct::OctantBWG{dim}) where dim = inside(oct,tree.b) + +""" + split_array(octree::OctreeBWG, a::OctantBWG) + split_array(octantarray, a::OctantBWG, b::Integer) +Algorithm 3.3 of [IBWG2015](@citet). Efficient binary search. +""" +function split_array(octantarray, a::OctantBWG{dim,N,T}, b::Integer) where {dim,N,T} + o = one(T) + 𝐤 = T[i==1 ? 1 : length(octantarray)+1 for i in 1:2^dim+1] + for i in 2:2^dim + m = 𝐤[i-1] + while m < 𝐤[i] + n = m + (𝐤[i] - m)÷2 + c = ancestor_id(octantarray[n], a.l+o, b) + if c < i + m = n+1 + else + for j in i:c + 𝐤[j] = n + end + end + end + end + #TODO non-allocating way? + return ntuple(i->view(octantarray,𝐤[i]:𝐤[i+1]-1),2^dim) +end + +split_array(tree::OctreeBWG, a::OctantBWG) = split_array(tree.leaves, a, tree.b) + +function search(octantarray, a::OctantBWG{dim,N,T1}, idxset::Vector{T2}, b::Integer, Match=match) where {dim,N,T1<:Integer,T2} + isempty(octantarray) && return + isleaf = (length(octantarray) == 1 && a ∈ octantarray) ? true : false + idxset_match = eltype(idxset)[] + for q in idxset + if Match(a,isleaf,q,b) + push!(idxset_match,q) + end + end + if isempty(idxset_match) && !isleaf + 𝐇 = split_array(octantarray,a,b) + _children = children(a,b) + for (child,h) in zip(_children,𝐇) + search(h,child,idxset_match,b) + end + end + return idxset_match +end + +search(tree::OctreeBWG, a::OctantBWG, idxset, Match=match) = search(tree.leaves, a, idxset, tree.b, match) + +""" + match(o::OctantBWG, isleaf::Bool, q) +from [IBWG2015](@citet) +> match returns true if there is a leaf r ∈ 𝒪 that is a descendant of o +> such that match_q(r) = true, and is allowed to return a false positive +> (i.e., true even if match_q(r) = false for all descendants leaves of o) +> if isleaf=true, then the return value of match is irrelevant +I don't understand what of a to check against index q +""" +function match(o::OctantBWG, isleaf::Bool, q, b) + isleaf && (return true) + println(q) + println(o) + return false +end + +""" + ForestBWG{dim, C<:AbstractAdaptiveCell, T<:Real} <: AbstractAdaptiveGrid{dim} +`p4est` adaptive grid implementation based on [BWG2011](@citet) +and [IBWG2015](@citet). + +## Constructor + ForestBWG(grid::AbstractGrid{dim}, b=_maxlevel[dim-1]) where dim +Builds an adaptive grid based on a non-adaptive one `grid` and a given max refinement level `b`. +""" +struct ForestBWG{dim, C<:OctreeBWG, T<:Real} <: AbstractAdaptiveGrid{dim} + cells::Vector{C} + nodes::Vector{Node{dim,T}} + # Sets + cellsets::Dict{String,OrderedSet{Int}} + nodesets::Dict{String,OrderedSet{Int}} + facetsets::Dict{String,OrderedSet{Ferrite.FacetIndex}} + vertexsets::Dict{String,OrderedSet{Ferrite.VertexIndex}} + #Topology + topology::ExclusiveTopology +end + +function ForestBWG(grid::Ferrite.AbstractGrid{dim},b=_maxlevel[dim-1]) where dim + cells = getcells(grid) + C = eltype(cells) + @assert isconcretetype(C) + @assert (C == Quadrilateral && dim == 2) || (C == Hexahedron && dim == 3) + topology = ExclusiveTopology(cells) + cells = OctreeBWG.(grid.cells,b) + nodes = getnodes(grid) + cellsets = Ferrite.getcellsets(grid) + nodesets = Ferrite.getnodesets(grid) + facetsets = Ferrite.getfacetsets(grid) + vertexsets = Ferrite.getvertexsets(grid) + return ForestBWG(cells,nodes,cellsets,nodesets,facetsets,vertexsets,topology) +end + +function Ferrite.get_facet_facet_neighborhood(g::ForestBWG{dim}) where dim + return Ferrite._get_facet_facet_neighborhood(g.topology,Val(dim)) +end + +function refine_all!(forest::ForestBWG,l) + for tree in forest.cells + for leaf in tree.leaves + if leaf.l != l-1 #maxlevel + continue + else + refine!(tree,leaf) + end + end + end +end + +function refine!(forest::ForestBWG, cellid::Integer) + nleaves_k = length(forest.cells[1].leaves) + prev_nleaves_k = 0 + k = 1 + while nleaves_k < cellid + k += 1 + prev_nleaves_k = nleaves_k + nleaves_k += length(forest.cells[k].leaves) + end + refine!(forest.cells[k],forest.cells[k].leaves[cellid-prev_nleaves_k]) +end + +function refine!(forest::ForestBWG, cellids::Vector{<:Integer}) + sort!(cellids) + ncells = getncells(forest) + shift = 0 + for cellid in cellids + refine!(forest,cellid+shift) + shift += getncells(forest) - ncells + ncells = getncells(forest) + end +end + +function coarsen_all!(forest::ForestBWG) + for tree in forest.cells + for leaf in tree.leaves + if child_id(leaf,tree.b) == 1 + coarsen!(tree,leaf) + end + end + end +end + +Ferrite.getneighborhood(forest::ForestBWG,idx) = getneighborhood(forest.topology,forest,idx) + +function Ferrite.getncells(grid::ForestBWG) + numcells = 0 + for tree in grid.cells + numcells += length(tree) + end + return numcells +end + +function Ferrite.getcells(forest::ForestBWG{dim,C}) where {dim,C} + treetype = C + ncells = getncells(forest) + nnodes = 2^dim + cellvector = Vector{OctantBWG{dim,nnodes,eltype(C)}}(undef,ncells) + o = one(eltype(C)) + cellid = o + for tree in forest.cells + for leaf in tree.leaves + cellvector[cellid] = leaf + cellid += o + end + end + return cellvector +end + +function Ferrite.getcells(forest::ForestBWG{dim}, cellid::Int) where dim + @warn "Slow dispatch, consider to call `getcells(forest)` once instead" maxlog=1 #TODO doc page for performance + #TODO should nleaves be saved by forest? + nleaves = length.(forest.cells) # cells=trees + #TODO remove that later by for loop or [IBWG2015](@citet) iterator approach + nleaves_cumsum = cumsum(nleaves) + k = findfirst(x->cellid<=x,nleaves_cumsum) + #TODO is this actually correct? + leafid = k == 1 ? cellid : cellid - (nleaves_cumsum[k] - nleaves[k]) + return forest.cells[k].leaves[leafid] +end + +Ferrite.getcelltype(grid::ForestBWG) = eltype(grid.cells) +Ferrite.getcelltype(grid::ForestBWG, i::Int) = eltype(grid.cells) # assume for now same cell type TODO + +""" + transform_pointBWG(forest, vertices) -> Vector{Vec{dim}} + transform_pointBWG(forest::ForestBWG{dim}, k::Integer, vertex::NTuple{dim,T}) where {dim,T} -> Vec{dim} + +Transformation of a octree coordinate system point `vertex` (or a collection `vertices`) to the corresponding physical coordinate system. +""" +function transform_pointBWG(forest::ForestBWG{dim}, k::Integer, vertex::NTuple{dim,T}) where {dim,T} + tree = forest.cells[k] + cellnodes = getnodes(forest,collect(tree.nodes)) .|> get_node_coordinate + vertex = vertex .* (2/(2^tree.b)) .- 1 + octant_physical_coordinates = sum(j-> cellnodes[j] * Ferrite.shape_value(Lagrange{Ferrite.RefHypercube{dim},1}(),Vec{dim}(vertex),j),1:length(cellnodes)) + return Vec{dim}(octant_physical_coordinates) +end + +transform_pointBWG(forest, vertices) = transform_pointBWG.((forest,), first.(vertices), last.(vertices)) + +""" + rotation_permutation(::Val{2},r,i) -> i′ +computes based on the rotation indicator `r` ∈ {0,1} and a given corner index `i` ∈ {1,2} the permuted corner index `i′` +""" +function rotation_permutation(r,i) + i′ = r == 0 ? i : 3-i + return i′ +end + +""" + rotation_permutation(f,f′,r,i) -> i′ +computes based on the rotation indicator `r` ∈ {0,...,3} and a given corner index `i` ∈ {1,...,4} the permuted corner index `i′` +See Table 3 and Theorem 2.2 [BWG2011](@citet). +""" +function rotation_permutation(f,f′,r,i) + return 𝒫[𝒬[ℛ[f,f′],r+1],i] +end + +p4est_opposite_face_index(f) = ((f - 1) ⊻ 0b1) + 1 +p4est_opposite_edge_index(e) = ((e - 1) ⊻ 0b11) + 1 + +#TODO: this function should wrap the LNodes Iterator of [IBWG2015](@citet) +""" + creategrid(forest::ForestBWG) -> NonConformingGrid +Materializes a p4est forest into a NonConformingGrid that can be used as usual to solve a finite element problem. +""" +function creategrid(forest::ForestBWG{dim,C,T}) where {dim,C,T} + nodes = Vector{Tuple{Int,NTuple{dim,Int32}}}() + sizehint!(nodes,getncells(forest)*2^dim) + _perm = dim == 2 ? 𝒱₂_perm : 𝒱₃_perm + _perminv = dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv + node_map = dim < 3 ? node_map₂ : node_map₃ + node_map_inv = dim < 3 ? node_map₂_inv : node_map₃_inv + nodeids = Dict{Tuple{Int,NTuple{dim,Int32}},Int}() + nodeowners = Dict{Tuple{Int,NTuple{dim,Int32}},Tuple{Int,NTuple{dim,Int32}}}() + facet_neighborhood = Ferrite.Ferrite.get_facet_facet_neighborhood(forest) + + # Phase 1: Assign node owners intra-octree + pivot_nodeid = 1 + for (k,tree) in enumerate(forest.cells) + for leaf in tree.leaves + _vertices = vertices(leaf,tree.b) + for v in _vertices + push!(nodes,(k,v)) + nodeids[(k,v)] = pivot_nodeid + pivot_nodeid += 1 + nodeowners[(k,v)] = (k,v) + end + end + end + + # Phase 2: Assign node owners inter-octree + for (k,tree) in enumerate(forest.cells) + _vertices = vertices(root(dim),tree.b) + # Vertex neighbors + @debug println("Setting vertex neighbors for octree $k") + for (v,vc) in enumerate(_vertices) + vertex_neighbor = forest.topology.vertex_vertex_neighbor[k,node_map[v]] + for (k′, v′) in vertex_neighbor + @debug println(" pair $v $v′") + if k > k′ + #delete!(nodes,(k,v)) + new_v = vertex(root(dim),node_map[v′],tree.b) + nodeids[(k,vc)] = nodeids[(k′,new_v)] + nodeowners[(k,vc)] = (k′,new_v) + @debug println(" Matching $vc (local) to $new_v (neighbor)") + end + end + # TODO check if we need to also update the face neighbors + end + if dim > 1 + _faces = faces(root(dim),tree.b) + # Face neighbors + @debug println("Updating face neighbors for octree $k") + for (f,fc) in enumerate(_faces) # f in p4est notation + # Skip boundary edges + facet_neighbor_ = facet_neighborhood[k,_perm[f]] + if length(facet_neighbor_) == 0 + continue + end + @debug @assert length(facet_neighbor_) == 1 + k′, f′_ferrite = facet_neighbor_[1] + f′ = _perminv[f′_ferrite] + @debug println(" Neighboring tree: $k′, face $f′_ferrite (Ferrite)/$f′ (p4est)") + if k > k′ # Owner + tree′ = forest.cells[k′] + for leaf in tree.leaves + fnodes = face(leaf, f , tree.b) + if !contains_facet(fc, fnodes) + @debug println(" Rejecting leaf $leaf because its facet $fnodes is not on the octant boundary") + continue + end + neighbor_candidate = transform_facet(forest,k′,f′,leaf) + # Candidate must be the face opposite to f' + f′candidate = p4est_opposite_face_index(f′) + fnodes_neighbor = face(neighbor_candidate, f′candidate, tree′.b) + r = compute_face_orientation(forest,k,f) + @debug println(" Trying to match $fnodes (local) to $fnodes_neighbor (neighbor $neighbor_candidate)") + if dim == 2 + for i ∈ 1:ncorners_face2D + i′ = rotation_permutation(r,i) + if haskey(nodeids, (k′,fnodes_neighbor[i′])) + @debug println(" Updating $((k,fnodes[i])) $(nodeids[(k,fnodes[i])]) -> $(nodeids[(k′,fnodes_neighbor[i′])])") + nodeids[(k,fnodes[i])] = nodeids[(k′,fnodes_neighbor[i′])] + nodeowners[(k,fnodes[i])] = (k′,fnodes_neighbor[i′]) + end + end + else + for i ∈ 1:ncorners_face3D + rotated_ξ = rotation_permutation(f′,f,r,i) + if haskey(nodeids, (k′,fnodes_neighbor[i])) + @debug println(" Updating $((k,fnodes[i])) $(nodeids[(k,fnodes[rotated_ξ])]) -> $(nodeids[(k′,fnodes_neighbor[i])])") + nodeids[(k,fnodes[rotated_ξ])] = nodeids[(k′,fnodes_neighbor[i])] + nodeowners[(k,fnodes[rotated_ξ])] = (k′,fnodes_neighbor[i]) + end + end + end + end + end + end + end + if dim > 2 + # edge neighbors + @debug println("Updating edge neighbors for octree $k") + for (e,ec) in enumerate(edges(root(dim),tree.b)) # e in p4est notation + # Skip boundary edges + edge_neighbor_ = forest.topology.edge_edge_neighbor[k,edge_perm[e]] + if length(edge_neighbor_) == 0 + continue + end + @debug @assert length(edge_neighbor_) == 1 + k′, e′_ferrite = edge_neighbor_[1] + e′ = edge_perm_inv[e′_ferrite] + @debug println(" Neighboring tree: $k′, edge $e′_ferrite (Ferrite)/$e′ (p4est)") + if k > k′ # Owner + tree′ = forest.cells[k′] + for leaf in tree.leaves + # First we skip edges which are not on the current edge of the root element + enodes = edge(leaf, e , tree.b) + if !contains_edge(ec, enodes) + @debug println(" Rejecting leaf $leaf because its edge $enodes is not on the octant boundary") + continue + end + neighbor_candidate = transform_edge(forest,k′,e′,leaf, false) + # Candidate must be the edge opposite to e' + e′candidate = p4est_opposite_edge_index(e′) + + enodes_neighbor = edge(neighbor_candidate, e′candidate, tree′.b) + r = compute_edge_orientation(forest,k,e) + @debug println(" Trying to match $enodes (local) to $enodes_neighbor (neighbor $neighbor_candidate)") + for i ∈ 1:ncorners_edge + i′ = rotation_permutation(r,i) + if haskey(nodeids, (k′,enodes_neighbor[i′])) + @debug println(" Updating $((k,enodes[i])) $(nodeids[(k,enodes[i])]) -> $(nodeids[(k′,enodes_neighbor[i′])])") + nodeids[(k,enodes[i])] = nodeids[(k′,enodes_neighbor[i′])] + nodeowners[(k,enodes[i])] = (k′,enodes_neighbor[i′]) + end + end + end + end + end + end + end + + # Phase 3: Compute unique physical nodes + nodeids_dedup = Dict{Int,Int}() + next_nodeid = 1 + for (kv,nodeid) in nodeids + if !haskey(nodeids_dedup, nodeid) + nodeids_dedup[nodeid] = next_nodeid + next_nodeid += 1 + end + end + nodes_physical_all = transform_pointBWG(forest,nodes) + nodes_physical = zeros(eltype(nodes_physical_all), next_nodeid-1) + for (ni, (kv,nodeid)) in enumerate(nodeids) + nodes_physical[nodeids_dedup[nodeid]] = nodes_physical_all[nodeid] + end + + # Phase 4: Generate cells + celltype = dim < 3 ? Quadrilateral : Hexahedron + cells = celltype[] + cellnodes = zeros(Int,2^dim) + for (k,tree) in enumerate(forest.cells) + for leaf in tree.leaves + _vertices = vertices(leaf,tree.b) + cellnodes = ntuple(i-> nodeids_dedup[nodeids[nodeowners[(k,_vertices[i])]]],length(_vertices)) + push!(cells,celltype(ntuple(i->cellnodes[node_map[i]],length(cellnodes)))) + end + end + + # Phase 5: Generate grid and haning nodes + facetsets = reconstruct_facetsets(forest) #TODO edge, node and cellsets + #facesets_ordered = OrderedSet(facetsets) + hnodes = hangingnodes(forest, nodeids, nodeowners) + hnodes_dedup = Dict{Int64, Vector{Int64}}() + for (constrained,constainers) in hnodes + hnodes_dedup[nodeids_dedup[constrained]] = [nodeids_dedup[constainer] for constainer in constainers] + end + return NonConformingGrid(cells, nodes_physical .|> Node, facetsets=facetsets, conformity_info=hnodes_dedup) +end + +function reconstruct_facetsets(forest::ForestBWG{dim}) where dim + _perm = dim == 2 ? 𝒱₂_perm : 𝒱₃_perm + _perm_inv = dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv + new_facesets = typeof(forest.facetsets)() + for (facetsetname, facetset) in forest.facetsets + new_facetset = typeof(facetset)() + for facetidx in facetset + pivot_tree = forest.cells[facetidx[1]] + last_cellid = facetidx[1] != 1 ? sum(length,@view(forest.cells[1:(facetidx[1]-1)])) : 0 + pivot_faceid = facetidx[2] + pivot_face = faces(root(dim),pivot_tree.b)[_perm_inv[pivot_faceid]] + for (leaf_idx,leaf) in enumerate(pivot_tree.leaves) + for (leaf_face_idx,leaf_face) in enumerate(faces(leaf,pivot_tree.b)) + if contains_facet(pivot_face,leaf_face) + ferrite_leaf_face_idx = _perm[leaf_face_idx] + push!(new_facetset,FacetIndex(last_cellid+leaf_idx,ferrite_leaf_face_idx)) + end + end + end + end + new_facesets[facetsetname] = new_facetset + end + return new_facesets +end + +""" + hangingnodes(forest,nodeids,nodeowners) +Constructs a map from constrained nodeids to the ones that constraint +""" +function hangingnodes(forest::ForestBWG{dim}, nodeids, nodeowners) where dim + _perm = dim == 2 ? 𝒱₂_perm : 𝒱₃_perm + _perminv = dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv + facetable = dim == 2 ? 𝒱₂ : 𝒱₃ + opposite_face = dim == 2 ? opposite_face_2 : opposite_face_3 + hnodes = Dict{Int,Vector{Int}}() + facet_neighborhood = Ferrite.Ferrite.get_facet_facet_neighborhood(forest) + for (k,tree) in enumerate(forest.cells) + rootfaces = faces(root(dim),tree.b) + for (l,leaf) in enumerate(tree.leaves) + c̃ = child_id(leaf,tree.b) + if leaf == root(dim) + continue + end + for (ci,c) in enumerate(vertices(leaf,tree.b)) + parent_ = parent(leaf,tree.b) + parentfaces = faces(parent_,tree.b) + for (pface_i, pface) in enumerate(parentfaces) + if iscenter(c,pface) #hanging node candidate + neighbor_candidate = facet_neighbor(parent_, pface_i, tree.b) + if inside(tree,neighbor_candidate) #intraoctree branch + neighbor_candidate_idx = findfirst(x->x==neighbor_candidate,tree.leaves) + if neighbor_candidate_idx !== nothing + neighbor_candidate_faces = faces(neighbor_candidate,tree.b) + nf = findfirst(x->x==pface,neighbor_candidate_faces) + hnodes[nodeids[nodeowners[(k,c)]]] = [nodeids[nodeowners[(k,nc)]] for nc in neighbor_candidate_faces[nf]] + if dim > 2 + vs = vertices(leaf,tree.b) + for ξ ∈ 1:ncorners_face3D + c′ = facetable[pface_i, ξ] + if c′ ∉ (c̃,ci) + neighbor_candidate_edges = edges(neighbor_candidate,tree.b) + ne = findfirst(x->iscenter(vs[c′],x),neighbor_candidate_edges) + if ne !== nothing + hnodes[nodeids[nodeowners[(k,vs[c′])]]] = [nodeids[nodeowners[(k,ne)]] for ne in neighbor_candidate_edges[ne]] + end + end + end + end + break + end + else #interoctree branch + for (ri,rf) in enumerate(rootfaces) + facet_neighbor_ = facet_neighborhood[k,_perm[ri]] + if length(facet_neighbor_) == 0 + continue + end + if contains_facet(rf, pface) + k′ = facet_neighbor_[1][1] + ri′ = _perminv[facet_neighbor_[1][2]] + interoctree_neighbor = transform_facet(forest, k′, ri′, neighbor_candidate) + interoctree_neighbor_candidate_idx = findfirst(x->x==interoctree_neighbor,forest.cells[k′].leaves) + if interoctree_neighbor_candidate_idx !== nothing + r = compute_face_orientation(forest,k,pface_i) + neighbor_candidate_faces = faces(neighbor_candidate,forest.cells[k′].b) + transformed_neighbor_faces = faces(interoctree_neighbor,forest.cells[k′].b) + + fnodes = transformed_neighbor_faces[ri′] + vs = vertices(leaf,tree.b) + if dim > 2 + rotated_ξ = ntuple(i->rotation_permutation(ri′,ri,r,i),ncorners_face3D) + else + rotated_ξ = ntuple(i->rotation_permutation(r,i),ncorners_face2D) + end + hnodes[nodeids[nodeowners[(k,c)]]] = [nodeids[nodeowners[(k′,fnodes[ξ])]] for ξ in rotated_ξ] + + if dim > 2 + for ξ in rotated_ξ + c′ = facetable[pface_i, ξ] + if c′ ∉ (c̃,c) + neighbor_candidate_edges = edges(interoctree_neighbor,tree.b) + ne = findfirst(x->iscenter(vs[c′],x),neighbor_candidate_edges) + if ne !== nothing + hnodes[nodeids[nodeowners[(k,vs[c′])]]] = [nodeids[nodeowners[(k,ne)]] for ne in neighbor_candidate_edges[ne]] + end + end + end + end + break + end + end + end + end + end + end + if dim > 2 + parentedges = edges(parent_,tree.b) + for (pedge_i, pedge) in enumerate(parentedges) + if iscenter(c,pedge) #hanging node candidate + neighbor_candidate = edge_neighbor(parent_, pedge_i, tree.b) + for (ri,re) in enumerate(edges(root(dim),tree.b)) + edge_neighbor_ = forest.topology.edge_edge_neighbor[k,edge_perm[ri]] + if length(edge_neighbor_) == 0 + continue + end + if contains_edge(re, pedge) + k′ = edge_neighbor_[1][1] + ri′ = edge_perm_inv[edge_neighbor_[1][2]] + interoctree_neighbor = transform_edge(forest, k′, ri′, neighbor_candidate, true) + interoctree_neighbor_candidate_idx = findfirst(x->x==interoctree_neighbor,forest.cells[k′].leaves) + if interoctree_neighbor_candidate_idx !== nothing + neighbor_candidate_edges = edges(neighbor_candidate,forest.cells[k′].b) + transformed_neighbor_edges = edges(interoctree_neighbor,forest.cells[k′].b) + ne = findfirst(x->iscenter(c,x),neighbor_candidate_edges) + if ne !== nothing + hnodes[nodeids[nodeowners[(k,c)]]] = [nodeids[nodeowners[(k′,nc)]] for nc in transformed_neighbor_edges[ne]] + else + @error "things are messed up for hanging edge constraint interoctree at octree $k edge $(_perm[ri])" + end + break + end + end + end + end + end + end + end + end + end + return hnodes +end + +function balance_corner(forest,k′,c′,o,s) + o.l == 1 && return # no balancing needed for pivot octant level == 1 + o′ = transform_corner(forest,k′,c′,o,false) + s′ = transform_corner(forest,k′,c′,s,true) #TODO verify the bool here; I think it's correct + neighbor_tree = forest.cells[k′] + if s′ ∉ neighbor_tree.leaves && parent(s′, neighbor_tree.b) ∉ neighbor_tree.leaves + if parent(parent(s′,neighbor_tree.b),neighbor_tree.b) ∈ neighbor_tree.leaves + refine!(neighbor_tree,parent(parent(s′,neighbor_tree.b),neighbor_tree.b)) + end + end +end + +function balance_face(forest,k′,f′,o,s) + o.l == 1 && return # no balancing needed for pivot octant level == 1 + o′ = transform_facet(forest,k′,f′,o) + s′ = transform_facet(forest,k′,f′,s) + neighbor_tree = forest.cells[k′] + if s′ ∉ neighbor_tree.leaves && parent(s′, neighbor_tree.b) ∉ neighbor_tree.leaves + if parent(parent(s′,neighbor_tree.b),neighbor_tree.b) ∈ neighbor_tree.leaves + refine!(neighbor_tree,parent(parent(s′,neighbor_tree.b),neighbor_tree.b)) + end + end +end + +function balance_edge(forest,k′,e′,o,s) + o.l == 1 && return # no balancing needed for pivot octant level == 1 + o′ = transform_edge(forest,k′,e′,o,false) + s′ = transform_edge(forest,k′,e′,s,true) + neighbor_tree = forest.cells[k′] + if s′ ∉ neighbor_tree.leaves && parent(s′, neighbor_tree.b) ∉ neighbor_tree.leaves + if parent(parent(s′,neighbor_tree.b),neighbor_tree.b) ∈ neighbor_tree.leaves + refine!(neighbor_tree,parent(parent(s′,neighbor_tree.b),neighbor_tree.b)) + end + end +end +""" + balanceforest!(forest) +Algorithm 17 of [BWG2011](@citet) +""" +function balanceforest!(forest::ForestBWG{dim}) where dim + perm_face = dim == 2 ? 𝒱₂_perm : 𝒱₃_perm + perm_face_inv = dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv + perm_corner = dim == 2 ? node_map₂ : node_map₃ + perm_corner_inv = dim == 2 ? node_map₂_inv : node_map₃_inv + root_ = root(dim) + nrefcells = 0 + facet_neighborhood = Ferrite.Ferrite.get_facet_facet_neighborhood(forest) + while nrefcells - getncells(forest) != 0 + nrefcells = getncells(forest) + for k in 1:length(forest.cells) + tree = forest.cells[k] + rootfaces = faces(root_,tree.b) + rootedges = dim == 3 ? edges(root_,tree.b) : nothing #TODO change + rootvertices = vertices(root_,tree.b) + balanced = balancetree(tree) + forest.cells[k] = balanced + for (o_i, o) in enumerate(forest.cells[k].leaves) + ss = possibleneighbors(o,o.l,tree.b,;insidetree=false) + isinside = inside.(ss,(tree.b,)) + notinsideidx = findall(.! isinside) + if !isempty(notinsideidx) + # s_i encodes the type of neighborhood, since it is the index of possibleneighbors + # see the docs of this function + for s_i in notinsideidx + s = ss[s_i] + if dim == 2 # need more clever s_i encoding + if s_i <= 4 #corner neighbor, only true for 2D see possibleneighbors + cc = forest.topology.vertex_vertex_neighbor[k,perm_corner[s_i]] + participating_faces_idx = findall(x->any(x .== s_i),𝒱₂) #TODO! optimize by using inverted table + pivot_faces = faces(o,tree.b) + if isempty(cc) + # the branch below checks if we are in a newly introduced topologic tree connection + # by checking if the corner neighbor is only accesible by transforming through a face + # TODO: enable a bool that either activates or deactivates the balancing over a corner + for face_idx in participating_faces_idx + face_idx = face_idx[1] + contained = contains_facet(rootfaces[face_idx],pivot_faces[face_idx]) + if contained + fc = facet_neighborhood[k,perm_face[face_idx]] + isempty(fc) && continue + @assert length(fc) == 1 + fc = fc[1] + k′, f′ = fc[1], perm_face_inv[fc[2]] + balance_face(forest,k′,f′,o,s) + end + end + continue + else + for corner_connection in cc + !(vertex(o,s_i,tree.b) == rootvertices[s_i]) && continue + k′, c′ = corner_connection[1], perm_corner_inv[corner_connection[2]] + balance_corner(forest,k′,c′,o,s) + end + end + else # face neighbor, only true for 2D + s_i -= 4 + fc = facet_neighborhood[k,perm_face[s_i]] + isempty(fc) && continue + @assert length(fc) == 1 + fc = fc[1] + k′, f′ = fc[1], perm_face_inv[fc[2]] + balance_face(forest,k′,f′,o,s) + end + else #TODO collapse this 3D branch with more clever s_i encoding into the 2D branch + if s_i <= 8 #corner neighbor, only true for 2D see possibleneighbors + #TODO a check of new introduced corner neighbors aka corner balancing, see 2D branch + cc = forest.topology.vertex_vertex_neighbor[k,perm_corner[s_i]] + isempty(cc) && continue + for corner_connection in cc + !(vertex(o,s_i,tree.b) == rootvertices[s_i]) && continue + k′, c′ = corner_connection[1], perm_corner_inv[corner_connection[2]] + balance_corner(forest,k′,c′,o,s) + end + elseif 8 < s_i <= 14 + s_i -= 8 + fc = facet_neighborhood[k,perm_face[s_i]] + isempty(fc) && continue + @assert length(fc) == 1 + fc = fc[1] + k′, f′ = fc[1], perm_face_inv[fc[2]] + balance_face(forest,k′,f′,o,s) + else + s_i -= 14 + ec = forest.topology.edge_edge_neighbor[k,edge_perm[s_i]] + pivot_edge = edge(o,s_i,tree.b) + contained_face = findall(x->face_contains_edge(x,pivot_edge),rootfaces) + if !isempty(contained_face) && !contains_edge(rootedges[s_i],pivot_edge) #check if pivot edge in interior of rootface and not octree edge + for face_idx in contained_face + fc = facet_neighborhood[k,perm_face[face_idx]] + isempty(fc) && continue + @assert length(fc) == 1 + fc = fc[1] + k′, f′ = fc[1], perm_face_inv[fc[2]] + balance_face(forest,k′,f′,o,s) + end + continue + end + isempty(ec) && continue + for edge_connection in ec + !contains_edge(rootedges[s_i],pivot_edge) && continue + k′, e′ = edge_connection[1], edge_perm_inv[edge_connection[2]] + balance_edge(forest,k′,e′,o,s) + end + end + end + end + end + end + end + end +end + +""" +Algorithm 7 of [SSB2008](@citet) + +TODO optimise the unnecessary allocations +""" +function balancetree(tree::OctreeBWG) + if length(tree.leaves) == 1 + return tree + end + W = copy(tree.leaves); P = eltype(tree.leaves)[]; R = eltype(tree.leaves)[] + for l in tree.b:-1:1 #TODO verify to do this until level 1 + Q = [o for o in W if o.l == l] + sort!(Q) + #construct T + T = eltype(Q)[] + for x in Q + if isempty(T) + push!(T,x) + continue + end + p = parent(x,tree.b) + if p ∉ parent.(T,(tree.b,)) + push!(T,x) + end + end + for t in T + push!(R,t,siblings(t,tree.b)...) + push!(P,possibleneighbors(parent(t,tree.b),l-1,tree.b)...) + end + append!(P,x for x in W if x.l == l-1) + filter!(x->!(x.l == l-1), W) #don't know why I have to negotiate like this, otherwise behaves weird + unique!(P) + append!(W,P) + empty!(P) + end + sort!(R) # be careful with sort, by=morton doesn't work due to ambuigity at max depth level + linearise!(R,tree.b) + return OctreeBWG(R,tree.b,tree.nodes) +end + +""" +Algorithm 8 of [SSB2008](@citet) + +Inverted the algorithm to delete! instead of add incrementally to a new array +""" +function linearise!(leaves::Vector{T},b) where T<:OctantBWG + inds = [i for i in 1:length(leaves)-1 if isancestor(leaves[i],leaves[i+1],b)] + deleteat!(leaves,inds) +end + +function siblings(o::OctantBWG,b;include_self=false) + siblings = children(parent(o,b),b) + if !include_self + siblings = filter(x-> x !== o, siblings) + end + return siblings +end + +""" + possibleneighbors(o::OctantBWG{2},l,b;insidetree=true) +Returns a list of possible neighbors, where the first four are corner neighbors that are exclusively connected via a corner. +The other four possible neighbors are face neighbors. +""" +function possibleneighbors(o::OctantBWG{2},l,b;insidetree=true) + neighbors = ntuple(8) do i + if i > 4 + j = i - 4 + facet_neighbor(o,j,b) + else + corner_neighbor(o,i,b) + end + end + if insidetree + neighbors = filter(x->inside(x,b),neighbors) + end + return neighbors +end + +""" + possibleneighbors(o::OctantBWG{3},l,b;insidetree=true) +Returns a list of possible neighbors, where the first eight are corner neighbors that are exclusively connected via a corner. +After the first eight corner neighbors, the 6 possible face neighbors follow and after them, the edge neighbors. +""" +function possibleneighbors(o::OctantBWG{3},l,b;insidetree=true) + neighbors = ntuple(26) do i + if 8 < i ≤ 14 + j = i - 8 + facet_neighbor(o,j,b) + elseif 14 < i ≤ 26 + j = i - 14 + edge_neighbor(o,j,b) + else + corner_neighbor(o,i,b) + end + end + if insidetree + neighbors = filter(x->inside(x,b),neighbors) + end + return neighbors +end + +""" + isancestor(o1,o2,b) -> Bool +Is o2 an ancestor of o1 +""" +function isancestor(o1,o2,b) + ancestor = false + l = o2.l - 1 + p = parent(o2,b) + while l > 0 + if p == o1 + ancestor = true + break + end + l -= 1 + p = parent(p,b) + end + return ancestor +end + +function contains_facet(mface::Tuple{Tuple{T1,T1},Tuple{T1,T1}},sface::Tuple{Tuple{T2,T2},Tuple{T2,T2}}) where {T1<:Integer,T2<:Integer} + if mface[1][1] == sface[1][1] && mface[2][1] == sface[2][1] # vertical + return mface[1][2] ≤ sface[1][2] ≤ sface[2][2] ≤ mface[2][2] + elseif mface[1][2] == sface[1][2] && mface[2][2] == sface[2][2] # horizontal + return mface[1][1] ≤ sface[1][1] ≤ sface[2][1] ≤ mface[2][1] + else + return false + end +end + +# currently checking if sface centroid lies in mface +# TODO should be checked if applicaple in general, I guess yes +function contains_facet(mface::NTuple{4,Tuple{T1,T1,T1}}, sface::NTuple{4,Tuple{T2,T2,T2}}) where {T1<:Integer,T2<:Integer} + sface_center = center(sface) + lower_left = ntuple(i->minimum(getindex.(mface,i)),3) + top_right = ntuple(i->maximum(getindex.(mface,i)),3) + if (lower_left[1] ≤ sface_center[1] ≤ top_right[1]) && (lower_left[2] ≤ sface_center[2] ≤ top_right[2]) && (lower_left[3] ≤ sface_center[3] ≤ top_right[3]) + return true + else + return false + end +end + +function face_contains_edge(f::NTuple{4,Tuple{T1,T1,T1}},e::Tuple{Tuple{T2,T2,T2},Tuple{T2,T2,T2}}) where {T1<:Integer,T2<:Integer} + edge_center = center(e) + lower_left = ntuple(i->minimum(getindex.(f,i)),3) + top_right = ntuple(i->maximum(getindex.(f,i)),3) + if (lower_left[1] ≤ edge_center[1] ≤ top_right[1]) && (lower_left[2] ≤ edge_center[2] ≤ top_right[2]) && (lower_left[3] ≤ edge_center[3] ≤ top_right[3]) + return true + else + return false + end +end + +function contains_edge(medge::Tuple{Tuple{T1,T1,T1},Tuple{T1,T1,T1}},sedge::Tuple{Tuple{T2,T2,T2},Tuple{T2,T2,T2}}) where {T1<:Integer,T2<:Integer} + if (medge[1][1] == sedge[1][1] && medge[2][1] == sedge[2][1]) && (medge[1][2] == sedge[1][2] && medge[2][2] == sedge[2][2]) # x1 & x2 aligned + return medge[1][3] ≤ sedge[1][3] ≤ sedge[2][3] ≤ medge[2][3] + elseif (medge[1][1] == sedge[1][1] && medge[2][1] == sedge[2][1]) && (medge[1][3] == sedge[1][3] && medge[2][3] == sedge[2][3])# x1 & x3 aligned + return medge[1][2] ≤ sedge[1][2] ≤ sedge[2][2] ≤ medge[2][2] + elseif (medge[1][2] == sedge[1][2] && medge[2][2] == sedge[2][2]) && (medge[1][3] == sedge[1][3] && medge[2][3] == sedge[2][3])# x2 & x3 aligned + return medge[1][1] ≤ sedge[1][1] ≤ sedge[2][1] ≤ medge[2][1] + else + return false + end +end + +function center(pivot_face) + centerpoint = ntuple(i->0,length(pivot_face[1])) + for c in pivot_face + centerpoint = c .+ centerpoint + end + return centerpoint .÷ length(pivot_face) +end + +iscenter(c,f) = c == center(f) + +function Base.show(io::IO, ::MIME"text/plain", agrid::ForestBWG) + println(io, "ForestBWG with ") + println(io, " $(getncells(agrid)) cells") + println(io, " $(length(agrid.cells)) trees") +end + +""" + child_id(octant::OctantBWG, b::Integer) +Given some OctantBWG `octant` and maximum refinement level `b`, compute the child_id of `octant` +note the following quote from Bursedde et al: + children are numbered from 0 for the front lower left child, + to 1 for the front lower right child, to 2 for the back lower left, and so on, with + 4, . . . , 7 being the four children on top of the children 0, . . . , 3. +shifted by 1 due to julia 1 based indexing +""" +function child_id(octant::OctantBWG{dim,N,T},b::Integer=_maxlevel[2]) where {dim,N,T<:Integer} + i = 0x00 + t = T(2) + z = zero(T) + h = T(_compute_size(b,octant.l)) + xyz = octant.xyz + for j in 0:(dim-1) + i = i | ((xyz[j+1] & h) != z ? t^j : z) + end + return i+0x01 +end + +""" + ancestor_id(octant::OctantBWG, l::Integer, b::Integer) +Algorithm 3.2 of [IBWG2015](@citet) that generalizes `child_id` for different queried levels. +Applied to a single octree, i.e. the array of leaves, yields a monotonic sequence +""" +function ancestor_id(octant::OctantBWG{dim,N,T}, l::Integer, b::Integer=_maxlevel[dim-1]) where {dim,N,T<:Integer} + @assert 0 < l ≤ octant.l + i = 0x00 + t = T(2) + z = zero(T) + h = T(_compute_size(b,l)) + for j in 0:(dim-1) + i = i | ((octant.xyz[j+1] & h) != z ? t^j : z) + end + return i+0x01 +end + +function parent(octant::OctantBWG{dim,N,T}, b::Integer=_maxlevel[dim-1]) where {dim,N,T} + if octant.l > zero(T) + h = T(_compute_size(b,octant.l)) + l = octant.l - one(T) + return OctantBWG(l,octant.xyz .& ~h) + else + root(dim) + end +end + +""" + descendants(octant::OctantBWG, b::Integer) +Given an `octant`, computes the two smallest possible octants that fit into the first and last corners +of `octant`, respectively. These computed octants are called first and last descendants of `octant` +since they are connected to `octant` by a path down the octree to the maximum level `b` +""" +function descendants(octant::OctantBWG{dim,N,T}, b::Integer=_maxlevel[dim-1]) where {dim,N,T} + l1 = b; l2 = b + h = T(_compute_size(b,octant.l)) + return OctantBWG(l1,octant.xyz), OctantBWG(l2,octant.xyz .+ (h-one(T))) +end + +""" + facet_neighbor(octant::OctantBWG{dim,N,T}, f::T, b::T=_maxlevel[2]) -> OctantBWG{dim,N,T} +Intraoctree face neighbor for a given faceindex `f` (in p4est, i.e. z order convention) and specified maximum refinement level `b`. +Implements Algorithm 5 of [BWG2011](@citet). + + x-------x-------x + | | | + | 3 | 4 | + | | | + x-------x-------x + | | | + o 1 * 2 | + | | | + x-------x-------x + +Consider octant 1 at `xyz=(0,0)`, a maximum refinement level of 1 and faceindex 2 (marked as `*`). +Then, the computed face neighbor will be octant 2 with `xyz=(1,0)`. +Note that the function is not sensitive in terms of leaving the octree boundaries. +For the above example, a query for face index 1 (marked as `o`) will return an octant outside of the octree with `xyz=(-1,0)`. +""" +function facet_neighbor(octant::OctantBWG{3,N,T}, f::T, b::T=_maxlevel[2]) where {N,T<:Integer} + l = octant.l + h = T(_compute_size(b,octant.l)) + x,y,z = octant.xyz + x += ((f == T(1)) ? -h : ((f == T(2)) ? h : zero(T))) + y += ((f == T(3)) ? -h : ((f == T(4)) ? h : zero(T))) + z += ((f == T(5)) ? -h : ((f == T(6)) ? h : zero(T))) + return OctantBWG(l,(x,y,z)) +end +function facet_neighbor(octant::OctantBWG{2,N,T}, f::T, b::T=_maxlevel[1]) where {N,T<:Integer} + l = octant.l + h = T(_compute_size(b,octant.l)) + x,y = octant.xyz + x += ((f == T(1)) ? -h : ((f == T(2)) ? h : zero(T))) + y += ((f == T(3)) ? -h : ((f == T(4)) ? h : zero(T))) + return OctantBWG(l,(x,y)) +end +facet_neighbor(o::OctantBWG{dim,N,T1}, f::T2, b::T3) where {dim,N,T1<:Integer,T2<:Integer,T3<:Integer} = facet_neighbor(o,T1(f),T1(b)) + +reference_faces_bwg(::Type{Ferrite.RefHypercube{2}}) = ((1,3) , (2,4), (1,2), (3,4)) +reference_faces_bwg(::Type{Ferrite.RefHypercube{3}}) = ((1,3,5,7) , (2,4,6,8), (1,2,5,6), (3,4,7,8), (1,2,3,4), (5,6,7,8)) # p4est consistent ordering +reference_edges_bwg(::Type{Ferrite.RefHypercube{3}}) = ((𝒰[1,1],𝒰[1,2]), (𝒰[2,1],𝒰[2,2]), (𝒰[3,1],𝒰[3,2]), + (𝒰[4,1],𝒰[4,2]), (𝒰[5,1],𝒰[5,2]), (𝒰[6,1],𝒰[6,2]), + (𝒰[7,1],𝒰[7,2]), (𝒰[8,1],𝒰[8,2]), (𝒰[9,1],𝒰[9,2]), + (𝒰[10,1],𝒰[10,2]), (𝒰[11,1],𝒰[11,2]), (𝒰[12,1],𝒰[12,2])) # TODO maybe remove, unnecessary, can directly use the table +# reference_faces_bwg(::Type{RefHypercube{3}}) = ((1,3,7,5) , (2,4,8,6), (1,2,6,5), (3,4,8,7), (1,2,4,4), (5,6,8,7)) # Note that this does NOT follow P4est order! + +""" + compute_face_orientation(forest::ForestBWG, k::Integer, f::Integer) +Slow implementation for the determination of the face orientation of face `f` from octree `k` following definition 2.1 from [BWG2011](@citet). + +TODO use table 3 for more vroom +""" +function compute_face_orientation(forest::ForestBWG{<:Any,<:OctreeBWG{dim,<:Any,T2}}, k::T1, f::T1) where {dim,T1,T2} + f_perm = (dim == 2 ? 𝒱₂_perm : 𝒱₃_perm) + f_perminv = (dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv) + n_perm = (dim == 2 ? node_map₂ : node_map₃) + n_perminv = (dim == 2 ? node_map₂_inv : node_map₃_inv) + + f_ferrite = f_perm[f] + facet_neighbor_table = Ferrite.get_facet_facet_neighborhood(forest) + k′, f′_ferrite = facet_neighbor_table[k,f_ferrite][1] + f′ = f_perminv[f′_ferrite] + reffacenodes = reference_faces_bwg(Ferrite.RefHypercube{dim}) + nodes_f = [forest.cells[k].nodes[n_perm[ni]] for ni in reffacenodes[f]] + nodes_f′ = [forest.cells[k′].nodes[n_perm[ni]] for ni in reffacenodes[f′]] + if f > f′ + return T2(findfirst(isequal(nodes_f′[1]), nodes_f)-1) + else + return T2(findfirst(isequal(nodes_f[1]), nodes_f′)-1) + end +end + +""" + compute_edge_orientation(forest::ForestBWG, k::Integer, e::Integer) +Slow implementation for the determination of the edge orientation of edge `e` from octree `k` following definition below Table 3 [BWG2011](@citet). + +TODO use some table? +""" +function compute_edge_orientation(forest::ForestBWG{<:Any,<:OctreeBWG{3,<:Any,T2}}, k::T1, e::T1) where {T1,T2} + n_perm = node_map₃ + n_perminv = node_map₃_inv + e_perm = edge_perm + e_perminv = edge_perm_inv + + e_ferrite = e_perm[e] + k′, e′_ferrite = forest.topology.edge_edge_neighbor[k,e_ferrite][1] + e′ = e_perminv[e′_ferrite] + refedgenodes = reference_edges_bwg(Ferrite.RefHypercube{3}) + nodes_e = ntuple(i->forest.cells[k].nodes[n_perm[refedgenodes[e][i]]],length(refedgenodes[e])) + nodes_e′ = ntuple(i->forest.cells[k′].nodes[n_perm[refedgenodes[e′][i]]],length(refedgenodes[e′])) + if nodes_e == nodes_e′ + s = T2(0) + else + s = T2(1) + end + return s +end + +""" + transform_facet_remote(forest::ForestBWG, k::T1, f::T1, o::OctantBWG{dim,N,T2}) -> OctantBWG{dim,N,T1,T2} + transform_facet_remote(forest::ForestBWG, f::FacetIndex, o::OctantBWG{dim,N,T2}) -> OctantBWG{dim,N,T2} +Interoctree coordinate transformation of an given octant `o` to the face-neighboring of octree `k` by virtually pushing `o`s coordinate system through `k`s face `f`. +Implements Algorithm 8 of [BWG2011](@citet). + + x-------x-------x + | | | + | 3 | 4 | + | | | + x-------x-------x + | | | + | 1 * 2 | + | | | + x-------x-------x + +Consider 4 octrees with a single leaf each and a maximum refinement level of 1 +This function transforms octant 1 into the coordinate system of octant 2 by specifying `k=2` and `f=1`. +While in the own octree coordinate system octant 1 is at `xyz=(0,0)`, the returned and transformed octant is located at `xyz=(-2,0)` +""" +function transform_facet_remote(forest::ForestBWG, k::T1, f::T1, o::OctantBWG{dim,N,T2}) where {dim,N,T1<:Integer,T2<:Integer} + _one = one(T2) + _two = T2(2) + _perm = (dim == 2 ? 𝒱₂_perm : 𝒱₃_perm) + _perminv = (dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv) + facet_neighbor_table = Ferrite.get_facet_facet_neighborhood(forest) + k′, f′ = facet_neighbor_table[k,_perm[f]][1] + f′ = _perminv[f′] + s′ = _one - (((f - _one) & _one) ⊻ ((f′ - _one) & _one)) + s = zeros(T2,dim-1) + a = zeros(T2,3) # Coordinate axes of f + b = zeros(T2,3) # Coordinate axes of f' + r = compute_face_orientation(forest,k,f) + a[3] = (f - _one) ÷ 2; b[3] = (f′ - _one) ÷ 2 # origin and target normal axis + if dim == 2 + a[1] = 1 - a[3]; b[1] = 1 - b[3]; s[1] = r + else + a[1] = (f < 3) ? 1 : 0; a[2] = (f < 5) ? 2 : 1 + u = (ℛ[1,f] - _one) ⊻ (ℛ[1,f′] - _one) ⊻ (((r == 0) | (r == 3))) + b[u+1] = (f′ < 3) ? 1 : 0; b[1-u+1] = (f′ < 5) ? 2 : 1 # r = 0 -> index 1 + if ℛ[f,f′] == 1+1 # R is one-based + s[2] = r & 1; s[1] = r & 2 + else + s[1] = r & 1; s[2] = r & 2 + end + end + maxlevel = forest.cells[1].b + l = o.l; g = 2^maxlevel - 2^(maxlevel-l) + xyz = zeros(T2,dim) + xyz[b[1] + _one] = T2((s[1] == 0) ? o.xyz[a[1] + _one] : g - o.xyz[a[1] + _one]) + xyz[b[3] + _one] = T2(((_two*((f′ - _one) & 1)) - _one)*2^maxlevel + s′*g + (1-2*s′)*o.xyz[a[3] + _one]) + if dim == 2 + return OctantBWG(l,(xyz[1],xyz[2])) + else + xyz[b[2] + _one] = T2((s[2] == 0) ? o.xyz[a[2] + _one] : g - o.xyz[a[2] + _one]) + return OctantBWG(l,(xyz[1],xyz[2],xyz[3])) + end +end + +transform_facet_remote(forest::ForestBWG,f::FacetIndex,oct::OctantBWG) = transform_facet_remote(forest,f[1],f[2],oct) + +""" + transform_facet(forest::ForestBWG, k', f', o::OctantBWG) -> OctantBWG + transform_facet(forest::ForestBWG, f'::FacetIndex, o::OctantBWG) -> OctantBWG +Interoctree coordinate transformation of an given octant `o` that lies outside of the pivot octree `k`, namely in neighbor octree `k'`. +However, the coordinate of `o` is given in octree coordinates of `k`. +Thus, this algorithm implements the transformation of the octree coordinates of `o` into the octree coordinates of `k'`. +Useful in order to check whether or not a possible neighbor exists in a neighboring octree. +Implements Algorithm 8 of [BWG2011](@citet). + + x-------x-------x + | | | + | 3 | 4 | + | | | + x-------x-------x + | | | + | 1 * 2 | + | | | + x-------x-------x + +Consider 4 octrees with a single leaf each and a maximum refinement level of 1 +This function transforms octant 1 into the coordinate system of octant 2 by specifying `k=1` and `f=2`. +While from the perspective of octree coordinates `k=2` octant 1 is at `xyz=(-2,0)`, the returned and transformed octant is located at `xyz=(0,0)` +""" +function transform_facet(forest::ForestBWG, k::T1, f::T1, o::OctantBWG{2,<:Any,T2}) where {T1<:Integer,T2<:Integer} + _one = one(T2) + _two = T2(2) + _perm = 𝒱₂_perm + _perminv = 𝒱₂_perm_inv + k′, f′ = forest.topology.edge_edge_neighbor[k,_perm[f]][1] + f′ = _perminv[f′] + + r = compute_face_orientation(forest,k,f) + # Coordinate axes of f + a = ( + f ≤ 2, # tangent + f > 2 # normal + ) + a_sign = _two*((f - _one) & 1) - _one + # Coordinate axes of f' + b = ( + f′ ≤ 2, # tangent + f′ > 2 # normal + ) + # b_sign = _two*(f′ & 1) - _one + + maxlevel = forest.cells[1].b + depth_offset = 2^maxlevel - 2^(maxlevel-o.l) + + s′ = _one - (((f - _one) & _one) ⊻ ((f′ - _one) & _one)) # arithmetic switch: TODO understand this. + + # xyz = zeros(T2, 2) + # xyz[a[1] + _one] = T2((r == 0) ? o.xyz[b[1] + _one] : depth_offset - o.xyz[b[1] + _one]) + # xyz[a[2] + _one] = T2(a_sign*2^maxlevel + s′*depth_offset + (1-2*s′)*o.xyz[b[2] + _one]) + # return OctantBWG(o.l,(xyz[1],xyz[2])) + + # We can do this because the permutation and inverse permutation are the same + xyz = ( + T2((r == 0) ? o.xyz[b[1] + _one] : depth_offset - o.xyz[b[1] + _one]), + T2(a_sign*2^maxlevel + s′*depth_offset + (1-2*s′)*o.xyz[b[2] + _one]) + ) + return OctantBWG(o.l,(xyz[a[1] + _one],xyz[a[2] + _one])) +end + +function transform_facet(forest::ForestBWG, k::T1, f::T1, o::OctantBWG{3,<:Any,T2}) where {T1<:Integer,T2<:Integer} + _one = one(T2) + _two = T2(2) + _perm = 𝒱₃_perm + _perminv = 𝒱₃_perm_inv + k′, f′ = forest.topology.face_face_neighbor[k,_perm[f]][1] + f′ = _perminv[f′] + s′ = _one - (((f - _one) & _one) ⊻ ((f′ - _one) & _one)) + r = compute_face_orientation(forest,k,f) + + # Coordinate axes of f + a = ( + (f ≤ 2) ? 1 : 0, + (f ≤ 4) ? 2 : 1, + (f - _one) ÷ 2 + ) + a_sign = _two*((f - _one) & 1) - _one + + # Coordinate axes of f' + b = if Bool(ℛ[1,f] - _one) ⊻ Bool(ℛ[1,f′] - _one) ⊻ (((r == 0) || (r == 3))) # What is this condition exactly? + ( + (f′ < 5) ? 2 : 1, + (f′ < 3) ? 1 : 0, + (f′ - _one) ÷ 2 + ) + else + ( + (f′ < 3) ? 1 : 0, + (f′ < 5) ? 2 : 1, + (f′ - _one) ÷ 2 + ) + end + # b_sign = _two*(f′ & 1) - _one + + s = if ℛ[f,f′] == 1+1 # R is one-based + (r & 2, r & 1) + else + (r & 1, r & 2) + end + maxlevel = forest.cells[1].b + depth_offset = 2^maxlevel - 2^(maxlevel-o.l) + xyz = zeros(T2,3) + xyz[a[1] + _one] = T2((s[1] == 0) ? o.xyz[b[1] + _one] : depth_offset - o.xyz[b[1] + _one]) + xyz[a[2] + _one] = T2((s[2] == 0) ? o.xyz[b[2] + _one] : depth_offset - o.xyz[b[2] + _one]) + xyz[a[3] + _one] = T2(a_sign*2^maxlevel + s′*depth_offset + (1-2*s′)*o.xyz[b[3] + _one]) + return OctantBWG(o.l,(xyz[1],xyz[2],xyz[3])) + + # xyz = ( + # T2((s[1] == 0) ? o.xyz[b[1] + _one] : depth_offset - o.xyz[b[1] + _one]), + # T2((s[2] == 0) ? o.xyz[b[2] + _one] : depth_offset - o.xyz[b[2] + _one]), + # T2(a_sign*2^maxlevel + s′*depth_offset + (1-2*s′)*o.xyz[b[3] + _one]) + # ) + # return OctantBWG(o.l,(xyz[a[1] + _one],xyz[a[2] + _one],xyz[a[3] + _one])) +end + +transform_facet(forest::ForestBWG,f::FacetIndex,oct::OctantBWG) = transform_facet(forest,f[1],f[2],oct) + +""" + transform_corner(forest,k,c',oct,inside::Bool) + transform_corner(forest,v::VertexIndex,oct,inside::Bool) + +Algorithm 12 but with flipped logic in [BWG2011](@citet) to transform corner into different octree coordinate system +Implements flipped logic in the sense of pushing the Octant `oct` through vertex v and stays within octree coordinate system `k`. +""" +function transform_corner(forest::ForestBWG,k::T1,c::T1,oct::OctantBWG{dim,N,T2},inside::Bool) where {dim,N,T1<:Integer,T2<:Integer} + _perm = dim == 2 ? node_map₂ : node_map₃ + _perminv = dim == 2 ? node_map₂_inv : node_map₃_inv + k′, c′ = forest.topology.vertex_vertex_neighbor[k,_perm[c]][1] + k′, c′ = forest.topology.vertex_vertex_neighbor[k′,c′][1] #get the corner connection of neighbor to pivot oct + c′ = _perminv[c′] + # make a dispatch that returns only the coordinates? + b = forest.cells[k].b + l = oct.l; g = 2^b - 2^(b-l) + h⁻ = inside ? 0 : -2^(b-l); h⁺ = inside ? g : 2^b + xyz = ntuple(i->((c′-1) & 2^(i-1) == 0) ? h⁻ : h⁺,dim) + return OctantBWG(l,xyz) +end + +transform_corner(forest::ForestBWG,v::VertexIndex,oct::OctantBWG,inside) = transform_corner(forest,v[1],v[2],oct,inside) + +""" + transform_corner_remote(forest,k,c',oct,inside::Bool) + transform_corner_remote(forest,v::VertexIndex,oct,inside::Bool) + +Algorithm 12 in [BWG2011](@citet) to transform corner into different octree coordinate system. +Follows exactly the version of the paper by taking `oct` and looking from the neighbor octree coordinate system (neighboring to `k`,`v`) at `oct`. +""" +function transform_corner_remote(forest::ForestBWG,k::T1,c::T1,oct::OctantBWG{dim,N,T2},inside::Bool) where {dim,N,T1<:Integer,T2<:Integer} + _perm = dim == 2 ? node_map₂ : node_map₃ + _perminv = dim == 2 ? node_map₂_inv : node_map₃_inv + k′, c′ = forest.topology.vertex_vertex_neighbor[k,_perm[c]][1] + c′ = _perminv[c′] + # make a dispatch that returns only the coordinates? + b = forest.cells[k].b + l = oct.l; g = 2^b - 2^(b-l) + h⁻ = inside ? 0 : -2^(b-l); h⁺ = inside ? g : 2^b + xyz = ntuple(i->((c′-1) & 2^(i-1) == 0) ? h⁻ : h⁺,dim) + return OctantBWG(l,xyz) +end + +transform_corner_remote(forest::ForestBWG,v::VertexIndex,oct::OctantBWG,inside) = transform_corner_remote(forest,v[1],v[2],oct,inside) + + +""" + transform_edge_remote(forest,k,e,oct,inside::Bool) + transform_edge_remote(forest,e::Edgeindex,oct,inside::Bool) + +Algorithm 10 in [BWG2011](@citet) to transform edge into different octree coordinate system. +This function looks at the octant from the octree coordinate system of the neighbor that can be found at (k,e) +""" +function transform_edge_remote(forest::ForestBWG,k::T1,e::T1,oct::OctantBWG{3,N,T2},inside::Bool) where {N,T1<:Integer,T2<:Integer} + _four = T2(4) + _one = T2(1) + _two = T2(2) + z = zero(T2) + e_perm = edge_perm + e_perminv = edge_perm_inv + + e_ferrite = e_perm[e] + k′, e′_ferrite = forest.topology.edge_edge_neighbor[k,e_ferrite][1] + e′ = e_perminv[e′_ferrite] + #see Algorithm 9, line 18 + 𝐛 = (((e′-_one) ÷ _four), + e′-_one < 4 ? 1 : 0, + e′-_one < 8 ? 2 : 1) + a₀ = ((e-_one) ÷ _four) #subtract 1 based index + a₀ += _one #add it again + b = forest.cells[k].b + l = oct.l; g = _two^b - _two^(b-l) + h⁻ = inside ? z : -_two^(b-l); h⁺ = inside ? g : _two^b + s = compute_edge_orientation(forest,k,e) + xyz = zeros(T2,3) + xyz[𝐛[1]+_one] = s*g+(_one-(_two*s))*oct.xyz[a₀] + xyz[𝐛[2]+_one] = ((e′-_one) & 1) == 0 ? h⁻ : h⁺ + xyz[𝐛[3]+_one] = ((e′-_one) & 2) == 0 ? h⁻ : h⁺ + return OctantBWG(l,(xyz[1],xyz[2],xyz[3])) +end + +transform_edge_remote(forest::ForestBWG,e::EdgeIndex,oct::OctantBWG,inside) = transform_edge_remote(forest,e[1],e[2],oct,inside) + +""" + transform_edge(forest,k,e,oct,inside::Bool) + transform_edge(forest,e::Edgeindex,oct,inside::Bool) + +Algorithm 10 in [BWG2011](@citet) to transform cedge into different octree coordinate system but reversed logic. +See `transform_edge_remote` with logic from paper. +In this function we stick to the coordinate system of the pivot tree k and transform an octant through edge e into this k-th octree coordinate system. +""" +function transform_edge(forest::ForestBWG,k::T1,e::T1,oct::OctantBWG{3,N,T2},inside::Bool) where {N,T1<:Integer,T2<:Integer} + _four = T2(4) + _one = T2(1) + _two = T2(2) + z = zero(T2) + e_perm = edge_perm + e_perminv = edge_perm_inv + + e_ferrite = e_perm[e] + k′, e′_ferrite = forest.topology.edge_edge_neighbor[k,e_ferrite][1] + k′, e′_ferrite = forest.topology.edge_edge_neighbor[k′,e′_ferrite][1] #get pivot connection from neighbor perspective + e′ = e_perminv[e′_ferrite] + #see Algorithm 9, line 18 + 𝐛 = (((e′-_one) ÷ _four), + e′-_one < 4 ? 1 : 0, + e′-_one < 8 ? 2 : 1) + a₀ = ((e-_one) ÷ _four) #subtract 1 based index + a₀ += _one #add it again + b = forest.cells[k].b + l = oct.l; g = _two^b - _two^(b-l) + h⁻ = inside ? z : -_two^(b-l); h⁺ = inside ? g : _two^b + s = compute_edge_orientation(forest,k′,e′) + xyz = zeros(T2,3) + xyz[𝐛[1]+_one] = s*g+(_one-(_two*s))*oct.xyz[a₀] + xyz[𝐛[2]+_one] = ((e′-_one) & 1) == 0 ? h⁻ : h⁺ + xyz[𝐛[3]+_one] = ((e′-_one) & 2) == 0 ? h⁻ : h⁺ + return OctantBWG(l,(xyz[1],xyz[2],xyz[3])) +end + +transform_edge(forest::ForestBWG,e::EdgeIndex,oct::OctantBWG,inside) = transform_edge(forest,e[1],e[2],oct,inside) + +""" + edge_neighbor(octant::OctantBWG, e::Integer, b::Integer) +Computes the edge neighbor octant which is only connected by the edge `e` to `octant`. +""" +function edge_neighbor(octant::OctantBWG{3,N,T}, e::T, b::T=_maxlevel[2]) where {N,T<:Integer} + @assert 1 ≤ e ≤ 12 + _one = one(T) + _two = T(2) + e -= _one + l = octant.l + h = T(_compute_size(b,octant.l)) + ox,oy,oz = octant.xyz + 𝐚 = (e ÷ 4, + e < 4 ? 1 : 0, + e < 8 ? 2 : 1) + xyz = zeros(T,3) + xyz[𝐚[1]+_one] = octant.xyz[𝐚[1]+_one] + xyz[𝐚[2]+_one] = octant.xyz[𝐚[2]+_one] + (_two*(e&_one)-_one)*h + xyz[𝐚[3]+_one] = octant.xyz[𝐚[3]+_one] + ((e & _two)-_one)*h + return OctantBWG(l,(xyz[1],xyz[2],xyz[3])) +end +edge_neighbor(o::OctantBWG{3,N,T1}, e::T2, b::T3) where {N,T1<:Integer,T2<:Integer,T3<:Integer} = edge_neighbor(o,T1(e),T1(b)) + +""" + corner_neighbor(octant::OctantBWG, c::Integer, b::Integer) +Computes the corner neighbor octant which is only connected by the corner `c` to `octant` +""" +function corner_neighbor(octant::OctantBWG{3,N,T}, c::T, b::T=_maxlevel[2]) where {N,T<:Integer} + c -= one(T) + l = octant.l + h = T(_compute_size(b,octant.l)) + ox,oy,oz = octant.xyz + _one = one(T) + _two = T(2) + x = ox + (_two*(c & _one) - _one)*h + y = oy + ((c & _two) - _one)*h + z = oz + ((c & T(4))÷_two - _one)*h + return OctantBWG(l,(x,y,z)) +end + +function corner_neighbor(octant::OctantBWG{2,N,T}, c::T, b::T=_maxlevel[1]) where {N,T<:Integer} + c -= one(T) + l = octant.l + h = _compute_size(b,octant.l) + ox,oy = octant.xyz + _one = one(T) + _two = T(2) + x = ox + (_two*(c & _one) - _one)*h + y = oy + ((c & _two) - _one)*h + return OctantBWG(l,(x,y)) +end +corner_neighbor(o::OctantBWG{dim,N,T1}, c::T2, b::T3) where {dim,N,T1<:Integer,T2<:Integer,T3<:Integer} = corner_neighbor(o,T1(c),T1(b)) + +function corner_face_participation(dim::T,c::T) where T<:Integer + if dim == 2 + return 𝒱₂_perm[findall(x->c ∈ x, eachrow(𝒱₂))] + else + return 𝒱₃_perm[findall(x->c ∈ x, eachrow(𝒱₃))] + end +end + +function Base.show(io::IO, ::MIME"text/plain", o::OctantBWG{3,N,M}) where {N,M} + x,y,z = o.xyz + println(io, "OctantBWG{3,$N,$M}") + println(io, " l = $(o.l)") + println(io, " xyz = $x,$y,$z") +end + +function Base.show(io::IO, ::MIME"text/plain", o::OctantBWG{2,N,M}) where {N,M} + x,y = o.xyz + println(io, "OctantBWG{2,$N,$M}") + println(io, " l = $(o.l)") + println(io, " xy = $x,$y") +end + +_compute_size(b::Integer,l::Integer) = 2^(b-l) +_maximum_size(b::Integer) = 2^(b) +# return the two adjacent faces $f_i$ adjacent to edge `edge` +_face(edge::Int) = 𝒮[edge, :] +# return the `i`-th adjacent face fᵢ to edge `edge` +_face(edge::Int, i::Int) = 𝒮[edge, i] +# return two face corners ξᵢ of the face `face` along edge `edge` +_face_edge_corners(edge::Int, face::Int) = 𝒯[edge,face] +# return the two `edge` corners cᵢ +_edge_corners(edge::Int) = 𝒰[edge,:] +# return the `i`-th edge corner of `edge` +_edge_corners(edge::Int,i::Int) = 𝒰[edge,i] +# finds face corner ξ′ in f′ for two associated faces f,f′ in {1,...,6} and their orientation r in {1,...,4}} +_neighbor_corner(f::Int,f′::Int,r::Int,ξ::Int) = 𝒫[𝒬[ℛ[f,f′],r],ξ] + +# map given `face` and `ξ` to corner `c`. Need to provide dim for different lookup +function _face_corners(dim::Int,face::Int,ξ::Int) + if dim == 2 + return 𝒱₂[face,ξ] + elseif dim == 3 + return 𝒱₃[face,ξ] + else + error("No corner-lookup table available") + end +end + +function _face_corners(dim::Int,face::Int) + if dim == 2 + return 𝒱₂[face,:] + elseif dim == 3 + return 𝒱₃[face,:] + else + error("No corner-lookup table available") + end +end + +##### OCTANT LOOK UP TABLES ###### +const 𝒮 = [3 5 + 4 5 + 3 6 + 4 6 + 1 5 + 2 5 + 1 6 + 2 6 + 1 3 + 2 3 + 1 4 + 2 4] + +# (0,0) non existing connections +const 𝒯 = [(0, 0) (0, 0) (1, 2) (0, 0) (1, 2) (0, 0) + (0, 0) (0, 0) (0, 0) (1, 2) (3, 4) (0, 0) + (0, 0) (0, 0) (3, 4) (0, 0) (0, 0) (1, 2) + (0, 0) (0, 0) (0, 0) (3, 4) (0, 0) (3, 4) + (1, 2) (0, 0) (0, 0) (0, 0) (1, 3) (0, 0) + (0, 0) (1, 2) (0, 0) (0, 0) (2, 4) (0, 0) + (3, 4) (0, 0) (0, 0) (0, 0) (0, 0) (1, 3) + (0, 0) (3, 4) (0, 0) (0, 0) (0, 0) (2, 4) + (1, 3) (0, 0) (1, 3) (0, 0) (0, 0) (0, 0) + (0, 0) (1, 3) (2, 4) (0, 0) (0, 0) (0, 0) + (2, 4) (0, 0) (0, 0) (1, 3) (0, 0) (0, 0) + (0, 0) (2, 4) (0, 0) (2, 4) (0, 0) (0, 0)] + +const 𝒰 = [1 2 + 3 4 + 5 6 + 7 8 + 1 3 + 2 4 + 5 7 + 6 8 + 1 5 + 2 6 + 3 7 + 4 8] + +const 𝒱₂ = [1 3 + 2 4 + 1 2 + 3 4] + +const 𝒱₃ = [1 3 5 7 + 2 4 6 8 + 1 2 5 6 + 3 4 7 8 + 1 2 3 4 + 5 6 7 8] + +# Face indices permutation from p4est idx to Ferrite idx +const 𝒱₂_perm = [4 + 2 + 1 + 3] + +# Face indices permutation from Ferrite idx to p4est idx +const 𝒱₂_perm_inv = [3 + 2 + 4 + 1] + +const 𝒱₃_perm = [5 + 3 + 2 + 4 + 1 + 6] + +const 𝒱₃_perm_inv = [5 + 3 + 2 + 4 + 1 + 6] + +# edge indices permutation from p4est idx to Ferrite idx +const edge_perm = [1 + 3 + 5 + 7 + 4 + 2 + 8 + 6 + 9 + 10 + 12 + 11] + +# edge indices permutation from Ferrite idx to p4est idx +const edge_perm_inv = [1 + 6 + 2 + 5 + 3 + 8 + 4 + 7 + 9 + 10 + 12 + 11] + +const ℛ = [1 2 2 1 1 2 + 3 1 1 2 2 1 + 3 1 1 2 2 1 + 1 3 3 1 1 2 + 1 3 3 1 1 2 + 3 1 1 3 3 1] + +const 𝒬 = [2 3 6 7 + 1 4 5 8 + 1 5 4 8] + +const 𝒫 = [1 2 3 4 + 1 3 2 4 + 2 1 4 3 + 2 4 1 3 + 3 1 4 2 + 3 4 1 2 + 4 2 3 1 + 4 3 2 1] + +const opposite_corner_2 = [4, + 3, + 2, + 1] + +const opposite_corner_3 = [8, + 7, + 6, + 5, + 4, + 3, + 2, + 1] + +const opposite_face_2 = [2, + 1, + 4, + 3] + +const opposite_face_3 = [2, + 1, + 4, + 3, + 6, + 5] + +# Node indices permutation from p4est idx to Ferrite idx +const node_map₂ = [1, + 2, + 4, + 3] + +# Node indices permutation from Ferrite idx to p4est idx +const node_map₂_inv = [1, + 2, + 4, + 3] + +const node_map₃ = [1, + 2, + 4, + 3, + 5, + 6, + 8, + 7] + +const node_map₃_inv = [1, + 2, + 4, + 3, + 5, + 6, + 8, + 7] diff --git a/src/Adaptivity/constraints.jl b/src/Adaptivity/constraints.jl new file mode 100644 index 0000000000..c86ed68e18 --- /dev/null +++ b/src/Adaptivity/constraints.jl @@ -0,0 +1,42 @@ +""" + This constraint can be passed to the constraint handler when working with non-conforming meshes to + add the affine constraints required to make the associated interpolation conforming. + + For a full example visit the AMR tutorial. +""" +struct ConformityConstraint + field_name::Symbol +end + +function Ferrite.add!(ch::ConstraintHandler{<:DofHandler{<:Any,<:Grid}}, cc::ConformityConstraint) + @warn "Trying to add conformity constraint to $(cc.field_name) on a conforming grid. Skipping." +end + +function Ferrite.add!(ch::ConstraintHandler{<:DofHandler{<:Any,<:NonConformingGrid}}, cc::ConformityConstraint) + @assert length(ch.dh.field_names) == 1 "Multiple fields not supported yet." + @assert cc.field_name ∈ ch.dh.field_names "Field $(cc.field_name) not found in provided dof handler. Available fields are $(ch.dh.field_names)." + # One set of linear contraints per hanging node + for sdh in ch.dh.subdofhandlers + field_idx = Ferrite._find_field(sdh, cc.field_name) + field_idx !== nothing && _add_conformity_constraint(ch, field_idx, sdh.field_interpolations[field_idx]) + end +end + +function _add_conformity_constraint(ch::ConstraintHandler, field_index::Int, interpolation::ScalarInterpolation) + for (hdof,mdof) in ch.dh.grid.conformity_info + @debug @assert length(mdof) == 2 + lc = AffineConstraint(ch.dh.vertexdicts[field_index][hdof],[ch.dh.vertexdicts[field_index][m] => 0.5 for m in mdof], 0.0) + add!(ch,lc) + end +end + +function _add_conformity_constraint(ch::ConstraintHandler, field_index::Int, interpolation::VectorizedInterpolation{vdim}) where vdim + for (hdof,mdof) in ch.dh.grid.conformity_info + @debug @assert length(mdof) == 2 + # One constraint per component + for vd in 1:vdim + lc = AffineConstraint(ch.dh.vertexdicts[field_index][hdof]+vd-1,[ch.dh.vertexdicts[field_index][m]+vd-1 => 0.5 for m in mdof], 0.0) # TODO change for other interpolation types than linear + add!(ch,lc) + end + end +end diff --git a/src/Adaptivity/ncgrid.jl b/src/Adaptivity/ncgrid.jl new file mode 100644 index 0000000000..fc53633bdb --- /dev/null +++ b/src/Adaptivity/ncgrid.jl @@ -0,0 +1,48 @@ +""" + NonConformingGrid{dim, C<:AbstractCell, T<:Real, CIT} <: AbstractGrid} + +A `NonConformingGrid` is a collection of `Cells` and `Node`s which covers the computational domain, together with Sets of cells, nodes, faces +and assocaited information about the conformity. +There are multiple helper structures to apply boundary conditions or define subdomains. They are gathered in the `cellsets`, `nodesets`, +`facesets`, `edgesets` and `vertexsets`. + +This grid serves as an entry point for non-intrusive adaptive grid libraries. + +# Fields +- `cells::Vector{C}`: stores all cells of the grid +- `nodes::Vector{Node{dim,T}}`: stores the `dim` dimensional nodes of the grid +- `cellsets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of cell ids +- `nodesets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of global node ids +- `facesets::Dict{String,Set{FaceIndex}}`: maps a `String` to a `Set` of `Set{FaceIndex} (global_cell_id, local_face_id)` +- `edgesets::Dict{String,Set{EdgeIndex}}`: maps a `String` to a `Set` of `Set{EdgeIndex} (global_cell_id, local_edge_id` +- `vertexsets::Dict{String,Set{VertexIndex}}`: maps a `String` key to a `Set` of local vertex ids +- `conformity_info::CIT`: a container for conformity information +- `boundary_matrix::SparseMatrixCSC{Bool,Int}`: optional, only needed by `onboundary` to check if a cell is on the boundary, see, e.g. Helmholtz example +""" +mutable struct NonConformingGrid{dim,C<:Ferrite.AbstractCell,T<:Real,CIT} <: Ferrite.AbstractGrid{dim} + cells::Vector{C} + nodes::Vector{Node{dim,T}} + # Sets + cellsets::Dict{String,OrderedSet{Int}} + nodesets::Dict{String,OrderedSet{Int}} + facetsets::Dict{String,OrderedSet{Ferrite.FacetIndex}} + vertexsets::Dict{String,OrderedSet{Ferrite.VertexIndex}} + conformity_info::CIT # TODO refine + # Boundary matrix (faces per cell × cell) + boundary_matrix::SparseMatrixCSC{Bool,Int} +end + +function NonConformingGrid( + cells::Vector{C}, + nodes::Vector{Node{dim,T}}; + cellsets::Dict{String,OrderedSet{Int}}=Dict{String,OrderedSet{Int}}(), + nodesets::Dict{String,OrderedSet{Int}}=Dict{String,OrderedSet{Int}}(), + facetsets::Dict{String,OrderedSet{Ferrite.FacetIndex}}=Dict{String,OrderedSet{Ferrite.FacetIndex}}(), + vertexsets::Dict{String,OrderedSet{Ferrite.VertexIndex}}=Dict{String,OrderedSet{Ferrite.VertexIndex}}(), + conformity_info, + boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0) + ) where {dim,C,T} + return NonConformingGrid(cells, nodes, cellsets, nodesets, facetsets, vertexsets, conformity_info, boundary_matrix) +end + +get_coordinate_type(::NonConformingGrid{dim,C,T}) where {dim,C,T} = Vec{dim,T} diff --git a/src/Adaptivity/visualization.jl b/src/Adaptivity/visualization.jl new file mode 100644 index 0000000000..e1b4a53cfe --- /dev/null +++ b/src/Adaptivity/visualization.jl @@ -0,0 +1,35 @@ +function visualize_grid(forest::ForestBWG{dim}) where dim + fig = GLMakie.Figure() + ax = dim < 3 ? GLMakie.Axis(fig[1,1]) : GLMakie.LScene(fig[1,1]) + for tree in forest.cells + for leaf in tree.leaves + cellnodes = getnodes(forest,collect(tree.nodes)) .|> get_node_coordinate |> collect + #request vertices and faces in octree coordinate system + _vertices = Ferrite.vertices(leaf,tree.b) + # transform from octree coordinate system to -1,1 by first shifting to 0,2 and later shift by -1 + _vertices = broadcast.(x->x .* 2/(2^tree.b) .- 1, _vertices) + octant_physical_coordinates = zeros(length(_vertices),dim) + for (i,v) in enumerate(_vertices) + octant_physical_coordinates[i,:] .= sum(j-> cellnodes[j] * Ferrite.shape_value(Lagrange{Ferrite.RefHypercube{dim},1}(),Vec{dim}(v),j),1:length(cellnodes)) + end + GLMakie.scatter!(ax,octant_physical_coordinates,color=:black,markersize=25) + center = sum(octant_physical_coordinates,dims=1) ./ 4 + #GLMakie.scatter!(ax,center,color=:black,markersize=25) + facetable = dim == 2 ? Ferrite.𝒱₂ : Ferrite.𝒱₃ + for faceids in eachrow(facetable) + if dim < 3 + x = octant_physical_coordinates[faceids,1] + (octant_physical_coordinates[faceids,1] .- center[1])*0.02 + y = octant_physical_coordinates[faceids,2] + (octant_physical_coordinates[faceids,2] .- center[2])*0.02 + GLMakie.lines!(ax,x,y,color=:black) + else + faceids = [faceids[1], faceids[2], faceids[4], faceids[3], faceids[1]] + x = octant_physical_coordinates[faceids,1] + (octant_physical_coordinates[faceids,1] .- center[1])*0.02 + y = octant_physical_coordinates[faceids,2] + (octant_physical_coordinates[faceids,2] .- center[2])*0.02 + z = octant_physical_coordinates[faceids,3] + (octant_physical_coordinates[faceids,3] .- center[3])*0.02 + GLMakie.lines!(ax,x,y,z,color=:black) + end + end + end + end + return fig, ax +end diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index d1eafbfa47..6c50f27366 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -238,6 +238,7 @@ function close!(ch::ConstraintHandler) i == 0 && continue icoeffs = ch.dofcoefficients[i] if !(icoeffs === nothing || isempty(icoeffs)) + @debug @info "Nested affine constraint detected for $coeffs and $i" error("nested affine constraints currently not supported") end end diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 8baf08d81b..85138a6689 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -104,6 +104,19 @@ struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler closed::ScalarWrapper{Bool} grid::G ndofs::ScalarWrapper{Int} + # 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}} end """ @@ -131,7 +144,7 @@ close!(dh) function DofHandler(grid::G) where {dim, G <: AbstractGrid{dim}} ncells = getncells(grid) sdhs = SubDofHandler{DofHandler{dim, G}}[] - DofHandler{dim, G}(sdhs, Symbol[], Int[], zeros(Int, ncells), zeros(Int, ncells), ScalarWrapper(false), grid, ScalarWrapper(-1)) + DofHandler{dim, G}(sdhs, Symbol[], Int[], zeros(Int, ncells), zeros(Int, ncells), ScalarWrapper(false), grid, ScalarWrapper(-1), [Int[]], Dict{Tuple{Int,Int}}[], Dict{NTuple{3,Int}}[]) end function Base.show(io::IO, mime::MIME"text/plain", dh::DofHandler) @@ -354,22 +367,16 @@ function __close!(dh::DofHandler{dim}) where {dim} end numfields = length(dh.field_names) - # NOTE: Maybe it makes sense to store *Index in the dicts instead. + resize!(dh.vertexdicts, numfields) + resize!(dh.edgedicts, numfields) + resize!(dh.facedicts, numfields) - # `vertexdict` keeps track of the visited vertices. The first dof added to vertex v is - # stored in vertexdict[v]. - # TODO: No need to allocate this vector for fields that don't have vertex dofs - vertexdicts = [zeros(Int, getnnodes(get_grid(dh))) for _ in 1:numfields] - - # `edgedict` keeps track of the visited edges. - # An edge is uniquely determined by two global vertices, with global direction going - # from low to high vertex node number, see sortedge - edgedicts = [Dict{NTuple{2, Int}, Int}() for _ in 1:numfields] - - # `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. - # A face is uniquely determined by 3 vertex nodes, see sortface - facedicts = [Dict{NTuple{3, Int}, Int}() for _ in 1:numfields] + nnodes = getnnodes(get_grid(dh)) + for fi in 1:numfields + dh.vertexdicts[fi] = zeros(Int, nnodes) + dh.edgedicts[fi] = Dict{Tuple{Int,Int}, Int}() + dh.facedicts[fi] = Dict{NTuple{3,Int}, Int}() + end # Set initial values nextdof = 1 # next free dof to distribute @@ -381,15 +388,15 @@ function __close!(dh::DofHandler{dim}) where {dim} sdh, sdhi, # TODO: Store in the SubDofHandler? nextdof, - vertexdicts, - edgedicts, - facedicts, + dh.vertexdicts, + dh.edgedicts, + dh.facedicts, ) end dh.ndofs[] = maximum(dh.cell_dofs; init=0) dh.closed[] = true - return dh, vertexdicts, edgedicts, facedicts + return dh, dh.vertexdicts, dh.edgedicts, dh.facedicts end diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 2b6931d723..8e91adcaa2 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -179,17 +179,17 @@ nodes_to_vtkorder(cell::QuadraticHexahedron) = [ cell.nodes[27], # interior ] -function create_vtk_griddata(grid::Grid{dim,C,T}) where {dim,C,T} +function create_vtk_griddata(grid::AbstractGrid{dim}) where dim cls = WriteVTK.MeshCell[] for cell in getcells(grid) celltype = Ferrite.cell_to_vtkcell(typeof(cell)) push!(cls, WriteVTK.MeshCell(celltype, nodes_to_vtkorder(cell))) end - coords = reshape(reinterpret(T, getnodes(grid)), (dim, getnnodes(grid))) + coords = reshape(reinterpret(get_coordinate_eltype(grid), getnodes(grid)), (dim, getnnodes(grid))) return coords, cls end -function create_vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; kwargs...) where {dim,C,T} +function create_vtk_grid(filename::AbstractString, grid::AbstractGrid{dim}; kwargs...) where dim coords, cls = create_vtk_griddata(grid) return WriteVTK.vtk_grid(filename, coords, cls; kwargs...) end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index d1b44276aa..7c6c369772 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -159,4 +159,8 @@ include("PointEvalHandler.jl") include("deprecations.jl") include("docs.jl") +# Adaptiviy +include("Adaptivity/AMR.jl") +using .AMR + end # module diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 4012d5675c..39b2d212cc 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -351,6 +351,7 @@ end Get the datatype for a single point in the grid. """ +get_coordinate_type(grid::AbstractGrid{dim}) where {dim} = get_coordinate_type(first(getnodes(grid))) get_coordinate_type(::Grid{dim,C,T}) where {dim,C,T} = Vec{dim,T} # Node is baked into the mesh type. """ diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl index 42eca02dc1..e39cd1bc37 100644 --- a/src/Grid/grid_generators.jl +++ b/src/Grid/grid_generators.jl @@ -540,3 +540,45 @@ function generate_grid(::Type{Tetrahedron}, cells_per_dim::NTuple{3,Int}, left:: return Grid(cells, nodes, facetsets=facetsets) end + +function generate_simple_disc_grid(::Type{Quadrilateral}, n; radius = 1.0) + nnodes = 2n + 1 + θ = deg2rad(360/2n) + + nodepos = Vec((0.0,radius)) + nodes = [rotate(nodepos, θ*i) for i ∈ 0:(2n-1)] + push!(nodes, Vec((0.0,0.0))) + + elements = [Quadrilateral((2i-1==0 ? nnodes-1 : 2i-1,2i,2i+1 == nnodes ? 1 : 2i+1,nnodes)) for i ∈ 1:n] + + facesets = Dict( + "boundary" => Set([FaceIndex(i,1) for i ∈ 1:n]) ∪ Set([FaceIndex(i,2) for i ∈ 1:n]), + ) + + return Grid(elements, Node.(nodes); facesets=facesets) +end + +function generate_simple_disc_grid(::Type{Hexahedron}, n; radius = 1.0, layers = 1, height = 1.0) + nnodes = 2n + 1 + θ = deg2rad(360/2n) + + nodepos_bottom = Vec((0.0,radius,0.0)) + nodes = [rotate(nodepos_bottom, Vec{3}((0,0,1)), θ*i) for i ∈ 0:(2n-1)] + push!(nodes, Vec((0.0,0.0,0.0))) + + # TODO generalize for n layers by looping over layers + nodepos_layer = Vec((0.0,radius,height)) + nodes_layer = [rotate(nodepos_layer, Vec{3}((0,0,1)), θ*i) for i ∈ 0:(2n-1)] + push!(nodes_layer, Vec((0.0,0.0,1.0))) + nodes = vcat(nodes,nodes_layer) + + elements = [Hexahedron((2i-1==0 ? nnodes-1 : 2i-1,2i,2i+1 == nnodes ? 1 : 2i+1,nnodes,2i-1==0 ? nnodes-1 : 2i-1+(2*n+1),2i+(2*n+1),2i+1 == nnodes ? (2*n+2) : 2i+1+(2*n+1),nnodes*(layers+1))) for i ∈ 1:n*layers] + + facesets = Dict( + "side" => Set([FaceIndex(i,1) for i ∈ 1:n]) ∪ Set([FaceIndex(i,2) for i ∈ 1:n]), + "top" => Set([FaceIndex(i,5) for i ∈ 1:n]), + "bottom" => Set([FaceIndex(i,6) for i ∈ 1:n]), + ) + + return Grid(elements, Node.(nodes))#; facesets=facesets) +end diff --git a/src/L2_projection.jl b/src/L2_projection.jl index b4be4e500f..080b49bcdb 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -6,6 +6,7 @@ struct L2Projector <: AbstractProjector geom_ip::Interpolation M_cholesky #::SuiteSparse.CHOLMOD.Factor{Float64} dh::DofHandler + ch::ConstraintHandler set::OrderedSet{Int} end @@ -56,10 +57,15 @@ function L2Projector( add!(sdh, :_, func_ip) # we need to create the field, but the interpolation is not used here close!(dh) - M = _assemble_L2_matrix(fe_values_mass, _set, dh) # the "mass" matrix - M_cholesky = cholesky(M) + ch_hanging = ConstraintHandler(dh) + if grid isa NonConformingGrid + add!(ch_hanging, ConformityConstraint(:_)) + end + close!(ch_hanging); + M = _assemble_L2_matrix(fe_values_mass, set, dh, ch_hanging) # the "mass" matrix + #apply!(M,ch_hanging) - return L2Projector(func_ip, geom_ip, M_cholesky, dh, _set) + return L2Projector(func_ip, geom_ip, M, dh, ch_hanging, _set) end # Quadrature sufficient for integrating a mass matrix @@ -78,10 +84,10 @@ function Base.show(io::IO, ::MIME"text/plain", proj::L2Projector) println(io, " geometric interpolation: ", proj.geom_ip) end -function _assemble_L2_matrix(fe_values, set, dh) +function _assemble_L2_matrix(fe_values, set, dh, ch) n = Ferrite.getnbasefunctions(fe_values) - M = create_symmetric_sparsity_pattern(dh) + M = create_sparsity_pattern(dh, ch) assembler = start_assemble(M) Me = zeros(n, n) @@ -190,6 +196,7 @@ function _project(vars, proj::L2Projector, fe_values::AbstractValues, M::Integer get_data(x::AbstractTensor, i) = x.data[i] get_data(x::Number, i) = x + ch = proj.ch ## Assemble contributions from each cell for (ic,cellnum) in enumerate(proj.set) @@ -216,8 +223,16 @@ function _project(vars, proj::L2Projector, fe_values::AbstractValues, M::Integer end end - # solve for the projected nodal values - projected_vals = proj.M_cholesky \ f + # Correctly apply affine constraints + projected_vals = similar(f) + rhsdata = get_rhs_data(ch,proj.M_cholesky) + apply!(proj.M_cholesky,ch) + for (i,col) in enumerate(eachcol(f)) + apply_rhs!(rhsdata, col, ch) + u = proj.M_cholesky \ col + apply!(u, ch) + projected_vals[:,i] = u + end # Recast to original input type make_T(vals) = T <: AbstractTensor ? T(Tuple(vals)) : vals[1] diff --git a/src/exports.jl b/src/exports.jl index 0e834f88a6..922a0a3aa6 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -101,6 +101,8 @@ export addcellset!, transform_coordinates!, generate_grid, +# AdaptiveGrid + ForestBWG, # Grid coloring create_coloring, @@ -130,6 +132,7 @@ export collect_periodic_facets!, PeriodicFacetPair, AffineConstraint, + ConformityConstraint, update!, apply!, apply_rhs!, diff --git a/test/runtests.jl b/test/runtests.jl index ba74624411..63980d7a0e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,37 +31,33 @@ else end end -include("test_utils.jl") - -# Unit tests -include("test_interpolations.jl") -include("test_cellvalues.jl") -include("test_facevalues.jl") -include("test_interfacevalues.jl") -include("test_quadrules.jl") -include("test_assemble.jl") -include("test_dofs.jl") -include("test_constraints.jl") -include("test_grid_dofhandler_vtk.jl") -include("test_vtk_export.jl") -include("test_abstractgrid.jl") -include("test_grid_addboundaryset.jl") -include("test_mixeddofhandler.jl") -include("test_l2_projection.jl") -include("test_pointevaluation.jl") -# include("test_notebooks.jl") -include("test_apply_rhs.jl") -include("test_apply_analytical.jl") -include("test_deprecations.jl") -HAS_EXTENSIONS && include("blockarrays.jl") -include("test_examples.jl") - -@test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined -#= See which is not defined if fails -for name in names(Ferrite) - isdefined(Ferrite, name) || @warn "Ferrite.$name is not defined but $name is exported" -end -=# - -# Integration tests -include("integration/test_simple_scalar_convergence.jl") +#include("test_utils.jl") +# +## Unit tests +#include("test_interpolations.jl") +#include("test_cellvalues.jl") +#include("test_facevalues.jl") +#include("test_interfacevalues.jl") +#include("test_quadrules.jl") +#include("test_assemble.jl") +#include("test_dofs.jl") +#include("test_constraints.jl") +#include("test_grid_dofhandler_vtk.jl") +#include("test_abstractgrid.jl") +#include("test_grid_addboundaryset.jl") +#include("test_mixeddofhandler.jl") +#include("test_l2_projection.jl") +#include("test_pointevaluation.jl") +## include("test_notebooks.jl") +#include("test_apply_rhs.jl") +#include("test_apply_analytical.jl") +#include("test_deprecations.jl") +#HAS_EXTENSIONS && include("blockarrays.jl") +#include("test_examples.jl") +#@test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined +# +## Integration tests +#include("integration/test_simple_scalar_convergence.jl") + +include("test_p4est.jl") +include("test_p4est_example.jl") diff --git a/test/test_p4est.jl b/test/test_p4est.jl new file mode 100644 index 0000000000..71f99757ee --- /dev/null +++ b/test/test_p4est.jl @@ -0,0 +1,1095 @@ +@testset "OctantBWG Lookup Tables" begin + @test Ferrite.AMR._face(1) == [3,5] + @test Ferrite.AMR._face(5) == [1,5] + @test Ferrite.AMR._face(12) == [2,4] + @test Ferrite.AMR._face(1,1) == 3 && Ferrite.AMR._face(1,2) == 5 + @test Ferrite.AMR._face(5,1) == 1 && Ferrite.AMR._face(5,2) == 5 + @test Ferrite.AMR._face(12,1) == 2 && Ferrite.AMR._face(12,2) == 4 + @test Ferrite.AMR._face(3,1) == 3 && Ferrite.AMR._face(3,2) == 6 + + @test Ferrite.AMR._face_edge_corners(1,1) == (0,0) + @test Ferrite.AMR._face_edge_corners(3,3) == (3,4) + @test Ferrite.AMR._face_edge_corners(8,6) == (2,4) + @test Ferrite.AMR._face_edge_corners(4,5) == (0,0) + @test Ferrite.AMR._face_edge_corners(5,4) == (0,0) + @test Ferrite.AMR._face_edge_corners(7,1) == (3,4) + @test Ferrite.AMR._face_edge_corners(11,1) == (2,4) + @test Ferrite.AMR._face_edge_corners(9,1) == (1,3) + @test Ferrite.AMR._face_edge_corners(10,2) == (1,3) + @test Ferrite.AMR._face_edge_corners(12,2) == (2,4) + + @test Ferrite.AMR.𝒱₃[1,:] == Ferrite.AMR.𝒰[1:4,1] == Ferrite.AMR._face_corners(3,1) + @test Ferrite.AMR.𝒱₃[2,:] == Ferrite.AMR.𝒰[1:4,2] == Ferrite.AMR._face_corners(3,2) + @test Ferrite.AMR.𝒱₃[3,:] == Ferrite.AMR.𝒰[5:8,1] == Ferrite.AMR._face_corners(3,3) + @test Ferrite.AMR.𝒱₃[4,:] == Ferrite.AMR.𝒰[5:8,2] == Ferrite.AMR._face_corners(3,4) + @test Ferrite.AMR.𝒱₃[5,:] == Ferrite.AMR.𝒰[9:12,1] == Ferrite.AMR._face_corners(3,5) + @test Ferrite.AMR.𝒱₃[6,:] == Ferrite.AMR.𝒰[9:12,2] == Ferrite.AMR._face_corners(3,6) + + @test Ferrite.AMR._edge_corners(1) == [1,2] + @test Ferrite.AMR._edge_corners(4) == [7,8] + @test Ferrite.AMR._edge_corners(12,2) == 8 + + #Test Figure 3a) of Burstedde, Wilcox, Ghattas [2011] + test_ξs = (1,2,3,4) + @test Ferrite.AMR._neighbor_corner.((1,),(2,),(1,),test_ξs) == test_ξs + #Test Figure 3b) + @test Ferrite.AMR._neighbor_corner.((3,),(5,),(3,),test_ξs) == (Ferrite.AMR.𝒫[5,:]...,) +end + +@testset "Index Permutation" begin + for i in 1:length(Ferrite.AMR.edge_perm) + @test i == Ferrite.AMR.edge_perm_inv[Ferrite.AMR.edge_perm[i]] + end + for i in 1:length(Ferrite.AMR.𝒱₂_perm) + @test i == Ferrite.AMR.𝒱₂_perm_inv[Ferrite.AMR.𝒱₂_perm[i]] + end + for i in 1:length(Ferrite.AMR.𝒱₃_perm) + @test i == Ferrite.AMR.𝒱₃_perm_inv[Ferrite.AMR.𝒱₃_perm[i]] + end + for i in 1:length(Ferrite.AMR.node_map₂) + @test i == Ferrite.AMR.node_map₂_inv[Ferrite.AMR.node_map₂[i]] + end + for i in 1:length(Ferrite.AMR.node_map₃) + @test i == Ferrite.AMR.node_map₃_inv[Ferrite.AMR.node_map₃[i]] + end +end + +@testset "OctantBWG Encoding" begin +# # Tests from Figure 3a) and 3b) of Burstedde et al + o = Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,21,3) + b = 3 + @test Ferrite.AMR.child_id(o,b) == 5 + @test Ferrite.AMR.child_id(Ferrite.AMR.parent(o,b),b) == 3 + @test Ferrite.AMR.parent(Ferrite.AMR.parent(o,b),b) == Ferrite.AMR.Ferrite.AMR.OctantBWG(3,0,1,b) + @test Ferrite.AMR.parent(Ferrite.AMR.parent(Ferrite.AMR.parent(o,b),b),b) == Ferrite.AMR.root(3) + o = Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,4,3) + @test Ferrite.AMR.child_id(o,b) == 4 + @test Ferrite.AMR.child_id(Ferrite.AMR.parent(o,b),b) == 1 + @test Ferrite.AMR.parent(Ferrite.AMR.parent(o,b),b) == Ferrite.AMR.Ferrite.AMR.OctantBWG(3,0,1,b) + @test Ferrite.AMR.parent(Ferrite.AMR.parent(Ferrite.AMR.parent(o,b),b),b) == Ferrite.AMR.root(3) + + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,1,1,3),3) == 1 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,1,2,3),3) == 2 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,1,3,3),3) == 3 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,1,4,3),3) == 4 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,2,1,3),3) == 1 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,1,3),3) == 1 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,2,3),3) == 2 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,3,3),3) == 3 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,4,3),3) == 4 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,16,3),3) == 8 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,24,3),3) == 8 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,64,3),3) == 8 + @test Ferrite.AMR.child_id(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,2,9,3),3) == 1 + #maxlevel = 10 takes too long + maxlevel = 6 + levels = collect(1:maxlevel) + morton_ids = [1:2^(2*l) for l in levels] + for (level,morton_range) in zip(levels,morton_ids) + for morton_id in morton_range + @test Int(Ferrite.AMR.morton(Ferrite.AMR.OctantBWG(2,level,morton_id,maxlevel),level,maxlevel)) == morton_id + end + end + morton_ids = [1:2^(3*l) for l in levels] + for (level,morton_range) in zip(levels,morton_ids) + for morton_id in morton_range + @test Int(Ferrite.AMR.morton(Ferrite.AMR.OctantBWG(3,level,morton_id,maxlevel),level,maxlevel)) == morton_id + end + end +end + +@testset "OctantBWG Operations" begin + o = Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(2,0,0)) + @test Ferrite.AMR.facet_neighbor(o,1,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,0,0)) + @test Ferrite.AMR.facet_neighbor(o,2,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(4,0,0)) + @test Ferrite.AMR.facet_neighbor(o,3,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(2,-2,0)) + @test Ferrite.AMR.facet_neighbor(o,4,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(2,2,0)) + @test Ferrite.AMR.facet_neighbor(o,5,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(2,0,-2)) + @test Ferrite.AMR.facet_neighbor(o,6,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(2,0,2)) + @test Ferrite.AMR.descendants(o,2) == (Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)), Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(3,1,1))) + @test Ferrite.AMR.descendants(o,3) == (Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(2,0,0)), Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(5,3,3))) + + o = Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,0,0)) + @test Ferrite.AMR.facet_neighbor(o,1,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(-2,0,0)) + @test Ferrite.AMR.facet_neighbor(o,2,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(2,0,0)) + @test Ferrite.AMR.facet_neighbor(o,3,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,-2,0)) + @test Ferrite.AMR.facet_neighbor(o,4,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,2,0)) + @test Ferrite.AMR.facet_neighbor(o,5,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,0,-2)) + @test Ferrite.AMR.facet_neighbor(o,6,2) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,0,2)) + o = Ferrite.AMR.Ferrite.AMR.OctantBWG(0,(0,0,0)) + @test Ferrite.AMR.descendants(o,2) == (Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)), Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(3,3,3))) + @test Ferrite.AMR.descendants(o,3) == (Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(0,0,0)), Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(7,7,7))) + + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),1,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,-2,-2)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),4,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,2,2)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),6,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,0,-2)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),9,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,-2,0)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),12,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,2,0)) + + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(0,0,0)),1,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(0,-2,-2)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(0,0,0)),12,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(3,(2,2,0)) + + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),1,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,-4,-4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),2,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,4,-4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),3,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,-4,4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),4,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,4,4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),5,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(-4,0,-4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),6,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,0,-4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),7,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(-4,0,4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),8,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,0,4)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),9,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(-4,-4,0)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),10,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,-4,0)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),11,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(-4,4,0)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,0,0)),12,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,4,0)) + + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,0,0)),1,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,-8,-8)) + @test Ferrite.AMR.edge_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(0,0,0)),12,4) == Ferrite.AMR.Ferrite.AMR.OctantBWG(1,(8,8,0)) + + @test Ferrite.AMR.corner_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),1,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,-2,-2)) + @test Ferrite.AMR.corner_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),4,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,2,-2)) + @test Ferrite.AMR.corner_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0,0)),8,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,2,2)) + + @test Ferrite.AMR.corner_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0)),1,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(0,-2)) + @test Ferrite.AMR.corner_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0)),2,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,-2)) + @test Ferrite.AMR.corner_neighbor(Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(2,0)),4,3) == Ferrite.AMR.Ferrite.AMR.OctantBWG(2,(4,2)) +end + +@testset "OctreeBWG Operations" begin + # maximum level == 3 + # Octant level 0 size == 2^3=8 + # Octant level 1 size == 2^3/2 = 4 + # Octant level 2 size == 2^3/2 = 2 + # Octant level 3 size == 2^3/2 = 1 + # test translation constructor + grid = generate_grid(Quadrilateral,(2,2)) + # Rotate face topologically + grid.cells[2] = Quadrilateral((grid.cells[2].nodes[2], grid.cells[2].nodes[3], grid.cells[2].nodes[4], grid.cells[2].nodes[1])) + # This is our root mesh + # x-----------x-----------x + # |4 4 3|4 4 3| + # | | | + # | ^ | ^ | + # |1 | 2|1 | 2| + # | +--> | +--> | + # | | | + # |1 3 2|1 3 2| + # x-----------x-----------x + # |4 4 3|3 2 2| + # | | | + # | ^ | ^ | + # |1 | 2|4 | 3| + # | +--> | <--+ | + # | | | + # |1 3 2|4 1 1| + # x-----------x-----------x + adaptive_grid = ForestBWG(grid,3) + for cell in adaptive_grid.cells + @test cell isa Ferrite.AMR.OctreeBWG + @test cell.leaves[1] == Ferrite.AMR.OctantBWG(2,0,1,cell.b) + end + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,2), adaptive_grid.cells[1].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(4,1), adaptive_grid.cells[3].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(3,2), adaptive_grid.cells[4].leaves[1]) == Ferrite.AMR.OctantBWG(0,(-8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,4), adaptive_grid.cells[3].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,-8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(4,3), adaptive_grid.cells[2].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(2,2), adaptive_grid.cells[4].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,-8)) + o = adaptive_grid.cells[1].leaves[1] + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,4), o) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), o) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), o) == Ferrite.AMR.OctantBWG(0,(0,-8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(4,1), o) == Ferrite.AMR.OctantBWG(0,(-8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(4,3), o) == Ferrite.AMR.OctantBWG(0,(0,-8)) + + grid_new = Ferrite.AMR.creategrid(adaptive_grid) + @test length(grid_new.nodes) == 9 + @test length(grid_new.conformity_info) == 0 + + grid.cells[4] = Quadrilateral((grid.cells[4].nodes[2], grid.cells[4].nodes[3], grid.cells[4].nodes[4], grid.cells[4].nodes[1])) + grid.cells[4] = Quadrilateral((grid.cells[4].nodes[2], grid.cells[4].nodes[3], grid.cells[4].nodes[4], grid.cells[4].nodes[1])) + # root mesh in Ferrite.AMR notation in p4est notation + # x-----------x-----------x x-----------x-----------x + # |4 3 3|2 1 1| |3 4 4|2 3 1| + # | | | | | | + # | ^ | <--+ | | ^ | <--+ | + # |4 | 2|2 | 4| |1 | 2|2 | 1| + # | +--> | v | | +--> | v | + # | | | | | | + # |1 1 2|3 3 4| |1 3 2|4 4 3| + # x-----------x-----------x x-----------x-----------x + # |4 3 3|3 2 2| |3 4 4|4 4 1| + # | | | | | | + # | ^ | ^ | | ^ | ^ | + # |4 | 2|3 | 1| |1 | 2|2 | 1| + # | +--> | <--+ | | +--> | <--+ | + # | | | | | | + # |1 1 2|4 4 1| |1 3 2|3 3 1| + # x-----------x-----------x x-----------x-----------x + adaptive_grid = ForestBWG(grid,3) + for cell in adaptive_grid.cells + @test cell isa Ferrite.AMR.OctreeBWG + @test cell.leaves[1] == Ferrite.AMR.OctantBWG(2,0,1,cell.b) + end + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,2), adaptive_grid.cells[1].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(4,2), adaptive_grid.cells[3].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(3,2), adaptive_grid.cells[4].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,4), adaptive_grid.cells[3].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,-8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(4,4), adaptive_grid.cells[2].leaves[1]) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(2,2), adaptive_grid.cells[4].leaves[1]) == Ferrite.AMR.OctantBWG(0,(0,8)) + + #@test Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(4,4), adaptive_grid.cells[1].leaves[1],false) == Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,4), adaptive_grid.cells[1].leaves[1], false) == Ferrite.AMR.OctantBWG(0,(8,8)) + #@test Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(3,2), adaptive_grid.cells[1].leaves[1], false) == Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(2,4), adaptive_grid.cells[1].leaves[1], false) == Ferrite.AMR.OctantBWG(0,(8,-8)) + #@test Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(4,4), adaptive_grid.cells[1].leaves[1],false) == Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,4), adaptive_grid.cells[1].leaves[1],false) == Ferrite.AMR.OctantBWG(0,(8,8)) + #@test Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(3,2), adaptive_grid.cells[1].leaves[1], false) == Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(2,4), adaptive_grid.cells[1].leaves[1], false) == Ferrite.AMR.OctantBWG(0,(8,-8)) + + o = adaptive_grid.cells[1].leaves[1] + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,4), o) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), o) == Ferrite.AMR.OctantBWG(0,(0,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), o) == Ferrite.AMR.OctantBWG(0,(0,-8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(4,2), o) == Ferrite.AMR.OctantBWG(0,(8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(4,4), o) == Ferrite.AMR.OctantBWG(0,(0,8)) + + + #simple first and second level refinement + # first case + # x-----------x-----------x + # | | | + # | | | + # | | | + # | | | + # | | | + # | | | + # | | | + # x-----x-----x-----------x + # | | | | + # | | | | + # | | | | + # x--x--x-----x | + # | | | | | + # x--x--x | | + # | | | | | + # x--x--x-----x-----------x + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + @test length(adaptive_grid.cells[1].leaves) == 4 + for (m,octant) in zip(1:4,adaptive_grid.cells[1].leaves) + @test octant == Ferrite.AMR.OctantBWG(2,1,m,adaptive_grid.cells[1].b) + end + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[5]) == Ferrite.AMR.OctantBWG(1,(0,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[7]) == Ferrite.AMR.OctantBWG(1,(4,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[6]) == Ferrite.AMR.OctantBWG(1,(0,-4)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[7]) == Ferrite.AMR.OctantBWG(1,(4,-4)) + + grid_new = Ferrite.AMR.creategrid(adaptive_grid) + @test length(grid_new.nodes) == 19 + @test length(grid_new.conformity_info) == 4 + + # octree holds now 3 first level and 4 second level + @test length(adaptive_grid.cells[1].leaves) == 7 + for (m,octant) in zip(1:4,adaptive_grid.cells[1].leaves) + @test octant == Ferrite.AMR.OctantBWG(2,2,m,adaptive_grid.cells[1].b) + end + + + # second case + # x-----------x-----------x + # | | | + # | | | + # | | | + # | | | + # | | | + # x-----x--x--x-----------x + # | | | | | + # | x--x--x | + # | | | | | + # x-----x--x--x | + # | | | | + # | | | | + # x-----x-----x-----------x + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + @test length(adaptive_grid.cells[1].leaves) == 7 + @test all(getproperty.(adaptive_grid.cells[1].leaves[1:3],:l) .== 1) + + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[2]) == Ferrite.AMR.OctantBWG(1,(0,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[5]) == Ferrite.AMR.OctantBWG(2,(4,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(2,4), adaptive_grid.cells[1].leaves[7]) == Ferrite.AMR.OctantBWG(2,(6,8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[3]) == Ferrite.AMR.OctantBWG(1,(0,-4)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[6]) == Ferrite.AMR.OctantBWG(2,(4,-2)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(3,3), adaptive_grid.cells[1].leaves[7]) == Ferrite.AMR.OctantBWG(2,(6,-2)) + + grid_new = Ferrite.AMR.creategrid(adaptive_grid) + @test length(grid_new.nodes) == 19 + @test length(grid_new.conformity_info) == 4 + + # more complex neighborhoods + grid = Ferrite.generate_simple_disc_grid(Quadrilateral, 6) + grid.cells[2] = Quadrilateral((grid.cells[2].nodes[2], grid.cells[2].nodes[3], grid.cells[2].nodes[4], grid.cells[2].nodes[1])) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[3],adaptive_grid.cells[3].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[5],adaptive_grid.cells[5].leaves[1]) + + grid_new = Ferrite.AMR.creategrid(adaptive_grid) + @test length(grid_new.nodes) == 23 + @test length(grid_new.conformity_info) == 4 + + ################################################################## + ####uniform refinement and coarsening for all cells and levels#### + ################################################################## + adaptive_grid = ForestBWG(grid,8) + for l in 1:8 + Ferrite.AMR.refine_all!(adaptive_grid,l) + for tree in adaptive_grid.cells + @test all(Ferrite.AMR.morton.(tree.leaves,l,8) == collect(1:2^(2*l))) + end + end + #check montonicity of ancestor_id + for tree in adaptive_grid.cells + ids = Ferrite.AMR.ancestor_id.(tree.leaves,(1,),(tree.b,)) + @test issorted(ids) + end + #now go back from finest to coarsest + for l in 7:-1:0 + Ferrite.AMR.coarsen_all!(adaptive_grid) + for tree in adaptive_grid.cells + @test all(Ferrite.AMR.morton.(tree.leaves,l,8) == collect(1:2^(2*l))) + end + end + ######################### + # now do the same with 3D + # some ascii picasso can insert here something beautiful + ######################### + # TODO add some test with higher refinement level which failed in my REPl (I think 8 should fail) + # TODO add some rotation and more elaborate case + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + o = adaptive_grid.cells[1].leaves[1] + + # faces + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,2), o) == Ferrite.AMR.OctantBWG(0,(8,0,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,2), o) == Ferrite.AMR.OctantBWG(0,(-8,0,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,4), o) == Ferrite.AMR.OctantBWG(0,(0,8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,4), o) == Ferrite.AMR.OctantBWG(0,(0,-8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,6), o) == Ferrite.AMR.OctantBWG(0,(0,0,8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,6), o) == Ferrite.AMR.OctantBWG(0,(0,0,-8)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(8,1), o) == Ferrite.AMR.OctantBWG(0,(-8,0,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(8,1), o) == Ferrite.AMR.OctantBWG(0,(8,0,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(8,3), o) == Ferrite.AMR.OctantBWG(0,(0,-8,0)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(8,3), o) == Ferrite.AMR.OctantBWG(0,(0,8,0)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(8,5), o) == Ferrite.AMR.OctantBWG(0,(0,0,-8)) + @test Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(8,5), o) == Ferrite.AMR.OctantBWG(0,(0,0,8)) + + @test_throws BoundsError Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,1), o) + @test_throws BoundsError Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,1), o) + @test_throws BoundsError Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,3), o) + @test_throws BoundsError Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,3), o) + @test_throws BoundsError Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,5), o) + @test_throws BoundsError Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(1,5), o) + @test_throws BoundsError Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(8,2), o) + @test_throws BoundsError Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(8,2), o) + @test_throws BoundsError Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(8,4), o) + @test_throws BoundsError Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(8,4), o) + @test_throws BoundsError Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(8,6), o) + @test_throws BoundsError Ferrite.AMR.transform_facet_remote(adaptive_grid, FacetIndex(8,6), o) + + #corners + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,1), o, false) == Ferrite.AMR.OctantBWG(0,(-8,-8,-8)) + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,2), o, false) == Ferrite.AMR.OctantBWG(0,(8,-8,-8)) + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,3), o, false) == Ferrite.AMR.OctantBWG(0,(-8,8,-8)) + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,4), o, false) == Ferrite.AMR.OctantBWG(0,(8,8,-8)) + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,5), o, false) == Ferrite.AMR.OctantBWG(0,(-8,-8,8)) + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,6), o, false) == Ferrite.AMR.OctantBWG(0,(8,-8,8)) + @test_throws BoundsError Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,7), o, false) == Ferrite.AMR.OctantBWG(0,(-8,8,8)) + @test Ferrite.AMR.transform_corner(adaptive_grid, VertexIndex(1,8), o, false) == Ferrite.AMR.OctantBWG(0,(8,8,8)) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,1), o, false) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,2), o, false) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,3), o, false) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,4), o, false) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,5), o, false) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,6), o, false) + @test_throws BoundsError Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,7), o, false) + Ferrite.AMR.transform_corner_remote(adaptive_grid, VertexIndex(1,8), o, false) == Ferrite.AMR.OctantBWG(0,(-8,-8,-8)) + + #edges + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,1), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,2), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,3), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,1), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,2), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,3), o, false) + @test Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,4), o, false) == Ferrite.AMR.OctantBWG(0,(0,8,8)) + @test Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,4), o, false) == Ferrite.AMR.OctantBWG(0,(0,-8,-8)) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,5), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,6), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,7), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,5), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,6), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,7), o, false) + @test Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,8), o, false) == Ferrite.AMR.OctantBWG(0,(8,0,8)) + @test Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,8), o, false) == Ferrite.AMR.OctantBWG(0,(-8,0,-8)) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,9), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,10), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,11), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,9), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,10), o, false) + @test_throws BoundsError Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,11), o, false) + @test Ferrite.AMR.transform_edge(adaptive_grid, EdgeIndex(1,12), o, false) == Ferrite.AMR.OctantBWG(0,(8,8,0)) + @test Ferrite.AMR.transform_edge_remote(adaptive_grid, EdgeIndex(1,12), o, false) == Ferrite.AMR.OctantBWG(0,(-8,-8,0)) + + # Rotate three dimensional case + # This is our root mesh top view + # x-----------x-----------x + # |6 3 5|8 4 7| + # | | | + # | ^ | ^ | + # |2 | 1|1 | 2| + # | <--+ | +--> | + # | | | + # |7 4 8|5 3 6| + # x-----------x-----------x + # |8 4 7|8 4 7| + # | | | + # | ^ | ^ | + # |1 | 2|1 | 2| + # | +--> | +--> | + # | | | + # |5 3 6|5 3 6| + # x-----------x-----------x + grid = generate_grid(Hexahedron,(2,2,2)) + # Rotate face topologically as decscribed in the ascii picture above + grid.cells[7] = Hexahedron((grid.cells[7].nodes[2], grid.cells[7].nodes[3], grid.cells[7].nodes[4], grid.cells[7].nodes[1], grid.cells[7].nodes[4+2], grid.cells[7].nodes[4+3], grid.cells[7].nodes[4+4], grid.cells[7].nodes[4+1])) + grid.cells[7] = Hexahedron((grid.cells[7].nodes[2], grid.cells[7].nodes[3], grid.cells[7].nodes[4], grid.cells[7].nodes[1], grid.cells[7].nodes[4+2], grid.cells[7].nodes[4+3], grid.cells[7].nodes[4+4], grid.cells[7].nodes[4+1])) + adaptive_grid = ForestBWG(grid,3) + @test Ferrite.AMR.transform_corner(adaptive_grid,7,3,Ferrite.AMR.OctantBWG(0,(0,0,0)),false) == Ferrite.AMR.OctantBWG(0,(-8,8,-8)) + + #refinement + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + @test length(adaptive_grid.cells[1].leaves) == 8 + for (m,octant) in zip(1:8,adaptive_grid.cells[1].leaves) + @test octant == Ferrite.AMR.OctantBWG(3,1,m,adaptive_grid.cells[1].b) + end + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + @test length(adaptive_grid.cells[1].leaves) == 15 + for (m,octant) in zip(1:8,adaptive_grid.cells[1].leaves) + @test octant == Ferrite.AMR.OctantBWG(3,2,m,adaptive_grid.cells[1].b) + end + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + @test length(adaptive_grid.cells[1].leaves) == 15 + @test all(getproperty.(adaptive_grid.cells[1].leaves[1:3],:l) .== 1) + @test all(getproperty.(adaptive_grid.cells[1].leaves[4:11],:l) .== 2) + @test all(getproperty.(adaptive_grid.cells[1].leaves[12:end],:l) .== 1) + adaptive_grid = ForestBWG(grid,5) + #go from coarsest to finest uniformly + for l in 1:5 + Ferrite.AMR.refine_all!(adaptive_grid,l) + for tree in adaptive_grid.cells + @test all(Ferrite.AMR.morton.(tree.leaves,l,5) == collect(1:2^(3*l))) + end + end + #now go back from finest to coarsest + for l in 4:-1:0 + Ferrite.AMR.coarsen_all!(adaptive_grid) + for tree in adaptive_grid.cells + @test all(Ferrite.AMR.morton.(tree.leaves,l,5) == collect(1:2^(3*l))) + end + end + + # Single + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[8]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + # Unrefined grid has 5 ^ dim nodes and the refined element introduces 6 face center, 12 edge center and 1 volume center nodes + @test length(transfered_grid.nodes) == 5^3 + (6 + 12 + 1) + # 6 faces and 12 edges of the single refined element induces one hanging node each + @test length(transfered_grid.conformity_info) == 6 + 12 + + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + # Unrefined grid has 5 ^ dim nodes and the refined element introduces 6 face center, 12 edge center and 1 volume center nodes + @test length(transfered_grid.nodes) == 5^3 + (6 + 12 + 1) + # 6 faces and 12 edges of the single refined element induces one hanging node each - minus 3 faces and 3 edges on the outer boundary + @test length(transfered_grid.conformity_info) == 6 + 12- 2*3 + + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[8],adaptive_grid.cells[8].leaves[8]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + # Unrefined grid has 5 ^ dim nodes and the refined element introduces 6 face center, 12 edge center and 1 volume center nodes + @test length(transfered_grid.nodes) == 5^3 + (6 + 12 + 1) + # 6 faces and 12 edges of the single refined element induces one hanging node each - minus 3 faces and 3 edges on the outer boundary + @test length(transfered_grid.conformity_info) == 6 + 12 - 2*3 + + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[8],adaptive_grid.cells[8].leaves[1]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + # Unrefined grid has 5 ^ dim nodes and the refined element introduces 6 face center, 12 edge center and 1 volume center nodes + @test length(transfered_grid.nodes) == 5^3 + (6 + 12 + 1) + # 6 faces and 12 edges of the single refined element induces one hanging node each + @test length(transfered_grid.conformity_info) == 6 + 12 + + # Combined + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[8]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + @test length(transfered_grid.nodes) == 5^3 + 2*(6 + 12 + 1) + @test length(transfered_grid.conformity_info) == 2*(6 + 12) - 2*3 + + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[8],adaptive_grid.cells[8].leaves[8]) + Ferrite.AMR.refine!(adaptive_grid.cells[8],adaptive_grid.cells[8].leaves[1]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + @test length(transfered_grid.nodes) == 5^3 + 2*(6 + 12 + 1) + @test length(transfered_grid.conformity_info) == 2*(6 + 12) - 2*3 + + # Combined + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[8]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[8],adaptive_grid.cells[8].leaves[8]) + Ferrite.AMR.refine!(adaptive_grid.cells[8],adaptive_grid.cells[8].leaves[1]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + @test length(transfered_grid.nodes) == 5^3 + 4*(6 + 12 + 1) + @test length(transfered_grid.conformity_info) == 4*(6 + 12) - 2*3 - 2*3 + + # Combined and not rotated + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[8]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[6],adaptive_grid.cells[6].leaves[6]) + Ferrite.AMR.refine!(adaptive_grid.cells[6],adaptive_grid.cells[6].leaves[3]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + # +5^3 on the coarse grid + # +4 refined elements a 6 face nodes, 12 edge nodes and 1 volume nodes + # -1 shared node between tree 1 and 6 + @test length(transfered_grid.nodes) == 5^3 + 4*(6 + 12 + 1) - 1 + # 4*(6 + 12) potential hanging nodes + # - 2 shared through common edge + # - 2* (2*3) outer boundary face and edge nodes + @test length(transfered_grid.conformity_info) == 4*(6 + 12) - 2 - 2*3 - 2*3 + + # Combined and rotated + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[8]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[6]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[3]) + transfered_grid = Ferrite.creategrid(adaptive_grid) + @test unique(transfered_grid.nodes) == transfered_grid.nodes + # +5^3 on the coarse grid + # +4 refined elements a 6 face nodes, 12 edge nodes and 1 volume nodes + # -1 shared node between tree 1 and 7 + @test length(transfered_grid.nodes) == 5^3 + 4*(6 + 12 + 1) - 1 + # 4*(6 + 12) potential hanging nodes + # - 2 shared through common edge + # - 2* (2*3) outer boundary face and edge nodes + @test length(transfered_grid.conformity_info) == 4*(6 + 12) - 2 - 2*3 - 2*3 + + # Reproducer test for Fig.3 BWG 11 + grid = generate_grid(Hexahedron,(2,1,1)) + # (a) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[2],adaptive_grid.cells[2].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[2],adaptive_grid.cells[2].leaves[3]) + @test adaptive_grid.cells[2].leaves[3+4] == Ferrite.AMR.OctantBWG(2,(0,4,2)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,2), adaptive_grid.cells[2].leaves[3+4]) == Ferrite.AMR.OctantBWG(2,(8,4,2)) + # (b) Rotate elements topologically + grid.cells[1] = Hexahedron((grid.cells[1].nodes[2], grid.cells[1].nodes[3], grid.cells[1].nodes[4], grid.cells[1].nodes[1], grid.cells[1].nodes[6], grid.cells[1].nodes[7], grid.cells[1].nodes[8], grid.cells[1].nodes[5])) + grid.cells[2] = Hexahedron((grid.cells[2].nodes[4], grid.cells[2].nodes[1], grid.cells[2].nodes[2], grid.cells[2].nodes[3], grid.cells[2].nodes[8], grid.cells[2].nodes[5], grid.cells[2].nodes[6], grid.cells[2].nodes[7])) + # grid.cells[2] = Hexahedron((grid.cells[2].nodes[1], grid.cells[2].nodes[3], grid.cells[2].nodes[4], grid.cells[2].nodes[8], grid.cells[2].nodes[6], grid.cells[2].nodes[2], grid.cells[2].nodes[7], grid.cells[2].nodes[5])) How to rotate along diagonal? :) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[2],adaptive_grid.cells[2].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[2],adaptive_grid.cells[2].leaves[1]) + @test adaptive_grid.cells[2].leaves[6] == Ferrite.AMR.OctantBWG(2,(2,0,2)) + @test Ferrite.AMR.transform_facet(adaptive_grid, FacetIndex(1,3), adaptive_grid.cells[2].leaves[6]) == Ferrite.AMR.OctantBWG(2,(4,-2,2)) +end + +@testset "ForestBWG AbstractGrid Interfacing" begin + maxlevel = 3 + grid = generate_grid(Quadrilateral,(2,2)) + adaptive_grid = ForestBWG(grid,maxlevel) + for l in 1:maxlevel + Ferrite.AMR.refine_all!(adaptive_grid,l) + @test getncells(adaptive_grid) == 2^(2*l) * 4 == length(getcells(adaptive_grid)) + end +end + +@testset "Balancing" begin + #2D cases + #simple one quad with one additional non-allowed non-conformity level + grid = generate_grid(Quadrilateral,(1,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[6]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[6]) + balanced = Ferrite.AMR.balancetree(adaptive_grid.cells[1]) + @test length(balanced.leaves) == 16 + + #more complex non-conformity level 3 and 4 that needs to be balanced + adaptive_grid = ForestBWG(grid,5) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[7]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[12]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[12]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[15]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[16]) + balanced = Ferrite.AMR.balancetree(adaptive_grid.cells[1]) + @test length(balanced.leaves) == 64 + + grid = generate_grid(Quadrilateral,(2,1)) + adaptive_grid = ForestBWG(grid,2) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 11 + + grid = generate_grid(Quadrilateral,(2,2)) + adaptive_grid = ForestBWG(grid,2) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 19 + + # 2D example with balancing over a corner connection that is not within the topology tables + grid = generate_grid(Quadrilateral,(2,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[5]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 23 + + #corner balance case but rotated + grid = generate_grid(Quadrilateral,(2,1)) + grid.cells[1] = Quadrilateral((grid.cells[1].nodes[2], grid.cells[1].nodes[3], grid.cells[1].nodes[4], grid.cells[1].nodes[1])) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 23 + + # 3D case intra treee simple test, non conformity level 2 + grid = generate_grid(Hexahedron,(1,1,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[6]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[6]) + balanced = Ferrite.AMR.balancetree(adaptive_grid.cells[1]) + @test length(balanced.leaves) == 43 + + #3D case intra tree non conformity level 3 at two different places + adaptive_grid = ForestBWG(grid,4) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[7]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[12]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[28]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[29]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[37]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[39]) + balanced = Ferrite.AMR.balancetree(adaptive_grid.cells[1]) + @test length(balanced.leaves) == 127 + + #3D case inter tree non conformity level 3 at two different places + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,4) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + #Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[7]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[1]) + #Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + transfered_grid_ref = Ferrite.AMR.creategrid(adaptive_grid) + + # Rotate three dimensional case + grid = generate_grid(Hexahedron,(2,2,2)) + # This is our root mesh top view + # x-----------x-----------x + # |7 2 6|8 4 7| + # | | | + # | ^ | ^ | + # |4 | 3|1 | 2| + # | <--+ | +--> | + # | | | + # |8 1 5|5 3 6| + # x-----------x-----------x + # |8 4 7|8 4 7| + # | | | + # | ^ | ^ | + # |1 | 2|1 | 2| + # | +--> | +--> | + # | | | + # |5 3 6|5 3 6| + # x-----------x-----------x + # Rotate face topologically + grid.cells[7] = Hexahedron((grid.cells[7].nodes[2], grid.cells[7].nodes[3], grid.cells[7].nodes[4], grid.cells[7].nodes[1], grid.cells[7].nodes[4+2], grid.cells[7].nodes[4+3], grid.cells[7].nodes[4+4], grid.cells[7].nodes[4+1])) + grid.cells[7] = Hexahedron((grid.cells[7].nodes[2], grid.cells[7].nodes[3], grid.cells[7].nodes[4], grid.cells[7].nodes[1], grid.cells[7].nodes[4+2], grid.cells[7].nodes[4+3], grid.cells[7].nodes[4+4], grid.cells[7].nodes[4+1])) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[4]) + #Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[7]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == length(transfered_grid_ref.cells) + @test length(transfered_grid.cells) == 92 + + # edge balancing for new introduced connection that is not within topology table + grid = generate_grid(Hexahedron, (2,1,1)); + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid, [1,2]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid, [4]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid, [5]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 51 + + #another edge balancing case + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid,1) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid,[2,4,6,8]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid,34) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 134 + + #yet another edge balancing case + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid,1) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid,[2,4,6,8]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid,30) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 120 + + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 15 + + #yet another edge balancing case + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 43 + + #yet another edge balancing case + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid.cells[3],adaptive_grid.cells[3].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[10]) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[3]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 71 + + #yet another edge balancing case + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[1]) + Ferrite.AMR.balanceforest!(adaptive_grid) + Ferrite.AMR.refine!(adaptive_grid.cells[4],adaptive_grid.cells[4].leaves[7]) + Ferrite.AMR.balanceforest!(adaptive_grid) + @test Ferrite.AMR.getncells(adaptive_grid) == 120 +end + +@testset "Materializing Grid" begin + ################################################# + ############ structured 2D examples ############# + ################################################# + + # 2D case with a single tree + grid = generate_grid(Quadrilateral,(1,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 10 + @test length(transfered_grid.nodes) == 19 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + + #2D case with four trees and somewhat refinement pattern + grid = generate_grid(Quadrilateral,(2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 22 + @test length(transfered_grid.nodes) == 35 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + + #more random refinement + grid = generate_grid(Quadrilateral,(3,3)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[3],adaptive_grid.cells[3].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[3],adaptive_grid.cells[3].leaves[2]) + Ferrite.AMR.refine!(adaptive_grid.cells[3],adaptive_grid.cells[3].leaves[3]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[3]) + Ferrite.AMR.refine!(adaptive_grid.cells[7],adaptive_grid.cells[7].leaves[5]) + Ferrite.AMR.refine!(adaptive_grid.cells[9],adaptive_grid.cells[9].leaves[end]) + Ferrite.AMR.refine!(adaptive_grid.cells[9],adaptive_grid.cells[9].leaves[end]) + Ferrite.AMR.refine!(adaptive_grid.cells[9],adaptive_grid.cells[9].leaves[end]) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 45 + @test length(transfered_grid.nodes) == 76 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + + ################################################# + ############ structured 3D examples ############# + ################################################# + + # 3D case with a single tree + grid = generate_grid(Hexahedron,(1,1,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 8+7+7 + @test length(transfered_grid.nodes) == 65 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + + # Test only Interoctree by face connection + grid = generate_grid(Hexahedron,(2,1,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 16 + @test length(transfered_grid.nodes) == 45 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + #rotate the case around + grid = generate_grid(Hexahedron,(1,2,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 16 + @test length(transfered_grid.nodes) == 45 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + grid = generate_grid(Hexahedron,(1,1,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 16 + @test length(transfered_grid.nodes) == 45 + @test unique(transfered_grid.nodes) == transfered_grid.nodes + + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 8^2 + @test length(transfered_grid.nodes) == 125 # 5 per edge + @test unique(transfered_grid.nodes) == transfered_grid.nodes + + # Rotate three dimensional case + grid = generate_grid(Hexahedron,(2,2,2)) + # Rotate face topologically + grid.cells[2] = Hexahedron((grid.cells[2].nodes[2], grid.cells[2].nodes[3], grid.cells[2].nodes[4], grid.cells[2].nodes[1], grid.cells[2].nodes[4+2], grid.cells[2].nodes[4+3], grid.cells[2].nodes[4+4], grid.cells[2].nodes[4+1])) + grid.cells[2] = Hexahedron((grid.cells[2].nodes[2], grid.cells[2].nodes[3], grid.cells[2].nodes[4], grid.cells[2].nodes[1], grid.cells[2].nodes[4+2], grid.cells[2].nodes[4+3], grid.cells[2].nodes[4+4], grid.cells[2].nodes[4+1])) + # This is our root mesh bottom view + # x-----------x-----------x + # |4 4 3|4 4 3| + # | | | + # | ^ | ^ | + # |1 | 2|1 | 2| + # | +--> | +--> | + # | | | + # |1 3 2|1 3 2| + # x-----------x-----------x + # |4 4 3|3 2 2| + # | | | + # | ^ | ^ | + # |1 | 2|4 | 3| + # | +--> | <--+ | + # | | | + # |1 3 2|4 1 1| + # x-----------x-----------x + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.cells) == 8^2 + @test length(transfered_grid.nodes) == 125 # 5 per edge + @test unique(transfered_grid.nodes) == transfered_grid.nodes + #TODO iterate over all rotated versions and check if det J > 0 +end + +@testset "hanging nodes" begin + #Easy Intraoctree + grid = generate_grid(Hexahedron,(1,1,1)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine_all!(adaptive_grid,1) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + # x-----------x-----------x + # | | | + # | | | + # | | | + # | | | + # | | | + # | | | + # | | | + # x-----x-----x-----------x + # | | | | + # | | | | + # | | | | + # x-----x-----x | + # | | | | + # | | | | + # | | | | + # x-----x-----x-----------x + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.conformity_info) == 12 + + # Easy Interoctree + grid = generate_grid(Hexahedron,(2,2,2)) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + # x-----------x-----------x + # | | | + # | | | + # | | | + # | | | + # | | | + # | | | + # | | | + # x-----x-----x-----------x + # | | | | + # | | | | + # | | | | + # x-----x-----x | + # | | | | + # | | | | + # | | | | + # x-----x-----x-----------x + transfered_grid = Ferrite.AMR.creategrid(adaptive_grid) + @test length(transfered_grid.conformity_info) == 12 + + #rotate the case from above in the first cell around + grid = generate_grid(Hexahedron,(2,2,2)) + # Rotate face topologically + grid.cells[1] = Hexahedron((grid.cells[1].nodes[2], grid.cells[1].nodes[3], grid.cells[1].nodes[4], grid.cells[1].nodes[1], grid.cells[1].nodes[4+2], grid.cells[1].nodes[4+3], grid.cells[1].nodes[4+4], grid.cells[1].nodes[4+1])) + grid.cells[1] = Hexahedron((grid.cells[1].nodes[2], grid.cells[1].nodes[3], grid.cells[1].nodes[4], grid.cells[1].nodes[1], grid.cells[1].nodes[4+2], grid.cells[1].nodes[4+3], grid.cells[1].nodes[4+4], grid.cells[1].nodes[4+1])) + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + transfered_grid_rotated = Ferrite.AMR.creategrid(adaptive_grid) + @test transfered_grid_rotated.conformity_info[5] == [28,25] + @test transfered_grid_rotated.conformity_info[20] == [10,15] + @test transfered_grid_rotated.conformity_info[30] == [15,18] + @test transfered_grid_rotated.conformity_info[1] == [2,18] + @test transfered_grid_rotated.conformity_info[19] == [16,28] + @test transfered_grid_rotated.conformity_info[22] == [10,15,2,18] + @test transfered_grid_rotated.conformity_info[41] == [10,16,2,28] + @test transfered_grid_rotated.conformity_info[43] == [18,25] + @test transfered_grid_rotated.conformity_info[11] == [10,2] + @test transfered_grid_rotated.conformity_info[36] == [2,28] + @test transfered_grid_rotated.conformity_info[40] == [10,16] + @test transfered_grid_rotated.conformity_info[38] == [2,18,28,25] + @test length(transfered_grid_rotated.conformity_info) == 12 + + #2D rotated case + grid = generate_grid(Quadrilateral,(2,2)) + # Rotate face topologically + grid.cells[2] = Quadrilateral((grid.cells[2].nodes[2], grid.cells[2].nodes[3], grid.cells[2].nodes[4], grid.cells[2].nodes[1])) + # This is our root mesh + # x-----------x-----------x + # |4 4 3|4 4 3| + # | | | + # | ^ | ^ | + # |1 | 2|1 | 2| + # | +--> | +--> | + # | | | + # |1 3 2|1 3 2| + # x-----------x-----------x + # |4 4 3|3 2 2| + # | | | + # | ^ | ^ | + # |1 | 2|4 | 3| + # | +--> | <--+ | + # | | | + # |1 3 2|4 1 1| + # x-----------x-----------x + adaptive_grid = ForestBWG(grid,3) + Ferrite.AMR.refine!(adaptive_grid.cells[2],adaptive_grid.cells[2].leaves[1]) + transfered_grid_rotated = Ferrite.AMR.creategrid(adaptive_grid) + @test transfered_grid_rotated.conformity_info[11] == [1,7] + @test transfered_grid_rotated.conformity_info[10] == [3,7] + + # multiple corner connections in 2D by disc discretization + grid = Ferrite.generate_simple_disc_grid(Quadrilateral,10) + adaptive_grid = ForestBWG(grid,3) + @test getncells(adaptive_grid) == 10 + Ferrite.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + Ferrite.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[3]) + @test getncells(adaptive_grid) == 16 + Ferrite.balanceforest!(adaptive_grid) + @test getncells(adaptive_grid) == 9*4 + 3 + 4 + + # multiple corner connections in 3D by cylinder discretization + grid = Ferrite.generate_simple_disc_grid(Hexahedron,10) + adaptive_grid = ForestBWG(grid,3) + @test getncells(adaptive_grid) == 10 + Ferrite.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[1]) + @test getncells(adaptive_grid) == 17 + Ferrite.refine!(adaptive_grid.cells[1],adaptive_grid.cells[1].leaves[3]) + @test getncells(adaptive_grid) == 24 + Ferrite.balanceforest!(adaptive_grid) + @test getncells(adaptive_grid) == 9*8 + 7 + 8 +end diff --git a/test/test_p4est_example.jl b/test/test_p4est_example.jl new file mode 100644 index 0000000000..7168ab9350 --- /dev/null +++ b/test/test_p4est_example.jl @@ -0,0 +1,197 @@ +using Ferrite, Test + +module ConvergenceTestHelper + +using Ferrite, SparseArrays, ForwardDiff, Test +import LinearAlgebra: diag + +get_geometry(::Ferrite.Interpolation{RefLine}) = Line +get_geometry(::Ferrite.Interpolation{RefQuadrilateral}) = Quadrilateral +get_geometry(::Ferrite.Interpolation{RefTriangle}) = Triangle +get_geometry(::Ferrite.Interpolation{RefPrism}) = Wedge +get_geometry(::Ferrite.Interpolation{RefHexahedron}) = Hexahedron +get_geometry(::Ferrite.Interpolation{RefTetrahedron}) = Tetrahedron + +get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = 2*order +get_quadrature_order(::Serendipity{shape, order}) where {shape, order} = 2*order +get_quadrature_order(::CrouzeixRaviart{shape, order}) where {shape, order} = 2*order+1 + +get_N(::Ferrite.Interpolation{shape, 1}) where {shape} = 19 +get_N(::Ferrite.Interpolation{shape, 2}) where {shape} = 12 +get_N(::Ferrite.Interpolation{shape, 3}) where {shape} = 8 +get_N(::Ferrite.Interpolation{shape, 4}) where {shape} = 5 +get_N(::Ferrite.Interpolation{shape, 5}) where {shape} = 3 + +analytical_solution(x) = cot(norm(x)) +analytical_rhs(x) = -sum(diag(ForwardDiff.hessian(analytical_solution,x))) + +# Standard assembly copy pasta for Poisson problem +function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues, coords) + n_basefuncs = getnbasefunctions(cellvalues) + ## Reset to 0 + fill!(Ke, 0) + fill!(fe, 0) + ## Loop over quadrature points + for q_point in 1:getnquadpoints(cellvalues) + ## Get the quadrature weight + dΩ = getdetJdV(cellvalues, q_point) + x = spatial_coordinate(cellvalues, q_point, coords) + ## Loop over test shape functions + for i in 1:n_basefuncs + δu = shape_value(cellvalues, q_point, i) + ∇δu = shape_gradient(cellvalues, q_point, i) + ## Add contribution to fe + fe[i] += analytical_rhs(x) * δu * dΩ + ## Loop over trial shape functions + for j in 1:n_basefuncs + ∇u = shape_gradient(cellvalues, q_point, j) + ## Add contribution to Ke + Ke[i, j] += (∇δu ⋅ ∇u) * dΩ + end + end + end + #@show norm(Ke), norm(fe) + return Ke, fe +end + +# Standard assembly copy pasta for Poisson problem +function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler) + ## Allocate the element stiffness matrix and element force vector + n_basefuncs = getnbasefunctions(cellvalues) + Ke = zeros(n_basefuncs, n_basefuncs) + fe = zeros(n_basefuncs) + ## Allocate global force vector f + f = zeros(ndofs(dh)) + ## Create an assembler + assembler = start_assemble(K, f) + ## Loop over all cels + for cell in CellIterator(dh) + ## Reinitialize cellvalues for this cell + reinit!(cellvalues, cell) + coords = getcoordinates(cell) + ## Compute element contribution + assemble_element!(Ke, fe, cellvalues, coords) + ## Assemble Ke and fe into K and f + assemble!(assembler, celldofs(cell), Ke, fe) + end + return K, f +end + +# Check L2 convergence +function check_and_compute_convergence(dh, u, cellvalues, testatol, forest) + L2norm = 0.0 + L∞norm = 0.0 + marked_cells = Int[] + for (cellid,cell) in enumerate(CellIterator(dh)) + L2loc = 0.0 + L∞loc = 0.0 + reinit!(cellvalues, cell) + n_basefuncs = getnbasefunctions(cellvalues) + coords = getcoordinates(cell) + uₑ = u[celldofs(cell)] + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + x = spatial_coordinate(cellvalues, q_point, coords) + uₐₙₐ = prod(cos, x*π/2) + uₐₚₚᵣₒₓ = function_value(cellvalues, q_point, uₑ) + L2norm += norm(uₐₙₐ-uₐₚₚᵣₒₓ)*dΩ + L∞norm = max(L∞norm, norm(uₐₙₐ-uₐₚₚᵣₒₓ)) + L2loc += norm(uₐₙₐ-uₐₚₚᵣₒₓ)*dΩ + L∞loc = max(L∞norm, norm(uₐₙₐ-uₐₚₚᵣₒₓ)) + #@test isapprox(uₐₙₐ, uₐₚₚᵣₒₓ; atol=testatol) + end + if L2loc > 2e-1 + push!(marked_cells,cellid) + end + end + Ferrite.refine!(forest,marked_cells) + Ferrite.balanceforest!(forest) + L2norm, L∞norm +end + +# Assemble and solve +function solve(dh, ch, cellvalues) + K, f = assemble_global(cellvalues, create_sparsity_pattern(dh,ch), dh); + apply!(K, f, ch) + u = K \ f; + apply!(u,ch) +end + +function setup_poisson_problem(grid, interpolation, interpolation_geo, qr, N) + # Construct Ferrite stuff + dh = DofHandler(grid) + add!(dh, :u, interpolation) + close!(dh); + + ch = ConstraintHandler(dh); + ∂Ω = union( + values(Ferrite.getfacetsets(grid))... + ); + dbc = Dirichlet(:u, ∂Ω, (x, t) -> analytical_solution(x)) + add!(ch, dbc); + add!(ch, ConformityConstraint(:u)) + close!(ch); + + cellvalues = CellValues(qr, interpolation, interpolation_geo); + + return dh, ch, cellvalues +end + +end # module ConvergenceTestHelper + +function compute_fluxes(cellvalues::CellValues{<:ScalarInterpolation}, dh::DofHandler, a::AbstractVector{T}) where T + + n = getnbasefunctions(cellvalues) + cell_dofs = zeros(Int, n) + nqp = getnquadpoints(cellvalues) + + # Allocate storage for the fluxes to store + q = [Vec{2,T}[] for _ in 1:getncells(dh.grid)] + + for (cell_num, cell) in enumerate(CellIterator(dh)) + q_cell = q[cell_num] + celldofs!(cell_dofs, dh, cell_num) + aᵉ = a[cell_dofs] + reinit!(cellvalues, cell) + + for q_point in 1:nqp + q_qp = - function_gradient(cellvalues, q_point, aᵉ) + push!(q_cell, q_qp) + end + end + return q +end + +@testset "convergence analysis" begin + @testset "$interpolation" for interpolation in ( + Lagrange{RefQuadrilateral, 1}(), + #Lagrange{RefHexahedron, 1}(), + ) + L2norm = Inf + # Generate a grid ... + geometry = ConvergenceTestHelper.get_geometry(interpolation) + interpolation_geo = interpolation + N = ConvergenceTestHelper.get_N(interpolation) + grid = generate_grid(geometry, ntuple(x->10, Ferrite.getrefdim(geometry))); + adaptive_grid = ForestBWG(grid,7) + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{Ferrite.getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + i = 0 + while L2norm > 1e-3 && i < 8 + grid_transfered = Ferrite.creategrid(adaptive_grid) + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid_transfered, interpolation, interpolation_geo, qr, N) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + # q_gp = compute_fluxes(cellvalues, dh, u); + # projector = L2Projector(interpolation, grid_transfered); + # q_projected = project(projector, q_gp, qr); + L2norm, _ = ConvergenceTestHelper.check_and_compute_convergence(dh, u, cellvalues, 1e-2, adaptive_grid) + # vtk_grid("p4est_test$(i).vtu",dh) do vtk + # vtk_point_data(vtk,dh,u) + # vtk_point_data(vtk, projector, q_projected, "q") + # end + i += 1 + end + end +end