mesh.cpp 19.6 KB
Newer Older
1
#include "common/mesh.h"
2

3
#include <chrono>
4

5
6
7
Mesh::Mesh( std::vector<Vector> nodes,
            std::vector<std::vector<unsigned>> cells,
            std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries )
8
9
    : _dim( nodes[0].size() ), _numCells( cells.size() ), _numNodes( nodes.size() ), _numNodesPerCell( cells[0].size() ),
      _numBoundaries( boundaries.size() ), _ghostCellID( _numCells ), _nodes( nodes ), _cells( cells ), _boundaries( boundaries ) {
10
11
12
13
14
15
16
    if( _dim == 2 ) {
        _numNodesPerBoundary = 2u;
    }
    else {
        ErrorMessages::Error( "Unsupported mesh dimension!", CURRENT_FUNCTION );
    }

17
    ComputeCellAreas();
Jonas Kusch's avatar
Jonas Kusch committed
18
    ComputeCellMidpoints();
19
    ComputeConnectivity();
20
    ComputeBounds();
21
}
22

23
Mesh::~Mesh() {}
24

25
void Mesh::ComputeConnectivity() {
26

27
    // get MPI info
28
29
30
    int comm_size, comm_rank;
    MPI_Comm_size( MPI_COMM_WORLD, &comm_size );
    MPI_Comm_rank( MPI_COMM_WORLD, &comm_rank );
31
32

    // determine number/chunk size and indices of cells treated by each mpi thread
33
34
35
    unsigned chunkSize    = std::ceil( static_cast<float>( _numCells ) / static_cast<float>( comm_size ) );
    unsigned mpiCellStart = comm_rank * chunkSize;
    unsigned mpiCellEnd   = std::min( ( comm_rank + 1 ) * chunkSize, _numCells );
36
37
38

    // 'flat' vectors are a flattened representation of the neighbors numCells<numNodesPerCell> nested vectors; easier for MPI
    // 'part' vectors store information for each single MPI thread
39
40
    std::vector<int> neighborsFlatPart( _numNodesPerCell * chunkSize, -1 );
    std::vector<Vector> normalsFlatPart( _numNodesPerCell * chunkSize, Vector( _dim, -1.0 ) );
41

42
    // pre sort cells and boundaries; sorting is needed for std::set_intersection
43
44
45
46
47
48
49
50
51
    auto sortedCells( _cells );
    for( unsigned i = 0; i < _numCells; ++i ) {
        std::sort( sortedCells[i].begin(), sortedCells[i].end() );
    }
    std::vector<std::vector<unsigned>> sortedBoundaries;
    for( unsigned i = 0; i < _numBoundaries; ++i ) {
        sortedBoundaries.push_back( _boundaries[i].second );
        std::sort( sortedBoundaries[i].begin(), sortedBoundaries[i].end() );
    }
52

53
    blaze::CompressedMatrix<bool> connMat( _numCells, _numNodes );
54
55
56
57
58
59
    for( unsigned i = mpiCellStart; i < mpiCellEnd; ++i ) {
        for( auto j : _cells[i] ) connMat.set( i, j, true );
    }

// determine neighbor cells and normals with MPI and OpenMP
#pragma omp parallel for
60
61
62
    for( unsigned i = mpiCellStart; i < mpiCellEnd; ++i ) {
        std::vector<unsigned>* cellsI = &sortedCells[i];
        for( unsigned j = 0; j < _numCells; ++j ) {
63
            if( i == j ) continue;
64
            if( static_cast<unsigned>( blaze::dot( blaze::row( connMat, i ), blaze::row( connMat, j ) ) ) ==
65
                _numNodesPerBoundary ) {    // in 2D cells are neighbors if they share two nodes std::vector<unsigned>* cellsJ = &sortedCells[j];
66
67
68
69
70
71
72
                std::vector<unsigned>* cellsJ = &sortedCells[j];
                std::vector<unsigned> commonElements;
                std::set_intersection( cellsI->begin(),
                                       cellsI->end(),
                                       cellsJ->begin(),
                                       cellsJ->end(),
                                       std::back_inserter( commonElements ) );    // find common nodes of two cells
73
                // determine unused index
74
75
                unsigned pos0 = _numNodesPerCell * ( i - mpiCellStart );
                unsigned pos  = pos0;
76
                while( neighborsFlatPart[pos] != -1 && pos < pos0 + _numNodesPerCell - 1 && pos < chunkSize * _numNodesPerCell - 1 ) pos++;
77
                neighborsFlatPart[pos] = j;
78
79
                // compute normal vector
                normalsFlatPart[pos] = ComputeOutwardFacingNormal( _nodes[commonElements[0]], _nodes[commonElements[1]], _cellMidPoints[i] );
80
81
            }
        }
82
        // boundaries are treated similarly to normal cells, but need a special treatment due to the absence of a neighboring cell
83
84
85
86
87
88
89
90
        for( unsigned k = 0; k < _boundaries.size(); ++k ) {
            std::vector<unsigned>* bNodes = &sortedBoundaries[k];
            for( unsigned j = 0; j < _boundaries[k].second.size(); ++j ) {
                std::vector<unsigned> commonElements;
                std::set_intersection( cellsI->begin(), cellsI->end(), bNodes->begin(), bNodes->end(), std::back_inserter( commonElements ) );
                if( commonElements.size() == _dim ) {
                    unsigned pos0 = _numNodesPerCell * ( i - mpiCellStart );
                    unsigned pos  = pos0;
91
                    while( neighborsFlatPart[pos] != -1 && pos < pos0 + _numNodesPerCell - 1 && pos < chunkSize * _numNodesPerCell - 1 ) pos++;
92
                    neighborsFlatPart[pos] = _ghostCellID;
93
                    normalsFlatPart[pos]   = ComputeOutwardFacingNormal( _nodes[commonElements[0]], _nodes[commonElements[1]], _cellMidPoints[i] );
94
95
96
97
                }
            }
        }
    }
98
99

    // gather distributed data on all MPI threads
Jannick Wolters's avatar
Jannick Wolters committed
100
    std::vector<int> neighborsFlat( _numNodesPerCell * chunkSize * comm_size, -1 );
101
    std::vector<Vector> normalsFlat( _numNodesPerCell * chunkSize * comm_size, Vector( _dim, 0.0 ) );
102
    if( comm_size == 1 ) {    // can be done directly if there is only one MPI thread
103
        neighborsFlat.assign( neighborsFlatPart.begin(), neighborsFlatPart.end() );
104
105
106
        normalsFlat.assign( normalsFlatPart.begin(), normalsFlatPart.end() );
    }
    else {
Jannick Wolters's avatar
Jannick Wolters committed
107
108
109
110
111
112
113
        MPI_Allgather( neighborsFlatPart.data(),
                       _numNodesPerCell * chunkSize,
                       MPI_INT,
                       neighborsFlat.data(),
                       _numNodesPerCell * chunkSize,
                       MPI_INT,
                       MPI_COMM_WORLD );
114
    }
Jannick Wolters's avatar
Jannick Wolters committed
115

116
    // check for any unassigned faces
117
    if( std::any_of( neighborsFlat.begin(), neighborsFlat.end(), []( int i ) { return i == -1; } ) ) {
118
        for( unsigned idx = 0; idx < neighborsFlat.size(); ++idx ) {
119
120
            if( neighborsFlat[idx] == -1 )
                ErrorMessages::Error( "Detected unassigned faces at index " + std::to_string( idx ) + " !", CURRENT_FUNCTION );
121
122
        }
    }
123

124
    // reorder neighbors and normals into nested structure
125
126
127
    _cellNeighbors.resize( _numCells );
    _cellNormals.resize( _numCells );
    for( unsigned i = 0; i < neighborsFlat.size(); ++i ) {
128
        unsigned IDi = static_cast<unsigned>( i / static_cast<double>( _numNodesPerCell ) );
129
        unsigned IDj = neighborsFlat[i];
130
131
        if( IDi == IDj ) continue;     // avoid self assignment
        if( IDj == _ghostCellID ) {    // cell is boundary cell
132
133
134
135
136
            if( std::find( _cellNeighbors[IDi].begin(), _cellNeighbors[IDi].end(), _ghostCellID ) == _cellNeighbors[IDi].end() ) {
                _cellNeighbors[IDi].push_back( _ghostCellID );
                _cellNormals[IDi].push_back( normalsFlat[i] );
            }
        }
137
        else {    // normal cell neighbor
138
139
140
141
142
143
144
            if( std::find( _cellNeighbors[IDi].begin(), _cellNeighbors[IDi].end(), IDj ) == _cellNeighbors[IDi].end() ) {
                _cellNeighbors[IDi].push_back( IDj );
                _cellNormals[IDi].push_back( normalsFlat[i] );
            }
        }
    }

145
    // assign boundary types to all cells
146
    _cellBoundaryTypes.resize( _numCells, BOUNDARY_TYPE::NONE );
Jannick Wolters's avatar
Jannick Wolters committed
147
    for( unsigned i = 0; i < _numCells; ++i ) {
148
149
150
151
        if( std::any_of( _cellNeighbors[i].begin(), _cellNeighbors[i].end(), [this]( unsigned i ) {
                return i == _ghostCellID;
            } ) ) {                           // cell is boundary cell
            for( auto bc : _boundaries ) {    // loop over all boundaries and boundary nodes and search for the respective nodeID
152
153
154
155
156
157
158
                for( auto cNodes : _cells[i] ) {
                    if( std::find( bc.second.begin(), bc.second.end(), cNodes ) != bc.second.end() ) {
                        _cellBoundaryTypes[i] = bc.first;
                    }
                }
            }
        }
Jannick Wolters's avatar
Jannick Wolters committed
159
    }
160
161
162
163
164
165
}

