All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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 <opm/common/Unused.hpp>
34 
35 #include <opm/common/Exceptions.hpp>
36 #include <opm/common/ErrorMacros.hpp>
37 
38 #include <dune/grid/common/intersectioniterator.hh>
39 #include <dune/grid/common/mcmgmapper.hh>
40 #include <dune/geometry/referenceelements.hh>
41 
42 #if HAVE_DUNE_LOCALFUNCTIONS
43 #include <dune/localfunctions/lagrange/pqkfactory.hh>
44 #endif // HAVE_DUNE_LOCALFUNCTIONS
45 
46 #include <dune/common/version.hh>
47 
48 #include <vector>
49 
50 namespace Ewoms {
51 
55 template <class Scalar, unsigned dim, unsigned basicGeomType>
56 class VcfvScvGeometries;
57 
59 // local geometries for 1D elements
61 template <class Scalar>
62 class VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>
63 {
64  enum { dim = 1 };
65  enum { numScv = 2 };
66  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
67 
68 public:
70 
71  static void init()
72  {
73  // 1D LINE SEGMENTS
74  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] = {
75  { // corners of the first sub control volume
76  { 0.0 },
77  { 0.5 }
78  },
79 
80  { // corners of the second sub control volume
81  { 0.5 },
82  { 1.0 }
83  }
84  };
85  for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
86  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
87  }
88 
89  static const ScvLocalGeometry& get(unsigned scvIdx)
90  { return scvGeoms_[scvIdx]; }
91 
92 private:
93  static ScvLocalGeometry scvGeoms_[numScv];
94 };
95 
96 template <class Scalar>
97 typename VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::ScvLocalGeometry
98 VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::scvGeoms_[
99  VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::numScv];
100 
101 template <class Scalar>
102 class VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::simplex>
103 {
104  enum { dim = 1 };
105  enum { numScv = 2 };
106  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
107 
108 public:
110 
111  static const ScvLocalGeometry& get(unsigned scvIdx OPM_UNUSED)
112  {
113  OPM_THROW(std::logic_error,
114  "Not implemented: VcfvScvGeometries<Scalar, 1, Dune::GeometryType::simplex>");
115  }
116 };
117 
119 // local geometries for 2D elements
121 template <class Scalar>
122 class VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>
123 {
124  enum { dim = 2 };
125  enum { numScv = 3 };
126  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
127 
128 public:
130 
131  static const ScvLocalGeometry& get(unsigned scvIdx)
132  { return scvGeoms_[scvIdx]; }
133 
134  static void init()
135  {
136  // 2D SIMPLEX
137  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
138  {
139  { // SCV 0 corners
140  { 0.0, 0.0 },
141  { 1.0/2, 0.0 },
142  { 1.0/3, 1.0/3 },
143  { 0.0, 1.0/2 },
144  },
145 
146  { // SCV 1 corners
147  { 1.0/2, 0.0 },
148  { 1.0, 0.0 },
149  { 1.0/3, 1.0/3 },
150  { 1.0/2, 1.0/2 },
151  },
152 
153  { // SCV 2 corners
154  { 0.0, 1.0/2 },
155  { 1.0/3, 1.0/3 },
156  { 0.0, 1.0 },
157  { 1.0/2, 1.0/2 },
158  }
159  };
160 
161  for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
162  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
163  }
164 
165 private:
166  static ScvLocalGeometry scvGeoms_[numScv];
167 };
168 
169 template <class Scalar>
170 typename VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::ScvLocalGeometry
171 VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::scvGeoms_[
172  VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::numScv];
173 
174 template <class Scalar>
175 class VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>
176 {
177  enum { dim = 2 };
178  enum { numScv = 4 };
179  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::cube;
180 
181 public:
183 
184  static const ScvLocalGeometry& get(unsigned scvIdx)
185  { return scvGeoms_[scvIdx]; }
186 
187  static void init()
188  {
189  // 2D CUBE
190  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
191  {
192  { // SCV 0 corners
193  { 0.0, 0.0 },
194  { 0.5, 0.0 },
195  { 0.0, 0.5 },
196  { 0.5, 0.5 }
197  },
198 
199  { // SCV 1 corners
200  { 0.5, 0.0 },
201  { 1.0, 0.0 },
202  { 0.5, 0.5 },
203  { 1.0, 0.5 }
204  },
205 
206  { // SCV 2 corners
207  { 0.0, 0.5 },
208  { 0.5, 0.5 },
209  { 0.0, 1.0 },
210  { 0.5, 1.0 }
211  },
212 
213  { // SCV 3 corners
214  { 0.5, 0.5 },
215  { 1.0, 0.5 },
216  { 0.5, 1.0 },
217  { 1.0, 1.0 }
218  }
219  };
220 
221  for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
222  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
223  }
224 
225  static ScvLocalGeometry scvGeoms_[numScv];
226 };
227 
228 template <class Scalar>
229 typename VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::ScvLocalGeometry
230 VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::scvGeoms_[
231  VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::numScv];
232 
234 // local geometries for 3D elements
236 template <class Scalar>
237 class VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::simplex>
238 {
239  enum { dim = 3 };
240  enum { numScv = 4 };
241  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
242 
243 public:
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 
307 private:
308  static ScvLocalGeometry scvGeoms_[numScv];
309 };
310 
311 template <class Scalar>
312 typename VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::simplex>::ScvLocalGeometry
313 VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::simplex>::scvGeoms_[
314  VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::simplex>::numScv];
315 
316 template <class Scalar>
317 class VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>
318 {
319  enum { dim = 3 };
320  enum { numScv = 8 };
321  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::cube;
322 
323 public:
325 
326  static const ScvLocalGeometry& get(unsigned scvIdx)
327  { return scvGeoms_[scvIdx]; }
328 
329  static void init()
330  {
331  // 3D CUBE
332  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
333  {
334  { // SCV 0 corners
335  { 0.0, 0.0, 0.0 },
336  { 1.0/2, 0.0, 0.0 },
337  { 0.0, 1.0/2, 0.0 },
338  { 1.0/2, 1.0/2, 0.0 },
339 
340  { 0.0, 0.0, 1.0/2 },
341  { 1.0/2, 0.0, 1.0/2 },
342  { 0.0, 1.0/2, 1.0/2 },
343  { 1.0/2, 1.0/2, 1.0/2 },
344  },
345 
346  { // SCV 1 corners
347  { 1.0/2, 0.0, 0.0 },
348  { 1.0, 0.0, 0.0 },
349  { 1.0/2, 1.0/2, 0.0 },
350  { 1.0, 1.0/2, 0.0 },
351 
352  { 1.0/2, 0.0, 1.0/2 },
353  { 1.0, 0.0, 1.0/2 },
354  { 1.0/2, 1.0/2, 1.0/2 },
355  { 1.0, 1.0/2, 1.0/2 },
356  },
357 
358  { // SCV 2 corners
359  { 0.0, 1.0/2, 0.0 },
360  { 1.0/2, 1.0/2, 0.0 },
361  { 0.0, 1.0, 0.0 },
362  { 1.0/2, 1.0, 0.0 },
363 
364  { 0.0, 1.0/2, 1.0/2 },
365  { 1.0/2, 1.0/2, 1.0/2 },
366  { 0.0, 1.0, 1.0/2 },
367  { 1.0/2, 1.0, 1.0/2 },
368  },
369 
370  { // SCV 3 corners
371  { 1.0/2, 1.0/2, 0.0 },
372  { 1.0, 1.0/2, 0.0 },
373  { 1.0/2, 1.0, 0.0 },
374  { 1.0, 1.0, 0.0 },
375 
376  { 1.0/2, 1.0/2, 1.0/2 },
377  { 1.0, 1.0/2, 1.0/2 },
378  { 1.0/2, 1.0, 1.0/2 },
379  { 1.0, 1.0, 1.0/2 },
380  },
381 
382  { // SCV 4 corners
383  { 0.0, 0.0, 1.0/2 },
384  { 1.0/2, 0.0, 1.0/2 },
385  { 0.0, 1.0/2, 1.0/2 },
386  { 1.0/2, 1.0/2, 1.0/2 },
387 
388  { 0.0, 0.0, 1.0 },
389  { 1.0/2, 0.0, 1.0 },
390  { 0.0, 1.0/2, 1.0 },
391  { 1.0/2, 1.0/2, 1.0 },
392  },
393 
394  { // SCV 5 corners
395  { 1.0/2, 0.0, 1.0/2 },
396  { 1.0, 0.0, 1.0/2 },
397  { 1.0/2, 1.0/2, 1.0/2 },
398  { 1.0, 1.0/2, 1.0/2 },
399 
400  { 1.0/2, 0.0, 1.0 },
401  { 1.0, 0.0, 1.0 },
402  { 1.0/2, 1.0/2, 1.0 },
403  { 1.0, 1.0/2, 1.0 },
404  },
405 
406  { // SCV 6 corners
407  { 0.0, 1.0/2, 1.0/2 },
408  { 1.0/2, 1.0/2, 1.0/2 },
409  { 0.0, 1.0, 1.0/2 },
410  { 1.0/2, 1.0, 1.0/2 },
411 
412  { 0.0, 1.0/2, 1.0 },
413  { 1.0/2, 1.0/2, 1.0 },
414  { 0.0, 1.0, 1.0 },
415  { 1.0/2, 1.0, 1.0 },
416  },
417 
418  { // SCV 7 corners
419  { 1.0/2, 1.0/2, 1.0/2 },
420  { 1.0, 1.0/2, 1.0/2 },
421  { 1.0/2, 1.0, 1.0/2 },
422  { 1.0, 1.0, 1.0/2 },
423 
424  { 1.0/2, 1.0/2, 1.0 },
425  { 1.0, 1.0/2, 1.0 },
426  { 1.0/2, 1.0, 1.0 },
427  { 1.0, 1.0, 1.0 },
428  },
429  };
430 
431  for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
432  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
433  }
434 private:
435  static ScvLocalGeometry scvGeoms_[numScv];
436 };
437 
438 template <class Scalar>
439 typename VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::ScvLocalGeometry
440 VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::scvGeoms_[
441  VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::numScv];
442 
466 template <class Scalar, class GridView>
468 {
469  enum{dim = GridView::dimension};
470  enum{dimWorld = GridView::dimensionworld};
471  enum{maxNC = (dim < 3 ? 4 : 8)};
472  enum{maxNE = (dim < 3 ? 4 : 12)};
473  enum{maxNF = (dim < 3 ? 1 : 6)};
474  enum{maxBF = (dim < 3 ? 8 : 24)};
475  typedef typename GridView::ctype CoordScalar;
476  typedef typename GridView::Traits::template Codim<0>::Entity Element;
477 public:
478  typedef typename GridView::Traits::template Codim<dim>::Entity Entity;
479 private:
480  typedef typename Element::Geometry Geometry;
481  typedef Dune::FieldVector<Scalar,dimWorld> DimVector;
482  typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
483  typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
484  typedef typename GridView::IntersectionIterator IntersectionIterator;
485 
487 
488 #if HAVE_DUNE_LOCALFUNCTIONS
489  typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
490  typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
491  typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
492  typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
493 #endif // HAVE_DUNE_LOCALFUNCTIONS
494 
495  Scalar quadrilateralArea(const GlobalPosition& p0,
496  const GlobalPosition& p1,
497  const GlobalPosition& p2,
498  const GlobalPosition& p3)
499  {
500  return 0.5*std::abs((p3[0] - p1[0])*(p2[1] - p0[1]) - (p3[1] - p1[1])*(p2[0] - p0[0]));
501  }
502 
503  void crossProduct(DimVector& c, const DimVector& a, const DimVector& b)
504  {
505  c[0] = a[1]*b[2] - a[2]*b[1];
506  c[1] = a[2]*b[0] - a[0]*b[2];
507  c[2] = a[0]*b[1] - a[1]*b[0];
508  }
509 
510  Scalar pyramidVolume(const GlobalPosition& p0,
511  const GlobalPosition& p1,
512  const GlobalPosition& p2,
513  const GlobalPosition& p3,
514  const GlobalPosition& p4)
515  {
516  DimVector a(p2); a -= p0;
517  DimVector b(p3); b -= p1;
518 
519  DimVector n;
520  crossProduct(n, a, b);
521 
522  a = p4; a -= p0;
523 
524  return 1.0/6.0*(n*a);
525  }
526 
527  Scalar prismVolume(const GlobalPosition& p0,
528  const GlobalPosition& p1,
529  const GlobalPosition& p2,
530  const GlobalPosition& p3,
531  const GlobalPosition& p4,
532  const GlobalPosition& p5)
533  {
534  DimVector a(p4);
535  for (unsigned k = 0; k < dimWorld; ++k)
536  a[k] -= p0[k];
537  DimVector b(p1);
538  for (unsigned k = 0; k < dimWorld; ++k)
539  b[k] -= p3[k];
540  DimVector m;
541  crossProduct(m, a, b);
542 
543  for (unsigned k = 0; k < dimWorld; ++k)
544  a[k] = p1[k] - p0[k];
545  for (unsigned k = 0; k < dimWorld; ++k)
546  b[k] = p2[k] - p0[k];
547  DimVector n;
548  crossProduct(n, a, b);
549  n += m;
550 
551  for (unsigned k = 0; k < dimWorld; ++k)
552  a[k] = p5[k] - p0[k];
553 
554  return std::abs(1.0/6.0*(n*a));
555  }
556 
557  Scalar hexahedronVolume(const GlobalPosition& p0,
558  const GlobalPosition& p1,
559  const GlobalPosition& p2,
560  const GlobalPosition& p3,
561  const GlobalPosition& p4,
562  const GlobalPosition& p5,
563  const GlobalPosition& p6,
564  const GlobalPosition& p7)
565  {
566  return
567  prismVolume(p0,p1,p2,p4,p5,p6)
568  + prismVolume(p0,p2,p3,p4,p6,p7);
569  }
570 
571  void normalOfQuadrilateral3D(DimVector& normal,
572  const GlobalPosition& p0,
573  const GlobalPosition& p1,
574  const GlobalPosition& p2,
575  const GlobalPosition& p3)
576  {
577  DimVector a(p2);
578  for (unsigned k = 0; k < dimWorld; ++k)
579  a[k] -= p0[k];
580  DimVector b(p3);
581  for (unsigned k = 0; k < dimWorld; ++k)
582  b[k] -= p1[k];
583 
584  crossProduct(normal, a, b);
585  normal *= 0.5;
586  }
587 
588  Scalar quadrilateralArea3D(const GlobalPosition& p0,
589  const GlobalPosition& p1,
590  const GlobalPosition& p2,
591  const GlobalPosition& p3)
592  {
593  DimVector normal;
594  normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
595  return normal.two_norm();
596  }
597 
598  void getFaceIndices(unsigned numElemVertices, unsigned k, unsigned& leftFace, unsigned& rightFace)
599  {
600  static const unsigned edgeToFaceTet[2][6] = {
601  {1, 0, 3, 2, 1, 3},
602  {0, 2, 0, 1, 3, 2}
603  };
604  static const unsigned edgeToFacePyramid[2][8] = {
605  {1, 2, 3, 4, 1, 3, 4, 2},
606  {0, 0, 0, 0, 3, 2, 1, 4}
607  };
608  static const unsigned edgeToFacePrism[2][9] = {
609  {1, 0, 2, 0, 1, 2, 4, 4, 4},
610  {0, 2, 1, 3, 3, 3, 0, 1, 2}
611  };
612  static const unsigned edgeToFaceHex[2][12] = {
613  {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
614  {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
615  };
616 
617  switch (numElemVertices) {
618  case 4:
619  leftFace = edgeToFaceTet[0][k];
620  rightFace = edgeToFaceTet[1][k];
621  break;
622  case 5:
623  leftFace = edgeToFacePyramid[0][k];
624  rightFace = edgeToFacePyramid[1][k];
625  break;
626  case 6:
627  leftFace = edgeToFacePrism[0][k];
628  rightFace = edgeToFacePrism[1][k];
629  break;
630  case 8:
631  leftFace = edgeToFaceHex[0][k];
632  rightFace = edgeToFaceHex[1][k];
633  break;
634  default:
635  OPM_THROW(std::logic_error,
636  "Not implemented: VcfvStencil::getFaceIndices for "
637  << numElemVertices << " vertices");
638  break;
639  }
640  }
641 
642  void getEdgeIndices(unsigned numElemVertices, unsigned face, unsigned vert, unsigned& leftEdge, unsigned& rightEdge)
643  {
644  static const int faceAndVertexToLeftEdgeTet[4][4] = {
645  { 0, 0, 2, -1},
646  { 0, 0, -1, 3},
647  { 1, -1, 1, 3},
648  {-1, 2, 2, 4}
649  };
650  static const int faceAndVertexToRightEdgeTet[4][4] = {
651  { 1, 2, 1, -1},
652  { 3, 4, -1, 4},
653  { 3, -1, 5, 5},
654  {-1, 4, 5, 5}
655  };
656  static const int faceAndVertexToLeftEdgePyramid[5][5] = {
657  { 0, 2, 3, 1, -1},
658  { 0, -1, 0, -1, 4},
659  {-1, 1, -1, 1, 5},
660  { 2, 2, -1, -1, 4},
661  {-1, -1, 3, 3, 7}
662  };
663  static const int faceAndVertexToRightEdgePyramid[5][5] = {
664  { 2, 1, 0, 3, -1},
665  { 4, -1, 6, -1, 6},
666  {-1, 5, -1, 7, 7},
667  { 4, 5, -1, -1, 5},
668  {-1, -1, 6, 7, 6}
669  };
670  static const int faceAndVertexToLeftEdgePrism[5][6] = {
671  { 3, 3, -1, 0, 1, -1},
672  { 4, -1, 4, 0, -1, 2},
673  {-1, 5, 5, -1, 1, 2},
674  { 3, 3, 5, -1, -1, -1},
675  {-1, -1, -1, 6, 6, 8}
676  };
677  static const int faceAndVertexToRightEdgePrism[5][6] = {
678  { 0, 1, -1, 6, 6, -1},
679  { 0, -1, 2, 7, -1, 7},
680  {-1, 1, 2, -1, 8, 8},
681  { 4, 5, 4, -1, -1, -1},
682  {-1, -1, -1, 7, 8, 7}
683  };
684  static const int faceAndVertexToLeftEdgeHex[6][8] = {
685  { 0, -1, 4, -1, 8, -1, 2, -1},
686  {-1, 5, -1, 3, -1, 1, -1, 9},
687  { 6, 1, -1, -1, 0, 10, -1, -1},
688  {-1, -1, 2, 7, -1, -1, 11, 3},
689  { 4, 6, 7, 5, -1, -1, -1, -1},
690  {-1, -1, -1, -1, 10, 9, 8, 11}
691  };
692  static const int faceAndVertexToRightEdgeHex[6][8] = {
693  { 4, -1, 2, -1, 0, -1, 8, -1},
694  {-1, 1, -1, 5, -1, 9, -1, 3},
695  { 0, 6, -1, -1, 10, 1, -1, -1},
696  {-1, -1, 7, 3, -1, -1, 2, 11},
697  { 6, 5, 4, 7, -1, -1, -1, -1},
698  {-1, -1, -1, -1, 8, 10, 11, 9}
699  };
700 
701  switch (numElemVertices) {
702  case 4:
703  leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgeTet[face][vert]);
704  rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgeTet[face][vert]);
705  break;
706  case 5:
707  leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgePyramid[face][vert]);
708  rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgePyramid[face][vert]);
709  break;
710  case 6:
711  leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgePrism[face][vert]);
712  rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgePrism[face][vert]);
713  break;
714  case 8:
715  leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgeHex[face][vert]);
716  rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgeHex[face][vert]);
717  break;
718  default:
719  OPM_THROW(std::logic_error,
720  "Not implemented: VcfvStencil::getFaceIndices for "
721  << numElemVertices << " vertices");
722  break;
723  }
724  }
725 
726 public:
727  typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,
728  Dune::MCMGVertexLayout > VertexMapper;
730  typedef VertexMapper Mapper;
731 
733  {
734  public:
735  const GlobalPosition center() const
736  { return global(localGeometry_->center()); }
737 
738  const GlobalPosition corner(unsigned cornerIdx) const
739  { return global(localGeometry_->corner(cornerIdx)); }
740 
741  const GlobalPosition global(const LocalPosition& localPos) const
742  { return element_->geometry().global(localPos); }
743 
744  const ScvLocalGeometry& localGeometry() const
745  { return *localGeometry_; }
746 
747  const ScvLocalGeometry *localGeometry_;
748  const Element *element_;
749  };
750 
752  {
753  const GlobalPosition& globalPos() const
754  { return global; }
755 
756  const GlobalPosition center() const
757  { return geometry_.center(); }
758 
759  Scalar volume() const
760  { return volume_; }
761 
762  const ScvLocalGeometry& localGeometry() const
763  { return geometry_.localGeometry(); }
764 
765  const ScvGeometry& geometry() const
766  { return geometry_; }
767 
769  LocalPosition local;
771  GlobalPosition global;
773  Scalar volume_;
776 
778  Dune::FieldVector<DimVector, maxNC> gradCenter;
779  };
780 
782  {
783  const DimVector& normal() const
784  { return normal_; }
785 
786  unsigned short interiorIndex() const
787  { return i; }
788 
789  unsigned short exteriorIndex() const
790  { return j; }
791 
792  Scalar area() const
793  { return area_; }
794 
795  const LocalPosition& localPos() const
796  { return ipLocal_; }
797 
798  const GlobalPosition& integrationPos() const
799  { return ipGlobal_; }
800 
802  unsigned short i,j;
804  LocalPosition ipLocal_;
806  GlobalPosition ipGlobal_;
808  DimVector normal_;
810  Scalar area_;
811  };
812 
815 
816  VcfvStencil(const GridView& gridView, const VertexMapper& vertexMapper)
817  : gridView_(gridView)
818  , vertexMapper_( vertexMapper )
819  , element_(*gridView.template begin</*codim=*/0>())
820  {
821  static bool localGeometriesInitialized = false;
822  if (!localGeometriesInitialized) {
823  localGeometriesInitialized = true;
824 
825  VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::init();
826  VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::init();
827  VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::init();
828  VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::init();
829  VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::simplex>::init();
830  }
831  }
832 
838  void updateTopology(const Element& e)
839  {
840  element_ = e;
841 
842  numVertices = e.subEntities(/*codim=*/dim);
843  numEdges = e.subEntities(/*codim=*/dim-1);
844  numFaces = (dim<3)?0:e.subEntities(/*codim=*/1);
845 
846  numBoundarySegments_ = 0; // TODO: really required here(?)
847 
848  // compute the local and global coordinates of the element
849  const Geometry& geometry = e.geometry();
850  geometryType_ = geometry.type();
851  const typename Dune::ReferenceElementContainer<CoordScalar,dim>::value_type&
852  referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
853  for (unsigned vertexIdx = 0; vertexIdx < numVertices; vertexIdx++) {
854  subContVol[vertexIdx].local = referenceElement.position(static_cast<int>(vertexIdx), dim);
855  subContVol[vertexIdx].global = geometry.corner(static_cast<int>(vertexIdx));
856  }
857  }
858 
859  void updatePrimaryTopology(const Element& element)
860  {
861  // since all degrees of freedom in a stencil are "primary" DOFs for the
862  // vertex-centered finite volume method, there's no difference to
863  // updateTopology()
864  updateTopology(element);
865  }
866 
867  void update(const Element& e)
868  {
869  updateTopology(e);
870 
871  const Geometry& geometry = e.geometry();
872  geometryType_ = geometry.type();
873 
874  const typename Dune::ReferenceElementContainer<CoordScalar,dim>::value_type&
875  referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
876 
877  elementVolume = geometry.volume();
878  elementLocal = referenceElement.position(0,0);
879  elementGlobal = geometry.global(elementLocal);
880 
881  // corners:
882  for (unsigned vert = 0; vert < numVertices; vert++) {
883  subContVol[vert].local = referenceElement.position(static_cast<int>(vert), dim);
884  subContVol[vert].global = geometry.global(subContVol[vert].local);
885  }
886 
887  // edges:
888  for (unsigned edge = 0; edge < numEdges; edge++) {
889  edgeCoord[edge] = geometry.global(referenceElement.position(static_cast<int>(edge), dim-1));
890  }
891 
892  // faces:
893  for (unsigned face = 0; face < numFaces; face++) {
894  faceCoord[face] = geometry.global(referenceElement.position(static_cast<int>(face), 1));
895  }
896 
897  // fill sub control volume data use specialization for this
898  // \todo maybe it would be a good idea to migrate everything
899  // which is dependend of the grid's dimension to
900  // _VcfvFVElemGeomHelper in order to benefit from more aggressive
901  // compiler optimizations...
902  fillSubContVolData_();
903 
904  // fill sub control volume face data:
905  for (unsigned k = 0; k < numEdges; k++) { // begin loop over edges / sub control volume faces
906  unsigned short i = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 0, dim));
907  unsigned short j = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 1, dim));
908  if (numEdges == 4 && (i == 2 || j == 2))
909  std::swap(i, j);
910  subContVolFace[k].i = i;
911  subContVolFace[k].j = j;
912 
913  // calculate the local integration point and
914  // the face normal. note that since dim is a
915  // constant which is known at compile time
916  // the compiler can optimize away all if
917  // cases which don't apply.
918  LocalPosition ipLocal_;
919  DimVector diffVec;
920  if (dim==1) {
921  subContVolFace[k].ipLocal_ = 0.5;
922  subContVolFace[k].normal_ = 1.0;
923  subContVolFace[k].area_ = 1.0;
924  ipLocal_ = subContVolFace[k].ipLocal_;
925  }
926  else if (dim==2) {
927  ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal;
928  ipLocal_ *= 0.5;
929  subContVolFace[k].ipLocal_ = ipLocal_;
930  for (unsigned m = 0; m < dimWorld; ++m)
931  diffVec[m] = elementGlobal[m] - edgeCoord[k][m];
932  subContVolFace[k].normal_[0] = diffVec[1];
933  subContVolFace[k].normal_[1] = -diffVec[0];
934 
935  for (unsigned m = 0; m < dimWorld; ++m)
936  diffVec[m] = subContVol[j].global[m] - subContVol[i].global[m];
937  // make sure the normal points to the right direction
938  if (subContVolFace[k].normal_ * diffVec < 0)
939  subContVolFace[k].normal_ *= -1;
940 
941  subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
942  subContVolFace[k].normal_ /= subContVolFace[k].area_;
943  }
944  else if (dim==3) {
945  unsigned leftFace;
946  unsigned rightFace;
947  getFaceIndices(numVertices, k, leftFace, rightFace);
948  ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal
949  + referenceElement.position(static_cast<int>(leftFace), 1)
950  + referenceElement.position(static_cast<int>(rightFace), 1);
951  ipLocal_ *= 0.25;
952  subContVolFace[k].ipLocal_ = ipLocal_;
953  normalOfQuadrilateral3D(subContVolFace[k].normal_,
954  edgeCoord[k], faceCoord[rightFace],
955  elementGlobal, faceCoord[leftFace]);
956  subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
957  subContVolFace[k].normal_ /= subContVolFace[k].area_;
958  }
959 
960  // get the global integration point and the Jacobian inverse
961  subContVolFace[k].ipGlobal_ = geometry.global(ipLocal_);
962  } // end loop over edges / sub control volume faces
963 
964  // fill boundary face data:
965  IntersectionIterator endit = gridView_.iend(e);
966  for (IntersectionIterator it = gridView_.ibegin(e); it != endit; ++it) {
967  const typename IntersectionIterator::Intersection& intersection = *it ;
968 
969  if ( ! intersection.boundary())
970  continue;
971 
972  unsigned face = static_cast<unsigned>(intersection.indexInInside());
973  unsigned numVerticesOfFace = static_cast<unsigned>(referenceElement.size(static_cast<int>(face), 1, dim));
974  for (unsigned vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
975  {
976  unsigned short vertInElement = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(face), 1, static_cast<int>(vertInFace), dim));
977  unsigned bfIdx = numBoundarySegments_;
978  ++numBoundarySegments_;
979 
980  if (dim == 1) {
981  boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim);
982  boundaryFace_[bfIdx].area_ = 1.0;
983  }
984  else if (dim == 2) {
985  boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
986  + referenceElement.position(static_cast<int>(face), 1);
987  boundaryFace_[bfIdx].ipLocal_ *= 0.5;
988  boundaryFace_[bfIdx].area_ = 0.5 * intersection.geometry().volume();
989  }
990  else if (dim == 3) {
991  unsigned leftEdge;
992  unsigned rightEdge;
993  getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
994  boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
995  + referenceElement.position(static_cast<int>(face), 1)
996  + referenceElement.position(static_cast<int>(leftEdge), dim-1)
997  + referenceElement.position(static_cast<int>(rightEdge), dim-1);
998  boundaryFace_[bfIdx].ipLocal_ *= 0.25;
999  boundaryFace_[bfIdx].area_ =
1000  quadrilateralArea3D(subContVol[vertInElement].global,
1001  edgeCoord[rightEdge],
1002  faceCoord[face],
1003  edgeCoord[leftEdge]);
1004  }
1005  else
1006  OPM_THROW(std::logic_error, "Not implemented:VcfvStencil for dim = " << dim);
1007 
1008  boundaryFace_[bfIdx].ipGlobal_ = geometry.global(boundaryFace_[bfIdx].ipLocal_);
1009  boundaryFace_[bfIdx].i = vertInElement;
1010  boundaryFace_[bfIdx].j = vertInElement;
1011 
1012  // ASSUME constant normal on the segment of the boundary face
1013  boundaryFace_[bfIdx].normal_ = intersection.centerUnitOuterNormal();
1014  }
1015  }
1016 
1017  updateScvGeometry(e);
1018  }
1019 
1020  void updateScvGeometry(const Element& element)
1021  {
1022  auto geomType = element.geometry().type();
1023 
1024  // get the local geometries of the sub control volumes
1025  if (geomType.isTriangle() || geomType.isTetrahedron()) {
1026  for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1027  subContVol[vertIdx].geometry_.element_ = &element;
1028  subContVol[vertIdx].geometry_.localGeometry_ =
1029  &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::simplex>::get(vertIdx);
1030  }
1031  }
1032  else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1033  for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1034  subContVol[vertIdx].geometry_.element_ = &element;
1035  subContVol[vertIdx].geometry_.localGeometry_ =
1036  &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::cube>::get(vertIdx);
1037  }
1038  }
1039  else
1040  OPM_THROW(std::logic_error,
1041  "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(0 <= scvIdx && 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(0 <= dofIdx && 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(0 <= dofIdx && dofIdx < numDof());
1110  return element_.template subEntity<dim>(static_cast<int>(dofIdx));
1111  }
1112 
1113 private:
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 (dim == 1) {
1122  // 1D
1123  subContVol[0].volume_ = 0.5*elementVolume;
1124  subContVol[1].volume_ = 0.5*elementVolume;
1125  }
1126  else if (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  OPM_THROW(std::logic_error,
1157  "Not implemented:VcfvStencil dim = " << dim
1158  << ", numVertices = " << numVertices);
1159  }
1160  }
1161  else if (dim == 3) {
1162  switch (numVertices) {
1163  case 4: // 3D, tetrahedron
1164  for (unsigned k = 0; k < numVertices; k++)
1165  subContVol[k].volume_ = elementVolume / 4.0;
1166  break;
1167  case 5: // 3D, pyramid
1168  subContVol[0].volume_ =
1169  hexahedronVolume(subContVol[0].global,
1170  edgeCoord[2],
1171  faceCoord[0],
1172  edgeCoord[0],
1173  edgeCoord[4],
1174  faceCoord[3],
1175  elementGlobal,
1176  faceCoord[1]);
1177  subContVol[1].volume_ =
1178  hexahedronVolume(subContVol[1].global,
1179  edgeCoord[1],
1180  faceCoord[0],
1181  edgeCoord[2],
1182  edgeCoord[5],
1183  faceCoord[2],
1184  elementGlobal,
1185  faceCoord[3]);
1186  subContVol[2].volume_ =
1187  hexahedronVolume(subContVol[2].global,
1188  edgeCoord[0],
1189  faceCoord[0],
1190  edgeCoord[3],
1191  edgeCoord[6],
1192  faceCoord[1],
1193  elementGlobal,
1194  faceCoord[4]);
1195  subContVol[3].volume_ =
1196  hexahedronVolume(subContVol[3].global,
1197  edgeCoord[3],
1198  faceCoord[0],
1199  edgeCoord[1],
1200  edgeCoord[7],
1201  faceCoord[4],
1202  elementGlobal,
1203  faceCoord[2]);
1204  subContVol[4].volume_ = elementVolume -
1205  subContVol[0].volume_ -
1206  subContVol[1].volume_ -
1207  subContVol[2].volume_ -
1208  subContVol[3].volume_;
1209  break;
1210  case 6: // 3D, prism
1211  subContVol[0].volume_ =
1212  hexahedronVolume(subContVol[0].global,
1213  edgeCoord[3],
1214  faceCoord[3],
1215  edgeCoord[4],
1216  edgeCoord[0],
1217  faceCoord[0],
1218  elementGlobal,
1219  faceCoord[1]);
1220  subContVol[1].volume_ =
1221  hexahedronVolume(subContVol[1].global,
1222  edgeCoord[5],
1223  faceCoord[3],
1224  edgeCoord[3],
1225  edgeCoord[1],
1226  faceCoord[2],
1227  elementGlobal,
1228  faceCoord[0]);
1229  subContVol[2].volume_ =
1230  hexahedronVolume(subContVol[2].global,
1231  edgeCoord[4],
1232  faceCoord[3],
1233  edgeCoord[5],
1234  edgeCoord[2],
1235  faceCoord[1],
1236  elementGlobal,
1237  faceCoord[2]);
1238  subContVol[3].volume_ =
1239  hexahedronVolume(edgeCoord[0],
1240  faceCoord[0],
1241  elementGlobal,
1242  faceCoord[1],
1243  subContVol[3].global,
1244  edgeCoord[6],
1245  faceCoord[4],
1246  edgeCoord[7]);
1247  subContVol[4].volume_ =
1248  hexahedronVolume(edgeCoord[1],
1249  faceCoord[2],
1250  elementGlobal,
1251  faceCoord[0],
1252  subContVol[4].global,
1253  edgeCoord[8],
1254  faceCoord[4],
1255  edgeCoord[6]);
1256  subContVol[5].volume_ =
1257  hexahedronVolume(edgeCoord[2],
1258  faceCoord[1],
1259  elementGlobal,
1260  faceCoord[2],
1261  subContVol[5].global,
1262  edgeCoord[7],
1263  faceCoord[4],
1264  edgeCoord[8]);
1265  break;
1266  case 8: // 3D, hexahedron
1267  subContVol[0].volume_ =
1268  hexahedronVolume(subContVol[0].global,
1269  edgeCoord[6],
1270  faceCoord[4],
1271  edgeCoord[4],
1272  edgeCoord[0],
1273  faceCoord[2],
1274  elementGlobal,
1275  faceCoord[0]);
1276  subContVol[1].volume_ =
1277  hexahedronVolume(subContVol[1].global,
1278  edgeCoord[5],
1279  faceCoord[4],
1280  edgeCoord[6],
1281  edgeCoord[1],
1282  faceCoord[1],
1283  elementGlobal,
1284  faceCoord[2]);
1285  subContVol[2].volume_ =
1286  hexahedronVolume(subContVol[2].global,
1287  edgeCoord[4],
1288  faceCoord[4],
1289  edgeCoord[7],
1290  edgeCoord[2],
1291  faceCoord[0],
1292  elementGlobal,
1293  faceCoord[3]);
1294  subContVol[3].volume_ =
1295  hexahedronVolume(subContVol[3].global,
1296  edgeCoord[7],
1297  faceCoord[4],
1298  edgeCoord[5],
1299  edgeCoord[3],
1300  faceCoord[3],
1301  elementGlobal,
1302  faceCoord[1]);
1303  subContVol[4].volume_ =
1304  hexahedronVolume(edgeCoord[0],
1305  faceCoord[2],
1306  elementGlobal,
1307  faceCoord[0],
1308  subContVol[4].global,
1309  edgeCoord[10],
1310  faceCoord[5],
1311  edgeCoord[8]);
1312  subContVol[5].volume_ =
1313  hexahedronVolume(edgeCoord[1],
1314  faceCoord[1],
1315  elementGlobal,
1316  faceCoord[2],
1317  subContVol[5].global,
1318  edgeCoord[9],
1319  faceCoord[5],
1320  edgeCoord[10]);
1321  subContVol[6].volume_ =
1322  hexahedronVolume(edgeCoord[2],
1323  faceCoord[0],
1324  elementGlobal,
1325  faceCoord[3],
1326  subContVol[6].global,
1327  edgeCoord[8],
1328  faceCoord[5],
1329  edgeCoord[11]);
1330  subContVol[7].volume_ =
1331  hexahedronVolume(edgeCoord[3],
1332  faceCoord[3],
1333  elementGlobal,
1334  faceCoord[1],
1335  subContVol[7].global,
1336  edgeCoord[11],
1337  faceCoord[5],
1338  edgeCoord[9]);
1339  break;
1340  default:
1341  OPM_THROW(std::logic_error,
1342  "Not implemented:VcfvStencil for dim = " << dim
1343  << ", numVertices = " << numVertices);
1344  }
1345  }
1346  else
1347  OPM_THROW(std::logic_error, "Not implemented:VcfvStencil for dim = " << dim);
1348  }
1349 #if __GNUC__ || __clang__
1350 #pragma GCC diagnostic pop
1351 #endif
1352 
1353  const GridView& gridView_;
1354  const VertexMapper& vertexMapper_;
1355 
1356  Element element_;
1357 
1358 #if HAVE_DUNE_LOCALFUNCTIONS
1359  static LocalFiniteElementCache feCache_;
1360 #endif // HAVE_DUNE_LOCALFUNCTIONS
1361 
1363  LocalPosition elementLocal;
1365  GlobalPosition elementGlobal;
1367  Scalar elementVolume;
1369  SubControlVolume subContVol[maxNC];
1371  SubControlVolumeFace subContVolFace[maxNE];
1373  BoundaryFace boundaryFace_[maxBF];
1374  unsigned numBoundarySegments_;
1376  GlobalPosition edgeCoord[maxNE];
1378  GlobalPosition faceCoord[maxNF];
1380  unsigned numVertices;
1382  unsigned numEdges;
1384  unsigned numFaces;
1385  Dune::GeometryType geometryType_;
1386 };
1387 
1388 #if HAVE_DUNE_LOCALFUNCTIONS
1389 template<class Scalar, class GridView>
1390 typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
1391 VcfvStencil<Scalar, GridView>::feCache_;
1392 #endif // HAVE_DUNE_LOCALFUNCTIONS
1393 
1394 } // namespace Ewoms
1395 
1396 #endif
1397 
DimVector normal_
normal on face pointing to CV j or outward of the domain with length equal to |scvf| ...
Definition: vcfvstencil.hh:808
LocalPosition local
local vertex position
Definition: vcfvstencil.hh:769
GlobalPosition ipGlobal_
integration point in global coords
Definition: vcfvstencil.hh:806
Quadrature geometry for quadrilaterals.
Definition: quadraturegeometries.hh:39
Entity entity(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1107
finite volume intersected with element
Definition: vcfvstencil.hh:751
Scalar area_
area of face
Definition: vcfvstencil.hh:810
GlobalPosition global
global vertex position
Definition: vcfvstencil.hh:771
unsigned short i
scvf seperates corner i and j of elem
Definition: vcfvstencil.hh:802
void updateTopology(const Element &e)
Update the non-geometric part of the stencil.
Definition: vcfvstencil.hh:838
interior face of a sub control volume
Definition: vcfvstencil.hh:781
ScvGeometry geometry_
The geometry of the sub-control volume in local coordinates.
Definition: vcfvstencil.hh:775
Definition: vcfvstencil.hh:732
LocalPosition ipLocal_
integration point in local coords
Definition: vcfvstencil.hh:804
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition: quadraturegeometries.hh:69
Quadrature geometry for quadrilaterals.
const GlobalPosition & corner(unsigned cornerIdx) const
Return the position of the corner with a given index.
Definition: quadraturegeometries.hh:127
VertexMapper Mapper
exported Mapper type
Definition: vcfvstencil.hh:730
SubControlVolumeFace BoundaryFace
compatibility typedef
Definition: vcfvstencil.hh:814
unsigned globalSpaceIndex(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1096
Represents the finite volume geometry of a single element in the VCFV discretization.
Definition: vcfvstencil.hh:467
Dune::FieldVector< DimVector, maxNC > gradCenter
derivative of shape function at the center of the sub control volume
Definition: vcfvstencil.hh:778
Scalar volume_
space occupied by the sub-control volume
Definition: vcfvstencil.hh:773