mesh.cpp 22.8 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
    // ComputePartitioning();
21
    ComputeBounds();
22
}
23

24
Mesh::~Mesh() {}
25

26
void Mesh::ComputeConnectivity() {
27

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

    // determine number/chunk size and indices of cells treated by each mpi thread
34
35
36
    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 );
37
38
39

    // '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
40
41
    std::vector<int> neighborsFlatPart( _numNodesPerCell * chunkSize, -1 );
    std::vector<Vector> normalsFlatPart( _numNodesPerCell * chunkSize, Vector( _dim, -1.0 ) );
42

43
    // pre sort cells and boundaries; sorting is needed for std::set_intersection
44
45
46
47
48
49
50
51
52
    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() );
    }
53

54
    blaze::CompressedMatrix<bool> connMat( _numCells, _numNodes );
55
56
57
58
59
60
    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
61
62
63
    for( unsigned i = mpiCellStart; i < mpiCellEnd; ++i ) {
        std::vector<unsigned>* cellsI = &sortedCells[i];
        for( unsigned j = 0; j < _numCells; ++j ) {
64
            if( i == j ) continue;
65
            if( static_cast<unsigned>( blaze::dot( blaze::row( connMat, i ), blaze::row( connMat, j ) ) ) ==
66
                _numNodesPerBoundary ) {    // in 2D cells are neighbors if they share two nodes std::vector<unsigned>* cellsJ = &sortedCells[j];
67
68
69
70
71
72
73
                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
74
                // determine unused index
75
76
                unsigned pos0 = _numNodesPerCell * ( i - mpiCellStart );
                unsigned pos  = pos0;
77
                while( neighborsFlatPart[pos] != -1 && pos < pos0 + _numNodesPerCell - 1 && pos < chunkSize * _numNodesPerCell - 1 ) pos++;
78
                neighborsFlatPart[pos] = j;
79
80
                // compute normal vector
                normalsFlatPart[pos] = ComputeOutwardFacingNormal( _nodes[commonElements[0]], _nodes[commonElements[1]], _cellMidPoints[i] );
81
82
            }
        }
83
        // boundaries are treated similarly to normal cells, but need a special treatment due to the absence of a neighboring cell
84
85
86
87
88
89
90
91
        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;
92
                    while( neighborsFlatPart[pos] != -1 && pos < pos0 + _numNodesPerCell - 1 && pos < chunkSize * _numNodesPerCell - 1 ) pos++;
93
                    neighborsFlatPart[pos] = _ghostCellID;
94
                    normalsFlatPart[pos]   = ComputeOutwardFacingNormal( _nodes[commonElements[0]], _nodes[commonElements[1]], _cellMidPoints[i] );
95
96
97
98
                }
            }
        }
    }
99
100

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

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

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

146
    // assign boundary types to all cells
147
    _cellBoundaryTypes.resize( _numCells, BOUNDARY_TYPE::NONE );
Jannick Wolters's avatar
Jannick Wolters committed
148
    for( unsigned i = 0; i < _numCells; ++i ) {
149
150
151
152
        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
153
154
155
156
157
158
159
                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
160
    }
161
162
163
164
165
166
}

void Mesh::ComputeCellAreas() {
    _cellAreas.resize( _numCells );
    for( unsigned i = 0; i < _numCells; ++i ) {
        switch( _numNodesPerCell ) {
167
            case 3: {    // triangular cells
168
169
170
171
172
173
                _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;
            }
174
            case 4: {    // quadrilateral cells
175
176
177
178
                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] };
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

                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: {
194
195
                ErrorMessages::Error( "Area computation for cells with " + std::to_string( _numNodesPerCell ) + " nodes is not implemented yet!",
                                      CURRENT_FUNCTION );
196
197
198
199
200
            }
        }
    }
}

Jonas Kusch's avatar
Jonas Kusch committed
201
202
203
204
205
206
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]];
        }
207
        _cellMidPoints[j] = _cellMidPoints[j] / static_cast<double>( _cells[j].size() );
