Template Numerical Library version\ main:481315e2
Loading...
Searching...
No Matches
Orthogonal grids

Table of Contents

Introduction

Grids are regular orthognal meshes. Similar to unstructured numerical meshes they provide indexing of mesh entites and express their adjacency. The difference, compared to the unstructured meshes, is that the adjacency of the mesh entities are not stored explicitly in the memory but the are computed on-the-fly. The interface of grids is as simillar as possible to the unstructured meshes but there are some differences. The main difference is that the mesh entities are given by their coordinates and orientation. The type and orientation of the entity is given by its basis and normals. Basis is a vector having one for axes, along which the entity has non-zero length, and zeros otherwise. Normals is a vector orthogonal to the basis vector, i.e. it has ones where basis vector has zeros and vice versa. The meaning of the normals vector is such that it is like a pack of all vectors of standart basis which are orthogonal to the grid entity. The following tables show all possible grid entities in 1D, 2D and 3D.

TODO: Bases and normals

Grid entities in 1D are as follows:

Entities in 1D Basis Normals Unpacked normal vectors
Cells ( 1 ) ( 0 ) N/a
Vertexes ( 0 ) ( 1 ) ( 1 )

Grid entities in 2D are as follows:

Entities in 2D Basis Normals Unpacked normal vectors
Cells ( 1, 1 ) ( 0, 0 ) N/A
Faces along x- axis ( 1, 0 ) ( 0, 1 ) ( 0, 1 )
Faces along y- axis ( 0, 1 ) ( 1, 0 ) ( 1, 0 )
Vertexes ( 0, 0 ) ( 1, 1 ) ( 1, 0 ), ( 0, 1 )

Grid entities in 3D are as follows:

Entities in 3D Basis Normals Unpacked normal vectors
Cells ( 1, 1, 1 ) ( 0, 0, 0 ) N/A
Faces along x- and y- axes ( 1, 1, 0 ) ( 0, 0, 1 ) ( 0, 0, 1 )
Faces along x- and z- axes ( 1, 0, 1 ) ( 0, 1, 0 ) ( 0, 1, 0 )
Faces along y- and z- axes ( 0, 1, 1 ) ( 1, 0, 0 ) ( 1, 0, 0 )
Edges along x-axis ( 1, 0, 0 ) ( 0, 1, 1 ) ( 0, 1, 0 ), ( 0, 0, 1 )
Edges along y-axis ( 0, 1, 0 ) ( 1, 0, 1 ) ( 1, 0, 0 ), ( 0, 0, 1 )
Edges along z-axis ( 0, 0, 1 ) ( 1, 1, 0 ) ( 1, 0, 0 ), ( 0, 1, 0 )
Vertexes ( 0, 0, 0 ) ( 1, 1, 1 ) ( 1, 0, 0 ), ( 0, 1, 0 ), ( 0, 0, 1 )

The grid entity stores the vector with packed normals, tha basis vector is always computed on the fly. So whenever it possible, using the the normals vector is preferred for better performance.

Remark: The entity orientation given by the normals or basis vector should be encoded staticaly in the type of the entity. This would make the implementation of the grid entities more efficient. Such implementation, however, requires suppport of the generic lambda function by the compiler. Since the CUDA compiler nvcc is not able to compile a code with the generic lambda functions we stay with current implementation which is not optimal. Therefore, in the future, the implementation of the grid entities may change.

The following figures show coordinates and indexing of the grid entities in 2D for demonstration. Indexing of cells looks as follows:

+-------+-------+-------+-------+-------+ +-------+-------+-------+-------+-------+
| (0,4) | (1,4) | (2,4) | (3,4) | (4,4) | | (20 ) | (21 ) | (22 ) | (23 ) | (24 ) |
+-------+-------+-------+-------+-------+ +-------+-------+-------+-------+-------+
| (0,3) | (1,3) | (2,3) | (3,3) | (4,3) | | (15 ) | (16 ) | (17 ) | (18 ) | (19 ) |
+-------+-------+-------+-------+-------+ +-------+-------+-------+-------+-------+
| (0,2) | (1,2) | (2,2) | (3,2) | (4,2) | | (10 ) | (11 ) | (12 ) | (13 ) | (14 ) |
+-------+-------+-------+-------+-------+ +-------+-------+-------+-------+-------+
| (0,1) | (1,1) | (2,1) | (3,1) | (4,1) | | ( 5 ) | ( 6 ) | ( 7 ) | ( 8 ) | ( 9 ) |
+-------+-------+-------+-------+-------+ +-------+-------+-------+-------+-------+
| (0,0) | (1,0) | (2,0) | (3,0) | (4,0) | | ( 0 ) | ( 1 ) | ( 2 ) | ( 3 ) | ( 4 ) |
+-------+-------+-------+-------+-------+ +-------+-------+-------+-------+-------+

Indexing of faces looks as:

