25#ifndef DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
26#define DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
29#include <dune/geometry/referenceelements.hh>
36template<
class Gr
idView,
int dim,
class ScvType,
class ScvfType>
40template <
class Gr
idView,
class ScvType,
class ScvfType>
44 using Scalar =
typename GridView::ctype;
45 using GlobalPosition =
typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
46 using ScvCornerStorage =
typename ScvType::Traits::CornerStorage;
47 using ScvfCornerStorage =
typename ScvfType::Traits::CornerStorage;
49 using Element =
typename GridView::template Codim<0>::Entity;
50 using Intersection =
typename GridView::Intersection;
54 static constexpr int maxPoints = 3;
59 : elementGeometry_(geometry), corners_(geometry.corners()),
60 p_({geometry.center(), geometry.corner(0), geometry.corner(1)}) {}
66 static const std::uint8_t map[2][2] =
72 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
73 p_[map[localScvIdx][1]]} };
79 return ScvfCornerStorage{{p_[0]}};
84 const typename Intersection::Geometry& geometry,
85 unsigned int indexInIntersection)
const
87 return ScvfCornerStorage{{geometry.corner(0)}};
91 GlobalPosition
normal(
const ScvfCornerStorage& scvfCorners,
92 const std::vector<unsigned int>& scvIndices)
const
94 auto normal = p_[2] - p_[1];
95 normal /= normal.two_norm();
100 Scalar
scvVolume(
const ScvCornerStorage& scvCorners)
const
102 return (scvCorners[1] - scvCorners[0]).two_norm();
106 Scalar
scvfArea(
const ScvfCornerStorage& scvfCorners)
const
114 std::array<GlobalPosition, maxPoints>
p_;
118template <
class Gr
idView,
class ScvType,
class ScvfType>
121 using Scalar =
typename GridView::ctype;
122 using GlobalPosition =
typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
123 using ScvCornerStorage =
typename ScvType::Traits::CornerStorage;
124 using ScvfCornerStorage =
typename ScvfType::Traits::CornerStorage;
126 using Element =
typename GridView::template Codim<0>::Entity;
127 using Intersection =
typename GridView::Intersection;
129 static constexpr auto dim = GridView::dimension;
130 static constexpr auto dimWorld = GridView::dimensionworld;
131 using ReferenceElements =
typename Dune::ReferenceElements<Scalar, dim>;
135 static constexpr int maxPoints = 9;
140 : elementGeometry_(geometry)
141 , corners_(geometry.corners())
143 const auto referenceElement = ReferenceElements::general(geometry.type());
146 p_[0] = geometry.center();
149 for (
int i = 0; i < corners_; ++i)
150 p_[i+1] = geometry.corner(i);
153 for (
int i = 0; i < referenceElement.size(1); ++i)
154 p_[i+corners_+1] = geometry.global(referenceElement.position(i, 1));
166 static const std::uint8_t vo = 1;
167 static const std::uint8_t fo = 4;
168 static const std::uint8_t map[3][4] =
170 {vo+0, fo+0, fo+1, 0},
171 {vo+1, fo+2, fo+0, 0},
172 {vo+2, fo+1, fo+2, 0}
175 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
176 p_[map[localScvIdx][1]],
177 p_[map[localScvIdx][2]],
178 p_[map[localScvIdx][3]]} };
183 static const std::uint8_t vo = 1;
184 static const std::uint8_t fo = 5;
185 static const std::uint8_t map[4][4] =
187 {vo+0, fo+2, fo+0, 0},
188 {vo+1, fo+1, fo+2, 0},
189 {vo+2, fo+0, fo+3, 0},
190 {vo+3, fo+3, fo+1, 0}
193 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
194 p_[map[localScvIdx][1]],
195 p_[map[localScvIdx][2]],
196 p_[map[localScvIdx][3]]} };
199 DUNE_THROW(Dune::NotImplemented,
"Box scv geometries for dim=" << dim
200 <<
" dimWorld=" << dimWorld
201 <<
" corners=" << corners_);
215 static const std::uint8_t fo = 4;
216 static const std::uint8_t map[3][2] =
223 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
224 p_[map[localScvfIdx][1]]} };
229 static const std::uint8_t fo = 5;
230 static const std::uint8_t map[4][2] =
238 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
239 p_[map[localScvfIdx][1]]} };
242 DUNE_THROW(Dune::NotImplemented,
"Box scvf geometries for dim=" << dim
243 <<
" dimWorld=" << dimWorld
244 <<
" corners=" << corners_);
250 const typename Intersection::Geometry& isGeom,
251 unsigned int indexInIntersection)
const
253 const auto referenceElement = ReferenceElements::general(elementGeometry_.type());
255 const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, indexInIntersection, dim);
256 if (indexInIntersection == 0)
257 return ScvfCornerStorage({p_[vIdxLocal+1], isGeom.center()});
258 else if (indexInIntersection == 1)
259 return ScvfCornerStorage({isGeom.center(), p_[vIdxLocal+1]});
261 DUNE_THROW(Dune::InvalidStateException,
"local index exceeds the number of corners of 2d intersections");
265 template <
int w = dimWorld>
266 typename std::enable_if<w == 3, GlobalPosition>::type
267 normal(
const ScvfCornerStorage& scvfCorners,
268 const std::vector<unsigned int>& scvIndices)
const
270 const auto v1 = elementGeometry_.corner(1) - elementGeometry_.corner(0);
271 const auto v2 = elementGeometry_.corner(2) - elementGeometry_.corner(0);
273 const auto t = scvfCorners[1] - scvfCorners[0];
275 normal /= normal.two_norm();
278 const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
279 const auto s = v*normal;
287 template <
int w = dimWorld>
288 typename std::enable_if<w == 2, GlobalPosition>::type
289 normal(
const ScvfCornerStorage& scvfCorners,
290 const std::vector<unsigned int>& scvIndices)
const
293 const auto t = scvfCorners[1] - scvfCorners[0];
294 GlobalPosition normal({-t[1], t[0]});
295 normal /= normal.two_norm();
298 const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
299 const auto s = v*normal;
307 template <
int w = dimWorld>
308 typename std::enable_if<w == 3, Scalar>::type
315 template <
int w = dimWorld>
316 typename std::enable_if<w == 2, Scalar>::type
328 return (p[1]-p[0]).two_norm();
334 std::array<GlobalPosition, maxPoints>
p_;
338template <
class Gr
idView,
class ScvType,
class ScvfType>
341 using Scalar =
typename GridView::ctype;
342 using GlobalPosition =
typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
343 using ScvCornerStorage =
typename ScvType::Traits::CornerStorage;
344 using ScvfCornerStorage =
typename ScvfType::Traits::CornerStorage;
346 using Element =
typename GridView::template Codim<0>::Entity;
347 using Intersection =
typename GridView::Intersection;
349 static constexpr auto dim = GridView::dimension;
350 static constexpr auto dimWorld = GridView::dimensionworld;
351 using ReferenceElements =
typename Dune::ReferenceElements<Scalar, dim>;
352 using FaceReferenceElements =
typename Dune::ReferenceElements<Scalar, dim-1>;
356 static constexpr int maxPoints = 27;
359 : elementGeometry_(geometry)
360 , corners_(geometry.corners())
362 const auto referenceElement = ReferenceElements::general(geometry.type());
365 p_[0] = geometry.center();
368 for (
int i = 0; i < corners_; ++i)
369 p_[i+1] = geometry.corner(i);
372 for (
int i = 0; i < referenceElement.size(dim-1); ++i)
373 p_[i+corners_+1] = geometry.global(referenceElement.position(i, dim-1));
376 for (
int i = 0; i < referenceElement.size(1); ++i)
377 p_[i+corners_+1+referenceElement.size(dim-1)] = geometry.global(referenceElement.position(i, 1));
389 static const std::uint8_t vo = 1;
390 static const std::uint8_t eo = 5;
391 static const std::uint8_t fo = 11;
392 static const std::uint8_t map[4][8] =
394 {vo+0, eo+0, eo+1, fo+0, eo+3, fo+1, fo+2, 0},
395 {vo+1, eo+2, eo+0, fo+0, eo+4, fo+3, fo+1, 0},
396 {vo+2, eo+1, eo+2, fo+0, eo+5, fo+2, fo+3, 0},
397 {vo+3, eo+3, eo+5, fo+2, eo+4, fo+1, fo+3, 0}
400 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
401 p_[map[localScvIdx][1]],
402 p_[map[localScvIdx][2]],
403 p_[map[localScvIdx][3]],
404 p_[map[localScvIdx][4]],
405 p_[map[localScvIdx][5]],
406 p_[map[localScvIdx][6]],
407 p_[map[localScvIdx][7]]} };
412 static const std::uint8_t vo = 1;
413 static const std::uint8_t eo = 7;
414 static const std::uint8_t fo = 16;
415 static const std::uint8_t map[6][8] =
417 {vo+0, eo+3, eo+4, fo+3, eo+0, fo+0, fo+1, 0},
418 {vo+1, eo+5, eo+3, fo+3, eo+1, fo+2, fo+0, 0},
419 {vo+2, eo+4, eo+5, fo+3, eo+2, fo+1, fo+2, 0},
420 {vo+3, eo+7, eo+6, fo+4, eo+0, fo+1, fo+0, 0},
421 {vo+4, eo+6, eo+8, fo+4, eo+1, fo+0, fo+2, 0},
422 {vo+5, eo+8, eo+7, fo+4, eo+2, fo+2, fo+1, 0}
425 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
426 p_[map[localScvIdx][1]],
427 p_[map[localScvIdx][2]],
428 p_[map[localScvIdx][3]],
429 p_[map[localScvIdx][4]],
430 p_[map[localScvIdx][5]],
431 p_[map[localScvIdx][6]],
432 p_[map[localScvIdx][7]]} };
437 static const std::uint8_t vo = 1;
438 static const std::uint8_t eo = 9;
439 static const std::uint8_t fo = 21;
440 static const std::uint8_t map[8][8] =
442 {vo+0, eo+6, eo+4, fo+4, eo+0, fo+2, fo+0, 0},
443 {vo+1, eo+5, eo+6, fo+4, eo+1, fo+1, fo+2, 0},
444 {vo+2, eo+4, eo+7, fo+4, eo+2, fo+0, fo+3, 0},
445 {vo+3, eo+7, eo+5, fo+4, eo+3, fo+3, fo+1, 0},
446 {vo+4, eo+8, eo+10, fo+5, eo+0, fo+0, fo+2, 0},
447 {vo+5, eo+10, eo+9, fo+5, eo+1, fo+2, fo+1, 0},
448 {vo+6, eo+11, eo+8, fo+5, eo+2, fo+3, fo+0, 0},
449 {vo+7, eo+9, eo+11, fo+5, eo+3, fo+1, fo+3, 0},
452 return ScvCornerStorage{ {p_[map[localScvIdx][0]],
453 p_[map[localScvIdx][1]],
454 p_[map[localScvIdx][2]],
455 p_[map[localScvIdx][3]],
456 p_[map[localScvIdx][4]],
457 p_[map[localScvIdx][5]],
458 p_[map[localScvIdx][6]],
459 p_[map[localScvIdx][7]]} };
462 DUNE_THROW(Dune::NotImplemented,
"Box scv geometries for dim=" << dim
463 <<
" dimWorld=" << dimWorld
464 <<
" corners=" << corners_);
477 static const std::uint8_t eo = 5;
478 static const std::uint8_t fo = 11;
479 static const std::uint8_t map[6][4] =
481 {eo+0, fo+0, fo+1, 0},
482 {fo+0, eo+1, 0, fo+2},
483 {eo+2, fo+0, fo+3, 0},
484 {fo+2, eo+3, 0, fo+1},
485 {fo+3, 0, eo+4, fo+1},
486 {eo+5, fo+2, fo+3, 0}
489 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
490 p_[map[localScvfIdx][1]],
491 p_[map[localScvfIdx][2]],
492 p_[map[localScvfIdx][3]]} };
497 static const std::uint8_t eo = 7;
498 static const std::uint8_t fo = 16;
499 static const std::uint8_t map[9][4] =
501 {eo+0, fo+0, fo+1, 0},
502 {eo+1, fo+2, fo+0, 0},
503 {eo+2, fo+1, fo+2, 0},
504 {eo+3, fo+0, fo+3, 0},
505 {eo+4, fo+3, fo+1, 0},
506 {eo+5, fo+2, fo+3, 0},
507 {eo+6, fo+4, fo+0, 0},
508 {eo+7, fo+1, fo+4, 0},
509 {eo+8, fo+4, fo+2, 0}
512 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
513 p_[map[localScvfIdx][1]],
514 p_[map[localScvfIdx][2]],
515 p_[map[localScvfIdx][3]]} };
520 static const std::uint8_t eo = 9;
521 static const std::uint8_t fo = 21;
522 static const std::uint8_t map[12][4] =
524 {fo+0, eo+0, 0, fo+2},
525 {fo+1, 0, eo+1, fo+2},
526 {fo+3, eo+2, 0, fo+0},
527 {eo+3, fo+3, fo+1, 0},
528 {fo+4, eo+4, 0, fo+0},
529 {eo+5, fo+4, fo+1, 0},
530 {eo+6, fo+4, fo+2, 0},
531 {fo+4, eo+7, 0, fo+3},
532 { 0, fo+0, fo+5, eo+8},
533 {eo+9, fo+1, fo+5, 0},
534 {eo+10, fo+2, fo+5, 0},
535 {eo+11, fo+5, fo+3, 0}
538 return ScvfCornerStorage{ {p_[map[localScvfIdx][0]],
539 p_[map[localScvfIdx][1]],
540 p_[map[localScvfIdx][2]],
541 p_[map[localScvfIdx][3]]} };
544 DUNE_THROW(Dune::NotImplemented,
"Box scv geometries for dim=" << dim
545 <<
" dimWorld=" << dimWorld
546 <<
" corners=" << corners_);
552 const typename Intersection::Geometry& geometry,
553 unsigned int indexInIntersection)
const
555 const auto referenceElement = ReferenceElements::general(elementGeometry_.type());
556 const auto faceRefElem = FaceReferenceElements::general(geometry.type());
558 GlobalPosition pi[9];
559 auto corners = geometry.corners();
562 pi[0] = geometry.center();
565 const auto idxInInside = is.indexInInside();
566 for (
int i = 0; i < corners; ++i)
568 const auto vIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim);
569 pi[i+1] = elementGeometry_.corner(vIdxLocal);
573 for (
int i = 0; i < faceRefElem.size(1); ++i)
575 const auto edgeIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim-1);
576 pi[i+corners+1] = p_[edgeIdxLocal+corners_+1];
585 static const std::uint8_t vo = 1;
586 static const std::uint8_t eo = 4;
587 static const std::uint8_t map[3][4] =
589 {vo+0, eo+0, eo+1, 0},
590 {vo+1, eo+2, eo+0, 0},
591 {vo+2, eo+1, eo+2, 0}
594 return ScvfCornerStorage{ {pi[map[indexInIntersection][0]],
595 pi[map[indexInIntersection][1]],
596 pi[map[indexInIntersection][2]],
597 pi[map[indexInIntersection][3]]} };
602 static const std::uint8_t vo = 1;
603 static const std::uint8_t eo = 5;
604 static const std::uint8_t map[4][4] =
606 {vo+0, eo+2, eo+0, 0},
607 {vo+1, eo+1, eo+2, 0},
608 {vo+2, eo+0, eo+3, 0},
609 {vo+3, eo+3, eo+1, 0}
612 return ScvfCornerStorage{ {pi[map[indexInIntersection][0]],
613 pi[map[indexInIntersection][1]],
614 pi[map[indexInIntersection][2]],
615 pi[map[indexInIntersection][3]]} };
618 DUNE_THROW(Dune::NotImplemented,
"Box scvf boundary geometries for dim=" << dim
619 <<
" dimWorld=" << dimWorld
620 <<
" corners=" << corners);
625 GlobalPosition
normal(
const ScvfCornerStorage& p,
626 const std::vector<unsigned int>& scvIndices)
const
629 normal /= normal.two_norm();
631 const auto v = elementGeometry_.corner(scvIndices[1]) - elementGeometry_.corner(scvIndices[0]);
632 const auto s = v*normal;
643 const auto v = p[7]-p[0];
659 std::array<GlobalPosition, maxPoints>
p_;
Define some often used mathematical functions.
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:631
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:660
Create sub control volumes and sub control volume face geometries.
Definition: boxgeometryhelper.hh:37
Scalar scvfArea(const ScvfCornerStorage &scvfCorners) const
get scvf area
Definition: boxgeometryhelper.hh:106
BoxGeometryHelper(const typename Element::Geometry &geometry)
Definition: boxgeometryhelper.hh:58
ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
Create a vector with the scv corners.
Definition: boxgeometryhelper.hh:63
std::array< GlobalPosition, maxPoints > p_
Definition: boxgeometryhelper.hh:114
GlobalPosition normal(const ScvfCornerStorage &scvfCorners, const std::vector< unsigned int > &scvIndices) const
get scvf normal vector
Definition: boxgeometryhelper.hh:91
ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
Create a vector with the corners of sub control volume faces.
Definition: boxgeometryhelper.hh:77
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:83
const Element::Geometry & elementGeometry_
Reference to the element geometry.
Definition: boxgeometryhelper.hh:112
Scalar scvVolume(const ScvCornerStorage &scvCorners) const
get scv volume
Definition: boxgeometryhelper.hh:100
std::size_t corners_
Definition: boxgeometryhelper.hh:113
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:289
Scalar scvfArea(const ScvfCornerStorage &p) const
get scvf area
Definition: boxgeometryhelper.hh:326
std::array< GlobalPosition, maxPoints > p_
Definition: boxgeometryhelper.hh:334
ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
Create a vector with the corners of sub control volume faces.
Definition: boxgeometryhelper.hh:207
const Element::Geometry & elementGeometry_
Reference to the element geometry.
Definition: boxgeometryhelper.hh:332
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:249
BoxGeometryHelper(const typename Element::Geometry &geometry)
Definition: boxgeometryhelper.hh:139
std::enable_if< w==3, Scalar >::type scvVolume(const ScvCornerStorage &p) const
get scv volume for dim == 2, dimworld == 3
Definition: boxgeometryhelper.hh:309
std::size_t corners_
Definition: boxgeometryhelper.hh:333
ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
Create a vector with the scv corners.
Definition: boxgeometryhelper.hh:158
std::enable_if< w==2, Scalar >::type scvVolume(const ScvCornerStorage &p) const
get scv volume for dim == 2, dimworld == 2
Definition: boxgeometryhelper.hh:317
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:267
Scalar scvfArea(const ScvfCornerStorage &p) const
get scvf area
Definition: boxgeometryhelper.hh:650
const Element::Geometry & elementGeometry_
Reference to the element geometry.
Definition: boxgeometryhelper.hh:657
GlobalPosition normal(const ScvfCornerStorage &p, const std::vector< unsigned int > &scvIndices) const
get scvf normal vector
Definition: boxgeometryhelper.hh:625
std::array< GlobalPosition, maxPoints > p_
Definition: boxgeometryhelper.hh:659
BoxGeometryHelper(const typename Element::Geometry &geometry)
Definition: boxgeometryhelper.hh:358
Scalar scvVolume(const ScvCornerStorage &p) const
get scv volume
Definition: boxgeometryhelper.hh:640
ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
Create a vector with the scv corners.
Definition: boxgeometryhelper.hh:381
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:551
std::size_t corners_
Definition: boxgeometryhelper.hh:658
ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
Create a vector with the scvf corners.
Definition: boxgeometryhelper.hh:469