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  Copyright (C) 2010-2013 by Andreas Lauser
5  Copyright (C) 2010 by Klaus Mosthaf
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 2 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
27 #ifndef EWOMS_VCFV_STENCIL_HH
28 #define EWOMS_VCFV_STENCIL_HH
29 
30 #include <ewoms/common/quadraturegeometries.hh>
31 
32 #include <opm/common/Exceptions.hpp>
33 #include <opm/common/ErrorMacros.hpp>
34 
35 #include <dune/grid/common/intersectioniterator.hh>
36 #include <dune/grid/common/mcmgmapper.hh>
37 #include <dune/geometry/referenceelements.hh>
38 #include <dune/localfunctions/lagrange/pqkfactory.hh>
39 #include <dune/common/version.hh>
40 
41 #include <vector>
42 
43 namespace Ewoms {
44 
48 template <class Scalar, int dim, int basicGeomType>
49 class VcfvScvGeometries;
50 
52 // local geometries for 1D elements
54 template <class Scalar>
55 class VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>
56 {
57  enum { dim = 1 };
58  enum { numScv = 2 };
59  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
60 
61 public:
62  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
63 
64  static int init()
65  {
66  // 1D LINE SEGMENTS
67  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] = {
68  { // corners of the first sub control volume
69  { 0.0 },
70  { 0.5 }
71  },
72 
73  { // corners of the second sub control volume
74  { 0.5 },
75  { 1.0 }
76  }
77  };
78  for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
79  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
80 
81  return 0;
82  }
83 
84  static const ScvLocalGeometry &get(int scvIdx)
85  { return scvGeoms_[scvIdx]; }
86 
87 private:
88  static ScvLocalGeometry scvGeoms_[numScv];
89 };
90 
91 template <class Scalar>
92 typename VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::ScvLocalGeometry
93 VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::scvGeoms_[
94  VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::cube>::numScv];
95 
96 template <class Scalar>
97 class VcfvScvGeometries<Scalar, /*dim=*/1, Dune::GeometryType::simplex>
98 {
99  enum { dim = 1 };
100  enum { numScv = 2 };
101  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
102 
103 public:
104  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
105 
106  static const ScvLocalGeometry &get(int scvIdx)
107  { OPM_THROW(std::logic_error,
108  "Not implemented: "
109  "VcfvScvGeometries<Scalar, 1, Dune::GeometryType::simplex>"); }
110 };
111 
113 // local geometries for 2D elements
115 template <class Scalar>
116 class VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>
117 {
118  enum { dim = 2 };
119  enum { numScv = 3 };
120  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
121 
122 public:
123  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
124 
125  static const ScvLocalGeometry &get(int scvIdx)
126  { return scvGeoms_[scvIdx]; }
127 
128  static int init()
129  {
130  // 2D SIMPLEX
131  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
132  {
133  { // SCV 0 corners
134  { 0.0, 0.0 },
135  { 1.0/2, 0.0 },
136  { 1.0/3, 1.0/3 },
137  { 0.0, 1.0/2 },
138  },
139 
140  { // SCV 1 corners
141  { 1.0/2, 0.0 },
142  { 1.0, 0.0 },
143  { 1.0/3, 1.0/3 },
144  { 1.0/2, 1.0/2 },
145  },
146 
147  { // SCV 2 corners
148  { 0.0, 1.0/2 },
149  { 1.0/3, 1.0/3 },
150  { 0.0, 1.0 },
151  { 1.0/2, 1.0/2 },
152  }
153  };
154 
155  for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
156  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
157 
158  return 0;
159  }
160 
161 private:
162  static ScvLocalGeometry scvGeoms_[numScv];
163 };
164 
165 template <class Scalar>
166 typename VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::ScvLocalGeometry
167 VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::scvGeoms_[
168  VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::simplex>::numScv];
169 
170 template <class Scalar>
171 class VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>
172 {
173  enum { dim = 2 };
174  enum { numScv = 4 };
175  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::cube;
176 
177 public:
178  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
179 
180  static const ScvLocalGeometry &get(int scvIdx)
181  { return scvGeoms_[scvIdx]; }
182 
183  static int init()
184  {
185  // 2D CUBE
186  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
187  {
188  { // SCV 0 corners
189  { 0.0, 0.0 },
190  { 0.5, 0.0 },
191  { 0.0, 0.5 },
192  { 0.5, 0.5 }
193  },
194 
195  { // SCV 1 corners
196  { 0.5, 0.0 },
197  { 1.0, 0.0 },
198  { 0.5, 0.5 },
199  { 1.0, 0.5 }
200  },
201 
202  { // SCV 2 corners
203  { 0.0, 0.5 },
204  { 0.5, 0.5 },
205  { 0.0, 1.0 },
206  { 0.5, 1.0 }
207  },
208 
209  { // SCV 3 corners
210  { 0.5, 0.5 },
211  { 1.0, 0.5 },
212  { 0.5, 1.0 },
213  { 1.0, 1.0 }
214  }
215  };
216 
217  for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
218  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
219 
220  return 0;
221  }
222 
223  static ScvLocalGeometry scvGeoms_[numScv];
224 };
225 
226 template <class Scalar>
227 typename VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::ScvLocalGeometry
228 VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::scvGeoms_[
229  VcfvScvGeometries<Scalar, /*dim=*/2, Dune::GeometryType::cube>::numScv];
230 
232 // local geometries for 3D elements
234 template <class Scalar>
235 class VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::simplex>
236 {
237  enum { dim = 3 };
238  enum { numScv = 4 };
239  static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
240 
241 public:
242  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
243 
244  static const ScvLocalGeometry &get(int scvIdx)
245  { return scvGeoms_[scvIdx]; }
246 
247  static int init()
248  {
249  // 3D SIMPLEX
250  Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
251  {
252  { // SCV 0 corners
253  { 0.0, 0.0, 0.0 },
254  { 1.0/2, 0.0, 0.0 },
255  { 0.0, 1.0/2, 0.0 },
256  { 1.0/3, 1.0/3, 0.0 },
257 
258  { 0.0, 0.0, 0.5 },
259  { 1.0/3, 0.0, 1.0/3 },
260  { 0.0, 1.0/3, 1.0/3 },
261  { 1.0/4, 1.0/4, 1.0/4 },
262  },
263 
264  { // SCV 1 corners
265  { 1.0/2, 0.0, 0.0 },
266  { 1.0, 0.0, 0.0 },
267  { 1.0/3, 1.0/3, 0.0 },
268  { 1.0/2, 1.0/2, 0.0 },
269 
270  { 1.0/3, 0.0, 1.0/3 },
271  { 1.0/2, 0.0, 1.0/2 },
272  { 1.0/4, 1.0/4, 1.0/4 },
273  { 1.0/3, 1.0/3, 1.0/3 },
274  },
275 
276  { // SCV 2 corners
277  { 0.0, 1.0/2, 0.0 },
278  { 1.0/3, 1.0/3, 0.0 },
279  { 0.0, 1.0, 0.0 },
280  { 1.0/2, 1.0/2, 0.0 },
281 
282  { 0.0, 1.0/3, 1.0/3 },
283  { 1.0/4, 1.0/4, 1.0/4 },
284  { 0.0, 1.0/2, 1.0/2 },
285  { 1.0/3, 1.0/3, 1.0/3 },
286  },
287 
288  { // SCV 3 corners
289  { 0.0, 0.0, 1.0/2 },
290  { 1.0/3, 0.0, 1.0/3 },
291  { 0.0, 1.0/3, 1.0/3 },
292  { 1.0/4, 1.0/4, 1.0/4 },
293 
294  { 0.0, 0.0, 1.0 },
295  { 1.0/2, 0.0, 1.0/2 },
296  { 0.0, 1.0/2, 1.0/2 },
297  { 1.0/3, 1.0/3, 1.0/3 },
298  }
299  };
300 
301  for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
302  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
303 
304  return 0;
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:
324  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
325 
326  static const ScvLocalGeometry &get(int scvIdx)
327  { return scvGeoms_[scvIdx]; }
328 
329  static int 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 (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
432  scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
433 
434  return 0;
435  }
436 private:
437  static ScvLocalGeometry scvGeoms_[numScv];
438 };
439 
440 template <class Scalar>
441 typename VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::ScvLocalGeometry
442 VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::scvGeoms_[
443  VcfvScvGeometries<Scalar, /*dim=*/3, Dune::GeometryType::cube>::numScv];
444 
468 template <class Scalar, class GridView>
469 class VcfvStencil
470 {
471  enum{dim = GridView::dimension};
472  enum{dimWorld = GridView::dimensionworld};
473  enum{maxNC = (dim < 3 ? 4 : 8)};
474  enum{maxNE = (dim < 3 ? 4 : 12)};
475  enum{maxNF = (dim < 3 ? 1 : 6)};
476  enum{maxBF = (dim < 3 ? 8 : 24)};
477  typedef typename GridView::ctype CoordScalar;
478  typedef typename GridView::Traits::template Codim<0>::Entity Element;
479  typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
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 
486  typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> ScvLocalGeometry;
487 
488  typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
489  typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
490  typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
491  typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
492 
493  Scalar quadrilateralArea(const GlobalPosition& p0,
494  const GlobalPosition& p1,
495  const GlobalPosition& p2,
496  const GlobalPosition& p3)
497  {
498  return 0.5*std::abs((p3[0] - p1[0])*(p2[1] - p0[1]) - (p3[1] - p1[1])*(p2[0] - p0[0]));
499  }
500 
501  void crossProduct(DimVector& c, const DimVector& a, const DimVector& b)
502  {
503  c[0] = a[1]*b[2] - a[2]*b[1];
504  c[1] = a[2]*b[0] - a[0]*b[2];
505  c[2] = a[0]*b[1] - a[1]*b[0];
506  }
507 
508  Scalar pyramidVolume(const GlobalPosition& p0,
509  const GlobalPosition& p1,
510  const GlobalPosition& p2,
511  const GlobalPosition& p3,
512  const GlobalPosition& p4)
513  {
514  DimVector a(p2); a -= p0;
515  DimVector b(p3); b -= p1;
516 
517  DimVector n;
518  crossProduct(n, a, b);
519 
520  a = p4; a -= p0;
521 
522  return 1.0/6.0*(n*a);
523  }
524 
525  Scalar prismVolume(const GlobalPosition& p0,
526  const GlobalPosition& p1,
527  const GlobalPosition& p2,
528  const GlobalPosition& p3,
529  const GlobalPosition& p4,
530  const GlobalPosition& p5)
531  {
532  DimVector a(p4);
533  for (int k = 0; k < dimWorld; ++k)
534  a[k] -= p0[k];
535  DimVector b(p1);
536  for (int k = 0; k < dimWorld; ++k)
537  b[k] -= p3[k];
538  DimVector m;
539  crossProduct(m, a, b);
540 
541  for (int k = 0; k < dimWorld; ++k)
542  a[k] = p1[k] - p0[k];
543  for (int k = 0; k < dimWorld; ++k)
544  b[k] = p2[k] - p0[k];
545  DimVector n;
546  crossProduct(n, a, b);
547  n += m;
548 
549  for (int k = 0; k < dimWorld; ++k)
550  a[k] = p5[k] - p0[k];
551 
552  return std::abs(1.0/6.0*(n*a));
553  }
554 
555  Scalar hexahedronVolume(const GlobalPosition& p0,
556  const GlobalPosition& p1,
557  const GlobalPosition& p2,
558  const GlobalPosition& p3,
559  const GlobalPosition& p4,
560  const GlobalPosition& p5,
561  const GlobalPosition& p6,
562  const GlobalPosition& p7)
563  {
564  return
565  prismVolume(p0,p1,p2,p4,p5,p6)
566  + prismVolume(p0,p2,p3,p4,p6,p7);
567  }
568 
569  void normalOfQuadrilateral3D(DimVector &normal,
570  const GlobalPosition& p0,
571  const GlobalPosition& p1,
572  const GlobalPosition& p2,
573  const GlobalPosition& p3)
574  {
575  DimVector a(p2);
576  for (int k = 0; k < dimWorld; ++k)
577  a[k] -= p0[k];
578  DimVector b(p3);
579  for (int k = 0; k < dimWorld; ++k)
580  b[k] -= p1[k];
581 
582  crossProduct(normal, a, b);
583  normal *= 0.5;
584  }
585 
586  Scalar quadrilateralArea3D(const GlobalPosition& p0,
587  const GlobalPosition& p1,
588  const GlobalPosition& p2,
589  const GlobalPosition& p3)
590  {
591  DimVector normal;
592  normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
593  return normal.two_norm();
594  }
595 
596  void getFaceIndices(int numElemVertices, int k, int& leftFace, int& rightFace)
597  {
598  static const int edgeToFaceTet[2][6] = {
599  {1, 0, 3, 2, 1, 3},
600  {0, 2, 0, 1, 3, 2}
601  };
602  static const int edgeToFacePyramid[2][8] = {
603  {1, 2, 3, 4, 1, 3, 4, 2},
604  {0, 0, 0, 0, 3, 2, 1, 4}
605  };
606  static const int edgeToFacePrism[2][9] = {
607  {1, 0, 2, 0, 1, 2, 4, 4, 4},
608  {0, 2, 1, 3, 3, 3, 0, 1, 2}
609  };
610  static const int edgeToFaceHex[2][12] = {
611  {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
612  {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
613  };
614 
615  switch (numElemVertices) {
616  case 4:
617  leftFace = edgeToFaceTet[0][k];
618  rightFace = edgeToFaceTet[1][k];
619  break;
620  case 5:
621  leftFace = edgeToFacePyramid[0][k];
622  rightFace = edgeToFacePyramid[1][k];
623  break;
624  case 6:
625  leftFace = edgeToFacePrism[0][k];
626  rightFace = edgeToFacePrism[1][k];
627  break;
628  case 8:
629  leftFace = edgeToFaceHex[0][k];
630  rightFace = edgeToFaceHex[1][k];
631  break;
632  default:
633  OPM_THROW(std::logic_error,
634  "Not implemented: VcfvStencil::getFaceIndices for "
635  << numElemVertices << " vertices");
636  break;
637  }
638  }
639 
640  void getEdgeIndices(int numElemVertices, int face, int vert, int& leftEdge, int& 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 = faceAndVertexToLeftEdgeTet[face][vert];
702  rightEdge = faceAndVertexToRightEdgeTet[face][vert];
703  break;
704  case 5:
705  leftEdge = faceAndVertexToLeftEdgePyramid[face][vert];
706  rightEdge = faceAndVertexToRightEdgePyramid[face][vert];
707  break;
708  case 6:
709  leftEdge = faceAndVertexToLeftEdgePrism[face][vert];
710  rightEdge = faceAndVertexToRightEdgePrism[face][vert];
711  break;
712  case 8:
713  leftEdge = faceAndVertexToLeftEdgeHex[face][vert];
714  rightEdge = faceAndVertexToRightEdgeHex[face][vert];
715  break;
716  default:
717  OPM_THROW(std::logic_error,
718  "Not implemented: VcfvStencil::getFaceIndices for "
719  << numElemVertices << " vertices");
720  break;
721  }
722  }
723 
724 public:
725  typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,
726  Dune::MCMGVertexLayout > VertexMapper;
727 
728  class ScvGeometry
729  {
730  public:
731  const GlobalPosition center() const
732  { return global(localGeometry_->center()); }
733 
734  const GlobalPosition corner(int cornerIdx) const
735  { return global(localGeometry_->corner(cornerIdx)); }
736 
737  const GlobalPosition global(const LocalPosition &localPos) const
738  { return element_->geometry().global(localPos); }
739 
740  const ScvLocalGeometry &localGeometry() const
741  { return *localGeometry_; }
742 
743  const ScvLocalGeometry *localGeometry_;
744  const Element *element_;
745  };
746 
747  struct SubControlVolume
748  {
749  const GlobalPosition &globalPos() const
750  { return global; }
751 
752  const GlobalPosition center() const
753  { return geometry_.center(); }
754 
755  Scalar volume() const
756  { return volume_; }
757 
758  const ScvLocalGeometry &localGeometry() const
759  { return geometry_.localGeometry(); }
760 
761  const ScvGeometry &geometry() const
762  { return geometry_; }
763 
765  LocalPosition local;
767  GlobalPosition global;
769  Scalar volume_;
771  ScvGeometry geometry_;
772 
774  Dune::FieldVector<DimVector, maxNC> gradCenter;
775  };
776 
777  struct SubControlVolumeFace
778  {
779  const DimVector &normal() const
780  { return normal_; }
781 
782  int interiorIndex() const
783  { return i; }
784 
785  int exteriorIndex() const
786  { return j; }
787 
788  Scalar area() const
789  { return area_; }
790 
791  const LocalPosition &localPos() const
792  { return ipLocal_; }
793 
794  const GlobalPosition &integrationPos() const
795  { return ipGlobal_; }
796 
798  int i,j;
800  LocalPosition ipLocal_;
802  GlobalPosition ipGlobal_;
804  DimVector normal_;
806  Scalar area_;
807  };
808 
810  typedef SubControlVolumeFace BoundaryFace;
811 
812  VcfvStencil(const GridView &gridView)
813  : gridView_(gridView)
814  , vertexMapper_(gridView)
815 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
816  , element_(*gridView.template begin</*codim=*/0>())
817 #else
818  , elementPtr_(gridView.template begin</*codim=*/0>())
819 #endif
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 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
841  element_ = e;
842 
843  numVertices = e.subEntities(/*codim=*/dim);
844  numEdges = e.subEntities(/*codim=*/dim-1);
845  numFaces = (dim<3)?0:e.subEntities(/*codim=*/1);
846 #else
847  elementPtr_ = ElementPointer(e);
848 
849  numVertices = e.template count</*codim=*/dim>();
850  numEdges = e.template count</*codim=*/dim-1>();
851  numFaces = (dim<3)?0:e.template count</*codim=*/1>();
852 #endif
853 
854  numBoundarySegments_ = 0; // TODO: really required here(?)
855  }
856 
857  void update(const Element& e)
858  {
859  updateTopology(e);
860 
861  const Geometry& geometry = e.geometry();
862  geometryType_ = geometry.type();
863 
864 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
865  const typename Dune::ReferenceElementContainer<CoordScalar,dim>::value_type&
866  referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
867 #else
868  const typename Dune::GenericReferenceElementContainer<CoordScalar,dim>::value_type&
869  referenceElement = Dune::GenericReferenceElements<CoordScalar,
870  dim>::general(geometryType_);
871 #endif
872 
873  elementVolume = geometry.volume();
874  elementLocal = referenceElement.position(0,0);
875  elementGlobal = geometry.global(elementLocal);
876 
877  // corners:
878  for (int vert = 0; vert < numVertices; vert++) {
879  subContVol[vert].local = referenceElement.position(vert, dim);
880  subContVol[vert].global = geometry.global(subContVol[vert].local);
881  }
882 
883  // edges:
884  for (int edge = 0; edge < numEdges; edge++) {
885  edgeCoord[edge] = geometry.global(referenceElement.position(edge, dim-1));
886  }
887 
888  // faces:
889  for (int face = 0; face < numFaces; face++) {
890  faceCoord[face] = geometry.global(referenceElement.position(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 (int k = 0; k < numEdges; k++) { // begin loop over edges / sub control volume faces
902  int i = referenceElement.subEntity(k, dim-1, 0, dim);
903  int j = referenceElement.subEntity(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(k, dim-1) + elementLocal;
924  ipLocal_ *= 0.5;
925  subContVolFace[k].ipLocal_ = ipLocal_;
926  for (int 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 (int 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  int leftFace;
942  int rightFace;
943  getFaceIndices(numVertices, k, leftFace, rightFace);
944  ipLocal_ = referenceElement.position(k, dim-1) + elementLocal
945  + referenceElement.position(leftFace, 1)
946  + referenceElement.position(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  int face = intersection.indexInInside();
969  int numVerticesOfFace = referenceElement.size(face, 1, dim);
970  for (int vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
971  {
972  int vertInElement = referenceElement.subEntity(face, 1, vertInFace, dim);
973  int bfIdx = numBoundarySegments_;
974  ++numBoundarySegments_;
975 
976  if (dim == 1) {
977  boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(vertInElement, dim);
978  boundaryFace_[bfIdx].area_ = 1.0;
979  }
980  else if (dim == 2) {
981  boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(vertInElement, dim)
982  + referenceElement.position(face, 1);
983  boundaryFace_[bfIdx].ipLocal_ *= 0.5;
984  boundaryFace_[bfIdx].area_ = 0.5 * intersection.geometry().volume();
985  }
986  else if (dim == 3) {
987  int leftEdge;
988  int rightEdge;
989  getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
990  boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(vertInElement, dim)
991  + referenceElement.position(face, 1)
992  + referenceElement.position(leftEdge, dim-1)
993  + referenceElement.position(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  OPM_THROW(std::logic_error, "Not implemented:VcfvStencil for dim = " << 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 (int vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1023  subContVol[vertIdx].geometry_.element_ = &element;
1024  subContVol[vertIdx].geometry_.localGeometry_ =
1025  &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::simplex>::get(vertIdx);
1026  }
1027  }
1028  else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1029  for (int vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1030  subContVol[vertIdx].geometry_.element_ = &element;
1031  subContVol[vertIdx].geometry_.localGeometry_ =
1032  &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::cube>::get(vertIdx);
1033  }
1034  }
1035  else
1036  OPM_THROW(std::logic_error,
1037  "Not implemented: SCV geometries for non hexahedron elements");
1038  }
1039 
1040  void updateCenterGradients()
1041  {
1042 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1043  const auto &localFiniteElement = feCache_.get(element_.type());
1044  const auto &geom = element_.geometry();
1045 #else
1046  const auto &localFiniteElement = feCache_.get(elementPtr_->type());
1047  const auto &geom = elementPtr_->geometry();
1048 #endif
1049  std::vector<ShapeJacobian> localJac;
1050 
1051  for (int scvIdx = 0; scvIdx < numVertices; ++ scvIdx) {
1052  const auto &localCenter = subContVol[scvIdx].localGeometry().center();
1053  localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
1054  const auto &globalPos = subContVol[scvIdx].geometry().center();
1055 
1056  const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
1057  for (int vert = 0; vert < numVertices; vert++) {
1058  jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
1059  }
1060  }
1061  }
1062 
1063  int numDof() const
1064  { return numVertices; }
1065 
1066  int numPrimaryDof() const
1067  { return numDof(); }
1068 
1069  Dune::PartitionType partitionType(int scvIdx) const
1070 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1071  { return element_.template subEntity</*codim=*/dim>(scvIdx)->partitionType(); }
1072 #else
1073  { return elementPtr_->template subEntity</*codim=*/dim>(scvIdx)->partitionType(); }
1074 #endif
1075 
1076  const SubControlVolume &subControlVolume(int scvIdx) const
1077  {
1078  assert(0 <= scvIdx && scvIdx < numDof());
1079  return subContVol[scvIdx];
1080  }
1081 
1082  int numInteriorFaces() const
1083  { return numEdges; }
1084 
1085  int numBoundaryFaces() const
1086  { return numBoundarySegments_; }
1087 
1088  const SubControlVolumeFace &interiorFace(int faceIdx) const
1089  { return subContVolFace[faceIdx]; }
1090 
1091  const BoundaryFace &boundaryFace(int bfIdx) const
1092  { return boundaryFace_[bfIdx]; }
1093 
1098  int globalSpaceIndex(int dofIdx) const
1099  {
1100  assert(0 <= dofIdx && dofIdx < numDof());
1101 
1102 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1103  return vertexMapper_.subIndex(element_, dofIdx, /*codim=*/dim);
1104 #else
1105  return vertexMapper_.map(*elementPtr_, dofIdx, /*codim=*/dim);
1106 #endif
1107  }
1108 
1109 private:
1110  void fillSubContVolData_()
1111  {
1112  if (dim == 1) {
1113  // 1D
1114  subContVol[0].volume_ = 0.5*elementVolume;
1115  subContVol[1].volume_ = 0.5*elementVolume;
1116  }
1117  else if (dim == 2) {
1118  switch (numVertices) {
1119  case 3: // 2D, triangle
1120  subContVol[0].volume_ = elementVolume/3;
1121  subContVol[1].volume_ = elementVolume/3;
1122  subContVol[2].volume_ = elementVolume/3;
1123  break;
1124  case 4: // 2D, quadrilinear
1125  subContVol[0].volume_ =
1126  quadrilateralArea(subContVol[0].global,
1127  edgeCoord[2],
1128  elementGlobal,
1129  edgeCoord[0]);
1130  subContVol[1].volume_ =
1131  quadrilateralArea(subContVol[1].global,
1132  edgeCoord[1],
1133  elementGlobal,
1134  edgeCoord[2]);
1135  subContVol[2].volume_ =
1136  quadrilateralArea(subContVol[2].global,
1137  edgeCoord[0],
1138  elementGlobal,
1139  edgeCoord[3]);
1140  subContVol[3].volume_ =
1141  quadrilateralArea(subContVol[3].global,
1142  edgeCoord[3],
1143  elementGlobal,
1144  edgeCoord[1]);
1145  break;
1146  default:
1147  OPM_THROW(std::logic_error,
1148  "Not implemented:VcfvStencil dim = " << dim
1149  << ", numVertices = " << numVertices);
1150  }
1151  }
1152  else if (dim == 3) {
1153  switch (numVertices) {
1154  case 4: // 3D, tetrahedron
1155  for (int k = 0; k < numVertices; k++)
1156  subContVol[k].volume_ = elementVolume / 4.0;
1157  break;
1158  case 5: // 3D, pyramid
1159  subContVol[0].volume_ = hexahedronVolume(subContVol[0].global,
1160  edgeCoord[2],
1161  faceCoord[0],
1162  edgeCoord[0],
1163  edgeCoord[4],
1164  faceCoord[3],
1165  elementGlobal,
1166  faceCoord[1]);
1167  subContVol[1].volume_ = hexahedronVolume(subContVol[1].global,
1168  edgeCoord[1],
1169  faceCoord[0],
1170  edgeCoord[2],
1171  edgeCoord[5],
1172  faceCoord[2],
1173  elementGlobal,
1174  faceCoord[3]);
1175  subContVol[2].volume_ = hexahedronVolume(subContVol[2].global,
1176  edgeCoord[0],
1177  faceCoord[0],
1178  edgeCoord[3],
1179  edgeCoord[6],
1180  faceCoord[1],
1181  elementGlobal,
1182  faceCoord[4]);
1183  subContVol[3].volume_ = hexahedronVolume(subContVol[3].global,
1184  edgeCoord[3],
1185  faceCoord[0],
1186  edgeCoord[1],
1187  edgeCoord[7],
1188  faceCoord[4],
1189  elementGlobal,
1190  faceCoord[2]);
1191  subContVol[4].volume_ = elementVolume -
1192  subContVol[0].volume_ -
1193  subContVol[1].volume_ -
1194  subContVol[2].volume_ -
1195  subContVol[3].volume_;
1196  break;
1197  case 6: // 3D, prism
1198  subContVol[0].volume_ = hexahedronVolume(subContVol[0].global,
1199  edgeCoord[3],
1200  faceCoord[3],
1201  edgeCoord[4],
1202  edgeCoord[0],
1203  faceCoord[0],
1204  elementGlobal,
1205  faceCoord[1]);
1206  subContVol[1].volume_ = hexahedronVolume(subContVol[1].global,
1207  edgeCoord[5],
1208  faceCoord[3],
1209  edgeCoord[3],
1210  edgeCoord[1],
1211  faceCoord[2],
1212  elementGlobal,
1213  faceCoord[0]);
1214  subContVol[2].volume_ = hexahedronVolume(subContVol[2].global,
1215  edgeCoord[4],
1216  faceCoord[3],
1217  edgeCoord[5],
1218  edgeCoord[2],
1219  faceCoord[1],
1220  elementGlobal,
1221  faceCoord[2]);
1222  subContVol[3].volume_ = hexahedronVolume(edgeCoord[0],
1223  faceCoord[0],
1224  elementGlobal,
1225  faceCoord[1],
1226  subContVol[3].global,
1227  edgeCoord[6],
1228  faceCoord[4],
1229  edgeCoord[7]);
1230  subContVol[4].volume_ = hexahedronVolume(edgeCoord[1],
1231  faceCoord[2],
1232  elementGlobal,
1233  faceCoord[0],
1234  subContVol[4].global,
1235  edgeCoord[8],
1236  faceCoord[4],
1237  edgeCoord[6]);
1238  subContVol[5].volume_ = hexahedronVolume(edgeCoord[2],
1239  faceCoord[1],
1240  elementGlobal,
1241  faceCoord[2],
1242  subContVol[5].global,
1243  edgeCoord[7],
1244  faceCoord[4],
1245  edgeCoord[8]);
1246  break;
1247  case 8: // 3D, hexahedron
1248  subContVol[0].volume_ = hexahedronVolume(subContVol[0].global,
1249  edgeCoord[6],
1250  faceCoord[4],
1251  edgeCoord[4],
1252  edgeCoord[0],
1253  faceCoord[2],
1254  elementGlobal,
1255  faceCoord[0]);
1256  subContVol[1].volume_ = hexahedronVolume(subContVol[1].global,
1257  edgeCoord[5],
1258  faceCoord[4],
1259  edgeCoord[6],
1260  edgeCoord[1],
1261  faceCoord[1],
1262  elementGlobal,
1263  faceCoord[2]);
1264  subContVol[2].volume_ = hexahedronVolume(subContVol[2].global,
1265  edgeCoord[4],
1266  faceCoord[4],
1267  edgeCoord[7],
1268  edgeCoord[2],
1269  faceCoord[0],
1270  elementGlobal,
1271  faceCoord[3]);
1272  subContVol[3].volume_ = hexahedronVolume(subContVol[3].global,
1273  edgeCoord[7],
1274  faceCoord[4],
1275  edgeCoord[5],
1276  edgeCoord[3],
1277  faceCoord[3],
1278  elementGlobal,
1279  faceCoord[1]);
1280  subContVol[4].volume_ = hexahedronVolume(edgeCoord[0],
1281  faceCoord[2],
1282  elementGlobal,
1283  faceCoord[0],
1284  subContVol[4].global,
1285  edgeCoord[10],
1286  faceCoord[5],
1287  edgeCoord[8]);
1288  subContVol[5].volume_ = hexahedronVolume(edgeCoord[1],
1289  faceCoord[1],
1290  elementGlobal,
1291  faceCoord[2],
1292  subContVol[5].global,
1293  edgeCoord[9],
1294  faceCoord[5],
1295  edgeCoord[10]);
1296  subContVol[6].volume_ = hexahedronVolume(edgeCoord[2],
1297  faceCoord[0],
1298  elementGlobal,
1299  faceCoord[3],
1300  subContVol[6].global,
1301  edgeCoord[8],
1302  faceCoord[5],
1303  edgeCoord[11]);
1304  subContVol[7].volume_ = hexahedronVolume(edgeCoord[3],
1305  faceCoord[3],
1306  elementGlobal,
1307  faceCoord[1],
1308  subContVol[7].global,
1309  edgeCoord[11],
1310  faceCoord[5],
1311  edgeCoord[9]);
1312  break;
1313  default:
1314  OPM_THROW(std::logic_error,
1315  "Not implemented:VcfvStencil for dim = " << dim
1316  << ", numVertices = " << numVertices);
1317  }
1318  }
1319  else
1320  OPM_THROW(std::logic_error, "Not implemented:VcfvStencil for dim = " << dim);
1321  }
1322 
1323  GridView gridView_;
1324  VertexMapper vertexMapper_;
1325 
1326 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1327  Element element_;
1328 #else
1329  ElementPointer elementPtr_;
1330 #endif
1331 
1332  static LocalFiniteElementCache feCache_;
1333 
1335  LocalPosition elementLocal;
1337  GlobalPosition elementGlobal;
1339  Scalar elementVolume;
1341  SubControlVolume subContVol[maxNC];
1343  SubControlVolumeFace subContVolFace[maxNE];
1345  BoundaryFace boundaryFace_[maxBF];
1346  int numBoundarySegments_;
1348  GlobalPosition edgeCoord[maxNE];
1350  GlobalPosition faceCoord[maxNF];
1352  int numVertices;
1354  int numEdges;
1356  int numFaces;
1357  Dune::GeometryType geometryType_;
1358 };
1359 
1360 template<class Scalar, class GridView>
1361 typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
1362 VcfvStencil<Scalar, GridView>::feCache_;
1363 
1364 } // namespace Ewoms
1365 
1366 #endif
1367 
Ewoms::VcfvStencil::SubControlVolumeFace::interiorIndex
int interiorIndex() const
Definition: vcfvstencil.hh:782
Ewoms::VcfvStencil::ScvGeometry::localGeometry
const ScvLocalGeometry & localGeometry() const
Definition: vcfvstencil.hh:740
Ewoms::VcfvStencil::SubControlVolumeFace::area_
Scalar area_
area of face
Definition: vcfvstencil.hh:806
Ewoms::VcfvStencil::subControlVolume
const SubControlVolume & subControlVolume(int scvIdx) const
Definition: vcfvstencil.hh:1076
quadraturegeometries.hh
Quadrature geometry for quadrilaterals.
Ewoms::VcfvStencil::updateCenterGradients
void updateCenterGradients()
Definition: vcfvstencil.hh:1040
Ewoms::VcfvStencil::SubControlVolumeFace::area
Scalar area() const
Definition: vcfvstencil.hh:788
Ewoms::QuadrialteralQuadratureGeometry::center
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition: quadraturegeometries.hh:67
Ewoms::VcfvStencil::ScvGeometry::global
const GlobalPosition global(const LocalPosition &localPos) const
Definition: vcfvstencil.hh:737
Ewoms::VcfvStencil::numInteriorFaces
int numInteriorFaces() const
Definition: vcfvstencil.hh:1082
Ewoms::VcfvStencil::BoundaryFace
SubControlVolumeFace BoundaryFace
compatibility typedef
Definition: vcfvstencil.hh:810
Ewoms::VcfvStencil::SubControlVolumeFace::exteriorIndex
int exteriorIndex() const
Definition: vcfvstencil.hh:785
Ewoms::VcfvStencil::SubControlVolume::global
GlobalPosition global
global vertex position
Definition: vcfvstencil.hh:767
Ewoms::VcfvStencil::SubControlVolume::local
LocalPosition local
local vertex position
Definition: vcfvstencil.hh:765
Ewoms::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:804
Ewoms::VcfvStencil::SubControlVolumeFace::i
int i
scvf seperates corner i and j of elem
Definition: vcfvstencil.hh:798
Ewoms::VcfvStencil::ScvGeometry::localGeometry_
const ScvLocalGeometry * localGeometry_
Definition: vcfvstencil.hh:743
Ewoms::VcfvStencil::SubControlVolume
finite volume intersected with element
Definition: vcfvstencil.hh:747
Ewoms::VcfvStencil::globalSpaceIndex
int globalSpaceIndex(int dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1098
Ewoms::VcfvStencil::SubControlVolumeFace::j
int j
Definition: vcfvstencil.hh:798
Ewoms::QuadrialteralQuadratureGeometry
Quadrature geometry for quadrilaterals.
Definition: quadraturegeometries.hh:37
Ewoms::VcfvStencil::SubControlVolume::gradCenter
Dune::FieldVector< DimVector, maxNC > gradCenter
derivative of shape function at the center of the sub control volume
Definition: vcfvstencil.hh:774
Ewoms::VcfvStencil::SubControlVolumeFace::ipGlobal_
GlobalPosition ipGlobal_
integration point in global coords
Definition: vcfvstencil.hh:802
Ewoms::VcfvStencil::ScvGeometry::corner
const GlobalPosition corner(int cornerIdx) const
Definition: vcfvstencil.hh:734
Dune
Definition: cartesianindexmapper.hh:31
Ewoms::VcfvStencil::SubControlVolume::volume_
Scalar volume_
space occupied by the sub-control volume
Definition: vcfvstencil.hh:769
Ewoms::VcfvStencil::SubControlVolume::geometry_
ScvGeometry geometry_
The geometry of the sub-control volume in local coordinates.
Definition: vcfvstencil.hh:771
Ewoms::VcfvStencil::partitionType
Dune::PartitionType partitionType(int scvIdx) const
Definition: vcfvstencil.hh:1069
Ewoms::VcfvStencil::SubControlVolumeFace
interior face of a sub control volume
Definition: vcfvstencil.hh:777
Ewoms::VcfvStencil::update
void update(const Element &e)
Definition: vcfvstencil.hh:857
Ewoms::VcfvStencil
Represents the finite volume geometry of a single element in the VCFV discretization.
Definition: vcfvstencil.hh:469
Ewoms::VcfvStencil::VcfvStencil
VcfvStencil(const GridView &gridView)
Definition: vcfvstencil.hh:812
Ewoms::VcfvStencil::SubControlVolumeFace::integrationPos
const GlobalPosition & integrationPos() const
Definition: vcfvstencil.hh:794
Ewoms::VcfvStencil::numDof
int numDof() const
Definition: vcfvstencil.hh:1063
Ewoms
Definition: baseauxiliarymodule.hh:35
Ewoms::VcfvStencil::SubControlVolumeFace::normal
const DimVector & normal() const
Definition: vcfvstencil.hh:779
Ewoms::VcfvStencil::ScvGeometry::element_
const Element * element_
Definition: vcfvstencil.hh:744
Ewoms::VcfvStencil::SubControlVolume::localGeometry
const ScvLocalGeometry & localGeometry() const
Definition: vcfvstencil.hh:758
Ewoms::VcfvStencil::SubControlVolume::globalPos
const GlobalPosition & globalPos() const
Definition: vcfvstencil.hh:749
Ewoms::VcfvStencil::updateTopology
void updateTopology(const Element &e)
Update the non-geometric part of the stencil.
Definition: vcfvstencil.hh:838
Ewoms::VcfvStencil::ScvGeometry::center
const GlobalPosition center() const
Definition: vcfvstencil.hh:731
Ewoms::VcfvStencil::VertexMapper
Dune::MultipleCodimMultipleGeomTypeMapper< GridView, Dune::MCMGVertexLayout > VertexMapper
Definition: vcfvstencil.hh:726
Ewoms::VcfvStencil::numPrimaryDof
int numPrimaryDof() const
Definition: vcfvstencil.hh:1066
Ewoms::VcfvStencil::numBoundaryFaces
int numBoundaryFaces() const
Definition: vcfvstencil.hh:1085
Ewoms::VcfvStencil::SubControlVolume::geometry
const ScvGeometry & geometry() const
Definition: vcfvstencil.hh:761
Ewoms::QuadrialteralQuadratureGeometry::corner
const GlobalPosition & corner(int cornerIdx) const
Return the position of the corner with a given index.
Definition: quadraturegeometries.hh:125
Ewoms::VcfvStencil::boundaryFace
const BoundaryFace & boundaryFace(int bfIdx) const
Definition: vcfvstencil.hh:1091
Ewoms::VcfvStencil::SubControlVolumeFace::ipLocal_
LocalPosition ipLocal_
integration point in local coords
Definition: vcfvstencil.hh:800
Ewoms::VcfvStencil::SubControlVolume::volume
Scalar volume() const
Definition: vcfvstencil.hh:755
Ewoms::VcfvStencil::SubControlVolumeFace::localPos
const LocalPosition & localPos() const
Definition: vcfvstencil.hh:791
Ewoms::VcfvStencil::SubControlVolume::center
const GlobalPosition center() const
Definition: vcfvstencil.hh:752
Ewoms::VcfvStencil::interiorFace
const SubControlVolumeFace & interiorFace(int faceIdx) const
Definition: vcfvstencil.hh:1088
Ewoms::VcfvStencil::ScvGeometry
Definition: vcfvstencil.hh:728
Ewoms::VcfvStencil::updateScvGeometry
void updateScvGeometry(const Element &element)
Definition: vcfvstencil.hh:1016