Interval and cartesian grid representation#
Introduction#
A cartesian grid is composed of several cells with a given edge length. The length can be the same for all the cells such as in the case of uniform grids. The core idea in Samurai is that a cartesian grid can be represented by intervals. To illustrate this purpose, we start with a simple 1D example. Consider
The whole domain is defined by the interval \([0, 13]\). It is made up of five cells which can also be interpreted as intervals (visiting the domain from left to right)
cell 1: \([0, 4]\) (green)
cell 2: \([4, 5]\) (red)
cell 3: \([5, 7]\) (blue)
cell 4: \([7, 9]\) (blue)
cell 5: \([9, 13]\) (green)
In this example, we also observe that several cells have the same width, so they can be collected in families of equal width (or resolution)
width of size 1: \([4, 5]\)
width of size 2: \([5, 9]\)
width of size 4: \([0, 4]\), \([9, 13]\)
Since we know the resolution, we can regroup the cells to construct contiguous intervals, still being able to “crop” as the were unmerged. Indeed, in the class of cells of width 2, we have two contiguous cells: \([5, 7]\) and \([7, 9]\) which form the interval \([5, 9]\). Therefore, since we know the resolution of this given interval, the two cells can be reconstructed.
If we plot the mesh over the initial domain using different levels of resolution, we obtain
In this example, we have chosen not to have intersections between cells at different levels. Still, this is not mandatory and we can imagine a mesh with overlapping regions on a given domain
Note
It is worthwhile observing that the construction of the cells in Samurai still has some constraints and it will not be possible to exactly construct the previously described domain. These constraints shall be explained in the next section.
Interval definition#
The data structure widely used by Samurai is the interval. An interval is described as follows
An interval is defined by its start and end values (left and right extrema values, in black). Moreover, we introduce two additional attributes
the index (in red), which is used to make the link (being an offset) with the actual structure where data are stored (more details in the sequel).
the step (in green), which represents the step utilized to navigate inside the interval it refers to.
Warning
In the previous introduction, we have not paid attention to the type of extrema values describing the interval. In Samurai, the start and the end of the interval are integers.
It is also important to notice that the end value of the interval is not included. Intervals are closed on the left and open on the right.
As in the introduction, we can handle multiple resolutions, which mean different cell sizes. The grid resolution is indexed and defined by a level. Thus, the size of a cell (edge length) is fixed by its level in the following way
At this stage of development, we take the size to be the same in all the directions, and, therefore, the cells are squares in 2D and cubes in 3D.
As emphasized in the previous remark, in Samurai, intervals are represented by integers instead that by real numbers as in the introduction. This is made possible by the fact that an interval is associated to some cells at a given level, whose size (a real number) can be simply obtained by (1). The figure below illustrates the idea
We have two cells \(0\) and \(1\) on the level \(l\). Therefore, we can write that the interval \([0, 2[\) describes this domain at the level \(l\).
Note
We could have chosen to describe closed intervals such as \([0, 1]\). Nevertheless, the fact of excluding the end of the interval has important consequences for the algebra of intervals (see Algebra of set).
Cell properties#
As hinted so far, given an interval and a level, we can reconstruct the correspoding cells. Indeed:
The number of cells contained in an interval is equal to the size of the interval.
The size of an interval is given by
The coordinates of the “minimal” corner (left in 1D, lower-left in 2D, lower-lower-left in 3D) of a cell is given by
where \(\Delta x\) is given by equation (1).
The coordinates of the center of a cell is given by
Constraints on the grid representation#
As stressed before, there are some constraints on the grid representation allowed by Samurai.
Since the extrema values for each interval (and thus the cells) are constrained to be integers and the corresponding real coordinates are reconstructed from these values via \(\Delta x\), the latter will take only certain values following this principle. Therefore, not every real coordinate can be represented by this system.
Suppose we are considering level 1 with associated step \(\Delta x = 0.5\) for a 1D problem. It is therefore impossible by definition to have a cell at this level with center at \(\frac{1}{3}\).
Note
Of course, we can imagine that in the future versions of Samurai, we could have some operator which transforming the grid represented by intervals onto cells of the real domain.
The remaining constraint is that a cell at the level \(l\) is included in a cell at a lower level.
Note
This property is important to perform mesh adaptation.
1D mesh example#
Let us now consider a 1D example with several levels and explain how the mesh is stored in Samurai.
For each level, the intervals are
level 0: \([0, 2[\), \([5, 6[\)
level 1: \([4, 7[\), \([8, 10[\)
level 2: \([14, 16[\)
The actual(real) intervals are given by the level and \(\Delta x\) defined by (1), thus yielding
level 0: \([0., 2.]\), \([5., 6.]\)
level 1: \([2., 3.5]\), \([4., 5.]\)
level 2: \([3.5, 4.]\)
Note
There are no overlapping regions in this example. The purpose is to make the example more understandable. Still, we shall see examples with overlaps in the following tutorials. Indeed, this latter case often happens when mesh adaptation is performed and ghost cells are needed to update the solution using a spatial operator with a stencil (such fluxes for Finite Volumes discretizations).
The following code implements this example using Samurai.
#include <iostream>
#include <samurai/cell_array.hpp>
#include <samurai/cell_list.hpp>
int main()
{
constexpr std::size_t dim = 1;
samurai::CellList<dim> cl;
cl[0][{}].add_interval({0, 2});
cl[0][{}].add_interval({5, 6});
cl[1][{}].add_interval({4, 7});
cl[1][{}].add_interval({8, 10});
cl[2][{}].add_interval({14, 16});
samurai::CellArray<dim> ca{cl};
std::cout << ca << std::endl;
return 0;
}
The associated output is
┌────────────────────┐
│ Level 0 │
└────────────────────┘
dim 0
cells = [0, 2[@0:1 [5, 6[@-3:1
┌────────────────────┐
│ Level 1 │
└────────────────────┘
dim 0
cells = [4, 7[@-1:1 [8, 10[@-2:1
┌────────────────────┐
│ Level 2 │
└────────────────────┘
dim 0
cells = [14, 16[@-6:1
The computation of the index values represented by the @ operator will be explained in a next section.
Two new data structures are used in this example samurai::CellList and samurai::CellArray which are C++ arrays of size max_level defined as a template parameter.
The default size is 16.
samurai::CellList is used to efficiently add new intervals when the mesh is constructed.
As its name suggests, samurai::CellList is nothing else than a list of intervals in the x-direction.
This list is stored in a map where the keys are the indices in the other dimensions (y, z, …) and the values are the list of intervals in the x-direction.
Note
We will give an example of a 2D case to better explain how the keys are constructed.
For 1D problems, the key obviously remains empty.
This is why we use {} in the construction of the samurai::CellList like in cl[0][{}].add_interval({ 0, 2});
Indeed, the samurai::CellList data structure is efficient at adding new elements (search, removal and insertion operations have logarithmic complexity for std::map).
On the other hand, when scientific computing algorithms, such as numerical schemes, come, it is crucial to efficiently loop over the cells without a search algorithm.
Moreover, we also want to apply algebraic operations on sets of intervals for a given dimension which means that a dimension should have its own representation by intervals in a compact writing.
To this end samurai::CellArray is precisely used to compress the representation of the mesh, obtain that each spatial dimension has its own interval list stored as an array.
In our 1D example, the samurai::CellList associated with this mesh is
level 0:
x: [0, 2[, [5, 6[
level 1:
x: [4, 7[, [8, 10[
level 2:
x: [14, 16[
and the samurai::CellArray is defined in the same fashion.
The only difference is that x is a std::forward_list of samurai::Interval in samurai::CellList and a std::vector of samurai::Interval in samurai::CellArray.
To understand the differences, we must address an example in 2D.
2D mesh example#
The example below is supposed to help the reader to better understand the subtle difference between samurai::CellList and samurai::CellArray.
Consider the following mesh
The samurai::CellList associated with this mesh is
level 0:
y: 0
x: [0, 4[
y: 1
x: [0, 1[, [3, 4[
y: 2
x: [0, 1[, [3, 4[
y: 3
x: [0, 3[
level 1:
y: 2
x: [2, 6[
y: 3
x: [2, 6[
y: 4
x: [2, 4[, [5, 6[
y: 5
x: [2, 6[
y: 6
x: [6, 8[
y: 7
x: [6, 7[
level 2:
y: 8
x: [8, 10[
y: 9
x: [8, 10[
y: 14
x: [14, 16[
y: 15
x: [14, 16[
The key of the map in samurai::CellList is the index in y and the value of the key is the list of intervals in the x-direction for this key (index).
On the other hand, the samurai::CellArray is
level 0:
x: [0, 4[, [0, 1[, [3, 4[, [0, 1[, [3, 4[, [0, 3[
y: [0, 4[@0
y-offset: [0, 1, 3, 5, 6]
level 1:
x: [2, 6[, [2, 6[, [2, 4[, [5, 6[, [2, 6[, [6, 8[, [6, 7[
y: [2, 8[@-2
y-offset: [0, 1, 3, 5, 6, 7, 8]
level 2:
x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
y: [8, 10[@-8, [14, 16[@-12
y-offset: [0, 1, 2, 3, 4]
How do we construct a samurai::CellArray from a samurai::CellList?
First, we concatenate the intervals met in the x-direction for each index y at a given level. Therefore, the x array at level 2 representing the intervals in the x-direction is
x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
Then, we try to construct intervals in the y-direction from the keys. For level 2, we have the keys y = 8, 9, 14, 15. Therefore, we can construct two intervals with consecutive elements: \([8, 10[\) and \([14, 16[\).
The compressed view of the samurai::CellList at level 2 is as follows
x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
y: [8, 10[, [14, 16[
Now, we have to connect each y entry to the corresponding intervals in the x-direction. To this end we utilize a new array called y-offset and the index of the y interval represented by the operator @.
Each y has one corresponding interval in the x-direction. The y-offset indicates for each y where are the corresponding intervals in the x-direction inside the array x.
for y = 8, there is one interval in x-direction,
for y = 9, there is one interval in x-direction,
for y = 14, there is one interval in x-direction,
for y = 15, there is one interval in x-direction.
Thus the y-offset is the array [0, 1, 2, 3, 4]. The size of this array is the number of y + 1 and indicates that \(y[i]\) corresponds to the intervals in the x-direction between \(y-offset[i]\) and \(y-offset[i+1]\) in the x array.
One point still has to be clarified: how many elements y have we already gone through to know where to look in the y-offsets? This is where the index comes into play. If we look at the y interval \([14, 16[\), we know that the corresponding index in y-offsets for y = 14 is y-offset[2]. To obtain the right index, we choose the index defined in the interval by the @ operator in order to have y + index is equal to the entry in the y-offset entry. Then for y = 14, if the index is equal to -12, we find y-offset[y + @index] = y-offset[14 - 12] = y-offset[2].
If the same operation is made to compute the y-offset and the index on each interval in the y-direction, we end up the corresponding samurai::CellArray.
level 0:
x: [0, 4[, [0, 1[, [3, 4[, [0, 1[, [3, 4[, [0, 3[
y: [0, 4[@0
y-offset: [0, 1, 3, 5, 6]
level 1:
x: [2, 6[, [2, 6[, [2, 4[, [5, 6[, [2, 6[, [6, 8[, [6, 7[
y: [2, 8[@-2
y-offset: [0, 1, 2, 4, 5, 6, 7]
level 2:
x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
y: [8, 10[@-8, [14, 16[@-12
y-offset: [0, 1, 2, 3, 4]
Note
The algorithm described here to compress the samurai::CellList to yield the samurai::CellArray for the 2D is a recursive algorithm and, thus, it can be easily used for other spatial dimensions.
Therefore, we are not obliged to stick with 1D, 2D, and 3D problems: we can construct a grid in higher dimensions.
The implementation of this example is
#include <iostream>
#include <samurai/cell_array.hpp>
#include <samurai/cell_list.hpp>
#include <samurai/hdf5.hpp>
int main()
{
constexpr std::size_t dim = 2;
samurai::CellList<dim> cl;
cl[0][{0}].add_interval({0, 4});
cl[0][{1}].add_interval({0, 1});
cl[0][{1}].add_interval({3, 4});
cl[0][{2}].add_interval({0, 1});
cl[0][{2}].add_interval({3, 4});
cl[0][{3}].add_interval({0, 3});
cl[1][{2}].add_interval({2, 6});
cl[1][{3}].add_interval({2, 6});
cl[1][{4}].add_interval({2, 4});
cl[1][{4}].add_interval({5, 6});
cl[1][{5}].add_interval({2, 6});
cl[1][{6}].add_interval({6, 8});
cl[1][{7}].add_interval({6, 7});
cl[2][{8}].add_interval({8, 10});
cl[2][{8}].add_interval({8, 10});
cl[2][{14}].add_interval({14, 16});
cl[2][{15}].add_interval({14, 16});
samurai::CellArray<dim> ca{cl};
std::cout << ca << std::endl;
samurai::save("2d_mesh_representation", ca);
return 0;
}
And the output is
┌────────────────────┐
│ Level 0 │
└────────────────────┘
dim 0
cells = [0, 4[@0:1 [0, 1[@4:1 [3, 4[@2:1 [0, 1[@6:1 [3, 4[@4:1 [0, 3[@8:1
dim 1
cells = [0, 4[@0:1
offsets = 0 1 3 5 6
┌────────────────────┐
│ Level 1 │
└────────────────────┘
dim 0
cells = [2, 6[@9:1 [2, 6[@13:1 [2, 4[@17:1 [5, 6[@16:1 [2, 6[@20:1 [6, 8[@20:1 [6, 7[@22:1
dim 1
cells = [2, 8[@-2:1
offsets = 0 1 2 4 5 6 7
┌────────────────────┐
│ Level 2 │
└────────────────────┘
dim 0
cells = [8, 10[@21:1 [8, 10[@23:1 [14, 16[@19:1 [14, 16[@21:1
dim 1
cells = [8, 10[@-8:1 [14, 16[@-12:1
offsets = 0 1 2 3 4
Build a grid from a box#
Since the beginning, we have used samurai::CellList to build samurai::CellArray.
Now we show that we can also easily initialize a samurai::CellArray at a given level with a uniform cartesian grid by defining a box.
The box can be 1D, 2D, or 3D.
The box can either define the cells involved in the grid using integers or in real coordinates.
The following example uses a box in real coordinates
#include <iostream>
#include <samurai/box.hpp>
#include <samurai/cell_array.hpp>
#include <samurai/hdf5.hpp>
int main()
{
constexpr std::size_t dim = 2;
std::size_t start_level = 3;
samurai::Box<double, dim> box({-1, -1}, {1, 1});
samurai::CellArray<dim> ca_box;
ca_box[start_level] = {start_level, box};
std::cout << ca_box << std::endl;
samurai::save("2d_mesh_box", ca_box);
return 0;
}
The box is defined by its “minimal” (lower-left in 2D) and “maximal” (upper-right in 2d) corners. In this example, the box is therefore \([-1, 1] \times [-1, 1]\). The space step is chosen from the given level which means \(\Delta x = 2^{-3} = 0.125\). The number of cells is defined by the length of the box and the space step.
┌────────────────────┐
│ Level 3 │
└────────────────────┘
dim 0
cells = [-8, 8[@8:1 [-8, 8[@24:1 [-8, 8[@40:1 [-8, 8[@56:1 [-8, 8[@72:1 [-8, 8[@88:1 [-8, 8[@104:1 [-8, 8[@120:1 [-8, 8[@136:1 [-8, 8[@152:1 [-8, 8[@168:1 [-8, 8[@184:1 [-8, 8[@200:1 [-8, 8[@216:1 [-8, 8[@232:1 [-8, 8[@248:1
dim 1
cells = [-8, 8[@8:1
offsets = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
We obtain the following mesh
Warning
Since the size of the cells is fixed at a given level and their coordinates are constrained by their integer representation, it is not always possible to build a box in real coordinates where the bounds of the box correspond to a corner point of a cell.