void Mesh::ComputeCellAreas() {
    _cellAreas.resize( _numCells );
    for( unsigned i = 0; i < _numCells; ++i ) {
        switch( _numNodesPerCell ) {
166
            case 3: {    // triangular cells
167
168
169
170
171
172
                _cellAreas[i] = std::abs( ( _nodes[_cells[i][0]][0] * ( _nodes[_cells[i][1]][1] - _nodes[_cells[i][2]][1] ) +
                                            _nodes[_cells[i][1]][0] * ( _nodes[_cells[i][2]][1] - _nodes[_cells[i][0]][1] ) +
                                            _nodes[_cells[i][2]][0] * ( _nodes[_cells[i][0]][1] - _nodes[_cells[i][1]][1] ) ) /
                                          2 );
                break;
            }
173
            case 4: {    // quadrilateral cells
174
175
176
177
                std::vector d1{ _nodes[_cells[i][0]][0] - _nodes[_cells[i][1]][0], _nodes[_cells[i][0]][1] - _nodes[_cells[i][1]][1] };
                std::vector d2{ _nodes[_cells[i][1]][0] - _nodes[_cells[i][2]][0], _nodes[_cells[i][1]][1] - _nodes[_cells[i][2]][1] };
                std::vector d3{ _nodes[_cells[i][2]][0] - _nodes[_cells[i][3]][0], _nodes[_cells[i][2]][1] - _nodes[_cells[i][3]][1] };
                std::vector d4{ _nodes[_cells[i][3]][0] - _nodes[_cells[i][0]][0], _nodes[_cells[i][3]][1] - _nodes[_cells[i][0]][1] };
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192

                double a = std::sqrt( d1[0] * d1[0] + d1[1] * d1[1] );
                double b = std::sqrt( d2[0] * d2[0] + d2[1] * d2[1] );
                double c = std::sqrt( d3[0] * d3[0] + d3[1] * d3[1] );
                double d = std::sqrt( d4[0] * d4[0] + d4[1] * d4[1] );
                double T = 0.5 * ( a + b + c + d );

                double alpha = std::acos( ( d4[0] * d1[0] + d4[1] * d1[1] ) / ( a * d ) );
                double beta  = std::acos( ( d2[0] * d3[0] + d2[1] * d3[1] ) / ( b * c ) );

                _cellAreas[i] = std::sqrt( ( T - a ) * ( T - b ) * ( T - c ) * ( T - d ) -
                                           a * b * c * d * std::cos( 0.5 * ( alpha + beta ) ) * std::cos( 0.5 * ( alpha + beta ) ) );
                break;
            }
            default: {
193
194
                ErrorMessages::Error( "Area computation for cells with " + std::to_string( _numNodesPerCell ) + " nodes is not implemented yet!",
                                      CURRENT_FUNCTION );
195
196
197
198
199
            }
        }
    }
}