Jonas Kusch's avatar
Jonas Kusch committed
208
209
210
    }
}

211
212
213
Vector Mesh::ComputeOutwardFacingNormal( const Vector& nodeA, const Vector& nodeB, const Vector& cellCenter ) {
    double dx = nodeA[0] - nodeB[0];
    double dy = nodeA[1] - nodeB[1];
214
215
216
217
    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
218
219
220
221
222
223
224
225
226
    }
    return n;
}

void Mesh::ComputePartitioning() {
    int comm_size, comm_rank;
    MPI_Comm comm = MPI_COMM_WORLD;
    MPI_Comm_size( comm, &comm_size );
    MPI_Comm_rank( comm, &comm_rank );
227
    unsigned ompNumThreads = omp_get_max_threads();
228

229
    // only run coloring if multiple OpenMP threads are present
230
    if( ompNumThreads > 1 ) {
231
        // setup adjacency relationship of nodes
Jannick Wolters's avatar
Jannick Wolters committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
        blaze::CompressedMatrix<bool> adjMatrix( _numNodes, _numNodes );
        for( unsigned i = 0; i < _numNodes; ++i ) {
            for( unsigned j = 0; j < _numCells; ++j ) {
                for( unsigned k = 0; k < _numNodesPerCell; ++k ) {
                    if( i == _cells[j][k] ) {
                        if( k == 0 ) {
                            adjMatrix.set( i, _cells[j][_numNodesPerCell - 1], true );
                            adjMatrix.set( i, _cells[j][1], true );
                        }
                        else if( k == _numNodesPerCell - 1 ) {
                            adjMatrix.set( i, _cells[j][_numNodesPerCell - 2], true );
                            adjMatrix.set( i, _cells[j][0], true );
                        }
                        else {
                            adjMatrix.set( i, _cells[j][k - 1], true );
                            adjMatrix.set( i, _cells[j][k + 1], true );
                        }
                    }
                }
            }
        }
253

254
        // for xadj and adjncy structure see parmetis documentation
Jannick Wolters's avatar
Jannick Wolters committed
255
        std::vector<int> xadj( _numNodes + 1 );
256
        int ctr = 0;
Jannick Wolters's avatar
Jannick Wolters committed
257
        std::vector<int> adjncy;
258
259
260
        for( unsigned i = 0; i < _numNodes; ++i ) {
            xadj[i] = ctr;
            for( unsigned j = 0; j < _numNodes; ++j ) {
Jannick Wolters's avatar
Jannick Wolters committed
261
262
                if( adjMatrix( i, j ) ) {
                    adjncy.push_back( static_cast<int>( j ) );
263
264
265
266
                    ctr++;
                }
            }
        }
Jannick Wolters's avatar
Jannick Wolters committed
267
268
        xadj[_numNodes] = ctr;

Jannick Wolters's avatar
Jannick Wolters committed
269
270
271
        std::vector<int> partitions;
        int edgecut  = 0;
        int ncon     = 1;
272
        int nparts   = ompNumThreads;
273
        real_t ubvec = static_cast<real_t>( 1.05 );    // parameter taken from SU2
Jannick Wolters's avatar
Jannick Wolters committed
274

275
276
        if( comm_size > 1 ) {    // if multiple MPI threads -> use parmetis
            // split adjacency information for each MPI thread -> see 'local' variables
277
            std::vector<int> local_xadj{ 0 };
Jannick Wolters's avatar
Jannick Wolters committed
278
279
280
281
282
283
284
285
286
287
288
289
            unsigned xadjChunk = std::ceil( static_cast<float>( _numNodes ) / static_cast<float>( comm_size ) );
            unsigned xadjStart = comm_rank * xadjChunk + 1;
            unsigned adjncyStart;
            for( unsigned i = 0; i < xadjChunk; ++i ) {
                if( i == 0 ) adjncyStart = xadj[i + xadjStart - 1];
                local_xadj.push_back( xadj[i + xadjStart] );
            }
            unsigned adjncyEnd = local_xadj.back();
            for( unsigned i = 1; i < local_xadj.size(); ++i ) {
                local_xadj[i] -= xadj[xadjStart - 1];
            }
            std::vector<int> local_adjncy( adjncy.begin() + adjncyStart, adjncy.begin() + adjncyEnd );
290

291
            // vtxdist describes what nodes are on which processor - see parmetis doc
292
            std::vector<int> vtxdist{ 0 };
Jannick Wolters's avatar
Jannick Wolters committed
293
294
            for( unsigned i = 1; i <= static_cast<unsigned>( comm_size ); i++ ) {
                vtxdist.push_back( std::min( i * xadjChunk, _numNodes ) );
295
296
            }

297
            // weights setup set as in SU2
Jannick Wolters's avatar
Jannick Wolters committed
298
299
300
301
            real_t* tpwgts = new real_t[nparts];
            for( unsigned i = 0; i < static_cast<unsigned>( nparts ); ++i ) {
                tpwgts[i] = 1.0 / static_cast<real_t>( nparts );
            }
302

303
            // setup as in SU2
Jannick Wolters's avatar
Jannick Wolters committed
304
305
            int wgtflag = 0;
            int numflag = 0;
Jannick Wolters's avatar
Jannick Wolters committed
306
307
            int options[1];
            options[0] = 0;
308

Jannick Wolters's avatar
Jannick Wolters committed
309
310
311
312
313
            unsigned chunkSize = local_xadj.size() - 1;
            int* part          = new int[chunkSize];
            ParMETIS_V3_PartKway( vtxdist.data(),
                                  local_xadj.data(),
                                  local_adjncy.data(),
Jannick Wolters's avatar
Jannick Wolters committed
314
315
316
317
318
319
                                  nullptr,
                                  nullptr,
                                  &wgtflag,
                                  &numflag,
                                  &ncon,
                                  &nparts,
Jannick Wolters's avatar
Jannick Wolters committed
320
                                  tpwgts,
Jannick Wolters's avatar
Jannick Wolters committed
321
322
323
324
325
                                  &ubvec,
                                  options,
                                  &edgecut,
                                  part,
                                  &comm );
Jannick Wolters's avatar
Jannick Wolters committed
326
            partitions.resize( chunkSize * comm_size );
327
            MPI_Allgather( part, chunkSize, MPI_INT, partitions.data(), chunkSize, MPI_INT, comm );    // gather all information in partitions
328
            partitions.resize( _numNodes );
329

330
            // free memory
Jannick Wolters's avatar
Jannick Wolters committed
331
            delete[] tpwgts;
332
        }
333
        else {    // with a singular MPI thread -> use metis
Jannick Wolters's avatar
Jannick Wolters committed
334
335
            int options[METIS_NOPTIONS];
            METIS_SetDefaultOptions( options );
Jannick Wolters's avatar
Jannick Wolters committed
336
337
            int nvtxs = _numNodes;
            int vsize;
Jannick Wolters's avatar
Jannick Wolters committed
338
            partitions.resize( _numNodes );
Jannick Wolters's avatar
Jannick Wolters committed
339
            METIS_PartGraphKway(
Jannick Wolters's avatar
Jannick Wolters committed
340
                &nvtxs, &ncon, xadj.data(), adjncy.data(), nullptr, &vsize, nullptr, &nparts, nullptr, &ubvec, options, &edgecut, partitions.data() );
Jannick Wolters's avatar
Jannick Wolters committed
341
        }
342
343

        // extract coloring information
Jannick Wolters's avatar
Jannick Wolters committed
344
345
346
347
        _colors.resize( _numCells );
        for( unsigned i = 0; i < _numCells; ++i ) {
            std::map<unsigned, int> occurances;
            for( unsigned j = 0; j < _numNodesPerCell; ++j ) {
Jannick Wolters's avatar
Jannick Wolters committed
348
                ++occurances[partitions[_cells[i][j]]];
Jannick Wolters's avatar
Jannick Wolters committed
349
350
            }
            _colors[i] = occurances.rbegin()->first;
351
352
353
        }
    }
    else {
Jannick Wolters's avatar
Jannick Wolters committed
354
        _colors.resize( _numCells, 0u );
355
356
357
    }
}

