Skip to content

Commit

Permalink
start some devdocs
Browse files Browse the repository at this point in the history
  • Loading branch information
koehlerson committed May 9, 2024
1 parent c7f3c6f commit 359c821
Show file tree
Hide file tree
Showing 2 changed files with 181 additions and 1 deletion.
180 changes: 180 additions & 0 deletions docs/src/devdocs/AMR.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
# AMR

## P4est

All of it 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, where 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.

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.
This coordinate system goes from [0,2^b]^{dim}. The size of an octant is always 1 at the lowest possible level `b`.

### 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> oct = OctantBWG(2,0,1,3)
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._compute_size(3,0)
8
```
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 4.
```julia
julia> Ferrite._compute_size(3,1)
4
```
Construct all level 1 octants based on mortonid:
```julia
# note the arguments are dim,level,mortonid,maximumlevel
julia> oct = OctantBWG(2,1,1,3)
OctantBWG{2,4,4}
l = 1
xy = 0,0

julia> o = OctantBWG(2,1,2,3)
OctantBWG{2,4,4}
l = 1
xy = 4,0

julia> o = OctantBWG(2,1,3,3)
OctantBWG{2,4,4}
l = 1
xy = 0,4

julia> o = OctantBWG(2,1,4,3)
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 = OctantBWG(2,1,8,3)
ERROR: AssertionError: m (one(T) + one(T)) ^ (dim * l)
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 = OctantBWG(2,1,2,3)
OctantBWG{2,4,4}
l = 1
xy = 4,0
```

### Intraoctree operation

```@docs
Ferrite.corner_neighbor
Ferrite.edge_neighbor
Ferrite.face_neighbor
```

### Interoctree operation

```@docs
Ferrite.transform_corner
Ferrite.transform_edge
Ferrite.transform_face
```
2 changes: 1 addition & 1 deletion docs/src/devdocs/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
```

0 comments on commit 359c821

Please sign in to comment.