Jonas Kusch's avatar
Jonas Kusch committed
200
201
202
203
204
205
void Mesh::ComputeCellMidpoints() {
    _cellMidPoints = std::vector( _numCells, Vector( _dim, 0.0 ) );
    for( unsigned j = 0; j < _numCells; ++j ) {
        for( unsigned l = 0; l < _cells[j].size(); ++l ) {
            _cellMidPoints[j] = _cellMidPoints[j] + _nodes[_cells[j][l]];
        }
206
        _cellMidPoints[j] = _cellMidPoints[j] / static_cast<double>( _cells[j].size() );
Jonas Kusch's avatar
Jonas Kusch committed
207
208
209
    }
}

210
211
212
Vector Mesh::ComputeOutwardFacingNormal( const Vector& nodeA, const Vector& nodeB, const Vector& cellCenter ) {
    double dx = nodeA[0] - nodeB[0];
    double dy = nodeA[1] - nodeB[1];
213
214
215
216
    Vector n{ -dy, dx };                    // normal vector
    Vector p{ nodeA[0], nodeA[1] };         // take arbitrary node
    if( dot( n, cellCenter - p ) > 0 ) {    // dot product is negative for outward facing normals
        n *= -1.0;                          // if dot product is positive -> flip normal
217
218
219
220
    }
    return n;
}