358
359
360
361
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 ) {
            // if( cell->IsBoundaryCell() ) continue; // skip ghost cells
vavrines's avatar
vavrines committed
362
            if( _cellBoundaryTypes[j] != 2 ) continue;    // skip ghost cells
363
364
365
366
367
368
369
370
371
            // 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
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
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
388
389
            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
390

vavrines's avatar
vavrines committed
391
392
            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
393
394
395
396
397
398
399
400

            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
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417

    double phi;
    VectorVector dPsiMax = std::vector( _numCells, Vector( nq, 0.0 ) );
    VectorVector dPsiMin = std::vector( _numCells, Vector( nq, 0.0 ) );
    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;

vavrines's avatar
vavrines committed
418
            // step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
vavrines's avatar
vavrines committed
419
420
421
422
423
424
425
426
427
428
            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
429
430
431
                // psiSample[j][k][l] = 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ); // interface central points
                psiSample[j][k][l] = psi[j][k] + psiDerX[j][k] * ( _nodes[_cells[j][l]][0] - _cellMidPoints[j][0] ) +
                                     psiDerY[j][k] * ( _nodes[_cells[j][l]][1] - _cellMidPoints[j][1] );    // vertex points
vavrines's avatar
vavrines committed
432
433
434

                // step 3: calculate Phi_ij at sample points
                if( psiSample[j][k][l] > psi[j][k] ) {
vavrines's avatar
vavrines committed
435
                    phiSample[j][k][l] = fmin( 0.5, dPsiMax[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
vavrines's avatar
vavrines committed
436
437
                }
                else if( psiSample[j][l][k] < psi[j][k] ) {
vavrines's avatar
vavrines committed
438
                    phiSample[j][k][l] = fmin( 0.5, dPsiMin[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
vavrines's avatar
vavrines committed
439
440
                }
                else {
vavrines's avatar
vavrines committed
441
                    phiSample[j][k][l] = 0.5;
vavrines's avatar
vavrines committed
442
443
444
                }
            }

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

vavrines's avatar
vavrines committed
448
            // step 5: limit the slope reconstructed from Gauss theorem
vavrines's avatar
vavrines committed
449
450
451
452
453
454
455
456
            for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
                psiDerX[j][k] *= phi;
                psiDerY[j][k] *= phi;
            }
        }
    }
}

457
458
459
460
461
462
463
464
465
466
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];
        }
    }
}

467
const std::vector<Vector>& Mesh::GetNodes() const { return _nodes; }
Jonas Kusch's avatar
Jonas Kusch committed
468
const std::vector<Vector>& Mesh::GetCellMidPoints() const { return _cellMidPoints; }
469
const std::vector<std::vector<unsigned>>& Mesh::GetCells() const { return _cells; }
Jannick Wolters's avatar
Jannick Wolters committed
470
471
const std::vector<double>& Mesh::GetCellAreas() const { return _cellAreas; }
const std::vector<unsigned>& Mesh::GetPartitionIDs() const { return _colors; }
Jonas Kusch's avatar
Jonas Kusch committed
472
473
const std::vector<std::vector<unsigned>>& Mesh::GetNeighbours() const { return _cellNeighbors; }
const std::vector<std::vector<Vector>>& Mesh::GetNormals() const { return _cellNormals; }
474
const std::vector<BOUNDARY_TYPE>& Mesh::GetBoundaryTypes() const { return _cellBoundaryTypes; }
475
const std::vector<std::pair<double, double>> Mesh::GetBounds() const { return _bounds; }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
476
477
478
479
480
481
482
483

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 );
}