My Project
Loading...
Searching...
No Matches
GridInterfaceEuler.hpp
1//===========================================================================
2//
3// File: GridInterfaceEuler.hpp
4//
5// Created: Mon Jun 15 12:53:38 2009
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Bård Skaflestad <bard.skaflestad@sintef.no>
9//
10// $Date$
11//
12// $Revision$
13//
14//===========================================================================
15
16/*
17 Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18 Copyright 2009, 2010 Statoil ASA.
19
20 This file is part of The Open Reservoir Simulator Project (OpenRS).
21
22 OpenRS is free software: you can redistribute it and/or modify
23 it under the terms of the GNU General Public License as published by
24 the Free Software Foundation, either version 3 of the License, or
25 (at your option) any later version.
26
27 OpenRS is distributed in the hope that it will be useful,
28 but WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 GNU General Public License for more details.
31
32 You should have received a copy of the GNU General Public License
33 along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34*/
35
36#ifndef OPENRS_GRIDINTERFACEEULER_HEADER
37#define OPENRS_GRIDINTERFACEEULER_HEADER
38
39#include <opm/grid/utility/SparseTable.hpp>
40#include <opm/grid/utility/StopWatch.hpp>
41#include <opm/grid/CpGrid.hpp> // How to avoid this? Needed for the explicit mapper specialization below.
42
43#include <opm/common/utility/platform_dependent/disable_warnings.h>
44
45#include <dune/common/fvector.hh>
46#include <dune/grid/common/mcmgmapper.hh>
47#include <dune/grid/common/defaultgridview.hh>
48
49#include <boost/iterator/iterator_facade.hpp>
50
51#include <opm/common/utility/platform_dependent/reenable_warnings.h>
52
53#include <climits>
54#include <iostream>
55#include <memory>
56
57
58
59namespace Opm
60{
61
63 template <class DuneGrid>
65 {
66 using Type = Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<DuneGrid>>>;
67 };
68
71 {
72 explicit CpGridCellMapper(const Dune::CpGrid&)
73 {
74 }
75 template<class EntityType>
76 int map (const EntityType& e) const
77 {
78 return e.index();
79 }
80 };
81
82 namespace GIE
83 {
84 template <class GridInterface, class EntityPointerType>
85 class Cell;
86
87 template <class GI>
88 class Face {
89 public:
90 typedef typename GI::DuneIntersectionIterator DuneIntersectionIter;
91 typedef typename GI::GridType::Traits::template Codim<0>::Entity CellPtr;
92 typedef typename GI::GridType::ctype Scalar;
93 typedef Dune::FieldVector<Scalar, GI::GridType::dimension> Vector;
94 typedef Dune::FieldVector<Scalar, GI::GridType::dimension-1> LocalVector;
95 typedef int Index;
97
98 enum { BoundaryMarkerIndex = -999,
99 LocalEndIndex = INT_MAX };
100
101 Face()
102 : pgrid_(0), iter_(), local_index_(-1)
103 {
104 }
105
106 Face(const GI& grid,
107 const DuneIntersectionIter& it,
108 const Index loc_ind)
109 : pgrid_(&grid), iter_(it), local_index_(loc_ind)
110 {
111 }
112
116 Scalar area() const
117 {
118 return iter_->geometry().volume();
119 }
120
124 Vector centroid() const
125 {
126 //return iter_->geometry().global(localCentroid());
127 return iter_->geometry().center();
128 }
129
133 Vector normal() const
134 {
135 //return iter_->unitOuterNormal(localCentroid());
136 return iter_->centerUnitOuterNormal();
137 }
138
142 bool boundary() const
143 {
144 return iter_->boundary();
145 }
146
150 int boundaryId() const
151 {
152 return iter_->boundaryId();
153 }
154
158 Cell cell() const
159 {
160 return Cell(*pgrid_, iter_->inside());
161 }
162
166 Index cellIndex() const
167 {
168 return pgrid_->mapper().index(iter_->inside());
169 }
170
175 {
176 return Cell(*pgrid_, iter_->outside());
177 }
178
182 Index neighbourCellIndex() const
183 {
184 if (iter_->boundary()) {
185 return BoundaryMarkerIndex;
186 } else {
187 return pgrid_->mapper().index(
188 iter_->outside());
189 }
190 }
191
195 Index index() const
196 {
197 return pgrid_->faceIndex(cellIndex(), localIndex());
198 }
199
203 Index localIndex() const
204 {
205 return local_index_;
206 }
207
211 Scalar neighbourCellVolume() const
212 {
213 return iter_->outside()->geometry().volume();
214 }
215
216 protected:
217 const GI* pgrid_;
218 DuneIntersectionIter iter_;
219 Index local_index_;
220
221 private:
222// LocalVector localCentroid() const
223// {
224// typedef Dune::ReferenceElements<Scalar, GI::GridType::dimension-1> RefElems;
225// return RefElems::general(iter_->type()).position(0,0);
226// }
227 };
228
229
235 template <class GridInterface>
237 public boost::iterator_facade<FaceIterator<GridInterface>,
238 const Face<GridInterface>,
239 boost::forward_traversal_tag>,
240 public Face<GridInterface>
241 {
242 private:
244
245 public:
249 typedef typename Face<GridInterface>::DuneIntersectionIter DuneIntersectionIter;
250
253 : FaceType()
254 {
255 }
256
268 FaceIterator(const GridInterface& grid,
269 const DuneIntersectionIter& it,
270 const int local_index)
271 : FaceType(grid, it, local_index)
272 {
273 }
274
277 {
278 return *this;
279 }
280
282 bool equal(const FaceIterator& other) const
283 {
284 // Note that we do not compare the local_index_ members,
285 // since they may or may not be equal for end iterators.
286 return FaceType::iter_ == other.FaceType::iter_;
287 }
288
291 {
292 ++FaceType::iter_;
293 ++FaceType::local_index_;
294 }
295
297 bool operator<(const FaceIterator& other) const
298 {
299 if (FaceType::cellIndex() == other.FaceType::cellIndex()) {
300 return FaceType::localIndex() < other.FaceType::localIndex();
301 } else {
302 return FaceType::cellIndex() < other.FaceType::cellIndex();
303 }
304 }
305 };
306
307
308 template <class GridInterface, class EntityType>
309 class Cell
310 {
311 public:
312 Cell()
313 : pgrid_(0), cell_()
314 {
315 }
316 Cell(const GridInterface& grid,
317 const EntityType& cell)
318 : pgrid_(&grid), cell_(cell)
319 {
320 }
322 typedef typename FaceIterator::Vector Vector;
323 typedef typename FaceIterator::Scalar Scalar;
324 typedef typename FaceIterator::Index Index;
325
326 FaceIterator facebegin() const
327 {
328 return FaceIterator(*pgrid_, cell_.ileafbegin(), 0);
329 }
330
331 FaceIterator faceend() const
332 {
333 return FaceIterator(*pgrid_, cell_.ileafend(),
334 FaceIterator::LocalEndIndex);
335 }
336
337 Scalar volume() const
338 {
339 return cell_.geometry().volume();
340 }
341
342 Vector centroid() const
343 {
344// typedef Dune::ReferenceElements<Scalar, GridInterface::GridType::dimension> RefElems;
345// Vector localpt
346// = RefElems::general(cell_->type()).position(0,0);
347// return cell_->geometry().global(localpt);
348 return cell_.geometry().center();
349 }
350
351 Index index() const
352 {
353 return pgrid_->mapper().index(cell_);
354 }
355 protected:
356 const GridInterface* pgrid_;
357 EntityType cell_;
358 };
359
360
361 template <class GridInterface>
363 : public boost::iterator_facade<CellIterator<GridInterface>,
364 const Cell<GridInterface,
365 typename GridInterface::GridType::template Codim<0>::LeafIterator>,
366 boost::forward_traversal_tag>,
367 public Cell<GridInterface, typename GridInterface::GridType::template Codim<0>::LeafIterator>
368 {
369 private:
370 typedef typename GridInterface::GridType::template Codim<0>::LeafIterator DuneCellIter;
372 public:
373 typedef typename CellType::Vector Vector;
374 typedef typename CellType::Scalar Scalar;
375 typedef typename CellType::Index Index;
376
377 CellIterator(const GridInterface& grid, DuneCellIter it)
378 : CellType(grid, it)
379 {
380 }
383 {
384 return *this;
385 }
387 bool equal(const CellIterator& other) const
388 {
389 return CellType::cell_ == other.CellType::cell_;
390 }
393 {
394 ++CellType::cell_;
395 }
396 };
397
398 } // namespace GIE
399
400
401 template <class DuneGrid>
403 {
404 public:
405 typedef typename DuneGrid::LeafIntersectionIterator DuneIntersectionIterator;
406 typedef DuneGrid GridType;
409 typedef typename CellIterator::Vector Vector;
410 typedef typename CellIterator::Scalar Scalar;
411 typedef typename CellIterator::Index Index;
412
413 typedef typename GICellMapper<DuneGrid>::Type Mapper;
414
415 enum { Dimension = DuneGrid::dimension };
416
418 : pgrid_(0), num_faces_(0), max_faces_per_cell_(0)
419 {
420 }
421 explicit GridInterfaceEuler(const DuneGrid& grid, bool build_facemap = true)
422 : pgrid_(&grid), pmapper_(new Mapper(grid.leafGridView(), Dune::mcmgElementLayout())), num_faces_(0), max_faces_per_cell_(0)
423 {
424 if (build_facemap) {
425 buildFaceIndices();
426 }
427 }
428 void init(const DuneGrid& grid, bool build_facemap = true)
429 {
430 pgrid_ = &grid;
431 pmapper_.reset(new Mapper(grid.leafGridView(), Dune::mcmgElementLayout()));
432 if (build_facemap) {
433 buildFaceIndices();
434 }
435 }
436 CellIterator cellbegin() const
437 {
438 return CellIterator(*this, grid().template leafbegin<0>());
439 }
440 CellIterator cellend() const
441 {
442 return CellIterator(*this, grid().template leafend<0>());
443 }
444 int numberOfCells() const
445 {
446 return grid().size(0);
447 }
448 int numberOfFaces() const
449 {
450 assert(num_faces_ != 0);
451 return num_faces_;
452 }
453 int maxFacesPerCell() const
454 {
455 assert(max_faces_per_cell_ != 0);
456 return max_faces_per_cell_;
457 }
458 const DuneGrid& grid() const
459 {
460 assert(pgrid_);
461 return *pgrid_;
462 }
463
464 // The following are primarily helpers for the implementation,
465 // perhaps they should be private?
466 const Mapper& mapper() const
467 {
468 assert (pmapper_);
469 return *pmapper_;
470 }
471 Index faceIndex(int cell_index, int local_face_index) const
472 {
473 assert(num_faces_ != 0);
474 return face_indices_[cell_index][local_face_index];
475 }
476 typedef Opm::SparseTable<int>::row_type Indices;
477 Indices faceIndices(int cell_index) const
478 {
479 assert(num_faces_ != 0);
480 return face_indices_[cell_index];
481 }
482 private:
483 const DuneGrid* pgrid_;
484 std::unique_ptr<Mapper> pmapper_;
485 int num_faces_;
486 int max_faces_per_cell_;
487 Opm::SparseTable<int> face_indices_;
488
489 void buildFaceIndices()
490 {
491#ifdef VERBOSE
492 std::cout << "Building unique face indices... " << std::flush;
493 Opm::time::StopWatch clock;
494 clock.start();
495#endif
496 typedef CellIterator CI;
497 typedef typename CI::FaceIterator FI;
498
499 // We build the actual cell to face mapping in two passes.
500 // [code mostly lifted from IncompFlowSolverHybrid::enumerateGridDof(),
501 // but with a twist: This code builds a mapping from cells in index
502 // order to unique face numbers, while the mapping built in the
503 // enumerateGridDof() method was ordered by cell iterator order]
504
505 // Allocate and reserve structures.
506 const int nc = numberOfCells();
507 std::vector<int> cell(nc, -1);
508 std::vector<int> num_faces(nc); // In index order.
509 std::vector<int> fpos; fpos .reserve(nc + 1);
510 std::vector<int> num_cf; num_cf.reserve(nc); // In iterator order.
511 std::vector<int> faces ;
512
513 // First pass: enumerate internal faces.
514 int cellno = 0; fpos.push_back(0);
515 int tot_ncf = 0, max_ncf = 0;
516 for (CI c = cellbegin(); c != cellend(); ++c, ++cellno) {
517 const int c0 = c->index();
518 assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
519 cell[c0] = cellno;
520 num_cf.push_back(0);
521 int& ncf = num_cf.back();
522 for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
523 if (!f->boundary()) {
524 const int c1 = f->neighbourCellIndex();
525 assert((0 <= c1) && (c1 < nc) && (c1 != c0));
526
527 if (cell[c1] == -1) {
528 // Previously undiscovered internal face.
529 faces.push_back(c1);
530 }
531 }
532 ++ncf;
533 }
534 num_faces[c0] = ncf;
535 fpos.push_back(int(faces.size()));
536 max_ncf = std::max(max_ncf, ncf);
537 tot_ncf += ncf;
538 }
539 assert(cellno == nc);
540
541 // Build cumulative face sizes enabling direct insertion of
542 // face indices into cfdata later.
543 std::vector<int> cumul_num_faces(numberOfCells() + 1);
544 cumul_num_faces[0] = 0;
545 std::partial_sum(num_faces.begin(), num_faces.end(), cumul_num_faces.begin() + 1);
546
547 // Avoid (most) allocation(s) inside 'c' loop.
548 std::vector<int> l2g;
549 l2g.reserve(max_ncf);
550 std::vector<double> cfdata(tot_ncf);
551 int total_num_faces = int(faces.size());
552
553 // Second pass: build cell-to-face mapping, including boundary.
554 typedef std::vector<int>::iterator VII;
555 for (CI c = cellbegin(); c != cellend(); ++c) {
556 const int c0 = c->index();
557 assert ((0 <= c0 ) && ( c0 < nc) &&
558 (0 <= cell[c0]) && (cell[c0] < nc));
559 const int ncf = num_cf[cell[c0]];
560 l2g.resize(ncf, 0);
561 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
562 if (f->boundary()) {
563 // External, not counted before. Add new face...
564 l2g[f->localIndex()] = total_num_faces++;
565 } else {
566 // Internal face. Need to determine during
567 // traversal of which cell we discovered this
568 // face first, and extract the face number
569 // from the 'faces' table range of that cell.
570
571 // Note: std::find() below is potentially
572 // *VERY* expensive (e.g., large number of
573 // seeks in moderately sized data in case of
574 // faulted cells).
575 const int c1 = f->neighbourCellIndex();
576 assert ((0 <= c1 ) && ( c1 < nc) &&
577 (0 <= cell[c1]) && (cell[c1] < nc));
578
579 int t = c0, seek = c1;
580 if (cell[seek] < cell[t])
581 std::swap(t, seek);
582 int s = fpos[cell[t]], e = fpos[cell[t] + 1];
583 VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
584 assert(p != faces.begin() + e);
585 l2g[f->localIndex()] = p - faces.begin();
586 }
587 }
588 assert(int(l2g.size()) == num_faces[c0]);
589 std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
590 }
591 num_faces_ = total_num_faces;
592 max_faces_per_cell_ = max_ncf;
593 face_indices_.assign(cfdata.begin(), cfdata.end(),
594 num_faces.begin(), num_faces.end());
595
596#ifdef VERBOSE
597 clock.stop();
598 double elapsed = clock.secsSinceStart();
599 std::cout << "done. Time elapsed: " << elapsed << std::endl;
600#endif
601 }
602
603 };
604
605
606
607} // namespace Opm
608
609
610#endif // OPENRS_GRIDINTERFACEEULER_HEADER
Definition GridInterfaceEuler.hpp:368
const CellIterator & dereference() const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:382
bool equal(const CellIterator &other) const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:387
void increment()
Used by iterator facade.
Definition GridInterfaceEuler.hpp:392
Definition GridInterfaceEuler.hpp:310
Intersection (face) iterator for solver-near grid interface.
Definition GridInterfaceEuler.hpp:241
bool equal(const FaceIterator &other) const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:282
void increment()
Used by iterator facade.
Definition GridInterfaceEuler.hpp:290
Face< GridInterface >::DuneIntersectionIter DuneIntersectionIter
Type of low-level intersection iterator.
Definition GridInterfaceEuler.hpp:249
bool operator<(const FaceIterator &other) const
Gives an ordering of intersectionIterators.
Definition GridInterfaceEuler.hpp:297
const FaceIterator & dereference() const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:276
FaceIterator(const GridInterface &grid, const DuneIntersectionIter &it, const int local_index)
Constructor.
Definition GridInterfaceEuler.hpp:268
FaceIterator()
Default constructor.
Definition GridInterfaceEuler.hpp:252
Definition GridInterfaceEuler.hpp:88
Vector normal() const
Definition GridInterfaceEuler.hpp:133
Vector centroid() const
Definition GridInterfaceEuler.hpp:124
Index neighbourCellIndex() const
Definition GridInterfaceEuler.hpp:182
bool boundary() const
Definition GridInterfaceEuler.hpp:142
Cell cell() const
Definition GridInterfaceEuler.hpp:158
Index cellIndex() const
Definition GridInterfaceEuler.hpp:166
Cell neighbourCell() const
Definition GridInterfaceEuler.hpp:174
int boundaryId() const
Definition GridInterfaceEuler.hpp:150
Index localIndex() const
Definition GridInterfaceEuler.hpp:203
Scalar neighbourCellVolume() const
Definition GridInterfaceEuler.hpp:211
Scalar area() const
Definition GridInterfaceEuler.hpp:116
Index index() const
Definition GridInterfaceEuler.hpp:195
Definition GridInterfaceEuler.hpp:403
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
A mapper for Dune::CpGrid cells only.
Definition GridInterfaceEuler.hpp:71
Mapper for general grids.
Definition GridInterfaceEuler.hpp:65