221
222
223
void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
    for( unsigned k = 0; k < nq; ++k ) {
        for( unsigned j = 0; j < _numCells; ++j ) {
224
225
226
            psiDerX[j][k] = 0.0;
            psiDerY[j][k] = 0.0;

227
            // if( cell->IsBoundaryCell() ) continue; // skip ghost cells
vavrines's avatar
vavrines committed
228
            if( _cellBoundaryTypes[j] != 2 ) continue;    // skip ghost cells
229
230
231
232
233
234
235
236
237
            // compute derivative by summing over cell boundary
            for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
                psiDerX[j][k] = psiDerX[j][k] + 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][0] / _cellAreas[j];
                psiDerY[j][k] = psiDerY[j][k] + 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][1] / _cellAreas[j];
            }
        }
    }
}

vavrines's avatar
vavrines committed
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
void Mesh::ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {

    double duxL;
    double duyL;
    double duxR;
    double duyR;

    for( unsigned k = 0; k < nq; ++k ) {
        for( unsigned j = 0; j < _numCells; ++j ) {
            // reset derivatives
            psiDerX[j][k] = 0.0;
            psiDerY[j][k] = 0.0;

            // skip boundary cells
            if( _cellBoundaryTypes[j] != 2 ) continue;

vavrines's avatar
vavrines committed
254
255
            duxL = ( psi[j][k] - psi[_cellNeighbors[j][0]][k] ) / ( _cellMidPoints[j][0] - _cellMidPoints[_cellNeighbors[j][0]][0] + 1e-8 );
            duyL = ( psi[j][k] - psi[_cellNeighbors[j][0]][k] ) / ( _cellMidPoints[j][1] - _cellMidPoints[_cellNeighbors[j][0]][1] + 1e-8 );
vavrines's avatar
vavrines committed
256

vavrines's avatar
vavrines committed
257
258
            duxR = ( psi[j][k] - psi[_cellNeighbors[j][2]][k] ) / ( _cellMidPoints[j][0] - _cellMidPoints[_cellNeighbors[j][2]][0] + 1e-8 );
            duyR = ( psi[j][k] - psi[_cellNeighbors[j][2]][k] ) / ( _cellMidPoints[j][1] - _cellMidPoints[_cellNeighbors[j][2]][1] + 1e-8 );
vavrines's avatar
vavrines committed
259
260
261
262
263
264
265
266

            psiDerX[j][k] = LMinMod( duxL, duxR ) * 0.5;
            psiDerY[j][k] = LMinMod( duyL, duyR ) * 0.5;
        }
    }
}