+-(0,6)-+-(1,6)-+-(2,6)-+-(3,6)-+-(4,6)-+ +-( 30)-+-( 31)-+-( 32)-+-( 33)-+-( 34)-+
| | | | | | | | | | | |
(0,5) (1,5) (2,5) (3,5) (4,5) (5,5) ( 65) ( 66) ( 67) ( 68) ( 69) ( 70)
| | | | | | | | | | | |
+-(0,5)-+-(1,5)-+-(2,5)-+-(3,5)-+-(4,5)-+ +-( 25)-+-( 26)-+-( 27)-+-( 28)-+-( 29)-+
| | | | | | | | | | | |
(0,4) (1,4) (2,4) (3,4) (4,4) (5,4) ( 59) ( 60) ( 61) ( 62) ( 63) ( 64)
| | | | | | | | | | | |
+-(0,4)-+-(1,4)-+-(2,4)-+-(3,4)-+-(4,4)-+ +-( 20)-+-( 21)-+-( 22)-+-( 23)-+-( 24)-+
| | | | | | | | | | | |
(0,3) (1,3) (2,3) (3,3) (4,3) (5,3) ( 53) ( 54) ( 55) ( 56) ( 57) ( 58)
| | | | | | | | | | | |
+-(0,3)-+-(1,3)-+-(2,3)-+-(3,3)-+-(4,3)-+ +-( 15)-+-( 16)-+-( 17)-+-( 18)-+-( 19)-+
| | | | | | | | | | | |
(0,2) (1,2) (2,2) (3,2) (4,2) (5,2) ( 47) ( 48) ( 49) ( 50) ( 51) ( 52)
| | | | | | | | | | | |
+-(0,2)-+-(1,2)-+-(2,2)-+-(3,2)-+-(4,2)-+ +-( 10)-+-( 11)-+-( 12)-+-( 13)-+-( 14)-+
| | | | | | | | | | | |
(0,1) (1,1) (2,1) (3,1) (4,1) (5,1) ( 41) ( 42) ( 43) ( 44) ( 45) ( 46)
| | | | | | | | | | | |
+-(0,1)-+-(1,1)-+-(2,1)-+-(3,1)-+-(4,1)-+ +-( 5 )-+-( 6 )-+-( 7 )-+-( 8 )-+-( 9 )-+
| | | | | | | | | | | |
(0,0) (1,0) (2,0) (3,0) (4,0) (5,0) ( 35) ( 36) ( 37) ( 38) ( 39) ( 40)
| | | | | | | | | | | |
+-(0,0)-+-(1,0)-+-(2,0)-+-(3,0)-+-(4,0)-+ +-( 0 )-+-( 1 )-+-( 2 )-+-( 3 )-+-( 4 )-+

And indexing of vertexes looks as follows:

(0,5)--(1,5)--(2,5)--(3,5)--(4,5)--(5,5) ( 30)--( 31)--( 32)--( 33)--( 34)--( 35)
| | | | | | | | | | | |
(0,4)--(1,4)--(2,4)--(3,4)--(4,4)--(5,4) ( 24)--( 25)--( 26)--( 27)--( 28)--( 29)
| | | | | | | | | | | |
(0,3)--(1,3)--(2,3)--(3,3)--(4,3)--(5,3) ( 18)--( 19)--( 20)--( 21)--( 22)--( 23)
| | | | | | | | | | | |
(0,2)--(1,2)--(2,2)--(3,2)--(4,2)--(5,2) ( 12)--( 13)--( 14)--( 15)--( 16)--( 17)
| | | | | | | | | | | |
(0,1)--(1,1)--(2,1)--(3,1)--(4,1)--(5,1) ( 6)--( 7)--( 8)--( 9)--( 10)--( 11)
| | | | | | | | | | | |
(0,0)--(1,0)--(2,0)--(3,0)--(4,0)--(5,0) ( 0)--( 1)--( 2)--( 3)--( 4)--( 5)

Grid may have arbitrary dimension i.e. even higher than 3D. It is represented by the templated class TNL::Meshes::Grid which has the wollowing template parameters:

  • Dimension is dimension of the grid. This can be any interger value greater than zero.
  • Real is a precision of the arithmetics used by the grid. It is double by default.
  • Device is the device where the grid shall be allocated. Currently it can be either TNL::Devices::Host for CPU or TNL::Devices::Cuda for CUDA supporting GPUs. It is TNL::Devices::Host by default.
  • Index is a type being used for indexing. It is int by default.

Grid creation

The grid is defined by its dimension, domain covered by the grid and its resolution. The following example shows how to create a grid:

1#include <iostream>
2#include <TNL/Meshes/Grid.h>
3
4int main( int argc, char* argv[] )
5{
6 TNL::Meshes::Grid< 1 > grid_1D( 10 );
7 grid_1D.setDomain( {0.0}, {1.0} );
8 std::cout << grid_1D << std::endl;
9
10 TNL::Meshes::Grid< 2 > grid_2D( 10, 20 );
11 grid_2D.setDomain( {0.0, 0.0}, {1.0,2.0});
12 std::cout << grid_2D << std::endl;
13
14 TNL::Meshes::Grid< 3 > grid_3D( 10, 20, 30 );
15 grid_3D.setDomain( {0.0, 0.0, 0.0}, {1.0, 2.0, 3.0} );
16 std::cout << grid_3D << std::endl;
17
18 return EXIT_SUCCESS;
19}
TNL::Meshes::Grid
Orthogonal n-dimensional grid.
Definition Grid.h:33
std::endl
T endl(T... args)

