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

Table of Contents

Introduction

Grids are regular orthogonal meshes. Similar to unstructured numerical meshes they provide indexing of mesh entities 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 similar 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 standard basis which are orthogonal to the grid entity. The following tables show all possible grid entities in 1D, 2D and 3D.

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, the basis vector is always computed on the fly. So whenever it possible, using the normals vector is preferred for better performance.

Remark: The entity orientation given by the normals or basis vector should be encoded statically in the type of the entity. This would make the implementation of the grid entities more efficient. Such implementation, however, requires support 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 following template parameters:

  • Dimension is dimension of the grid. This can be any integer 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
5main( int argc, char* argv[] )
6{
7 TNL::Meshes::Grid< 1 > grid_1D( 10 );
8 grid_1D.setDomain( { 0.0 }, { 1.0 } );
9 std::cout << grid_1D << std::endl;
10
11 TNL::Meshes::Grid< 2 > grid_2D( 10, 20 );
12 grid_2D.setDomain( { 0.0, 0.0 }, { 1.0, 2.0 } );
13 std::cout << grid_2D << std::endl;
14
15 TNL::Meshes::Grid< 3 > grid_3D( 10, 20, 30 );
16 grid_3D.setDomain( { 0.0, 0.0, 0.0 }, { 1.0, 2.0, 3.0 } );
17 std::cout << grid_3D << std::endl;
18
19 return EXIT_SUCCESS;
20}
TNL::Meshes::Grid
Orthogonal n-dimensional grid.
Definition Grid.h:30
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 is the same along each axis. The code looks as follows:

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

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

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