void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
vavrines's avatar
vavrines committed
267
268

    double phi;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
269
270
    // VectorVector dPsiMax = std::vector( _numCells, Vector( nq, 0.0 ) );
    // VectorVector dPsiMin = std::vector( _numCells, Vector( nq, 0.0 ) );
271
272
273
    VectorVector dPsiMax( _numCells, Vector( nq, 0.0 ) );
    VectorVector dPsiMin( _numCells, Vector( nq, 0.0 ) );

vavrines's avatar
vavrines committed
274
275
276
277
278
279
280
281
282
283
284
285
    std::vector<std::vector<Vector>> psiSample( _numCells, std::vector<Vector>( nq, Vector( _numNodesPerCell, 0.0 ) ) );
    std::vector<std::vector<Vector>> phiSample( _numCells, std::vector<Vector>( nq, Vector( _numNodesPerCell, 0.0 ) ) );

    for( unsigned k = 0; k < nq; ++k ) {
        for( unsigned j = 0; j < _numCells; ++j ) {
            // reset derivatives
            psiDerX[j][k] = 0.0;
            psiDerY[j][k] = 0.0;

            // skip boundary cells
            if( _cellBoundaryTypes[j] != 2 ) continue;

286
287
            /*
            // version 1: original taste
vavrines's avatar
vavrines committed
288
            // step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
vavrines's avatar
vavrines committed
289
290
291
292
293
294
295
296
297
298
            for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
                if( psi[_cellNeighbors[j][l]][k] - psi[j][k] > dPsiMax[j][k] ) dPsiMax[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
                if( psi[_cellNeighbors[j][l]][k] - psi[j][k] < dPsiMin[j][k] ) dPsiMin[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];

                psiDerX[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][0] / _cellAreas[j];
                psiDerY[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][1] / _cellAreas[j];
            }

            for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
                // step 2: choose sample points
vavrines's avatar
vavrines committed
299
                // psiSample[j][k][l] = 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ); // interface central points
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
300
                psiSample[j][k][l] = psi[j][k] +
301
                                     psiDerX[j][k] * ( _nodes[_cells[j][l]][0] - _cellMidPoints[j][0] ) +
vavrines's avatar
vavrines committed
302
                                     psiDerY[j][k] * ( _nodes[_cells[j][l]][1] - _cellMidPoints[j][1] );    // vertex points
vavrines's avatar
vavrines committed
303
304
305

                // step 3: calculate Phi_ij at sample points
                if( psiSample[j][k][l] > psi[j][k] ) {
306
                    phiSample[j][k][l] = fmin( 1.0, dPsiMax[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
vavrines's avatar
vavrines committed
307
308
                }
                else if( psiSample[j][l][k] < psi[j][k] ) {
309
                    phiSample[j][k][l] = fmin( 1.0, dPsiMin[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
vavrines's avatar
vavrines committed
310
311
                }
                else {
312
                    phiSample[j][k][l] = 1.0;
vavrines's avatar
vavrines committed
313
314
315
                }
            }

vavrines's avatar
vavrines committed
316
            // step 4: find minimum limiter function phi
vavrines's avatar
vavrines committed
317
318
            phi = min( phiSample[j][k] );

vavrines's avatar
vavrines committed
319
            // step 5: limit the slope reconstructed from Gauss theorem
320
321
322
            psiDerX[j][k] *= phi;
            psiDerY[j][k] *= phi;
            */
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
323

324
325
326
327
328
329
330
331
332
333
334
            double eps = 1e-6;
            // version 2: Venkatakrishnan limiter
            // step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
            for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
                if( psi[_cellNeighbors[j][l]][k] - psi[j][k] > dPsiMax[j][k] ) dPsiMax[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
                if( psi[_cellNeighbors[j][l]][k] - psi[j][k] < dPsiMin[j][k] ) dPsiMin[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];

                psiDerX[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][0] / _cellAreas[j];
                psiDerY[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][1] / _cellAreas[j];
            }

vavrines's avatar
vavrines committed
335
            for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
336
                // step 2: choose sample points
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
337
338
                psiSample[j][k][l] = 10.0 * psiDerX[j][k] * ( _cellMidPoints[_cellNeighbors[j][l]][0] - _cellMidPoints[j][0] ) +
                                     10.0 * psiDerY[j][k] * ( _cellMidPoints[_cellNeighbors[j][l]][1] - _cellMidPoints[j][1] );
339
340
341

                // step 3: calculate Phi_ij at sample points
                if( psiSample[j][k][l] > 0.0 ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
342
343
344
                    phiSample[j][k][l] =
                        ( dPsiMax[j][k] * dPsiMax[j][k] + 2.0 * dPsiMax[j][k] * psiSample[j][k][l] + eps ) /
                        ( dPsiMax[j][k] * dPsiMax[j][k] + dPsiMax[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + eps );
345
346
                }
                else if( psiSample[j][k][l] < 0.0 ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
347
348
349
350
                    phiSample[j][k][l] =
                        ( dPsiMin[j][k] * dPsiMin[j][k] + 2.0 * dPsiMin[j][k] * psiSample[j][k][l] + eps ) /
                        ( dPsiMin[j][k] * dPsiMin[j][k] + dPsiMin[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + eps );
                    ;
351
352
353
354
                }
                else {
                    phiSample[j][k][l] = 1.0;
                }
vavrines's avatar
vavrines committed
355
            }
356
357
358

            // step 4: find minimum limiter function phi
            phi = min( phiSample[j][k] );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
359
            // phi = fmin( 0.5, abs(phi) );
360
361
362
363

            // step 5: limit the slope reconstructed from Gauss theorem
            psiDerX[j][k] *= phi;
            psiDerY[j][k] *= phi;
vavrines's avatar
vavrines committed
364
365
366
367
        }
    }
}