Here we create set of grids with different dimension. For each grid we set different resolution along each axis (using the constructor of the grid) and different length along each axis (by calling method TNL::Meshes::Grid::setDomain).

The result looks as follows:

| Dimensions: [ 10 ] |
| Origin: [ 0 ] |
| Proportions: [ 1 ] |
| Space steps: [ 0.1 ] |
| Entities count with basis [ 0 ]: 11 |
| Entities count with basis [ 1 ]: 10 |
| Dimensions: [ 10, 20 ] |
| Origin: [ 0, 0 ] |
| Proportions: [ 1, 2 ] |
| Space steps: [ 0.1, 0.1 ] |
| Entities count with basis [ 0, 0 ]: 231 |
| Entities count with basis [ 1, 0 ]: 210 |
| Entities count with basis [ 0, 1 ]: 220 |
| Entities count with basis [ 1, 1 ]: 200 |
| Dimensions: [ 10, 20, 30 ] |
| Origin: [ 0, 0, 0 ] |
| Proportions: [ 1, 2, 3 ] |
| Space steps: [ 0.1, 0.1, 0.1 ] |
| Entities count with basis [ 0, 0, 0 ]: 7161 |
| Entities count with basis [ 1, 0, 0 ]: 6510 |
| Entities count with basis [ 0, 1, 0 ]: 6820 |
| Entities count with basis [ 0, 0, 1 ]: 6930 |
| Entities count with basis [ 1, 1, 0 ]: 6200 |
| Entities count with basis [ 1, 0, 1 ]: 6300 |
| Entities count with basis [ 0, 1, 1 ]: 6600 |
| Entities count with basis [ 1, 1, 1 ]: 6000 |

The following example shows creation of a grid independently on the grid dimension. The domain covered by the grid is \( [0,1]^d\) where \( d \) is the grid dimension. The resolution os the same along each axis. Tho code looks as follows:

1#include <iostream>
2#include <TNL/Meshes/Grid.h>
3
4template< int Dimension >
5void createGrid()
6{
7 using GridType = TNL::Meshes::Grid< Dimension >;
8 using CoordinatesType = typename GridType::CoordinatesType;
9 GridType grid( 10 );
10 CoordinatesType origin( 0.0 ), proportions( 1.0 );
11 grid.setDomain( origin, proportions );
12
13 std::cout << grid << std::endl;
14}
15
16int main( int argc, char* argv[] )
17{
18 std::cout << "Creating grid with dimension equal one." << std::endl;
19 createGrid< 1 >();
20
21 std::cout << "Creating grid with dimension equal two." << std::endl;
22 createGrid< 2 >();
23
24 std::cout << "Creating grid with dimension equal three." << std::endl;
25 createGrid< 3 >();
26
27 std::cout << "Creating grid with dimension equal four." << std::endl;
28 createGrid< 4 >();
29
30 return EXIT_SUCCESS;
31}

The result looks as follows:

Creating grid with dimension equal one.
| Dimensions: [ 10 ] |
| Origin: [ 0 ] |
| Proportions: [ 1 ] |
| Space steps: [ 0.1 ] |
| Entities count with basis [ 0 ]: 11 |
| Entities count with basis [ 1 ]: 10 |
Creating grid with dimension equal two.
| Dimensions: [ 10, 10 ] |
| Origin: [ 0, 0 ] |
| Proportions: [ 1, 1 ] |
| Space steps: [ 0.1, 0.1 ] |
| Entities count with basis [ 0, 0 ]: 121 |
| Entities count with basis [ 1, 0 ]: 110 |
| Entities count with basis [ 0, 1 ]: 110 |
| Entities count with basis [ 1, 1 ]: 100 |
Creating grid with dimension equal three.
| Dimensions: [ 10, 10, 10 ] |
| Origin: [ 0, 0, 0 ] |
| Proportions: [ 1, 1, 1 ] |
| Space steps: [ 0.1, 0.1, 0.1 ] |
| Entities count with basis [ 0, 0, 0 ]: 1331 |
| Entities count with basis [ 1, 0, 0 ]: 1210 |
| Entities count with basis [ 0, 1, 0 ]: 1210 |
| Entities count with basis [ 0, 0, 1 ]: 1210 |
| Entities count with basis [ 1, 1, 0 ]: 1100 |
| Entities count with basis [ 1, 0, 1 ]: 1100 |
| Entities count with basis [ 0, 1, 1 ]: 1100 |
| Entities count with basis [ 1, 1, 1 ]: 1000 |
Creating grid with dimension equal four.
| Dimensions: [ 10, 10, 10, 10 ] |
| Origin: [ 0, 0, 0, 0 ] |
| Proportions: [ 1, 1, 1, 1 ] |
| Space steps: [ 0.1, 0.1, 0.1, 0.1 ] |
| Entities count with basis [ 0, 0, 0, 0 ]: 14641 |
| Entities count with basis [ 1, 0, 0, 0 ]: 13310 |
| Entities count with basis [ 0, 1, 0, 0 ]: 13310 |
| Entities count with basis [ 0, 0, 1, 0 ]: 13310 |
| Entities count with basis [ 0, 0, 0, 1 ]: 13310 |
| Entities count with basis [ 1, 1, 0, 0 ]: 12100 |
| Entities count with basis [ 1, 0, 1, 0 ]: 12100 |
| Entities count with basis [ 1, 0, 0, 1 ]: 12100 |
| Entities count with basis [ 0, 1, 1, 0 ]: 12100 |
| Entities count with basis [ 0, 1, 0, 1 ]: 12100 |
| Entities count with basis [ 0, 0, 1, 1 ]: 12100 |
| Entities count with basis [ 1, 1, 1, 0 ]: 11000 |
| Entities count with basis [ 1, 1, 0, 1 ]: 11000 |
| Entities count with basis [ 1, 0, 1, 1 ]: 11000 |
| Entities count with basis [ 0, 1, 1, 1 ]: 11000 |
| Entities count with basis [ 1, 1, 1, 1 ]: 10000 |