The whole example consists of the following steps:

  1. Setting dimension and resolution of the grid, defining necessary types that will be used later (grid type GridType, type for coordinates of the grid entities CoordinatesType, type for the real-world coordinates PointType, and type of container for storing the values of particular grid entities VectorType), and definition of types for grid entities (GridCell, GridFace, and GridVertex):

    // Define grid dimension and size.
    static constexpr int Dimension = 2;
    const int grid_size = 5;
    // Setup necessary types.
    using GridType = TNL::Meshes::Grid< Dimension, double, Device >;
    using CoordinatesType = typename GridType::CoordinatesType;
    using PointType = typename GridType::PointType;
    using VectorType = TNL::Containers::Vector< double, Device >;
    // Setup types of grid entities.
    using GridCell = typename GridType::Cell;
    using GridFace = typename GridType::Face;
    using GridVertex = typename GridType::Vertex;

  2. We create an instance of the grid:

    // Create an instance of a grid.
    GridType grid( grid_size );
    PointType origin( 0.0 );
    PointType proportions( 1.0 );
    grid.setDomain( origin, proportions );

    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 which are passed to the method TNL::Meshes::Grid::setDomain.

  3. Next we allocate vectors for storing of values in particular grid entities:

    // Allocate vectors for values stored in particular grid entities.
    VectorType cells( grid.template getEntitiesCount< Dimension >(), 0.0 );
    VectorType faces( grid.template getEntitiesCount< Dimension - 1 >(), 0.0 );
    VectorType vertexes( grid.template getEntitiesCount< 0 >(), 0.0 );

    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.

  4. In the next step, we prepare vector views (cells_view, faces_view and vertexes_view) which we will need later in lambda functions:

    // Prepare views for the data at the grid entities so that we can
    // manipulate them in lambda functions runnig eventually on GPU.
    auto cells_view = cells.getView();
    auto faces_view = faces.getView();
    auto vertexes_view = vertexes.getView();

  5. We iterate over all grid cells using templated method TNL::Meshes::Grid::forAllEntities. The template parameter again says the 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 lambda 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.

    // Setup value of each cell to its index in the grid.
    grid.template forAllEntities< Dimension >(
    [ = ] __cuda_callable__( const GridCell& cell ) mutable
    {
    cells_view[ cell.getIndex() ] = cell.getIndex();
    } );

  6. Next we print the values of all cells. This must be done sequentially on the CPU and therefore we do not use parallel for (TNL::Meshes::Grid::forAllEntities) anymore. Now we iterate over all cells sequentially row by row:

    // Print values of all cells in the grid.
    std::cout << "Values of cells .... " << std::endl;
    for( int i = grid_size - 1; i >= 0; i-- ) {
    for( int j = 0; j < grid_size; j++ ) {
    GridCell cell( grid, { j, i } );
    auto idx = cell.getIndex();
    std::cout << std::right << std::setw( 12 ) << cells.getElement( idx );
    }
    std::cout << std::endl;
    }
    std::cout << std::endl;

    In the inner loop we create a grid cell with given coordinates and we ask for its global index which we use for access to the vector with cell values (cell). The result of this is a matrix of cell values printed to the standard output.

  7. We iterate over all faces of the grid. Again we use the method TNL::Meshes::Grid::forAllEntities but now the template parameter specifying the dimension of the entities is Dimension-1. The lambda function being performed for each face now gets one parameter which is a grid entity (TNL::Meshes::GridEntity) with dimension equal to one:

    // Setup values of all faces to an average value of its neighbour cells.
    grid.template forAllEntities< Dimension - 1 >(
    [ = ] __cuda_callable__( const GridFace& face ) mutable
    {
    const CoordinatesType& normal = face.getNormals();
    double sum = 0.0;
    double count = 0.0;
    if( TNL::all( greaterEqual( face.getCoordinates() - normal, 0 ) ) ) {
    auto neighbour = face.template getNeighbourEntity< Dimension >( -normal );
    sum += cells_view[ neighbour.getIndex() ];
    count++;
    }
    if( TNL::all( less( face.getCoordinates(), face.getGrid().getDimensions() ) ) ) {
    auto neighbour = face.template getNeighbourEntity< Dimension >( { 0, 0 } );
    sum += cells_view[ neighbour.getIndex() ];
    count++;
    }
    faces_view[ face.getIndex() ] = sum / count;
    } );

    First of all we fetch the normal vector of the face by calling the method TNL::Meshes::GridEntity::getNormals which in case of faces is really 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 neighbor 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 (first if statement in the previous snippet) and if it is the case, we fetch its global index, add its value to the sum and then we increment the variable count which is the number of values that contributed to the average. Next we have to get the value from the cell that has the same coordinates as the face itself. We first check if there is such a cell (second if statement in the previous snippet) - 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 neighbor cell, we fetch its index, add its value to the sum and increment the variable count. Finally, we compute the average value of the neighbor cells and store in the vector with values of faces at the position given by the global index of the current face.

  8. Next we print the values of the faces on the console. Line by line, we print first horizontal faces with row index equal to five, followed by vertical faces with row index equal to five, next horizontal faces with row index equal to 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):

    // Print values of all faces in the grid.
    std::cout << "Values of faces ..." << std::endl;
    for( int i = grid_size; i >= 0; i-- ) {
    std::cout << std::right << std::setw( 6 ) << " ";
    for( int j = 0; j < grid_size; j++ ) {
    GridFace face( grid, { j, i }, { 0, 1 } );
    auto idx = face.getIndex();
    std::cout << std::right << std::setw( 12 ) << faces.getElement( idx );
    }
    std::cout << std::endl;
    if( i > 0 )
    for( int j = 0; j <= grid_size; j++ ) {
    GridFace face( grid, { j, i - 1 }, { 1, 0 } );
    auto idx = face.getIndex();
    std::cout << std::right << std::setw( 12 ) << faces.getElement( idx );
    }
    std::cout << std::endl;
    }

    So we iterate over all rows of the grid and for each row we first print the horizontal faces and then the vertical faces. Note that in both cases we create an instance of a grid face where the last parameter of the constructor determines the orientation of the face - {0, 1} is the normal of horizontal faces and {1, 0} is the normal of vertical faces.

  9. Finally, we iterate over all vertexes and compute average value of all neighboring cells. The vertexes have no orientation so we do not need to care about the normals. We only check the number of neighbor 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:

    // Setup values of all vertexes to an average value of its neighboring cells
    grid.template forAllEntities< 0 >(
    [ = ] __cuda_callable__( const GridVertex& vertex ) mutable
    {
    double sum = 0.0;
    double count = 0.0;
    auto grid_dimensions = vertex.getGrid().getDimensions();
    if( vertex.getCoordinates().x() > 0 && vertex.getCoordinates().y() > 0 ) {
    auto neighbour = vertex.template getNeighbourEntity< Dimension >( { -1, -1 } );
    sum += cells_view[ neighbour.getIndex() ];
    count++;
    }
    if( vertex.getCoordinates().x() > 0 && vertex.getCoordinates().y() < grid_dimensions.y() ) {
    auto neighbour = vertex.template getNeighbourEntity< Dimension >( { -1, 0 } );
    sum += cells_view[ neighbour.getIndex() ];
    count++;
    }
    if( vertex.getCoordinates().x() < grid_dimensions.x() && vertex.getCoordinates().y() > 0 ) {
    auto neighbour = vertex.template getNeighbourEntity< Dimension >( { 0, -1 } );
    sum += cells_view[ neighbour.getIndex() ];
    count++;
    }
    if( TNL::all( less( vertex.getCoordinates(), vertex.getGrid().getDimensions() ) ) ) {
    auto neighbour = vertex.template getNeighbourEntity< Dimension >( { 0, 0 } );
    sum += cells_view[ neighbour.getIndex() ];
    count++;
    }
    vertexes_view[ vertex.getIndex() ] = sum / count;
    } );

    The lambda function we perform on each vertex, has the parameter vertex which represents the vertex that we currently operate on. We just fetch the coordinates of the vertex and then we check what are the neighbor cells (four if statements in the previous snippet). For each such a cell we add its value to the sum. Finally, we compute the average value and store it in the vector.

  10. At the end we print values of all vertexes the same way as we did it with cells:

    // Print values of all vertexes in the grid.
    std::cout << "Values of vertexes .... " << std::endl;
    for( int i = grid_size; i >= 0; i-- ) {
    for( int j = 0; j <= grid_size; j++ ) {
    GridVertex vertex( grid, { j, i } );
    auto idx = vertex.getIndex();
    std::cout << std::right << std::setw( 12 ) << vertexes.getElement( idx );
    }
    std::cout << std::endl;
    }

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

Generated on Mon Apr 29 2024 10:28:40 for Template Numerical Library by doxygen 1.10.0