368
369
370
371
372
373
374
375
376
377
void Mesh::ComputeBounds() {
    _bounds = std::vector( _dim, std::make_pair( std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity() ) );
    for( unsigned i = 0; i < _numNodes; ++i ) {
        for( unsigned j = 0; j < _dim; ++j ) {
            if( _nodes[i][j] < _bounds[j].first ) _bounds[j].first = _nodes[i][j];
            if( _nodes[i][j] > _bounds[j].second ) _bounds[j].second = _nodes[i][j];
        }
    }
}

378
const std::vector<Vector>& Mesh::GetNodes() const { return _nodes; }
Jonas Kusch's avatar
Jonas Kusch committed
379
const std::vector<Vector>& Mesh::GetCellMidPoints() const { return _cellMidPoints; }
380
const std::vector<std::vector<unsigned>>& Mesh::GetCells() const { return _cells; }
Jannick Wolters's avatar
Jannick Wolters committed
381
const std::vector<double>& Mesh::GetCellAreas() const { return _cellAreas; }
Jonas Kusch's avatar
Jonas Kusch committed
382
383
const std::vector<std::vector<unsigned>>& Mesh::GetNeighbours() const { return _cellNeighbors; }
const std::vector<std::vector<Vector>>& Mesh::GetNormals() const { return _cellNormals; }
384
const std::vector<BOUNDARY_TYPE>& Mesh::GetBoundaryTypes() const { return _cellBoundaryTypes; }
385
const std::vector<std::pair<double, double>> Mesh::GetBounds() const { return _bounds; }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
386
387
388
389
390
391
392
393

double Mesh::GetDistanceToOrigin( unsigned idx_cell ) const {
    double distance = 0.0;
    for( unsigned idx_dim = 0; idx_dim < _dim; idx_dim++ ) {
        distance += _cellMidPoints[idx_cell][idx_dim] * _cellMidPoints[idx_cell][idx_dim];
    }
    return sqrt( distance );
}