My Project
Loading...
Searching...
No Matches
vcfvstencil.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_VCFV_STENCIL_HH
29#define EWOMS_VCFV_STENCIL_HH
30
32
33#include <dune/grid/common/intersectioniterator.hh>
34#include <dune/grid/common/mcmgmapper.hh>
35#include <dune/geometry/referenceelements.hh>
36
37#if HAVE_DUNE_LOCALFUNCTIONS
38#include <dune/localfunctions/lagrange/pqkfactory.hh>
39#endif // HAVE_DUNE_LOCALFUNCTIONS
40
41#include <dune/common/version.hh>
42
43#include <stdexcept>
44#include <vector>
45
46namespace Opm {
47
52{
53 none,
54 simplex,
55 cube,
56};
57
61template <class Scalar, unsigned dim, unsigned basicGeomType>
63
65// local geometries for 1D elements
67template <class Scalar>
68class VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>
69{
70 enum { dim = 1 };
71 enum { numScv = 2 };
72
73public:
74 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
75
76 static void init()
77 {
78 // 1D LINE SEGMENTS
79 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] = {
80 { // corners of the first sub control volume
81 { 0.0 },
82 { 0.5 }
83 },
84
85 { // corners of the second sub control volume
86 { 0.5 },
87 { 1.0 }
88 }
89 };
90 for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
91 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
92 }
93
94 static const ScvLocalGeometry& get(unsigned scvIdx)
95 { return scvGeoms_[scvIdx]; }
96
97private:
98 static ScvLocalGeometry scvGeoms_[numScv];
99};
100
101template <class Scalar>
102typename VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::ScvLocalGeometry
103VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::scvGeoms_[
104 VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::numScv];
105
106template <class Scalar>
107class VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::simplex>
108{
109 enum { dim = 1 };
110 enum { numScv = 2 };
111
112public:
113 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
114
115 static const ScvLocalGeometry& get(unsigned)
116 {
117 throw std::logic_error("Not implemented: VcfvScvGeometries<Scalar, 1, ElementType::simplex>");
118 }
119};
120
122// local geometries for 2D elements
124template <class Scalar>
125class VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>
126{
127 enum { dim = 2 };
128 enum { numScv = 3 };
129
130public:
131 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
132
133 static const ScvLocalGeometry& get(unsigned scvIdx)
134 { return scvGeoms_[scvIdx]; }
135
136 static void init()
137 {
138 // 2D SIMPLEX
139 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
140 {
141 { // SCV 0 corners
142 { 0.0, 0.0 },
143 { 1.0/2, 0.0 },
144 { 1.0/3, 1.0/3 },
145 { 0.0, 1.0/2 },
146 },
147
148 { // SCV 1 corners
149 { 1.0/2, 0.0 },
150 { 1.0, 0.0 },
151 { 1.0/3, 1.0/3 },
152 { 1.0/2, 1.0/2 },
153 },
154
155 { // SCV 2 corners
156 { 0.0, 1.0/2 },
157 { 1.0/3, 1.0/3 },
158 { 0.0, 1.0 },
159 { 1.0/2, 1.0/2 },
160 }
161 };
162
163 for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
164 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
165 }
166
167private:
168 static ScvLocalGeometry scvGeoms_[numScv];
169};
170
171template <class Scalar>
172typename VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::ScvLocalGeometry
173VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::scvGeoms_[
174 VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::numScv];
175
176template <class Scalar>
177class VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>
178{
179 enum { dim = 2 };
180 enum { numScv = 4 };
181
182public:
183 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
184
185 static const ScvLocalGeometry& get(unsigned scvIdx)
186 { return scvGeoms_[scvIdx]; }
187
188 static void init()
189 {
190 // 2D CUBE
191 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
192 {
193 { // SCV 0 corners
194 { 0.0, 0.0 },
195 { 0.5, 0.0 },
196 { 0.0, 0.5 },
197 { 0.5, 0.5 }
198 },
199
200 { // SCV 1 corners
201 { 0.5, 0.0 },
202 { 1.0, 0.0 },
203 { 0.5, 0.5 },
204 { 1.0, 0.5 }
205 },
206
207 { // SCV 2 corners
208 { 0.0, 0.5 },
209 { 0.5, 0.5 },
210 { 0.0, 1.0 },
211 { 0.5, 1.0 }
212 },
213
214 { // SCV 3 corners
215 { 0.5, 0.5 },
216 { 1.0, 0.5 },
217 { 0.5, 1.0 },
218 { 1.0, 1.0 }
219 }
220 };
221
222 for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
223 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
224 }
225
226 static ScvLocalGeometry scvGeoms_[numScv];
227};
228
229template <class Scalar>
230typename VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::ScvLocalGeometry
231VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::scvGeoms_[
232 VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::numScv];
233
235// local geometries for 3D elements
237template <class Scalar>
238class VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>
239{
240 enum { dim = 3 };
241 enum { numScv = 4 };
242
243public:
244 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
245
246 static const ScvLocalGeometry& get(unsigned scvIdx)
247 { return scvGeoms_[scvIdx]; }
248
249 static void init()
250 {
251 // 3D SIMPLEX
252 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
253 {
254 { // SCV 0 corners
255 { 0.0, 0.0, 0.0 },
256 { 1.0/2, 0.0, 0.0 },
257 { 0.0, 1.0/2, 0.0 },
258 { 1.0/3, 1.0/3, 0.0 },
259
260 { 0.0, 0.0, 0.5 },
261 { 1.0/3, 0.0, 1.0/3 },
262 { 0.0, 1.0/3, 1.0/3 },
263 { 1.0/4, 1.0/4, 1.0/4 },
264 },
265
266 { // SCV 1 corners
267 { 1.0/2, 0.0, 0.0 },
268 { 1.0, 0.0, 0.0 },
269 { 1.0/3, 1.0/3, 0.0 },
270 { 1.0/2, 1.0/2, 0.0 },
271
272 { 1.0/3, 0.0, 1.0/3 },
273 { 1.0/2, 0.0, 1.0/2 },
274 { 1.0/4, 1.0/4, 1.0/4 },
275 { 1.0/3, 1.0/3, 1.0/3 },
276 },
277
278 { // SCV 2 corners
279 { 0.0, 1.0/2, 0.0 },
280 { 1.0/3, 1.0/3, 0.0 },
281 { 0.0, 1.0, 0.0 },
282 { 1.0/2, 1.0/2, 0.0 },
283
284 { 0.0, 1.0/3, 1.0/3 },
285 { 1.0/4, 1.0/4, 1.0/4 },
286 { 0.0, 1.0/2, 1.0/2 },
287 { 1.0/3, 1.0/3, 1.0/3 },
288 },
289
290 { // SCV 3 corners
291 { 0.0, 0.0, 1.0/2 },
292 { 1.0/3, 0.0, 1.0/3 },
293 { 0.0, 1.0/3, 1.0/3 },
294 { 1.0/4, 1.0/4, 1.0/4 },
295
296 { 0.0, 0.0, 1.0 },
297 { 1.0/2, 0.0, 1.0/2 },
298 { 0.0, 1.0/2, 1.0/2 },
299 { 1.0/3, 1.0/3, 1.0/3 },
300 }
301 };
302
303 for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
304 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
305 }
306
307private:
308 static ScvLocalGeometry scvGeoms_[numScv];
309};
310
311template <class Scalar>
312typename VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::ScvLocalGeometry
313VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::scvGeoms_[
314 VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::numScv];
315
316template <class Scalar>
317class VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>
318{
319 enum { dim = 3 };
320 enum { numScv = 8 };
321
322public:
323 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
324
325 static const ScvLocalGeometry& get(unsigned scvIdx)
326 { return scvGeoms_[scvIdx]; }
327
328 static void init()
329 {
330 // 3D CUBE
331 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
332 {
333 { // SCV 0 corners
334 { 0.0, 0.0, 0.0 },
335 { 1.0/2, 0.0, 0.0 },
336 { 0.0, 1.0/2, 0.0 },
337 { 1.0/2, 1.0/2, 0.0 },
338
339 { 0.0, 0.0, 1.0/2 },
340 { 1.0/2, 0.0, 1.0/2 },
341 { 0.0, 1.0/2, 1.0/2 },
342 { 1.0/2, 1.0/2, 1.0/2 },
343 },
344
345 { // SCV 1 corners
346 { 1.0/2, 0.0, 0.0 },
347 { 1.0, 0.0, 0.0 },
348 { 1.0/2, 1.0/2, 0.0 },
349 { 1.0, 1.0/2, 0.0 },
350
351 { 1.0/2, 0.0, 1.0/2 },
352 { 1.0, 0.0, 1.0/2 },
353 { 1.0/2, 1.0/2, 1.0/2 },
354 { 1.0, 1.0/2, 1.0/2 },
355 },
356
357 { // SCV 2 corners
358 { 0.0, 1.0/2, 0.0 },
359 { 1.0/2, 1.0/2, 0.0 },
360 { 0.0, 1.0, 0.0 },
361 { 1.0/2, 1.0, 0.0 },
362
363 { 0.0, 1.0/2, 1.0/2 },
364 { 1.0/2, 1.0/2, 1.0/2 },
365 { 0.0, 1.0, 1.0/2 },
366 { 1.0/2, 1.0, 1.0/2 },
367 },
368
369 { // SCV 3 corners
370 { 1.0/2, 1.0/2, 0.0 },
371 { 1.0, 1.0/2, 0.0 },
372 { 1.0/2, 1.0, 0.0 },
373 { 1.0, 1.0, 0.0 },
374
375 { 1.0/2, 1.0/2, 1.0/2 },
376 { 1.0, 1.0/2, 1.0/2 },
377 { 1.0/2, 1.0, 1.0/2 },
378 { 1.0, 1.0, 1.0/2 },
379 },
380
381 { // SCV 4 corners
382 { 0.0, 0.0, 1.0/2 },
383 { 1.0/2, 0.0, 1.0/2 },
384 { 0.0, 1.0/2, 1.0/2 },
385 { 1.0/2, 1.0/2, 1.0/2 },
386
387 { 0.0, 0.0, 1.0 },
388 { 1.0/2, 0.0, 1.0 },
389 { 0.0, 1.0/2, 1.0 },
390 { 1.0/2, 1.0/2, 1.0 },
391 },
392
393 { // SCV 5 corners
394 { 1.0/2, 0.0, 1.0/2 },
395 { 1.0, 0.0, 1.0/2 },
396 { 1.0/2, 1.0/2, 1.0/2 },
397 { 1.0, 1.0/2, 1.0/2 },
398
399 { 1.0/2, 0.0, 1.0 },
400 { 1.0, 0.0, 1.0 },
401 { 1.0/2, 1.0/2, 1.0 },
402 { 1.0, 1.0/2, 1.0 },
403 },
404
405 { // SCV 6 corners
406 { 0.0, 1.0/2, 1.0/2 },
407 { 1.0/2, 1.0/2, 1.0/2 },
408 { 0.0, 1.0, 1.0/2 },
409 { 1.0/2, 1.0, 1.0/2 },
410
411 { 0.0, 1.0/2, 1.0 },
412 { 1.0/2, 1.0/2, 1.0 },
413 { 0.0, 1.0, 1.0 },
414 { 1.0/2, 1.0, 1.0 },
415 },
416
417 { // SCV 7 corners
418 { 1.0/2, 1.0/2, 1.0/2 },
419 { 1.0, 1.0/2, 1.0/2 },
420 { 1.0/2, 1.0, 1.0/2 },
421 { 1.0, 1.0, 1.0/2 },
422
423 { 1.0/2, 1.0/2, 1.0 },
424 { 1.0, 1.0/2, 1.0 },
425 { 1.0/2, 1.0, 1.0 },
426 { 1.0, 1.0, 1.0 },
427 },
428 };
429
430 for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
431 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
432 }
433private:
434 static ScvLocalGeometry scvGeoms_[numScv];
435};
436
437template <class Scalar>
438typename VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::ScvLocalGeometry
439VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::scvGeoms_[
440 VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::numScv];
441
465template <class Scalar, class GridView>
467{
468 enum{dim = GridView::dimension};
469 enum{dimWorld = GridView::dimensionworld};
470 enum{maxNC = (dim < 3 ? 4 : 8)};
471 enum{maxNE = (dim < 3 ? 4 : 12)};
472 enum{maxNF = (dim < 3 ? 1 : 6)};
473 enum{maxBF = (dim < 3 ? 8 : 24)};
474 using CoordScalar = typename GridView::ctype;
475 using Element = typename GridView::Traits::template Codim<0>::Entity ;
476public:
477 using Entity = typename GridView::Traits::template Codim<dim>::Entity ;
478private:
479 using Geometry = typename Element::Geometry;
480 using DimVector = Dune::FieldVector<Scalar,dimWorld>;
481 using GlobalPosition = Dune::FieldVector<CoordScalar,dimWorld>;
482 using LocalPosition = Dune::FieldVector<CoordScalar,dim>;
483 using IntersectionIterator = typename GridView::IntersectionIterator;
484
486
487#if HAVE_DUNE_LOCALFUNCTIONS
488 using LocalFiniteElementCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
489 using LocalFiniteElement = typename LocalFiniteElementCache::FiniteElementType;
490 using LocalBasisTraits = typename LocalFiniteElement::Traits::LocalBasisType::Traits;
491 using ShapeJacobian = typename LocalBasisTraits::JacobianType;
492#endif // HAVE_DUNE_LOCALFUNCTIONS
493
494 Scalar quadrilateralArea(const GlobalPosition& p0,
495 const GlobalPosition& p1,
496 const GlobalPosition& p2,
497 const GlobalPosition& p3)
498 {
499 return 0.5*std::abs((p3[0] - p1[0])*(p2[1] - p0[1]) - (p3[1] - p1[1])*(p2[0] - p0[0]));
500 }
501
502 void crossProduct(DimVector& c, const DimVector& a, const DimVector& b)
503 {
504 c[0] = a[1]*b[2] - a[2]*b[1];
505 c[1] = a[2]*b[0] - a[0]*b[2];
506 c[2] = a[0]*b[1] - a[1]*b[0];
507 }
508
509 Scalar pyramidVolume(const GlobalPosition& p0,
510 const GlobalPosition& p1,
511 const GlobalPosition& p2,
512 const GlobalPosition& p3,
513 const GlobalPosition& p4)
514 {
515 DimVector a(p2); a -= p0;
516 DimVector b(p3); b -= p1;
517
518 DimVector n;
519 crossProduct(n, a, b);
520
521 a = p4; a -= p0;
522
523 return 1.0/6.0*(n*a);
524 }
525
526 Scalar prismVolume(const GlobalPosition& p0,
527 const GlobalPosition& p1,
528 const GlobalPosition& p2,
529 const GlobalPosition& p3,
530 const GlobalPosition& p4,
531 const GlobalPosition& p5)
532 {
533 DimVector a(p4);
534 for (unsigned k = 0; k < dimWorld; ++k)
535 a[k] -= p0[k];
536 DimVector b(p1);
537 for (unsigned k = 0; k < dimWorld; ++k)
538 b[k] -= p3[k];
539 DimVector m;
540 crossProduct(m, a, b);
541
542 for (unsigned k = 0; k < dimWorld; ++k)
543 a[k] = p1[k] - p0[k];
544 for (unsigned k = 0; k < dimWorld; ++k)
545 b[k] = p2[k] - p0[k];
546 DimVector n;
547 crossProduct(n, a, b);
548 n += m;
549
550 for (unsigned k = 0; k < dimWorld; ++k)
551 a[k] = p5[k] - p0[k];
552
553 return std::abs(1.0/6.0*(n*a));
554 }
555
556 Scalar hexahedronVolume(const GlobalPosition& p0,
557 const GlobalPosition& p1,
558 const GlobalPosition& p2,
559 const GlobalPosition& p3,
560 const GlobalPosition& p4,
561 const GlobalPosition& p5,
562 const GlobalPosition& p6,
563 const GlobalPosition& p7)
564 {
565 return
566 prismVolume(p0,p1,p2,p4,p5,p6)
567 + prismVolume(p0,p2,p3,p4,p6,p7);
568 }
569
570 void normalOfQuadrilateral3D(DimVector& normal,
571 const GlobalPosition& p0,
572 const GlobalPosition& p1,
573 const GlobalPosition& p2,
574 const GlobalPosition& p3)
575 {
576 DimVector a(p2);
577 for (unsigned k = 0; k < dimWorld; ++k)
578 a[k] -= p0[k];
579 DimVector b(p3);
580 for (unsigned k = 0; k < dimWorld; ++k)
581 b[k] -= p1[k];
582
583 crossProduct(normal, a, b);
584 normal *= 0.5;
585 }
586
587 Scalar quadrilateralArea3D(const GlobalPosition& p0,
588 const GlobalPosition& p1,
589 const GlobalPosition& p2,
590 const GlobalPosition& p3)
591 {
592 DimVector normal;
593 normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
594 return normal.two_norm();
595 }
596
597 void getFaceIndices(unsigned numElemVertices, unsigned k, unsigned& leftFace, unsigned& rightFace)
598 {
599 static const unsigned edgeToFaceTet[2][6] = {
600 {1, 0, 3, 2, 1, 3},
601 {0, 2, 0, 1, 3, 2}
602 };
603 static const unsigned edgeToFacePyramid[2][8] = {
604 {1, 2, 3, 4, 1, 3, 4, 2},
605 {0, 0, 0, 0, 3, 2, 1, 4}
606 };
607 static const unsigned edgeToFacePrism[2][9] = {
608 {1, 0, 2, 0, 1, 2, 4, 4, 4},
609 {0, 2, 1, 3, 3, 3, 0, 1, 2}
610 };
611 static const unsigned edgeToFaceHex[2][12] = {
612 {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
613 {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
614 };
615
616 switch (numElemVertices) {
617 case 4:
620 break;
621 case 5:
624 break;
625 case 6:
628 break;
629 case 8:
632 break;
633 default:
634 throw std::logic_error("Not implemented: VcfvStencil::getFaceIndices for "
635 +std::to_string(numElemVertices)+" vertices");
636 break;
637 }
638 }
639
640 void getEdgeIndices(unsigned numElemVertices, unsigned face, unsigned vert, unsigned& leftEdge, unsigned& rightEdge)
641 {
642 static const int faceAndVertexToLeftEdgeTet[4][4] = {
643 { 0, 0, 2, -1},
644 { 0, 0, -1, 3},
645 { 1, -1, 1, 3},
646 {-1, 2, 2, 4}
647 };
648 static const int faceAndVertexToRightEdgeTet[4][4] = {
649 { 1, 2, 1, -1},
650 { 3, 4, -1, 4},
651 { 3, -1, 5, 5},
652 {-1, 4, 5, 5}
653 };
654 static const int faceAndVertexToLeftEdgePyramid[5][5] = {
655 { 0, 2, 3, 1, -1},
656 { 0, -1, 0, -1, 4},
657 {-1, 1, -1, 1, 5},
658 { 2, 2, -1, -1, 4},
659 {-1, -1, 3, 3, 7}
660 };
661 static const int faceAndVertexToRightEdgePyramid[5][5] = {
662 { 2, 1, 0, 3, -1},
663 { 4, -1, 6, -1, 6},
664 {-1, 5, -1, 7, 7},
665 { 4, 5, -1, -1, 5},
666 {-1, -1, 6, 7, 6}
667 };
668 static const int faceAndVertexToLeftEdgePrism[5][6] = {
669 { 3, 3, -1, 0, 1, -1},
670 { 4, -1, 4, 0, -1, 2},
671 {-1, 5, 5, -1, 1, 2},
672 { 3, 3, 5, -1, -1, -1},
673 {-1, -1, -1, 6, 6, 8}
674 };
675 static const int faceAndVertexToRightEdgePrism[5][6] = {
676 { 0, 1, -1, 6, 6, -1},
677 { 0, -1, 2, 7, -1, 7},
678 {-1, 1, 2, -1, 8, 8},
679 { 4, 5, 4, -1, -1, -1},
680 {-1, -1, -1, 7, 8, 7}
681 };
682 static const int faceAndVertexToLeftEdgeHex[6][8] = {
683 { 0, -1, 4, -1, 8, -1, 2, -1},
684 {-1, 5, -1, 3, -1, 1, -1, 9},
685 { 6, 1, -1, -1, 0, 10, -1, -1},
686 {-1, -1, 2, 7, -1, -1, 11, 3},
687 { 4, 6, 7, 5, -1, -1, -1, -1},
688 {-1, -1, -1, -1, 10, 9, 8, 11}
689 };
690 static const int faceAndVertexToRightEdgeHex[6][8] = {
691 { 4, -1, 2, -1, 0, -1, 8, -1},
692 {-1, 1, -1, 5, -1, 9, -1, 3},
693 { 0, 6, -1, -1, 10, 1, -1, -1},
694 {-1, -1, 7, 3, -1, -1, 2, 11},
695 { 6, 5, 4, 7, -1, -1, -1, -1},
696 {-1, -1, -1, -1, 8, 10, 11, 9}
697 };
698
699 switch (numElemVertices) {
700 case 4:
701 leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgeTet[face][vert]);
702 rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgeTet[face][vert]);
703 break;
704 case 5:
705 leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgePyramid[face][vert]);
706 rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgePyramid[face][vert]);
707 break;
708 case 6:
709 leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgePrism[face][vert]);
710 rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgePrism[face][vert]);
711 break;
712 case 8:
713 leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgeHex[face][vert]);
714 rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgeHex[face][vert]);
715 break;
716 default:
717 throw std::logic_error("Not implemented: VcfvStencil::getFaceIndices for "
718 +std::to_string(numElemVertices)+" vertices");
719 break;
720 }
721 }
722
723public:
725 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
726
728 {
729 public:
730 const GlobalPosition center() const
731 { return global(localGeometry_->center()); }
732
733 const GlobalPosition corner(unsigned cornerIdx) const
734 { return global(localGeometry_->corner(cornerIdx)); }
735
736 const GlobalPosition global(const LocalPosition& localPos) const
737 { return element_->geometry().global(localPos); }
738
739 const ScvLocalGeometry& localGeometry() const
740 { return *localGeometry_; }
741
742 const ScvLocalGeometry *localGeometry_;
743 const Element *element_;
744 };
745
747 {
748 const GlobalPosition& globalPos() const
749 { return global; }
750
751 const GlobalPosition center() const
752 { return geometry_.center(); }
753
754 Scalar volume() const
755 { return volume_; }
756
757 const ScvLocalGeometry& localGeometry() const
758 { return geometry_.localGeometry(); }
759
760 const ScvGeometry& geometry() const
761 { return geometry_; }
762
764 LocalPosition local;
766 GlobalPosition global;
768 Scalar volume_;
771
773 Dune::FieldVector<DimVector, maxNC> gradCenter;
774 };
775
777 {
778 const DimVector& normal() const
779 { return normal_; }
780
781 unsigned short interiorIndex() const
782 { return i; }
783
784 unsigned short exteriorIndex() const
785 { return j; }
786
787 Scalar area() const
788 { return area_; }
789
790 const LocalPosition& localPos() const
791 { return ipLocal_; }
792
793 const GlobalPosition& integrationPos() const
794 { return ipGlobal_; }
795
797 unsigned short i,j;
799 LocalPosition ipLocal_;
801 GlobalPosition ipGlobal_;
803 DimVector normal_;
805 Scalar area_;
806 };
807
810
811 VcfvStencil(const GridView& gridView, const Mapper& mapper)
812 : gridView_(gridView)
813 , vertexMapper_(mapper )
814 , element_(*gridView.template begin</*codim=*/0>())
815 {
816 // try to check if the mapper really maps the vertices
817 assert(static_cast<int>(gridView.size(/*codim=*/dimWorld)) == static_cast<int>(mapper.size()));
818
819 static bool localGeometriesInitialized = false;
822
823 VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::init();
824 VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::init();
825 VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::init();
826 VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::init();
827 VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::init();
828 }
829 }
830
836 void updateTopology(const Element& e)
837 {
838 element_ = e;
839
840 numVertices = e.subEntities(/*codim=*/dim);
841 numEdges = e.subEntities(/*codim=*/dim-1);
842 if constexpr (dim == 3) {
843 numFaces = e.subEntities(/*codim=*/1);
844 }
845 else {
846 numFaces = 0;
847 }
848
849 numBoundarySegments_ = 0; // TODO: really required here(?)
850
851 // compute the local and global coordinates of the element
852 const Geometry& geometry = e.geometry();
853 geometryType_ = geometry.type();
854 const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
855 for (unsigned vertexIdx = 0; vertexIdx < numVertices; vertexIdx++) {
856 subContVol[vertexIdx].local = referenceElement.position(static_cast<int>(vertexIdx), dim);
857 subContVol[vertexIdx].global = geometry.corner(static_cast<int>(vertexIdx));
858 }
859 }
860
861 void updatePrimaryTopology(const Element& element)
862 {
863 // since all degrees of freedom in a stencil are "primary" DOFs for the
864 // vertex-centered finite volume method, there's no difference to
865 // updateTopology()
866 updateTopology(element);
867 }
868
869 void update(const Element& e)
870 {
872
873 const Geometry& geometry = e.geometry();
874 geometryType_ = geometry.type();
875
876 const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
877
878 elementVolume = geometry.volume();
879 elementLocal = referenceElement.position(0,0);
880 elementGlobal = geometry.global(elementLocal);
881
882 // corners:
883 for (unsigned vert = 0; vert < numVertices; vert++) {
884 subContVol[vert].local = referenceElement.position(static_cast<int>(vert), dim);
885 subContVol[vert].global = geometry.global(subContVol[vert].local);
886 }
887
888 // edges:
889 for (unsigned edge = 0; edge < numEdges; edge++) {
890 edgeCoord[edge] = geometry.global(referenceElement.position(static_cast<int>(edge), dim-1));
891 }
892
893 // faces:
894 for (unsigned face = 0; face < numFaces; face++) {
895 faceCoord[face] = geometry.global(referenceElement.position(static_cast<int>(face), 1));
896 }
897
898 // fill sub control volume data use specialization for this
899 // \todo maybe it would be a good idea to migrate everything
900 // which is dependend of the grid's dimension to
901 // _VcfvFVElemGeomHelper in order to benefit from more aggressive
902 // compiler optimizations...
903 fillSubContVolData_();
904
905 // fill sub control volume face data:
906 for (unsigned k = 0; k < numEdges; k++) { // begin loop over edges / sub control volume faces
907 unsigned short i = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 0, dim));
908 unsigned short j = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 1, dim));
909 if (numEdges == 4 && (i == 2 || j == 2))
910 std::swap(i, j);
911 subContVolFace[k].i = i;
912 subContVolFace[k].j = j;
913
914 // calculate the local integration point and
915 // the face normal. note that since dim is a
916 // constant which is known at compile time
917 // the compiler can optimize away all if
918 // cases which don't apply.
919 LocalPosition ipLocal_;
920 DimVector diffVec;
921 if constexpr (dim == 1) {
922 subContVolFace[k].ipLocal_ = 0.5;
923 subContVolFace[k].normal_ = 1.0;
924 subContVolFace[k].area_ = 1.0;
925 ipLocal_ = subContVolFace[k].ipLocal_;
926 }
927 else if constexpr (dim == 2) {
928 ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal;
929 ipLocal_ *= 0.5;
930 subContVolFace[k].ipLocal_ = ipLocal_;
931 for (unsigned m = 0; m < dimWorld; ++m)
932 diffVec[m] = elementGlobal[m] - edgeCoord[k][m];
933 subContVolFace[k].normal_[0] = diffVec[1];
934 subContVolFace[k].normal_[1] = -diffVec[0];
935
936 for (unsigned m = 0; m < dimWorld; ++m)
937 diffVec[m] = subContVol[j].global[m] - subContVol[i].global[m];
938 // make sure the normal points to the right direction
939 if (subContVolFace[k].normal_ * diffVec < 0)
940 subContVolFace[k].normal_ *= -1;
941
942 subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
943 subContVolFace[k].normal_ /= subContVolFace[k].area_;
944 }
945 else if constexpr (dim == 3) {
946 unsigned leftFace;
947 unsigned rightFace;
948 getFaceIndices(numVertices, k, leftFace, rightFace);
949 ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal
950 + referenceElement.position(static_cast<int>(leftFace), 1)
951 + referenceElement.position(static_cast<int>(rightFace), 1);
952 ipLocal_ *= 0.25;
953 subContVolFace[k].ipLocal_ = ipLocal_;
954 normalOfQuadrilateral3D(subContVolFace[k].normal_,
955 edgeCoord[k], faceCoord[rightFace],
956 elementGlobal, faceCoord[leftFace]);
957 subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
958 subContVolFace[k].normal_ /= subContVolFace[k].area_;
959 }
960
961 // get the global integration point and the Jacobian inverse
962 subContVolFace[k].ipGlobal_ = geometry.global(ipLocal_);
963 } // end loop over edges / sub control volume faces
964
965 // fill boundary face data:
966 IntersectionIterator endit = gridView_.iend(e);
967 for (IntersectionIterator it = gridView_.ibegin(e); it != endit; ++it) {
968 const typename IntersectionIterator::Intersection& intersection = *it ;
969
970 if ( ! intersection.boundary())
971 continue;
972
973 unsigned face = static_cast<unsigned>(intersection.indexInInside());
974 unsigned numVerticesOfFace = static_cast<unsigned>(referenceElement.size(static_cast<int>(face), 1, dim));
975 for (unsigned vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
976 {
977 unsigned short vertInElement = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(face), 1, static_cast<int>(vertInFace), dim));
978 unsigned bfIdx = numBoundarySegments_;
979 ++numBoundarySegments_;
980
981 if constexpr (dim == 1) {
982 boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim);
983 boundaryFace_[bfIdx].area_ = 1.0;
984 }
985 else if constexpr (dim == 2) {
986 boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
987 + referenceElement.position(static_cast<int>(face), 1);
988 boundaryFace_[bfIdx].ipLocal_ *= 0.5;
989 boundaryFace_[bfIdx].area_ = 0.5 * intersection.geometry().volume();
990 }
991 else if constexpr (dim == 3) {
992 unsigned leftEdge;
993 unsigned rightEdge;
994 getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
995 boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
996 + referenceElement.position(static_cast<int>(face), 1)
997 + referenceElement.position(static_cast<int>(leftEdge), dim-1)
998 + referenceElement.position(static_cast<int>(rightEdge), dim-1);
999 boundaryFace_[bfIdx].ipLocal_ *= 0.25;
1000 boundaryFace_[bfIdx].area_ =
1001 quadrilateralArea3D(subContVol[vertInElement].global,
1002 edgeCoord[rightEdge],
1003 faceCoord[face],
1004 edgeCoord[leftEdge]);
1005 }
1006 else
1007 throw std::logic_error("Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1008
1009 boundaryFace_[bfIdx].ipGlobal_ = geometry.global(boundaryFace_[bfIdx].ipLocal_);
1010 boundaryFace_[bfIdx].i = vertInElement;
1011 boundaryFace_[bfIdx].j = vertInElement;
1012
1013 // ASSUME constant normal on the segment of the boundary face
1014 boundaryFace_[bfIdx].normal_ = intersection.centerUnitOuterNormal();
1015 }
1016 }
1017
1018 updateScvGeometry(e);
1019 }
1020
1021 void updateScvGeometry(const Element& element)
1022 {
1023 auto geomType = element.geometry().type();
1024
1025 // get the local geometries of the sub control volumes
1026 if (geomType.isTriangle() || geomType.isTetrahedron()) {
1027 for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1028 subContVol[vertIdx].geometry_.element_ = &element;
1029 subContVol[vertIdx].geometry_.localGeometry_ =
1030 &VcfvScvGeometries<Scalar, dim, ElementType::simplex>::get(vertIdx);
1031 }
1032 }
1033 else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1034 for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1035 subContVol[vertIdx].geometry_.element_ = &element;
1036 subContVol[vertIdx].geometry_.localGeometry_ =
1037 &VcfvScvGeometries<Scalar, dim, ElementType::cube>::get(vertIdx);
1038 }
1039 }
1040 else
1041 throw std::logic_error("Not implemented: SCV geometries for non hexahedron elements");
1042 }
1043
1044#if HAVE_DUNE_LOCALFUNCTIONS
1045 void updateCenterGradients()
1046 {
1047 const auto& localFiniteElement = feCache_.get(element_.type());
1048 const auto& geom = element_.geometry();
1049
1050 std::vector<ShapeJacobian> localJac;
1051
1052 for (unsigned scvIdx = 0; scvIdx < numVertices; ++ scvIdx) {
1053 const auto& localCenter = subContVol[scvIdx].localGeometry().center();
1054 localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
1055 const auto& globalPos = subContVol[scvIdx].geometry().center();
1056
1057 const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
1058 for (unsigned vert = 0; vert < numVertices; vert++) {
1059 jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
1060 }
1061 }
1062 }
1063#endif
1064
1065 unsigned numDof() const
1066 { return numVertices; }
1067
1068 unsigned numPrimaryDof() const
1069 { return numDof(); }
1070
1071 Dune::PartitionType partitionType(unsigned scvIdx) const
1072 { return element_.template subEntity</*codim=*/dim>(scvIdx)->partitionType(); }
1073
1074 const SubControlVolume& subControlVolume(unsigned scvIdx) const
1075 {
1076 assert(scvIdx < numDof());
1077 return subContVol[scvIdx];
1078 }
1079
1080 unsigned numInteriorFaces() const
1081 { return numEdges; }
1082
1083 unsigned numBoundaryFaces() const
1084 { return numBoundarySegments_; }
1085
1086 const SubControlVolumeFace& interiorFace(unsigned faceIdx) const
1087 { return subContVolFace[faceIdx]; }
1088
1089 const BoundaryFace& boundaryFace(unsigned bfIdx) const
1090 { return boundaryFace_[bfIdx]; }
1091
1096 unsigned globalSpaceIndex(unsigned dofIdx) const
1097 {
1098 assert(dofIdx < numDof());
1099
1100 return static_cast<unsigned>(vertexMapper_.subIndex(element_, static_cast<int>(dofIdx), /*codim=*/dim));
1101 }
1102
1107 Entity entity(unsigned dofIdx) const
1108 {
1109 assert(dofIdx < numDof());
1110 return element_.template subEntity<dim>(static_cast<int>(dofIdx));
1111 }
1112
1113private:
1114#if __GNUC__ || __clang__
1115#pragma GCC diagnostic push
1116#pragma GCC diagnostic ignored "-Wpragmas"
1117#pragma GCC diagnostic ignored "-Warray-bounds"
1118#endif
1119 void fillSubContVolData_()
1120 {
1121 if constexpr (dim == 1) {
1122 // 1D
1123 subContVol[0].volume_ = 0.5*elementVolume;
1124 subContVol[1].volume_ = 0.5*elementVolume;
1125 }
1126 else if constexpr (dim == 2) {
1127 switch (numVertices) {
1128 case 3: // 2D, triangle
1129 subContVol[0].volume_ = elementVolume/3;
1130 subContVol[1].volume_ = elementVolume/3;
1131 subContVol[2].volume_ = elementVolume/3;
1132 break;
1133 case 4: // 2D, quadrilinear
1134 subContVol[0].volume_ =
1135 quadrilateralArea(subContVol[0].global,
1136 edgeCoord[2],
1137 elementGlobal,
1138 edgeCoord[0]);
1139 subContVol[1].volume_ =
1140 quadrilateralArea(subContVol[1].global,
1141 edgeCoord[1],
1142 elementGlobal,
1143 edgeCoord[2]);
1144 subContVol[2].volume_ =
1145 quadrilateralArea(subContVol[2].global,
1146 edgeCoord[0],
1147 elementGlobal,
1148 edgeCoord[3]);
1149 subContVol[3].volume_ =
1150 quadrilateralArea(subContVol[3].global,
1151 edgeCoord[3],
1152 elementGlobal,
1153 edgeCoord[1]);
1154 break;
1155 default:
1156 throw std::logic_error("Not implemented:VcfvStencil dim = "+std::to_string(dim)
1157 +", numVertices = "+std::to_string(numVertices));
1158 }
1159 }
1160 else if constexpr (dim == 3) {
1161 switch (numVertices) {
1162 case 4: // 3D, tetrahedron
1163 for (unsigned k = 0; k < numVertices; k++)
1164 subContVol[k].volume_ = elementVolume / 4.0;
1165 break;
1166 case 5: // 3D, pyramid
1167 subContVol[0].volume_ =
1168 hexahedronVolume(subContVol[0].global,
1169 edgeCoord[2],
1170 faceCoord[0],
1171 edgeCoord[0],
1172 edgeCoord[4],
1173 faceCoord[3],
1174 elementGlobal,
1175 faceCoord[1]);
1176 subContVol[1].volume_ =
1177 hexahedronVolume(subContVol[1].global,
1178 edgeCoord[1],
1179 faceCoord[0],
1180 edgeCoord[2],
1181 edgeCoord[5],
1182 faceCoord[2],
1183 elementGlobal,
1184 faceCoord[3]);
1185 subContVol[2].volume_ =
1186 hexahedronVolume(subContVol[2].global,
1187 edgeCoord[0],
1188 faceCoord[0],
1189 edgeCoord[3],
1190 edgeCoord[6],
1191 faceCoord[1],
1192 elementGlobal,
1193 faceCoord[4]);
1194 subContVol[3].volume_ =
1195 hexahedronVolume(subContVol[3].global,
1196 edgeCoord[3],
1197 faceCoord[0],
1198 edgeCoord[1],
1199 edgeCoord[7],
1200 faceCoord[4],
1201 elementGlobal,
1202 faceCoord[2]);
1203 subContVol[4].volume_ = elementVolume -
1204 subContVol[0].volume_ -
1205 subContVol[1].volume_ -
1206 subContVol[2].volume_ -
1207 subContVol[3].volume_;
1208 break;
1209 case 6: // 3D, prism
1210 subContVol[0].volume_ =
1211 hexahedronVolume(subContVol[0].global,
1212 edgeCoord[3],
1213 faceCoord[3],
1214 edgeCoord[4],
1215 edgeCoord[0],
1216 faceCoord[0],
1217 elementGlobal,
1218 faceCoord[1]);
1219 subContVol[1].volume_ =
1220 hexahedronVolume(subContVol[1].global,
1221 edgeCoord[5],
1222 faceCoord[3],
1223 edgeCoord[3],
1224 edgeCoord[1],
1225 faceCoord[2],
1226 elementGlobal,
1227 faceCoord[0]);
1228 subContVol[2].volume_ =
1229 hexahedronVolume(subContVol[2].global,
1230 edgeCoord[4],
1231 faceCoord[3],
1232 edgeCoord[5],
1233 edgeCoord[2],
1234 faceCoord[1],
1235 elementGlobal,
1236 faceCoord[2]);
1237 subContVol[3].volume_ =
1238 hexahedronVolume(edgeCoord[0],
1239 faceCoord[0],
1240 elementGlobal,
1241 faceCoord[1],
1242 subContVol[3].global,
1243 edgeCoord[6],
1244 faceCoord[4],
1245 edgeCoord[7]);
1246 subContVol[4].volume_ =
1247 hexahedronVolume(edgeCoord[1],
1248 faceCoord[2],
1249 elementGlobal,
1250 faceCoord[0],
1251 subContVol[4].global,
1252 edgeCoord[8],
1253 faceCoord[4],
1254 edgeCoord[6]);
1255 subContVol[5].volume_ =
1256 hexahedronVolume(edgeCoord[2],
1257 faceCoord[1],
1258 elementGlobal,
1259 faceCoord[2],
1260 subContVol[5].global,
1261 edgeCoord[7],
1262 faceCoord[4],
1263 edgeCoord[8]);
1264 break;
1265 case 8: // 3D, hexahedron
1266 subContVol[0].volume_ =
1267 hexahedronVolume(subContVol[0].global,
1268 edgeCoord[6],
1269 faceCoord[4],
1270 edgeCoord[4],
1271 edgeCoord[0],
1272 faceCoord[2],
1273 elementGlobal,
1274 faceCoord[0]);
1275 subContVol[1].volume_ =
1276 hexahedronVolume(subContVol[1].global,
1277 edgeCoord[5],
1278 faceCoord[4],
1279 edgeCoord[6],
1280 edgeCoord[1],
1281 faceCoord[1],
1282 elementGlobal,
1283 faceCoord[2]);
1284 subContVol[2].volume_ =
1285 hexahedronVolume(subContVol[2].global,
1286 edgeCoord[4],
1287 faceCoord[4],
1288 edgeCoord[7],
1289 edgeCoord[2],
1290 faceCoord[0],
1291 elementGlobal,
1292 faceCoord[3]);
1293 subContVol[3].volume_ =
1294 hexahedronVolume(subContVol[3].global,
1295 edgeCoord[7],
1296 faceCoord[4],
1297 edgeCoord[5],
1298 edgeCoord[3],
1299 faceCoord[3],
1300 elementGlobal,
1301 faceCoord[1]);
1302 subContVol[4].volume_ =
1303 hexahedronVolume(edgeCoord[0],
1304 faceCoord[2],
1305 elementGlobal,
1306 faceCoord[0],
1307 subContVol[4].global,
1308 edgeCoord[10],
1309 faceCoord[5],
1310 edgeCoord[8]);
1311 subContVol[5].volume_ =
1312 hexahedronVolume(edgeCoord[1],
1313 faceCoord[1],
1314 elementGlobal,
1315 faceCoord[2],
1316 subContVol[5].global,
1317 edgeCoord[9],
1318 faceCoord[5],
1319 edgeCoord[10]);
1320 subContVol[6].volume_ =
1321 hexahedronVolume(edgeCoord[2],
1322 faceCoord[0],
1323 elementGlobal,
1324 faceCoord[3],
1325 subContVol[6].global,
1326 edgeCoord[8],
1327 faceCoord[5],
1328 edgeCoord[11]);
1329 subContVol[7].volume_ =
1330 hexahedronVolume(edgeCoord[3],
1331 faceCoord[3],
1332 elementGlobal,
1333 faceCoord[1],
1334 subContVol[7].global,
1335 edgeCoord[11],
1336 faceCoord[5],
1337 edgeCoord[9]);
1338 break;
1339 default:
1340 throw std::logic_error("Not implemented:VcfvStencil for dim = "+std::to_string(dim)
1341 +", numVertices = "+std::to_string(numVertices));
1342 }
1343 }
1344 else
1345 throw std::logic_error("Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1346 }
1347#if __GNUC__ || __clang__
1348#pragma GCC diagnostic pop
1349#endif
1350
1351 const GridView& gridView_;
1352 const Mapper& vertexMapper_;
1353
1354 Element element_;
1355
1356#if HAVE_DUNE_LOCALFUNCTIONS
1358#endif // HAVE_DUNE_LOCALFUNCTIONS
1359
1361 LocalPosition elementLocal;
1363 GlobalPosition elementGlobal;
1365 Scalar elementVolume;
1367 SubControlVolume subContVol[maxNC];
1369 SubControlVolumeFace subContVolFace[maxNE];
1371 BoundaryFace boundaryFace_[maxBF];
1372 unsigned numBoundarySegments_;
1374 GlobalPosition edgeCoord[maxNE];
1376 GlobalPosition faceCoord[maxNF];
1378 unsigned numVertices;
1380 unsigned numEdges;
1382 unsigned numFaces;
1383 Dune::GeometryType geometryType_;
1384};
1385
1386#if HAVE_DUNE_LOCALFUNCTIONS
1387template<class Scalar, class GridView>
1388typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
1389VcfvStencil<Scalar, GridView>::feCache_;
1390#endif // HAVE_DUNE_LOCALFUNCTIONS
1391
1392} // namespace Opm
1393
1394#endif
1395
Quadrature geometry for quadrilaterals.
Definition quadraturegeometries.hh:40
const GlobalPosition & corner(unsigned cornerIdx) const
Return the position of the corner with a given index.
Definition quadraturegeometries.hh:127
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition quadraturegeometries.hh:69
Definition vcfvstencil.hh:728
Represents the finite volume geometry of a single element in the VCFV discretization.
Definition vcfvstencil.hh:467
unsigned globalSpaceIndex(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition vcfvstencil.hh:1096
SubControlVolumeFace BoundaryFace
compatibility alias
Definition vcfvstencil.hh:809
Dune::MultipleCodimMultipleGeomTypeMapper< GridView > Mapper
exported Mapper type
Definition vcfvstencil.hh:725
Entity entity(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition vcfvstencil.hh:1107
void updateTopology(const Element &e)
Update the non-geometric part of the stencil.
Definition vcfvstencil.hh:836
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
ElementType
The types of reference elements available.
Definition vcfvstencil.hh:52
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Quadrature geometry for quadrilaterals.
interior face of a sub control volume
Definition vcfvstencil.hh:777
Scalar area_
area of face
Definition vcfvstencil.hh:805
GlobalPosition ipGlobal_
integration point in global coords
Definition vcfvstencil.hh:801
DimVector normal_
normal on face pointing to CV j or outward of the domain with length equal to |scvf|
Definition vcfvstencil.hh:803
unsigned short i
scvf seperates corner i and j of elem
Definition vcfvstencil.hh:797
LocalPosition ipLocal_
integration point in local coords
Definition vcfvstencil.hh:799
finite volume intersected with element
Definition vcfvstencil.hh:747
ScvGeometry geometry_
The geometry of the sub-control volume in local coordinates.
Definition vcfvstencil.hh:770
LocalPosition local
local vertex position
Definition vcfvstencil.hh:764
GlobalPosition global
global vertex position
Definition vcfvstencil.hh:766
Scalar volume_
space occupied by the sub-control volume
Definition vcfvstencil.hh:768
Dune::FieldVector< DimVector, maxNC > gradCenter
derivative of shape function at the center of the sub control volume
Definition vcfvstencil.hh:773