Traversing the grid

The grid does not store any data it only provides only indexing of the grid entities. The indexes then serve for accessing data stored in an array or vector. The grid entities may be traversed in parallel as we show in the following example:

1#include <iostream>
2#include <iomanip>
3#include <TNL/Meshes/Grid.h>
4#include <TNL/Containers/Vector.h>
5
6template< typename Device >
7void traverseGrid()
8{
9 /***
10 * Define grid dimension and size.
11 */
12 static constexpr int Dimension = 2;
13 const int grid_size = 5;
14
15 /***
16 * Setup necessary type.
17 */
18 using GridType = TNL::Meshes::Grid< Dimension, double, Device >;
19 using CoordinatesType = typename GridType::CoordinatesType;
20 using PointType = typename GridType::PointType;
21 using VectorType = TNL::Containers::Vector< double, Device >;
22
23 /***
24 * Setup types of grid entities.
25 */
26 using GridCell = typename GridType::Cell;
27 using GridFace = typename GridType::Face;
28 using GridVertex = typename GridType::Vertex;
29
30 /***
31 * Create an instance of a grid.
32 */
33 GridType grid( grid_size );
34 PointType origin( 0.0 ), proportions( 1.0 );
35 grid.setDomain( origin, proportions );
36
37 /***
38 * Allocate vectors for values stored in particular grid entities.
39 */
40 VectorType cells( grid.template getEntitiesCount< Dimension >(), 0.0 );
41 VectorType faces( grid.template getEntitiesCount< Dimension - 1 >(), 0.0 );
42 VectorType vertexes( grid.template getEntitiesCount< 0 >(), 0.0 );
43
44 /***
45 * Prepare views for the data at the grid entities so that we can
46 * manipulate them in lambda functions runnig eventually on GPU.
47 */
48 auto cells_view = cells.getView();
49 auto faces_view = faces.getView();
50 auto vertexes_view = vertexes.getView();
51
52 /***
53 * Setup value of each cell to its index in the grid.
54 */
55 grid.template forAllEntities< Dimension >( [=] __cuda_callable__ ( const GridCell& cell ) mutable {
56 cells_view[ cell.getIndex() ] = cell.getIndex();
57 } );
58
59 /***
60 * Print values of all cells in the grid.
61 */
62 std::cout << "Values of cells .... " << std::endl;
63 for( int i = grid_size-1; i>= 0; i-- ) {
64 for( int j = 0; j < grid_size; j++ ) {
65 GridCell cell( grid, {j, i} );
66 auto idx = cell.getIndex();
67 std::cout << std::right << std::setw( 12 ) << cells.getElement( idx );
68 }
69 std::cout << std::endl;
70 }
71 std::cout << std::endl;
72
73 /***
74 * Setup values of all faces to an average value of its neighbour cells.
75 */
76 grid.template forAllEntities< Dimension - 1 >( [=] __cuda_callable__ ( const GridFace& face ) mutable {
77 const CoordinatesType normal = face.getNormals();
78 double sum = 0.0;
79 double count = 0.0;
80 if( TNL::all(greaterEqual( face.getCoordinates() - normal, 0 )) ) {
81 auto neighbour = face.template getNeighbourEntity< Dimension >( -normal );
82 sum += cells_view[ neighbour.getIndex() ];
83 count++;
84 }
85 if( TNL::all(less( face.getCoordinates(), face.getGrid().getDimensions() )) ) {
86 auto neighbour = face.template getNeighbourEntity< Dimension >( { 0, 0 } );
87 sum += cells_view[ neighbour.getIndex() ];
88 count++;
89 }
90 faces_view[ face.getIndex() ] = sum / count;
91 } );
92
93 /***
94 * Print values of all faces in the grid.
95 */
96 std::cout << "Values of faces ..." << std::endl;
97 for( int i = grid_size; i>= 0; i-- ) {
98 std::cout << std::right << std::setw( 6 ) << " ";
99 for( int j = 0; j < grid_size; j++ ) {
100 GridFace face( grid, {j, i}, {0,1} );
101 auto idx = face.getIndex();
102 std::cout << std::right << std::setw( 12 ) << faces.getElement( idx );
103 }
104 std::cout << std::endl;
105 if( i > 0 )
106 for( int j = 0; j <= grid_size; j++ ) {
107 GridFace face( grid, {j, i - 1}, { 1,0 } );
108 auto idx = face.getIndex();
109 std::cout << std::right << std::setw( 12 ) << faces.getElement( idx );
110 }
111 std::cout << std::endl;
112 }
113
114 /***
115 * Setup values of all vertexes to an average value of its neighbouring cells.
116 */
117 grid.template forAllEntities< 0 >( [=] __cuda_callable__ ( const GridVertex& vertex ) mutable {
118 double sum = 0.0;
119 double count = 0.0;
120 auto grid_dimensions = vertex.getGrid().getDimensions();
121 if( vertex.getCoordinates().x() > 0 && vertex.getCoordinates().y() > 0 ) {
122 auto neighbour = vertex.template getNeighbourEntity< Dimension >( { -1,-1 } );
123 sum += cells_view[ neighbour.getIndex() ];
124 count++;
125 }
126 if( vertex.getCoordinates().x() > 0 && vertex.getCoordinates().y() < grid_dimensions.y() ) {
127 auto neighbour = vertex.template getNeighbourEntity< Dimension >( { -1,0 } );
128 sum += cells_view[ neighbour.getIndex() ];
129 count++;
130 }
131 if( vertex.getCoordinates().x() < grid_dimensions.x() && vertex.getCoordinates().y() > 0 ) {
132 auto neighbour = vertex.template getNeighbourEntity< Dimension >( { 0,-1 } );
133 sum += cells_view[ neighbour.getIndex() ];
134 count++;
135 }
136 if( TNL::all(less( vertex.getCoordinates(), vertex.getGrid().getDimensions() )) ) {
137 auto neighbour = vertex.template getNeighbourEntity< Dimension >( {0,0} );
138 sum += cells_view[ neighbour.getIndex() ];
139 count++;
140 }
141 vertexes_view[ vertex.getIndex() ] = sum / count;
142 } );
143
144 /***
145 * Print values of all vertexes in the grid.
146 */
147 std::cout << "Values of vertexes .... " << std::endl;
148 for( int i = grid_size; i>= 0; i-- ) {
149 for( int j = 0; j <= grid_size; j++ ) {
150 GridVertex vertex( grid, {j, i} );
151 auto idx = vertex.getIndex();
152 std::cout << std::right << std::setw( 12 ) << vertexes.getElement( idx );
153 }
154 std::cout << std::endl;
155 }
156}
157
158int main( int argc, char* argv[] )
159{
160 std::cout << "Traversing grid on CPU..." << std::endl;
161 traverseGrid< TNL::Devices::Host >();
162
163#ifdef __CUDACC__
164 std::cout << "Traversing grid on CUDA GPU..." << std::endl;
165 traverseGrid< TNL::Devices::Cuda >();
166#endif
167 return EXIT_SUCCESS;
168}
__cuda_callable__
#define __cuda_callable__
Definition CudaCallable.h:22
TNL::Containers::Vector
Vector extends Array with algebraic operations.
Definition Vector.h:39
std::count
T count(T... args)
std::right
T right(T... args)
std::setw
T setw(T... args)

