Algebra of set#
The second core ingredient of Samurai is the possibility to construct a subset from sets of intervals. It is useful when you want to navigate between levels and implement operators used when two or more sets have interaction.
Imagine the following 1D mesh
And we want to make the projection of the cells from level 1 to level 0 by computing the mean of the values as illustrate by this figure
Using an algebra of set, we can describe the subset as the intersection of the cells of level 1 with the cells of level 0.
The following code performs the projection on this intersection
#include <iostream>
#include <samurai/cell_array.hpp>
#include <samurai/cell_list.hpp>
#include <samurai/field.hpp>
int main()
{
constexpr std::size_t dim = 1;
// Mesh creation
samurai::CellList<dim> cl;
cl[0][{}].add_interval({0, 4});
cl[1][{}].add_interval({0, 4});
cl[1][{}].add_interval({6, 8});
samurai::CellArray<dim> ca{cl};
// Initialize field u on this mesh
auto u = samurai::make_field<double, 1>("u", ca);
samurai::for_each_cell(ca,
[&](auto cell)
{
u[cell] = cell.indices[0];
});
std::cout << "before projection" << std::endl;
std::cout << u << std::endl;
// Make projection on the intersection
auto subset = samurai::intersection(ca[0], ca[1]).on(0);
subset(
[&](const auto& i, auto)
{
u(0, i) = 0.5 * (u(1, 2 * i) + u(1, 2 * i + 1));
});
std::cout << "after projection" << std::endl;
std::cout << u << std::endl;
return 0;
}
And the output is
before projection
Field u
level: 0 coords: { 0.5} value: 0.
level: 0 coords: { 1.5} value: 1.
level: 0 coords: { 2.5} value: 2.
level: 0 coords: { 3.5} value: 3.
level: 1 coords: { 0.25} value: 0.
level: 1 coords: { 0.75} value: 1.
level: 1 coords: { 1.25} value: 2.
level: 1 coords: { 1.75} value: 3.
level: 1 coords: { 3.25} value: 6.
level: 1 coords: { 3.75} value: 7.
after projection
Field u
level: 0 coords: { 0.5} value: 0.5
level: 0 coords: { 1.5} value: 2.5
level: 0 coords: { 2.5} value: 2.
level: 0 coords: { 3.5} value: 6.5
level: 1 coords: { 0.25} value: 0.
level: 1 coords: { 0.75} value: 1.
level: 1 coords: { 1.25} value: 2.
level: 1 coords: { 1.75} value: 3.
level: 1 coords: { 3.25} value: 6.
level: 1 coords: { 3.75} value: 7.
The relevant part in this code is these few lines
// Make projection on the intersection
auto subset = samurai::intersection(ca[0], ca[1]).on(0);
subset([&](const auto& i, auto)
{
u(0, i) = 0.5*(u(1, 2 * i) + u(1, 2 * i + 1));
});
There are two different parts:
the construction of the subset
auto subset = samurai::intersection(ca[0], ca[1]).on(0);
In Samurai, you already have an algebra of set operators implemented (intersection, difference, union, translation, …). And you can easily create your operator as we will see in the following. The subset can be made up of sets of different levels. Since the intervals of these sets are represented by integers, it is important to be able to compare these integers at the same level. By default, the common level is the finest (here level = 1). But we want to apply our operator on another level. This is why we specify it using .on(0).
the numerical computation on this subset
subset([&](const auto& i, auto)
{
u(0, i) = 0.5*(u(1, 2 * i) + u(1, 2 * i + 1));
});
If the subset exists, we can then apply a numerical kernel to it. There exist different ways: operator(), apply_op where you describe an operator which can be different depending on the problem dimension. The method using operator() as in our example takes a lambda function with two parameters: an x-interval and an index array with the values for the other dimensions. Since our example is a 1d problem, we omit the index.
How to find the subset ?#
The 1d algorithm with intervals at the same level is really simple to understand and can be illustrated by the following steps.
Initialization step
We want to create the subset of two sets A and B which is their intersection. Each set starts at the beginning of their first interval. We define two other variables: scan which will navigate in the sets and sentinel which will tell us that we are at the end of the set. At the beginning of the algorithm, scan is the minimum value of the starts of the first interval of each set and sentinel the maximum value of the ends of the last interval of each set + 1.
Step 1
We begin the algorithm to see if scan is in the current interval of each set. Here, scan is only in A. We want the intersection of A and B which means that scan must be in A and B. This is not the case here so we update scan: if the current value of a set is equal to scan, we move the current value to the next value (the end of the current interval or the beginning of a new one). Then, we update scan with the minimum value of the current value of each set.
Note that out is the possible result of this subset and needs to be populated to be valid.
Step 2
This time, scan is in A and B and thus the intersection expression is valid. We can start the possible result of the subset materialized by out. Therefore, out starts by 4. And we update the current value of each set and scan accordingly.
Step 3
scan is no more in A which means that we are outside the interval result. Since it is the first outside value encountered, we can close out with this one. We now have our first result: the intersection of A and B contains the interval [4, 5[. We can apply a function to it if defined like in your first example with the operator().
And we update again the current value of each interval and scan accordingly.
Step 4
We can continue the algorithm. Since out was populated previously, we have to begin a new interval result. We can see in this example that there are no more intervals that can be in the intersection of A and B. The current value of each interval will be updated until the sentinel value is reached. Then, scan will be equal to sentinel and the algorithm will stop.
We can observe that it is easy to construct more complex subset since we just have to modify the formula \(x_{A \cap B} = scan \in A \; \text{and} \; scan \in B\). For example, if we want the union, the formula becomes \(x_{A \cup B} = scan \in A \; \text{or} \; scan \in B\).
How to manage the level ?#
The sets used to define the subset can be at different levels. Since the mesh is described using integers, we need to project each set on the same level to be able to compare the intervals and the subset makes sense. This figure illustrates the purpose
On level 0, we have the interval \([1, 3[\) and on level 1 the interval \([0, 5[\). We want to compute the intersection of these two sets. Without projection, the result is \([1, 3[\) and we can see that it is wrong.
In Samurai, we have two ways to define the level of the result:
if we don’t mention it, then the biggest level is chosen.
we can mention it using the method on(level).
Example:
auto subset = intersection(setA, setB);
// or
auto subset = intersection(setA, setB).on(level);
With the first version, our example gives the result \([2, 4[\) since the default level is the biggest one (here 1). And the second version with on(0) gives the result \([1, 2[\).
Internally, it means that we have to shift the value of the intervals of a set before to compute scan and sentinel values. But when we have to know if scan is in the interval, we need to shift scan to go back to the level of the set.
The other important thing is if we use on(level) and the level is biggest than the biggest level of the sets composing the subset, then it is not necessary to project the intervals on level since we can reconstruct the result from the result interval founded from the biggest level of the sets. In our example, if we use on(3), the result is \([8, 16[\) which is nothing less than the result on the default level with a shift \([2 << s, 4 << s[\) with \(s=on_{level} - default_{level}\).
Recursion on the dimension#
Since the beginning, the algorithm is made on 1d problem. The problems manage by Samurai are not only 1d problem. Then the algorithm must work for the largest dimensions (2, 3, 4, …). Recall that the data structure used to store the mesh is LevelCellArray where we have an array of size dimension composed for each dimension by an array of intervals and an array of size dimension - 1 for the offsets.
x: [0, 4[@0, [0, 1[@4, [3, 4[@2, [0, 1[@6, [3, 4[@4, [0, 3[@8
y: [0, 4[@0
y-offset: [0, 1, 3, 5, 6]
z: [0, 1]@0
z-offset: [0, 1]
If we take independently each dimension, we can observe it is once again 1d problem with a list of intervals. It means that we can apply the algorithm previously defined beginning by the largest dimension and decrement the dimension d to d-1 when we have found a result interval on d until the dimension 0 is reached.