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
31#include <opm/models/utils/quadraturegeometries.hh>
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
51enum ElementType
52{
53 none,
54 simplex,
55 cube,
56};
57
61template <class Scalar, unsigned dim, unsigned basicGeomType>
62class VcfvScvGeometries;
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>
466class VcfvStencil
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
485 using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;
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:
618 leftFace = edgeToFaceTet[0][k];
619 rightFace = edgeToFaceTet[1][k];
620 break;
621 case 5:
622 leftFace = edgeToFacePyramid[0][k];
623 rightFace = edgeToFacePyramid[1][k];
624 break;
625 case 6:
626 leftFace = edgeToFacePrism[0][k];
627 rightFace = edgeToFacePrism[1][k];
628 break;
629 case 8:
630 leftFace = edgeToFaceHex[0][k];
631 rightFace = edgeToFaceHex[1][k];
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
727 class ScvGeometry
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
746 struct SubControlVolume
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_;
770 ScvGeometry geometry_;
771
773 Dune::FieldVector<DimVector, maxNC> gradCenter;
774 };
775
776 struct SubControlVolumeFace
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
809 using BoundaryFace = SubControlVolumeFace;
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;
820 if (!localGeometriesInitialized) {
821 localGeometriesInitialized = true;
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 numFaces = (dim<3)?0:e.subEntities(/*codim=*/1);
843
844 numBoundarySegments_ = 0; // TODO: really required here(?)
845
846 // compute the local and global coordinates of the element
847 const Geometry& geometry = e.geometry();
848 geometryType_ = geometry.type();
849 const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
850 for (unsigned vertexIdx = 0; vertexIdx < numVertices; vertexIdx++) {
851 subContVol[vertexIdx].local = referenceElement.position(static_cast<int>(vertexIdx), dim);
852 subContVol[vertexIdx].global = geometry.corner(static_cast<int>(vertexIdx));
853 }
854 }
855
856 void updatePrimaryTopology(const Element& element)
857 {
858 // since all degrees of freedom in a stencil are "primary" DOFs for the
859 // vertex-centered finite volume method, there's no difference to
860 // updateTopology()
861 updateTopology(element);
862 }
863
864 void update(const Element& e)
865 {
866 updateTopology(e);
867
868 const Geometry& geometry = e.geometry();
869 geometryType_ = geometry.type();
870
871 const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
872
873 elementVolume = geometry.volume();
874 elementLocal = referenceElement.position(0,0);
875 elementGlobal = geometry.global(elementLocal);
876
877 // corners:
878 for (unsigned vert = 0; vert < numVertices; vert++) {
879 subContVol[vert].local = referenceElement.position(static_cast<int>(vert), dim);
880 subContVol[vert].global = geometry.global(subContVol[vert].local);
881 }
882
883 // edges:
884 for (unsigned edge = 0; edge < numEdges; edge++) {
885 edgeCoord[edge] = geometry.global(referenceElement.position(static_cast<int>(edge), dim-1));
886 }
887
888 // faces:
889 for (unsigned face = 0; face < numFaces; face++) {
890 faceCoord[face] = geometry.global(referenceElement.position(static_cast<int>(face), 1));
891 }
892
893 // fill sub control volume data use specialization for this
894 // \todo maybe it would be a good idea to migrate everything
895 // which is dependend of the grid's dimension to
896 // _VcfvFVElemGeomHelper in order to benefit from more aggressive
897 // compiler optimizations...
898 fillSubContVolData_();
899
900 // fill sub control volume face data:
901 for (unsigned k = 0; k < numEdges; k++) { // begin loop over edges / sub control volume faces
902 unsigned short i = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 0, dim));
903 unsigned short j = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 1, dim));
904 if (numEdges == 4 && (i == 2 || j == 2))
905 std::swap(i, j);
906 subContVolFace[k].i = i;
907 subContVolFace[k].j = j;
908
909 // calculate the local integration point and
910 // the face normal. note that since dim is a
911 // constant which is known at compile time
912 // the compiler can optimize away all if
913 // cases which don't apply.
914 LocalPosition ipLocal_;
915 DimVector diffVec;
916 if (dim==1) {
917 subContVolFace[k].ipLocal_ = 0.5;
918 subContVolFace[k].normal_ = 1.0;
919 subContVolFace[k].area_ = 1.0;
920 ipLocal_ = subContVolFace[k].ipLocal_;
921 }
922 else if (dim==2) {
923 ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal;
924 ipLocal_ *= 0.5;
925 subContVolFace[k].ipLocal_ = ipLocal_;
926 for (unsigned m = 0; m < dimWorld; ++m)
927 diffVec[m] = elementGlobal[m] - edgeCoord[k][m];
928 subContVolFace[k].normal_[0] = diffVec[1];
929 subContVolFace[k].normal_[1] = -diffVec[0];
930
931 for (unsigned m = 0; m < dimWorld; ++m)
932 diffVec[m] = subContVol[j].global[m] - subContVol[i].global[m];
933 // make sure the normal points to the right direction
934 if (subContVolFace[k].normal_ * diffVec < 0)
935 subContVolFace[k].normal_ *= -1;
936
937 subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
938 subContVolFace[k].normal_ /= subContVolFace[k].area_;
939 }
940 else if (dim==3) {
941 unsigned leftFace;
942 unsigned rightFace;
943 getFaceIndices(numVertices, k, leftFace, rightFace);
944 ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal
945 + referenceElement.position(static_cast<int>(leftFace), 1)
946 + referenceElement.position(static_cast<int>(rightFace), 1);
947 ipLocal_ *= 0.25;
948 subContVolFace[k].ipLocal_ = ipLocal_;
949 normalOfQuadrilateral3D(subContVolFace[k].normal_,
950 edgeCoord[k], faceCoord[rightFace],
951 elementGlobal, faceCoord[leftFace]);
952 subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
953 subContVolFace[k].normal_ /= subContVolFace[k].area_;
954 }
955
956 // get the global integration point and the Jacobian inverse
957 subContVolFace[k].ipGlobal_ = geometry.global(ipLocal_);
958 } // end loop over edges / sub control volume faces
959
960 // fill boundary face data:
961 IntersectionIterator endit = gridView_.iend(e);
962 for (IntersectionIterator it = gridView_.ibegin(e); it != endit; ++it) {
963 const typename IntersectionIterator::Intersection& intersection = *it ;
964
965 if ( ! intersection.boundary())
966 continue;
967
968 unsigned face = static_cast<unsigned>(intersection.indexInInside());
969 unsigned numVerticesOfFace = static_cast<unsigned>(referenceElement.size(static_cast<int>(face), 1, dim));
970 for (unsigned vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
971 {
972 unsigned short vertInElement = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(face), 1, static_cast<int>(vertInFace), dim));
973 unsigned bfIdx = numBoundarySegments_;
974 ++numBoundarySegments_;
975
976 if (dim == 1) {
977 boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim);
978 boundaryFace_[bfIdx].area_ = 1.0;
979 }
980 else if (dim == 2) {
981 boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
982 + referenceElement.position(static_cast<int>(face), 1);
983 boundaryFace_[bfIdx].ipLocal_ *= 0.5;
984 boundaryFace_[bfIdx].area_ = 0.5 * intersection.geometry().volume();
985 }
986 else if (dim == 3) {
987 unsigned leftEdge;
988 unsigned rightEdge;
989 getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
990 boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
991 + referenceElement.position(static_cast<int>(face), 1)
992 + referenceElement.position(static_cast<int>(leftEdge), dim-1)
993 + referenceElement.position(static_cast<int>(rightEdge), dim-1);
994 boundaryFace_[bfIdx].ipLocal_ *= 0.25;
995 boundaryFace_[bfIdx].area_ =
996 quadrilateralArea3D(subContVol[vertInElement].global,
997 edgeCoord[rightEdge],
998 faceCoord[face],
999 edgeCoord[leftEdge]);
1000 }
1001 else
1002 throw std::logic_error("Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1003
1004 boundaryFace_[bfIdx].ipGlobal_ = geometry.global(boundaryFace_[bfIdx].ipLocal_);
1005 boundaryFace_[bfIdx].i = vertInElement;
1006 boundaryFace_[bfIdx].j = vertInElement;
1007
1008 // ASSUME constant normal on the segment of the boundary face
1009 boundaryFace_[bfIdx].normal_ = intersection.centerUnitOuterNormal();
1010 }
1011 }
1012
1013 updateScvGeometry(e);
1014 }
1015
1016 void updateScvGeometry(const Element& element)
1017 {
1018 auto geomType = element.geometry().type();
1019
1020 // get the local geometries of the sub control volumes
1021 if (geomType.isTriangle() || geomType.isTetrahedron()) {
1022 for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1023 subContVol[vertIdx].geometry_.element_ = &element;
1024 subContVol[vertIdx].geometry_.localGeometry_ =
1025 &VcfvScvGeometries<Scalar, dim, ElementType::simplex>::get(vertIdx);
1026 }
1027 }
1028 else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1029 for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1030 subContVol[vertIdx].geometry_.element_ = &element;
1031 subContVol[vertIdx].geometry_.localGeometry_ =
1032 &VcfvScvGeometries<Scalar, dim, ElementType::cube>::get(vertIdx);
1033 }
1034 }
1035 else
1036 throw std::logic_error("Not implemented: SCV geometries for non hexahedron elements");
1037 }
1038
1039#if HAVE_DUNE_LOCALFUNCTIONS
1040 void updateCenterGradients()
1041 {
1042 const auto& localFiniteElement = feCache_.get(element_.type());
1043 const auto& geom = element_.geometry();
1044
1045 std::vector<ShapeJacobian> localJac;
1046
1047 for (unsigned scvIdx = 0; scvIdx < numVertices; ++ scvIdx) {
1048 const auto& localCenter = subContVol[scvIdx].localGeometry().center();
1049 localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
1050 const auto& globalPos = subContVol[scvIdx].geometry().center();
1051
1052 const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
1053 for (unsigned vert = 0; vert < numVertices; vert++) {
1054 jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
1055 }
1056 }
1057 }
1058#endif
1059
1060 unsigned numDof() const
1061 { return numVertices; }
1062
1063 unsigned numPrimaryDof() const
1064 { return numDof(); }
1065
1066 Dune::PartitionType partitionType(unsigned scvIdx) const
1067 { return element_.template subEntity</*codim=*/dim>(scvIdx)->partitionType(); }
1068
1069 const SubControlVolume& subControlVolume(unsigned scvIdx) const
1070 {
1071 assert(scvIdx < numDof());
1072 return subContVol[scvIdx];
1073 }
1074
1075 unsigned numInteriorFaces() const
1076 { return numEdges; }
1077
1078 unsigned numBoundaryFaces() const
1079 { return numBoundarySegments_; }
1080
1081 const SubControlVolumeFace& interiorFace(unsigned faceIdx) const
1082 { return subContVolFace[faceIdx]; }
1083
1084 const BoundaryFace& boundaryFace(unsigned bfIdx) const
1085 { return boundaryFace_[bfIdx]; }
1086
1091 unsigned globalSpaceIndex(unsigned dofIdx) const
1092 {
1093 assert(dofIdx < numDof());
1094
1095 return static_cast<unsigned>(vertexMapper_.subIndex(element_, static_cast<int>(dofIdx), /*codim=*/dim));
1096 }
1097
1102 Entity entity(unsigned dofIdx) const
1103 {
1104 assert(dofIdx < numDof());
1105 return element_.template subEntity<dim>(static_cast<int>(dofIdx));
1106 }
1107
1108private:
1109#if __GNUC__ || __clang__
1110#pragma GCC diagnostic push
1111#pragma GCC diagnostic ignored "-Wpragmas"
1112#pragma GCC diagnostic ignored "-Warray-bounds"
1113#endif
1114 void fillSubContVolData_()
1115 {
1116 if (dim == 1) {
1117 // 1D
1118 subContVol[0].volume_ = 0.5*elementVolume;
1119 subContVol[1].volume_ = 0.5*elementVolume;
1120 }
1121 else if (dim == 2) {
1122 switch (numVertices) {
1123 case 3: // 2D, triangle
1124 subContVol[0].volume_ = elementVolume/3;
1125 subContVol[1].volume_ = elementVolume/3;
1126 subContVol[2].volume_ = elementVolume/3;
1127 break;
1128 case 4: // 2D, quadrilinear
1129 subContVol[0].volume_ =
1130 quadrilateralArea(subContVol[0].global,
1131 edgeCoord[2],
1132 elementGlobal,
1133 edgeCoord[0]);
1134 subContVol[1].volume_ =
1135 quadrilateralArea(subContVol[1].global,
1136 edgeCoord[1],
1137 elementGlobal,
1138 edgeCoord[2]);
1139 subContVol[2].volume_ =
1140 quadrilateralArea(subContVol[2].global,
1141 edgeCoord[0],
1142 elementGlobal,
1143 edgeCoord[3]);
1144 subContVol[3].volume_ =
1145 quadrilateralArea(subContVol[3].global,
1146 edgeCoord[3],
1147 elementGlobal,
1148 edgeCoord[1]);
1149 break;
1150 default:
1151 throw std::logic_error("Not implemented:VcfvStencil dim = "+std::to_string(dim)
1152 +", numVertices = "+std::to_string(numVertices));
1153 }
1154 }
1155 else if (dim == 3) {
1156 switch (numVertices) {
1157 case 4: // 3D, tetrahedron
1158 for (unsigned k = 0; k < numVertices; k++)
1159 subContVol[k].volume_ = elementVolume / 4.0;
1160 break;
1161 case 5: // 3D, pyramid
1162 subContVol[0].volume_ =
1163 hexahedronVolume(subContVol[0].global,
1164 edgeCoord[2],
1165 faceCoord[0],
1166 edgeCoord[0],
1167 edgeCoord[4],
1168 faceCoord[3],
1169 elementGlobal,
1170 faceCoord[1]);
1171 subContVol[1].volume_ =
1172 hexahedronVolume(subContVol[1].global,
1173 edgeCoord[1],
1174 faceCoord[0],
1175 edgeCoord[2],
1176 edgeCoord[5],
1177 faceCoord[2],
1178 elementGlobal,
1179 faceCoord[3]);
1180 subContVol[2].volume_ =
1181 hexahedronVolume(subContVol[2].global,
1182 edgeCoord[0],
1183 faceCoord[0],
1184 edgeCoord[3],
1185 edgeCoord[6],
1186 faceCoord[1],
1187 elementGlobal,
1188 faceCoord[4]);
1189 subContVol[3].volume_ =
1190 hexahedronVolume(subContVol[3].global,
1191 edgeCoord[3],
1192 faceCoord[0],
1193 edgeCoord[1],
1194 edgeCoord[7],
1195 faceCoord[4],
1196 elementGlobal,
1197 faceCoord[2]);
1198 subContVol[4].volume_ = elementVolume -
1199 subContVol[0].volume_ -
1200 subContVol[1].volume_ -
1201 subContVol[2].volume_ -
1202 subContVol[3].volume_;
1203 break;
1204 case 6: // 3D, prism
1205 subContVol[0].volume_ =
1206 hexahedronVolume(subContVol[0].global,
1207 edgeCoord[3],
1208 faceCoord[3],
1209 edgeCoord[4],
1210 edgeCoord[0],
1211 faceCoord[0],
1212 elementGlobal,
1213 faceCoord[1]);
1214 subContVol[1].volume_ =
1215 hexahedronVolume(subContVol[1].global,
1216 edgeCoord[5],
1217 faceCoord[3],
1218 edgeCoord[3],
1219 edgeCoord[1],
1220 faceCoord[2],
1221 elementGlobal,
1222 faceCoord[0]);
1223 subContVol[2].volume_ =
1224 hexahedronVolume(subContVol[2].global,
1225 edgeCoord[4],
1226 faceCoord[3],
1227 edgeCoord[5],
1228 edgeCoord[2],
1229 faceCoord[1],
1230 elementGlobal,
1231 faceCoord[2]);
1232 subContVol[3].volume_ =
1233 hexahedronVolume(edgeCoord[0],
1234 faceCoord[0],
1235 elementGlobal,
1236 faceCoord[1],
1237 subContVol[3].global,
1238 edgeCoord[6],
1239 faceCoord[4],
1240 edgeCoord[7]);
1241 subContVol[4].volume_ =
1242 hexahedronVolume(edgeCoord[1],
1243 faceCoord[2],
1244 elementGlobal,
1245 faceCoord[0],
1246 subContVol[4].global,
1247 edgeCoord[8],
1248 faceCoord[4],
1249 edgeCoord[6]);
1250 subContVol[5].volume_ =
1251 hexahedronVolume(edgeCoord[2],
1252 faceCoord[1],
1253 elementGlobal,
1254 faceCoord[2],
1255 subContVol[5].global,
1256 edgeCoord[7],
1257 faceCoord[4],
1258 edgeCoord[8]);
1259 break;
1260 case 8: // 3D, hexahedron
1261 subContVol[0].volume_ =
1262 hexahedronVolume(subContVol[0].global,
1263 edgeCoord[6],
1264 faceCoord[4],
1265 edgeCoord[4],
1266 edgeCoord[0],
1267 faceCoord[2],
1268 elementGlobal,
1269 faceCoord[0]);
1270 subContVol[1].volume_ =
1271 hexahedronVolume(subContVol[1].global,
1272 edgeCoord[5],
1273 faceCoord[4],
1274 edgeCoord[6],
1275 edgeCoord[1],
1276 faceCoord[1],
1277 elementGlobal,
1278 faceCoord[2]);
1279 subContVol[2].volume_ =
1280 hexahedronVolume(subContVol[2].global,
1281 edgeCoord[4],
1282 faceCoord[4],
1283 edgeCoord[7],
1284 edgeCoord[2],
1285 faceCoord[0],
1286 elementGlobal,
1287 faceCoord[3]);
1288 subContVol[3].volume_ =
1289 hexahedronVolume(subContVol[3].global,
1290 edgeCoord[7],
1291 faceCoord[4],
1292 edgeCoord[5],
1293 edgeCoord[3],
1294 faceCoord[3],
1295 elementGlobal,
1296 faceCoord[1]);
1297 subContVol[4].volume_ =
1298 hexahedronVolume(edgeCoord[0],
1299 faceCoord[2],
1300 elementGlobal,
1301 faceCoord[0],
1302 subContVol[4].global,
1303 edgeCoord[10],
1304 faceCoord[5],
1305 edgeCoord[8]);
1306 subContVol[5].volume_ =
1307 hexahedronVolume(edgeCoord[1],
1308 faceCoord[1],
1309 elementGlobal,
1310 faceCoord[2],
1311 subContVol[5].global,
1312 edgeCoord[9],
1313 faceCoord[5],
1314 edgeCoord[10]);
1315 subContVol[6].volume_ =
1316 hexahedronVolume(edgeCoord[2],
1317 faceCoord[0],
1318 elementGlobal,
1319 faceCoord[3],
1320 subContVol[6].global,
1321 edgeCoord[8],
1322 faceCoord[5],
1323 edgeCoord[11]);
1324 subContVol[7].volume_ =
1325 hexahedronVolume(edgeCoord[3],
1326 faceCoord[3],
1327 elementGlobal,
1328 faceCoord[1],
1329 subContVol[7].global,
1330 edgeCoord[11],
1331 faceCoord[5],
1332 edgeCoord[9]);
1333 break;
1334 default:
1335 throw std::logic_error("Not implemented:VcfvStencil for dim = "+std::to_string(dim)
1336 +", numVertices = "+std::to_string(numVertices));
1337 }
1338 }
1339 else
1340 throw std::logic_error("Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1341 }
1342#if __GNUC__ || __clang__
1343#pragma GCC diagnostic pop
1344#endif
1345
1346 const GridView& gridView_;
1347 const Mapper& vertexMapper_;
1348
1349 Element element_;
1350
1351#if HAVE_DUNE_LOCALFUNCTIONS
1352 static LocalFiniteElementCache feCache_;
1353#endif // HAVE_DUNE_LOCALFUNCTIONS
1354
1356 LocalPosition elementLocal;
1358 GlobalPosition elementGlobal;
1360 Scalar elementVolume;
1362 SubControlVolume subContVol[maxNC];
1364 SubControlVolumeFace subContVolFace[maxNE];
1366 BoundaryFace boundaryFace_[maxBF];
1367 unsigned numBoundarySegments_;
1369 GlobalPosition edgeCoord[maxNE];
1371 GlobalPosition faceCoord[maxNF];
1373 unsigned numVertices;
1375 unsigned numEdges;
1377 unsigned numFaces;
1378 Dune::GeometryType geometryType_;
1379};
1380
1381#if HAVE_DUNE_LOCALFUNCTIONS
1382template<class Scalar, class GridView>
1383typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
1384VcfvStencil<Scalar, GridView>::feCache_;
1385#endif // HAVE_DUNE_LOCALFUNCTIONS
1386
1387} // namespace Opm
1388
1389#endif
1390
Opm::QuadrialteralQuadratureGeometry
Quadrature geometry for quadrilaterals.
Definition: quadraturegeometries.hh:40
Opm::QuadrialteralQuadratureGeometry::corner
const GlobalPosition & corner(unsigned cornerIdx) const
Return the position of the corner with a given index.
Definition: quadraturegeometries.hh:127
Opm::QuadrialteralQuadratureGeometry::center
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition: quadraturegeometries.hh:69
Opm::VcfvStencil::ScvGeometry
Definition: vcfvstencil.hh:728
Opm::VcfvStencil::ScvGeometry::center
const GlobalPosition center() const
Definition: vcfvstencil.hh:730
Opm::VcfvStencil::ScvGeometry::localGeometry_
const ScvLocalGeometry * localGeometry_
Definition: vcfvstencil.hh:742
Opm::VcfvStencil::ScvGeometry::element_
const Element * element_
Definition: vcfvstencil.hh:743
Opm::VcfvStencil::ScvGeometry::corner
const GlobalPosition corner(unsigned cornerIdx) const
Definition: vcfvstencil.hh:733
Opm::VcfvStencil::ScvGeometry::global
const GlobalPosition global(const LocalPosition &localPos) const
Definition: vcfvstencil.hh:736
Opm::VcfvStencil::ScvGeometry::localGeometry
const ScvLocalGeometry & localGeometry() const
Definition: vcfvstencil.hh:739
Opm::VcfvStencil
Represents the finite volume geometry of a single element in the VCFV discretization.
Definition: vcfvstencil.hh:467
Opm::VcfvStencil::numPrimaryDof
unsigned numPrimaryDof() const
Definition: vcfvstencil.hh:1063
Opm::VcfvStencil::subControlVolume
const SubControlVolume & subControlVolume(unsigned scvIdx) const
Definition: vcfvstencil.hh:1069
Opm::VcfvStencil::VcfvStencil
VcfvStencil(const GridView &gridView, const Mapper &mapper)
Definition: vcfvstencil.hh:811
Opm::VcfvStencil::numBoundaryFaces
unsigned numBoundaryFaces() const
Definition: vcfvstencil.hh:1078
Opm::VcfvStencil::globalSpaceIndex
unsigned globalSpaceIndex(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1091
Opm::VcfvStencil::BoundaryFace
SubControlVolumeFace BoundaryFace
compatibility alias
Definition: vcfvstencil.hh:809
Opm::VcfvStencil::numDof
unsigned numDof() const
Definition: vcfvstencil.hh:1060
Opm::VcfvStencil::Mapper
Dune::MultipleCodimMultipleGeomTypeMapper< GridView > Mapper
exported Mapper type
Definition: vcfvstencil.hh:725
Opm::VcfvStencil::update
void update(const Element &e)
Definition: vcfvstencil.hh:864
Opm::VcfvStencil::numInteriorFaces
unsigned numInteriorFaces() const
Definition: vcfvstencil.hh:1075
Opm::VcfvStencil::entity
Entity entity(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1102
Opm::VcfvStencil::updatePrimaryTopology
void updatePrimaryTopology(const Element &element)
Definition: vcfvstencil.hh:856
Opm::VcfvStencil::Entity
typename GridView::Traits::template Codim< dim >::Entity Entity
Definition: vcfvstencil.hh:477
Opm::VcfvStencil::interiorFace
const SubControlVolumeFace & interiorFace(unsigned faceIdx) const
Definition: vcfvstencil.hh:1081
Opm::VcfvStencil::boundaryFace
const BoundaryFace & boundaryFace(unsigned bfIdx) const
Definition: vcfvstencil.hh:1084
Opm::VcfvStencil::partitionType
Dune::PartitionType partitionType(unsigned scvIdx) const
Definition: vcfvstencil.hh:1066
Opm::VcfvStencil::updateTopology
void updateTopology(const Element &e)
Update the non-geometric part of the stencil.
Definition: vcfvstencil.hh:836
Opm::VcfvStencil::updateScvGeometry
void updateScvGeometry(const Element &element)
Definition: vcfvstencil.hh:1016
Opm::Properties::dim
static const int dim
Definition: structuredgridvanguard.hh:67
Opm
Definition: blackoilboundaryratevector.hh:37
Opm::ElementType
ElementType
The types of reference elements available.
Definition: vcfvstencil.hh:52
Opm::simplex
@ simplex
Definition: vcfvstencil.hh:54
Opm::cube
@ cube
Definition: vcfvstencil.hh:55
Opm::none
@ none
Definition: vcfvstencil.hh:53
Opm::VcfvStencil::SubControlVolumeFace
interior face of a sub control volume
Definition: vcfvstencil.hh:777
Opm::VcfvStencil::SubControlVolumeFace::area
Scalar area() const
Definition: vcfvstencil.hh:787
Opm::VcfvStencil::SubControlVolumeFace::area_
Scalar area_
area of face
Definition: vcfvstencil.hh:805
Opm::VcfvStencil::SubControlVolumeFace::j
unsigned short j
Definition: vcfvstencil.hh:797
Opm::VcfvStencil::SubControlVolumeFace::exteriorIndex
unsigned short exteriorIndex() const
Definition: vcfvstencil.hh:784
Opm::VcfvStencil::SubControlVolumeFace::ipGlobal_
GlobalPosition ipGlobal_
integration point in global coords
Definition: vcfvstencil.hh:801
Opm::VcfvStencil::SubControlVolumeFace::integrationPos
const GlobalPosition & integrationPos() const
Definition: vcfvstencil.hh:793
Opm::VcfvStencil::SubControlVolumeFace::interiorIndex
unsigned short interiorIndex() const
Definition: vcfvstencil.hh:781
Opm::VcfvStencil::SubControlVolumeFace::normal_
DimVector normal_
normal on face pointing to CV j or outward of the domain with length equal to |scvf|
Definition: vcfvstencil.hh:803
Opm::VcfvStencil::SubControlVolumeFace::i
unsigned short i
scvf seperates corner i and j of elem
Definition: vcfvstencil.hh:797
Opm::VcfvStencil::SubControlVolumeFace::normal
const DimVector & normal() const
Definition: vcfvstencil.hh:778
Opm::VcfvStencil::SubControlVolumeFace::localPos
const LocalPosition & localPos() const
Definition: vcfvstencil.hh:790
Opm::VcfvStencil::SubControlVolumeFace::ipLocal_
LocalPosition ipLocal_
integration point in local coords
Definition: vcfvstencil.hh:799
Opm::VcfvStencil::SubControlVolume
finite volume intersected with element
Definition: vcfvstencil.hh:747
Opm::VcfvStencil::SubControlVolume::geometry_
ScvGeometry geometry_
The geometry of the sub-control volume in local coordinates.
Definition: vcfvstencil.hh:770
Opm::VcfvStencil::SubControlVolume::geometry
const ScvGeometry & geometry() const
Definition: vcfvstencil.hh:760
Opm::VcfvStencil::SubControlVolume::local
LocalPosition local
local vertex position
Definition: vcfvstencil.hh:764
Opm::VcfvStencil::SubControlVolume::global
GlobalPosition global
global vertex position
Definition: vcfvstencil.hh:766
Opm::VcfvStencil::SubControlVolume::volume
Scalar volume() const
Definition: vcfvstencil.hh:754
Opm::VcfvStencil::SubControlVolume::center
const GlobalPosition center() const
Definition: vcfvstencil.hh:751
Opm::VcfvStencil::SubControlVolume::volume_
Scalar volume_
space occupied by the sub-control volume
Definition: vcfvstencil.hh:768
Opm::VcfvStencil::SubControlVolume::globalPos
const GlobalPosition & globalPos() const
Definition: vcfvstencil.hh:748
Opm::VcfvStencil::SubControlVolume::localGeometry
const ScvLocalGeometry & localGeometry() const
Definition: vcfvstencil.hh:757
Opm::VcfvStencil::SubControlVolume::gradCenter
Dune::FieldVector< DimVector, maxNC > gradCenter
derivative of shape function at the center of the sub control volume
Definition: vcfvstencil.hh:773