In this example we start with writting an index of each cell into the cell. Next, we fill each face with average number computed from the neighbour cells and at the end, we do the same with vertices. We also write values stored in particular grid entities to the console.

The whoel example consits of the following steps:

  1. Setting dimension and resoltion of the grid (lines 12 and 13).
  2. We define necessary types which we will need later. It is the grid type ( GridType on line 18), type for coordinates of the grid entities ( CoordinatesType line 19), type for the real-world coordinates (PointType on lines 20) and type of container for storing the values of particular grid entities (VectorType on line 21).
  3. We defiene types of the following grid entities - cell (GridCell, line 26), face (GridFace, line 27) and vertex (GridVertex, line 28).
  4. We create an instance of the grid (lines 33-35). The resolution of the grid, which equals 5x5, is set by the grid constructor. The domain covered by the grid is given by points origin and proportions (defined on the line 34) which are passed to the method TNL::Meshes::Grid::setDomain (line 35).
  5. Next we allocate vectors for storing of values in particular grid entities (lines 40-42). The number of grid entities are given by templated method TNL::Meshes::Grid::getEntitiesCount. The template parameter defines dimension of the grid entities the number of which we are asking for.
  6. In the next step, we prepare vector views (cells_view, faces_view and vertexes_view, lines 48-50) which we will need later in lambda functions.
  7. On the lines 55-57, we iterate over all grid cells using templated method TNL::Meshes::Grid::forAllEntities. The template parameter again says dimension of the grid entities we want to iterate over. The method takes one parameter which is a lambda function that is supposed to be evaluated for each cell. The lamda function receives one parameter, type of which is TNL::Meshes::GridEntity. This grid entity represents particular cells over which we iterate. The index of the cell is obtained by method TNL::Meshes::GridEntity::getIndex.
  8. Next we print the values of all cells. This must be done sequentially on the CPU and tehrefore we do not use parallal for (TNL::Meshes::Grid::forAllEntities) anymore. Now we iterate over all cells sequentialy row by row (for loops on lines 63 and 64). We create grid cell with given coordinates (line 65) and we ask for its global index (line 66) which we use for acces to the vector with cell values (cell). The result of this is a matrix of cell values printed on a console.
  9. On the lines (76-91), we iterate over all faces of the grid. Again we use the method TNL::Meshes::Grid::forAllEntities but now the template parameter telling the dimension of the entites is dimensions of the grid minus one (Dimension-1). The lambda function being performed for each face now gets one parameter which is a grid entity (TNL::Meshes::GridEntity) with dimensions equal to one. First of all we fetch the normal vector of the face (line 77) by calling a method TNL::Meshes::GridEntity::getNormals which in case of faces is realy a normal to the face (there is no more than one orthogonal normal to the face in the space having the same dimension as the grid). For horizontal faces, the normal is \(n_h=(0,1)\), and for vertical faces, the normal is \(n_v=(1,0)\). Note that if the resolution of the grid is \((N_x,N_y)\) then the number of horizontal faces along \(x\) and \(y\) axes is \((N_x,N_y+1)=(N_x,N_y)+n_h\) and the number of vertical faces along \(x\) and \(y\) axes is \((N_x+1,N_y)=(N_x,N_y)+n_v\). To compute the average number of values of the neighbour cells we need to get value from the cell in the direction opposite to the face normal. We first check if there is such a cell (if statement on the line 80) and if it is the case we fetch its global index (line 81), add its value to the sum (line 82) and then we increase variable count which is the number off values that contrtibuted to the average (line 83). Next we have to get value from the cell having the same coordinates as the face itself. We first check if there is such a cell (line 85) - note that the number of faces can be higher than the number of cells along the \(x\) axis for vertical faces and along \(y\) axis for horizontal faces. If there is such neighbour cell we fatch its index (line 86), add its value to variable sum (line 87) and increase variable count (line 88). Finally we compute the average value of the neighbout cells and store in the vector with values of faces at the position given by the global index of the current face (line 90).
  10. Next we print the values of the faces on the console. Line by line, we print first horizontal faces with row index equal five, followed by vertical faces with row index equal five, next horizontal faces with row index four and so on. As before, this must be done sequentially on the CPU so we do not use parallel for (TNL::Meshes::Grid::forAllEntities). So we iterate over all rows of the grid (line 97) and for each row we first print the horizontal faces (lines 99-103) and then the vertical ones (lines 106-110). Note that in both cases we create instance of grid face (lines 100 and 107) where the last parameter of the constructor defines the orientatio of the face - {0,1} is normal of horizontal faces (line 100) and {1,0} is normal of vertical faces.
  11. Finally, we iterate over all vertexes and compute average value of all neghbouring cells. The vertexes have no orientation so we do not need to care the normals. We only check the number of neighbour cells based on the vertex coordinates. To iterate over all vertexes in parallel we use the method TNL::Meshes::Grid::forAllEntities with entity dimension set to zero (line 117). The lambda function, we perform on each vertex, provides us parameter vertex which represents the vertex that we currently operate on. We just fetch the coordinates of the vertex (line 120) and then we check what are the neighbour cells (if statements on lines 121, 126, 131 and 136). For each such a cell we add its value to variable sum. Finaly we compute the average value and store it in the vertex (line 141).
  12. At the end we print values of all vertexes the same way as we did it with cells (lines 148-155).

