25#ifndef DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
26#define DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
35template<
class Gr
idView,
int dim,
class ScvType,
class ScvfType>
39template <
class Gr
idView,
class ScvType,
class ScvfType>
43 using Scalar =
typename GridView::ctype;
44 using GlobalPosition =
typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
45 using ScvCornerStorage =
typename ScvType::Traits::CornerStorage;
46 using ScvfCornerStorage =
typename ScvfType::Traits::CornerStorage;
48 using Element =
typename GridView::template Codim<0>::Entity;
49 using Intersection =
typename GridView::Intersection;
53 static constexpr int maxPoints = 3;
58 : elementGeometry_(geometry), corners_(geometry.corners()),
59 p_({geometry.center(), geometry.corner(0), geometry.corner(1)}) {}
65 static const std::uint8_t map[2][2] =
71 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
72 p_[map[localScvIdx][1]]} };
78 return ScvfCornerStorage{{p_[0]}};
83 const typename Intersection::Geometry& geometry,
84 unsigned int indexInIntersection)
const
86 return ScvfCornerStorage{{geometry.corner(0)}};
90 GlobalPosition
normal(
const ScvfCornerStorage& scvfCorners,
91 const std::vector<unsigned int>& scvIndices)
const
93 auto normal = p_[2] - p_[1];
99 Scalar
scvVolume(
const ScvCornerStorage& scvCorners)
const
101 return (scvCorners[1] - scvCorners[0]).two_norm();
105 Scalar
scvfArea(
const ScvfCornerStorage& scvfCorners)
const
113 std::array<GlobalPosition, maxPoints>
p_;
117template <
class Gr
idView,
class ScvType,
class ScvfType>
120 using Scalar =
typename GridView::ctype;
121 using GlobalPosition =
typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
122 using ScvCornerStorage =
typename ScvType::Traits::CornerStorage;
123 using ScvfCornerStorage =
typename ScvfType::Traits::CornerStorage;
125 using Element =
typename GridView::template Codim<0>::Entity;
126 using Intersection =
typename GridView::Intersection;
128 static constexpr auto dim = GridView::dimension;
129 static constexpr auto dimWorld = GridView::dimensionworld;
133 static constexpr int maxPoints = 9;
138 : elementGeometry_(geometry)
139 , corners_(geometry.corners())
141 const auto refElement = referenceElement(geometry);
144 p_[0] = geometry.center();
147 for (
int i = 0; i < corners_; ++i)
148 p_[i+1] = geometry.corner(i);
151 for (
int i = 0; i < refElement.size(1); ++i)
152 p_[i+corners_+1] = geometry.global(refElement.position(i, 1));
164 static const std::uint8_t vo = 1;
165 static const std::uint8_t fo = 4;
166 static const std::uint8_t map[3][4] =
168 {vo+0, fo+0, fo+1, 0},
169 {vo+1, fo+2, fo+0, 0},
170 {vo+2, fo+1, fo+2, 0}
173 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
174 p_[map[localScvIdx][1]],
175 p_[map[localScvIdx][2]],
176 p_[map[localScvIdx][3]]} };
181 static const std::uint8_t vo = 1;
182 static const std::uint8_t fo = 5;
183 static const std::uint8_t map[4][4] =
185 {vo+0, fo+2, fo+0, 0},
186 {vo+1, fo+1, fo+2, 0},
187 {vo+2, fo+0, fo+3, 0},
188 {vo+3, fo+3, fo+1, 0}
191 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
192 p_[map[localScvIdx][1]],
193 p_[map[localScvIdx][2]],
194 p_[map[localScvIdx][3]]} };
197 DUNE_THROW(Dune::NotImplemented,
"Box scv geometries for dim=" << dim
198 <<
" dimWorld=" << dimWorld
199 <<
" corners=" << corners_);
213 static const std::uint8_t fo = 4;
214 static const std::uint8_t map[3][2] =
221 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
222 p_[map[localScvfIdx][1]]} };
227 static const std::uint8_t fo = 5;
228 static const std::uint8_t map[4][2] =
236 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
237 p_[map[localScvfIdx][1]]} };
240 DUNE_THROW(Dune::NotImplemented,
"Box scvf geometries for dim=" << dim
241 <<
" dimWorld=" << dimWorld
242 <<
" corners=" << corners_);
248 const typename Intersection::Geometry& isGeom,
249 unsigned int indexInIntersection)
const
251 const auto refElement = referenceElement(elementGeometry_);
253 const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, indexInIntersection, dim);
254 if (indexInIntersection == 0)
255 return ScvfCornerStorage({p_[vIdxLocal+1], isGeom.center()});
256 else if (indexInIntersection == 1)
257 return ScvfCornerStorage({isGeom.center(), p_[vIdxLocal+1]});
259 DUNE_THROW(Dune::InvalidStateException,
"local index exceeds the number of corners of 2d intersections");
263 template <
int w = dimWorld>
264 typename std::enable_if<w == 3, GlobalPosition>::type
265 normal(
const ScvfCornerStorage& scvfCorners,
266 const std::vector<unsigned int>& scvIndices)
const
268 const auto v1 = elementGeometry_.corner(1) - elementGeometry_.corner(0);
269 const auto v2 = elementGeometry_.corner(2) - elementGeometry_.corner(0);
271 const auto t = scvfCorners[1] - scvfCorners[0];
276 const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
285 template <
int w = dimWorld>
286 typename std::enable_if<w == 2, GlobalPosition>::type
287 normal(
const ScvfCornerStorage& scvfCorners,
288 const std::vector<unsigned int>& scvIndices)
const
291 const auto t = scvfCorners[1] - scvfCorners[0];
292 GlobalPosition
normal({-t[1], t[0]});
296 const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
305 template <
int w = dimWorld>
306 typename std::enable_if<w == 3, Scalar>::type
313 template <
int w = dimWorld>
314 typename std::enable_if<w == 2, Scalar>::type
326 return (p[1]-p[0]).two_norm();
332 std::array<GlobalPosition, maxPoints>
p_;
336template <
class Gr
idView,
class ScvType,
class ScvfType>
339 using Scalar =
typename GridView::ctype;
340 using GlobalPosition =
typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
341 using ScvCornerStorage =
typename ScvType::Traits::CornerStorage;
342 using ScvfCornerStorage =
typename ScvfType::Traits::CornerStorage;
344 using Element =
typename GridView::template Codim<0>::Entity;
345 using Intersection =
typename GridView::Intersection;
347 static constexpr auto dim = GridView::dimension;
348 static constexpr auto dimWorld = GridView::dimensionworld;
352 static constexpr int maxPoints = 27;
355 : elementGeometry_(geometry)
356 , corners_(geometry.corners())
358 const auto refElement = referenceElement(geometry);
361 p_[0] = geometry.center();
364 for (
int i = 0; i < corners_; ++i)
365 p_[i+1] = geometry.corner(i);
368 for (
int i = 0; i < refElement.size(dim-1); ++i)
369 p_[i+corners_+1] = geometry.global(refElement.position(i, dim-1));
372 for (
int i = 0; i < refElement.size(1); ++i)
373 p_[i+corners_+1+refElement.size(dim-1)] = geometry.global(refElement.position(i, 1));
385 static const std::uint8_t vo = 1;
386 static const std::uint8_t eo = 5;
387 static const std::uint8_t fo = 11;
388 static const std::uint8_t map[4][8] =
390 {vo+0, eo+0, eo+1, fo+0, eo+3, fo+1, fo+2, 0},
391 {vo+1, eo+2, eo+0, fo+0, eo+4, fo+3, fo+1, 0},
392 {vo+2, eo+1, eo+2, fo+0, eo+5, fo+2, fo+3, 0},
393 {vo+3, eo+3, eo+5, fo+2, eo+4, fo+1, fo+3, 0}
396 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
397 p_[map[localScvIdx][1]],
398 p_[map[localScvIdx][2]],
399 p_[map[localScvIdx][3]],
400 p_[map[localScvIdx][4]],
401 p_[map[localScvIdx][5]],
402 p_[map[localScvIdx][6]],
403 p_[map[localScvIdx][7]]} };
408 static const std::uint8_t vo = 1;
409 static const std::uint8_t eo = 7;
410 static const std::uint8_t fo = 16;
411 static const std::uint8_t map[6][8] =
413 {vo+0, eo+3, eo+4, fo+3, eo+0, fo+0, fo+1, 0},
414 {vo+1, eo+5, eo+3, fo+3, eo+1, fo+2, fo+0, 0},
415 {vo+2, eo+4, eo+5, fo+3, eo+2, fo+1, fo+2, 0},
416 {vo+3, eo+7, eo+6, fo+4, eo+0, fo+1, fo+0, 0},
417 {vo+4, eo+6, eo+8, fo+4, eo+1, fo+0, fo+2, 0},
418 {vo+5, eo+8, eo+7, fo+4, eo+2, fo+2, fo+1, 0}
421 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
422 p_[map[localScvIdx][1]],
423 p_[map[localScvIdx][2]],
424 p_[map[localScvIdx][3]],
425 p_[map[localScvIdx][4]],
426 p_[map[localScvIdx][5]],
427 p_[map[localScvIdx][6]],
428 p_[map[localScvIdx][7]]} };
433 static const std::uint8_t vo = 1;
434 static const std::uint8_t eo = 9;
435 static const std::uint8_t fo = 21;
436 static const std::uint8_t map[8][8] =
438 {vo+0, eo+6, eo+4, fo+4, eo+0, fo+2, fo+0, 0},
439 {vo+1, eo+5, eo+6, fo+4, eo+1, fo+1, fo+2, 0},
440 {vo+2, eo+4, eo+7, fo+4, eo+2, fo+0, fo+3, 0},
441 {vo+3, eo+7, eo+5, fo+4, eo+3, fo+3, fo+1, 0},
442 {vo+4, eo+8, eo+10, fo+5, eo+0, fo+0, fo+2, 0},
443 {vo+5, eo+10, eo+9, fo+5, eo+1, fo+2, fo+1, 0},
444 {vo+6, eo+11, eo+8, fo+5, eo+2, fo+3, fo+0, 0},
445 {vo+7, eo+9, eo+11, fo+5, eo+3, fo+1, fo+3, 0},
448 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
449 p_[map[localScvIdx][1]],
450 p_[map[localScvIdx][2]],
451 p_[map[localScvIdx][3]],
452 p_[map[localScvIdx][4]],
453 p_[map[localScvIdx][5]],
454 p_[map[localScvIdx][6]],
455 p_[map[localScvIdx][7]]} };
458 DUNE_THROW(Dune::NotImplemented,
"Box scv geometries for dim=" << dim
459 <<
" dimWorld=" << dimWorld
460 <<
" corners=" << corners_);
473 static const std::uint8_t eo = 5;
474 static const std::uint8_t fo = 11;
475 static const std::uint8_t map[6][4] =
477 {eo+0, fo+0, fo+1, 0},
478 {fo+0, eo+1, 0, fo+2},
479 {eo+2, fo+0, fo+3, 0},
480 {fo+2, eo+3, 0, fo+1},
481 {fo+3, 0, eo+4, fo+1},
482 {eo+5, fo+2, fo+3, 0}
485 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
486 p_[map[localScvfIdx][1]],
487 p_[map[localScvfIdx][2]],
488 p_[map[localScvfIdx][3]]} };
493 static const std::uint8_t eo = 7;
494 static const std::uint8_t fo = 16;
495 static const std::uint8_t map[9][4] =
497 {eo+0, fo+0, fo+1, 0},
498 {eo+1, fo+2, fo+0, 0},
499 {eo+2, fo+1, fo+2, 0},
500 {eo+3, fo+0, fo+3, 0},
501 {eo+4, fo+3, fo+1, 0},
502 {eo+5, fo+2, fo+3, 0},
503 {eo+6, fo+4, fo+0, 0},
504 {eo+7, fo+1, fo+4, 0},
505 {eo+8, fo+4, fo+2, 0}
508 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
509 p_[map[localScvfIdx][1]],
510 p_[map[localScvfIdx][2]],
511 p_[map[localScvfIdx][3]]} };
516 static const std::uint8_t eo = 9;
517 static const std::uint8_t fo = 21;
518 static const std::uint8_t map[12][4] =
520 {fo+0, eo+0, 0, fo+2},
521 {fo+1, 0, eo+1, fo+2},
522 {fo+3, eo+2, 0, fo+0},
523 {eo+3, fo+3, fo+1, 0},
524 {fo+4, eo+4, 0, fo+0},
525 {eo+5, fo+4, fo+1, 0},
526 {eo+6, fo+4, fo+2, 0},
527 {fo+4, eo+7, 0, fo+3},
528 { 0, fo+0, fo+5, eo+8},
529 {eo+9, fo+1, fo+5, 0},
530 {eo+10, fo+2, fo+5, 0},
531 {eo+11, fo+5, fo+3, 0}
534 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
535 p_[map[localScvfIdx][1]],
536 p_[map[localScvfIdx][2]],
537 p_[map[localScvfIdx][3]]} };
540 DUNE_THROW(Dune::NotImplemented,
"Box scv geometries for dim=" << dim
541 <<
" dimWorld=" << dimWorld
542 <<
" corners=" << corners_);
548 const typename Intersection::Geometry& geometry,
549 unsigned int indexInIntersection)
const
551 const auto refElement = referenceElement(elementGeometry_);
552 const auto faceRefElem = referenceElement(geometry);
554 GlobalPosition pi[9];
555 auto corners = geometry.corners();
558 pi[0] = geometry.center();
561 const auto idxInInside = is.indexInInside();
562 for (
int i = 0; i < corners; ++i)
564 const auto vIdxLocal = refElement.subEntity(idxInInside, 1, i, dim);
565 pi[i+1] = elementGeometry_.corner(vIdxLocal);
569 for (
int i = 0; i < faceRefElem.size(1); ++i)
571 const auto edgeIdxLocal = refElement.subEntity(idxInInside, 1, i, dim-1);
572 pi[i+corners+1] = p_[edgeIdxLocal+corners_+1];
581 static const std::uint8_t vo = 1;
582 static const std::uint8_t eo = 4;
583 static const std::uint8_t map[3][4] =
585 {vo+0, eo+0, eo+1, 0},
586 {vo+1, eo+2, eo+0, 0},
587 {vo+2, eo+1, eo+2, 0}
590 return ScvfCornerStorage{ {pi[map[indexInIntersection][0]],
591 pi[map[indexInIntersection][1]],
592 pi[map[indexInIntersection][2]],
593 pi[map[indexInIntersection][3]]} };
598 static const std::uint8_t vo = 1;
599 static const std::uint8_t eo = 5;
600 static const std::uint8_t map[4][4] =
602 {vo+0, eo+2, eo+0, 0},
603 {vo+1, eo+1, eo+2, 0},
604 {vo+2, eo+0, eo+3, 0},
605 {vo+3, eo+3, eo+1, 0}
608 return ScvfCornerStorage{ {pi[map[indexInIntersection][0]],
609 pi[map[indexInIntersection][1]],
610 pi[map[indexInIntersection][2]],
611 pi[map[indexInIntersection][3]]} };
614 DUNE_THROW(Dune::NotImplemented,
"Box scvf boundary geometries for dim=" << dim
615 <<
" dimWorld=" << dimWorld
616 <<
" corners=" << corners);
621 GlobalPosition
normal(
const ScvfCornerStorage& p,
622 const std::vector<unsigned int>& scvIndices)
const
627 const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
639 const auto v = p[7]-p[0];
655 std::array<GlobalPosition, maxPoints>
p_;
Define some often used mathematical functions.
Vector normal(const Vector &v)
Create a vector normal to the given one (v is expected to be non-zero)
Definition: normal.hh:36
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:654
Scalar tripleProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2, const Dune::FieldVector< Scalar, 3 > &vec3)
Triple product of three vectors in three-dimensional Euclidean space retuning scalar.
Definition: math.hh:683
Create sub control volumes and sub control volume face geometries.
Definition: boxgeometryhelper.hh:36
Scalar scvfArea(const ScvfCornerStorage &scvfCorners) const
get scvf area
Definition: boxgeometryhelper.hh:105
BoxGeometryHelper(const typename Element::Geometry &geometry)
Definition: boxgeometryhelper.hh:57
ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
Create a vector with the scv corners.
Definition: boxgeometryhelper.hh:62
std::array< GlobalPosition, maxPoints > p_
Definition: boxgeometryhelper.hh:113
GlobalPosition normal(const ScvfCornerStorage &scvfCorners, const std::vector< unsigned int > &scvIndices) const
get scvf normal vector
Definition: boxgeometryhelper.hh:90
ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
Create a vector with the corners of sub control volume faces.
Definition: boxgeometryhelper.hh:76
ScvfCornerStorage getBoundaryScvfCorners(const Intersection &is, const typename Intersection::Geometry &geometry, unsigned int indexInIntersection) const
Create the sub control volume face geometries on the boundary.
Definition: boxgeometryhelper.hh:82
const Element::Geometry & elementGeometry_
Reference to the element geometry.
Definition: boxgeometryhelper.hh:111
Scalar scvVolume(const ScvCornerStorage &scvCorners) const
get scv volume
Definition: boxgeometryhelper.hh:99
std::size_t corners_
Definition: boxgeometryhelper.hh:112
std::enable_if< w==2, GlobalPosition >::type normal(const ScvfCornerStorage &scvfCorners, const std::vector< unsigned int > &scvIndices) const
get scvf normal vector for dim == 2, dimworld == 2
Definition: boxgeometryhelper.hh:287
Scalar scvfArea(const ScvfCornerStorage &p) const
get scvf area
Definition: boxgeometryhelper.hh:324
std::array< GlobalPosition, maxPoints > p_
Definition: boxgeometryhelper.hh:332
ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
Create a vector with the corners of sub control volume faces.
Definition: boxgeometryhelper.hh:205
const Element::Geometry & elementGeometry_
Reference to the element geometry.
Definition: boxgeometryhelper.hh:330
ScvfCornerStorage getBoundaryScvfCorners(const Intersection &is, const typename Intersection::Geometry &isGeom, unsigned int indexInIntersection) const
Create the sub control volume face geometries on the boundary.
Definition: boxgeometryhelper.hh:247
BoxGeometryHelper(const typename Element::Geometry &geometry)
Definition: boxgeometryhelper.hh:137
std::enable_if< w==3, Scalar >::type scvVolume(const ScvCornerStorage &p) const
get scv volume for dim == 2, dimworld == 3
Definition: boxgeometryhelper.hh:307
std::size_t corners_
Definition: boxgeometryhelper.hh:331
ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
Create a vector with the scv corners.
Definition: boxgeometryhelper.hh:156
std::enable_if< w==2, Scalar >::type scvVolume(const ScvCornerStorage &p) const
get scv volume for dim == 2, dimworld == 2
Definition: boxgeometryhelper.hh:315
std::enable_if< w==3, GlobalPosition >::type normal(const ScvfCornerStorage &scvfCorners, const std::vector< unsigned int > &scvIndices) const
get scvf normal vector for dim == 2, dimworld == 3
Definition: boxgeometryhelper.hh:265
Scalar scvfArea(const ScvfCornerStorage &p) const
get scvf area
Definition: boxgeometryhelper.hh:646
const Element::Geometry & elementGeometry_
Reference to the element geometry.
Definition: boxgeometryhelper.hh:653
GlobalPosition normal(const ScvfCornerStorage &p, const std::vector< unsigned int > &scvIndices) const
get scvf normal vector
Definition: boxgeometryhelper.hh:621
std::array< GlobalPosition, maxPoints > p_
Definition: boxgeometryhelper.hh:655
BoxGeometryHelper(const typename Element::Geometry &geometry)
Definition: boxgeometryhelper.hh:354
Scalar scvVolume(const ScvCornerStorage &p) const
get scv volume
Definition: boxgeometryhelper.hh:636
ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
Create a vector with the scv corners.
Definition: boxgeometryhelper.hh:377
ScvfCornerStorage getBoundaryScvfCorners(const Intersection &is, const typename Intersection::Geometry &geometry, unsigned int indexInIntersection) const
Create the sub control volume face geometries on the boundary.
Definition: boxgeometryhelper.hh:547
std::size_t corners_
Definition: boxgeometryhelper.hh:654
ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
Create a vector with the scvf corners.
Definition: boxgeometryhelper.hh:465