The result looks as follows:

Traversing grid on CPU...
Values of cells ....
20 21 22 23 24
15 16 17 18 19
10 11 12 13 14
5 6 7 8 9
0 1 2 3 4
Values of faces ...
20 21 22 23 24
20 20.5 21.5 22.5 23.5 24
17.5 18.5 19.5 20.5 21.5
15 15.5 16.5 17.5 18.5 19
12.5 13.5 14.5 15.5 16.5
10 10.5 11.5 12.5 13.5 14
7.5 8.5 9.5 10.5 11.5
5 5.5 6.5 7.5 8.5 9
2.5 3.5 4.5 5.5 6.5
0 0.5 1.5 2.5 3.5 4
0 1 2 3 4
Values of vertexes ....
20 20.5 21.5 22.5 23.5 24
17.5 18 19 20 21 21.5
12.5 13 14 15 16 16.5
7.5 8 9 10 11 11.5
2.5 3 4 5 6 6.5
0 0.5 1.5 2.5 3.5 4
Traversing grid on CUDA GPU...
Values of cells ....
20 21 22 23 24
15 16 17 18 19
10 11 12 13 14
5 6 7 8 9
0 1 2 3 4
Values of faces ...
20 21 22 23 24
20 20.5 21.5 22.5 23.5 24
17.5 18.5 19.5 20.5 21.5
15 15.5 16.5 17.5 18.5 19
12.5 13.5 14.5 15.5 16.5
10 10.5 11.5 12.5 13.5 14
7.5 8.5 9.5 10.5 11.5
5 5.5 6.5 7.5 8.5 9
2.5 3.5 4.5 5.5 6.5
0 0.5 1.5 2.5 3.5 4
0 1 2 3 4
Values of vertexes ....
20 20.5 21.5 22.5 23.5 24
17.5 18 19 20 21 21.5
12.5 13 14 15 16 16.5
7.5 8 9 10 11 11.5
2.5 3 4 5 6 6.5
0 0.5 1.5 2.5 3.5 4

Writers

Writers help to export data linked with a grid to some of standart formats like VTK, VTI or Gnuplot. The following example shows the use of the grid writers:

1#include <iostream>
2#include <iomanip>
3#include <TNL/Meshes/Grid.h>
4#include <TNL/Meshes/Writers/VTIWriter.h>
5#include <TNL/Meshes/Writers/VTKWriter.h>
6#include <TNL/Meshes/Writers/GnuplotWriter.h>
7#include <TNL/Containers/Vector.h>
8
9template< typename Device >
10void writeGrid()
11{
12 /***
13 * Define grid dimension and size.
14 */
15 static constexpr int Dimension = 2;
16 const int grid_size = 5;
17
18 /***
19 * Setup necessary type.
20 */
21 using GridType = TNL::Meshes::Grid< Dimension, double, Device >;
22 using PointType = typename GridType::PointType;
23 using VectorType = TNL::Containers::Vector< double, Device >;
24
25 /***
26 * Setup types of grid entities.
27 */
28 using GridCell = typename GridType::Cell;
29 using GridVertex = typename GridType::Vertex;
30
31 /***
32 * Create an instance of a grid.
33 */
34 GridType grid( grid_size );
35 PointType origin( 0.0 ), proportions( 1.0 );
36 grid.setDomain( origin, proportions );
37
38 /***
39 * Allocate vectors for values stored in particular grid entities.
40 */
41 VectorType cells( grid.template getEntitiesCount< Dimension >(), 0.0 );
42 VectorType vertexes( grid.template getEntitiesCount< 0 >(), 0.0 );
43
44 /***
45 * Prepare views for the data at the grid entities so that we can
46 * manipulate them in lambda functions runnig eventually on GPU.
47 */
48 auto cells_view = cells.getView();
49 auto vertexes_view = vertexes.getView();
50
51 /***
52 * Setup value of each cell to its index in the grid.
53 */
54 grid.template forAllEntities< Dimension >( [=] __cuda_callable__ ( const GridCell& cell ) mutable {
55 cells_view[ cell.getIndex() ] = cell.getIndex();
56 } );
57
58 /***
59 * Write values of all cells in the grid into a file in VTI format.
60 */
61 TNL::String cells_file_name_vti( "GridExample-cells-values-" + TNL::getType( Device{}) + ".vti" );
62 std::cout << "Writing a file " << cells_file_name_vti << " ..." << std::endl;
63 std::fstream cells_file_vti;
64 cells_file_vti.open( cells_file_name_vti.getString(), std::ios::out );
65 TNL::Meshes::Writers::VTIWriter< GridType > cells_vti_writer( cells_file_vti );
66 cells_vti_writer.writeImageData( grid );
67 cells_vti_writer.writeCellData( cells, "cell-values");
68
69 /***
70 * Write values of all cells in the grid into a file in VTK format.
71 */
72 TNL::String cells_file_name_vtk( "GridExample-cells-values-" + TNL::getType( Device{}) + ".vtk" );
73 std::cout << "Writing a file " << cells_file_name_vtk << " ..." << std::endl;
74 std::fstream cells_file_vtk;
75 cells_file_vtk.open( cells_file_name_vtk.getString(), std::ios::out );
76 TNL::Meshes::Writers::VTKWriter< GridType > cells_vtk_writer( cells_file_vtk );
77 cells_vtk_writer.writeEntities( grid );
78 cells_vtk_writer.writeCellData( cells, "cell-values");
79
80 /***
81 * Write values of all cells in the grid into a file in Gnuplot format.
82 */
83 TNL::String cells_file_name_gplt( "GridExample-cells-values-" + TNL::getType( Device{}) + ".gplt" );
84 std::cout << "Writing a file " << cells_file_name_gplt << " ..." << std::endl;
85 std::fstream cells_file_gplt;
86 cells_file_gplt.open( cells_file_name_gplt.getString(), std::ios::out );
87 TNL::Meshes::Writers::GnuplotWriter< GridType > cells_gplt_writer( cells_file_gplt );
88 cells_gplt_writer.writeEntities( grid );
89 cells_gplt_writer.writeCellData( grid, cells, "cell-values");
90
91 /***
92 * Setup values of all vertexes to an average value of its neighbouring cells.
93 */
94 grid.template forAllEntities< 0 >( [=] __cuda_callable__ ( const GridVertex& vertex ) mutable {
95 double sum = 0.0;
96 double count = 0.0;
97 auto grid_dimensions = vertex.getGrid().getDimensions();
98 if( vertex.getCoordinates().x() > 0 && vertex.getCoordinates().y() > 0 ) {
99 auto neighbour = vertex.template getNeighbourEntity< Dimension >( { -1,-1 } );
100 sum += cells_view[ neighbour.getIndex() ];
101 count++;
102 }
103 if( vertex.getCoordinates().x() > 0 && vertex.getCoordinates().y() < grid_dimensions.y() ) {
104 auto neighbour = vertex.template getNeighbourEntity< Dimension >( { -1,0 } );
105 sum += cells_view[ neighbour.getIndex() ];
106 count++;
107 }
108 if( vertex.getCoordinates().x() < grid_dimensions.x() && vertex.getCoordinates().y() > 0 ) {
109 auto neighbour = vertex.template getNeighbourEntity< Dimension >( { 0,-1 } );
110 sum += cells_view[ neighbour.getIndex() ];
111 count++;
112 }
113 if( TNL::all(less( vertex.getCoordinates(), vertex.getGrid().getDimensions() )) ) {
114 auto neighbour = vertex.template getNeighbourEntity< Dimension >( {0,0} );
115 sum += cells_view[ neighbour.getIndex() ];
116 count++;
117 }
118 vertexes_view[ vertex.getIndex() ] = sum / count;
119 } );
120
121 /***
122 * Write values of all vertexes in the grid to a file in VTI format
123 */
124 TNL::String vertexes_file_name_vti( "GridExample-vertexes-values-" + TNL::getType( Device{} ) + ".vti" );
125 std::cout << "Writing a file " << vertexes_file_name_vti << " ..." << std::endl;
126 std::fstream vertexes_file_vti;
127 vertexes_file_vti.open( vertexes_file_name_vti.getString(), std::ios::out );
128 TNL::Meshes::Writers::VTIWriter< GridType > vertexes_vti_writer( vertexes_file_vti );
129 vertexes_vti_writer.writeImageData( grid );
130 vertexes_vti_writer.writePointData( vertexes, "vertexes-values" );
131
132 /***
133 * Write values of all vertexes in the grid to a file in VTK format
134 */
135 TNL::String vertexes_file_name_vtk( "GridExample-vertexes-values-" + TNL::getType( Device{} ) + ".vtk" );
136 std::cout << "Writing a file " << vertexes_file_name_vtk << " ..." << std::endl;
137 std::fstream vertexes_file_vtk;
138 vertexes_file_vtk.open( vertexes_file_name_vtk.getString(), std::ios::out );
139 TNL::Meshes::Writers::VTIWriter< GridType > vertexes_vtk_writer( vertexes_file_vtk );
140 vertexes_vtk_writer.writeEntities( grid );
141 vertexes_vtk_writer.writePointData( vertexes, "vertexes-values" );
142
143 /***
144 * Write values of all vertexes in the grid to a file in Gnuplot format
145 */
146 TNL::String vertexes_file_name_gplt( "GridExample-vertexes-values-" + TNL::getType( Device{} ) + ".gplt" );
147 std::cout << "Writing a file " << vertexes_file_name_gplt << " ..." << std::endl;
148 std::fstream vertexes_file_gplt;
149 vertexes_file_gplt.open( vertexes_file_name_gplt.getString(), std::ios::out );
150 TNL::Meshes::Writers::GnuplotWriter< GridType > vertexes_gplt_writer( vertexes_file_gplt );
151 vertexes_gplt_writer.writeEntities( grid );
152 vertexes_gplt_writer.writePointData( grid, vertexes, "vertexes-values" );
153}
154
155int main( int argc, char* argv[] )
156{
157 std::cout << "Traversing grid on CPU..." << std::endl;
158 writeGrid< TNL::Devices::Host >();
159
160#ifdef __CUDACC__
161 std::cout << "Traversing grid on CUDA GPU..." << std::endl;
162 writeGrid< TNL::Devices::Cuda >();
163#endif
164 return EXIT_SUCCESS;
165}
TNL::Meshes::Writers::GnuplotWriter
Writer of data linked with meshes into Gnuplot format.
Definition GnuplotWriter.h:21
TNL::Meshes::Writers::VTIWriter
Writer of data linked with meshes into VTI format.
Definition VTIWriter.h:24
TNL::Meshes::Writers::VTKWriter
Writer of data linked with meshes into VTK format.
Definition VTKWriter.h:23
TNL::String
Class for managing strings.
Definition String.h:33
TNL::getType
std::string getType()
Returns a human-readable string representation of given type.
Definition TypeInfo.h:75
std::fstream::open
T open(T... args)

It is a modification of the previous example. As before, we first set values linked with the grid cells and vertexes and then we write the into output files. The values linked with cells are exported to VTI format on the lines 63-69, to VTK format on the lines 74-80 and the Gnuplot format on lines. The writers have the same interface. They have constructors which require output stream as a parameter (TNL::Meshes::Writers::VTIWriter::VTIWriter, TNL::Meshes::Writers::VTKWriter::VTKWriter ). Next we call a method writeEntities (TNL::Meshes::Writers::VTIWriter::writeEntities, TNL::Meshes::Writers::VTKWriter::writeEntities) which exports the grid entities to the file. In the case of the Gnuplot format, this method just writes simple header to the file and does not need to be called (the method is mainly for the compatibility with other writers). Finaly we may export data linked with the grid cells using a method writeCellData (TNL::Meshes::Writers::VTIWriter::writeCellData, TNL::Meshes::Writers::VTKWriter::writeCellData).

In the same way, we can export data linked with vertexes (line 125-152).

Generated on Thu Sep 7 2023 14:21:33 for Template Numerical Library by doxygen 1.9.8