3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
geometryhelper.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program 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 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_POROUSMEDIUMFLOW_BOXDFM_GEOMETRY_HELPER_HH
26#define DUMUX_POROUSMEDIUMFLOW_BOXDFM_GEOMETRY_HELPER_HH
27
29
30namespace Dumux
31{
32
34template<class GridView, int dim, class ScvType, class ScvfType>
36
38template <class GridView, class ScvType, class ScvfType>
39class BoxDfmGeometryHelper<GridView, 2, ScvType, ScvfType> : public BoxGeometryHelper<GridView, 2, ScvType, ScvfType>
40{
42
43 using Intersection = typename GridView::Intersection;
44 using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
45
46 static constexpr auto dim = GridView::dimension;
47 using Scalar = typename GridView::ctype;
48
49public:
50
52 using ParentType::ParentType;
53
56 ScvfCornerStorage getFractureScvfCorners(const Intersection& is,
57 const typename Intersection::Geometry& isGeom,
58 unsigned int idxOnIntersection = 0) const
59 { return ScvfCornerStorage({isGeom.center()}); }
60
61
64 typename ScvfType::Traits::GlobalPosition
65 fractureNormal(const ScvfCornerStorage& p,
66 const Intersection& is,
67 unsigned int edgeIndexInIntersection = 0) const
68 {
69 const auto refElement = referenceElement(this->elementGeometry_);
70 const auto vIdxLocal0 = refElement.subEntity(is.indexInInside(), 1, 0, dim);
71 const auto vIdxLocal1 = refElement.subEntity(is.indexInInside(), 1, 1, dim);
72 auto n = this->elementGeometry_.corner(vIdxLocal1) - this->elementGeometry_.corner(vIdxLocal0);
73 n /= n.two_norm();
74 return n;
75 }
76};
77
79template <class GridView, class ScvType, class ScvfType>
80class BoxDfmGeometryHelper<GridView, 3, ScvType, ScvfType> : public BoxGeometryHelper<GridView, 3, ScvType, ScvfType>
81{
83
84 using Intersection = typename GridView::Intersection;
85 using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
86
87 static constexpr auto dim = GridView::dimension;
88 static constexpr auto dimWorld = GridView::dimensionworld;
89 using Scalar = typename GridView::ctype;
90
91public:
92
94 using ParentType::ParentType;
95
97 ScvfCornerStorage getFractureScvfCorners(const Intersection& is,
98 const typename Intersection::Geometry& isGeom,
99 unsigned int edgeIndexInIntersection) const
100 {
101 const auto refElement = referenceElement(this->elementGeometry_);
102 const auto faceRefElem = referenceElement(isGeom);
103
104 // create point vector for this geometry
105 typename ScvfType::Traits::GlobalPosition pi[9];
106
107 // the facet center
108 pi[0] = isGeom.center();
109
110 // facet edge midpoints
111 const auto idxInInside = is.indexInInside();
112 for (int i = 0; i < faceRefElem.size(1); ++i)
113 {
114 const auto edgeIdxLocal = refElement.subEntity(idxInInside, 1, i, dim-1);
115 pi[i+1] = this->p_[edgeIdxLocal+this->corners_+1];
116 }
117
118 // proceed according to number of corners
119 const auto corners = isGeom.corners();
120 switch (corners)
121 {
122 case 3: // triangle
123 {
125 static const std::uint8_t fo = 1;
126 static const std::uint8_t map[3][2] =
127 {
128 {fo+0, 0},
129 {0, fo+1},
130 {fo+2, 0}
131 };
132
133 return ScvfCornerStorage{ {pi[map[edgeIndexInIntersection][0]],
134 pi[map[edgeIndexInIntersection][1]]} };
135 }
136 case 4: // quadrilateral
137 {
139 static const std::uint8_t fo = 1;
140 static const std::uint8_t map[4][2] =
141 {
142 {0, fo+0},
143 {fo+1, 0},
144 {fo+2, 0},
145 {0, fo+3}
146 };
147
148 return ScvfCornerStorage{ {pi[map[edgeIndexInIntersection][0]],
149 pi[map[edgeIndexInIntersection][1]]} };
150 }
151 default:
152 DUNE_THROW(Dune::NotImplemented, "Box fracture scvf geometries for dim=" << dim
153 << " dimWorld=" << dimWorld
154 << " corners=" << corners);
155 }
156 }
157
159 typename ScvfType::Traits::GlobalPosition
160 fractureNormal(const ScvfCornerStorage& p,
161 const Intersection& is,
162 unsigned int edgeIndexInIntersection) const
163 {
164 const auto refElement = referenceElement(this->elementGeometry_);
165
166 // first get the intersection corners (maximum "4" is for quadrilateral face)
167 typename ScvfType::Traits::GlobalPosition c[4];
168
169 const auto corners = refElement.size(is.indexInInside(), 1, dim);
170 for (int i = 0; i < corners; ++i)
171 {
172 const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, i, dim);
173 c[i] = this->elementGeometry_.corner(vIdxLocal);
174 }
175
176 // compute edge vector depending on number of corners
177 const auto gridEdge = [&] ()
178 {
179 // triangles
180 if (corners == 3)
181 {
182 if (edgeIndexInIntersection == 0) return c[1]-c[0];
183 else if (edgeIndexInIntersection == 1) return c[2]-c[0];
184 else if (edgeIndexInIntersection == 2) return c[2]-c[1];
185 else DUNE_THROW(Dune::InvalidStateException, "Invalid edge index");
186 }
187 else if (corners == 4)
188 {
189 if (edgeIndexInIntersection == 0) return c[2]-c[0];
190 else if (edgeIndexInIntersection == 1) return c[3]-c[1];
191 else if (edgeIndexInIntersection == 2) return c[1]-c[0];
192 else if (edgeIndexInIntersection == 3) return c[3]-c[2];
193 else DUNE_THROW(Dune::InvalidStateException, "Invalid edge index");
194 }
195 else
196 DUNE_THROW(Dune::InvalidStateException, "Invalid face geometry");
197 } ();
198
199 // compute lower edge of the scvf
200 assert(p.size() == 2);
201 const auto scvfEdge = p[1]-p[0];
202
203 // compute scvf normal via 2 cross products
204 const auto faceN = crossProduct(gridEdge, scvfEdge);
205 auto n = crossProduct(scvfEdge, faceN);
206 n /= n.two_norm();
207 return n;
208 }
209};
210
211} // end namespace Dumux
212
213#endif
Helper class constructing the dual grid finite volume geometries for the box discretizazion method.
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:640
Definition: adapt.hh:29
Create sub control volumes and sub control volume face geometries.
Definition: boxgeometryhelper.hh:36
A class to create sub control volume and sub control volume face geometries per element.
Definition: boxgeometryhelper.hh:119
A class to create sub control volume and sub control volume face geometries per element.
Definition: boxgeometryhelper.hh:338
Create sub control volumes and sub control volume face geometries.
Definition: geometryhelper.hh:35
ScvfType::Traits::GlobalPosition fractureNormal(const ScvfCornerStorage &p, const Intersection &is, unsigned int edgeIndexInIntersection=0) const
Definition: geometryhelper.hh:65
ScvfCornerStorage getFractureScvfCorners(const Intersection &is, const typename Intersection::Geometry &isGeom, unsigned int idxOnIntersection=0) const
Definition: geometryhelper.hh:56
ScvfCornerStorage getFractureScvfCorners(const Intersection &is, const typename Intersection::Geometry &isGeom, unsigned int edgeIndexInIntersection) const
Create the sub control volume face geometries on an intersection marked as fracture.
Definition: geometryhelper.hh:97
ScvfType::Traits::GlobalPosition fractureNormal(const ScvfCornerStorage &p, const Intersection &is, unsigned int edgeIndexInIntersection) const
get fracture scvf normal vector
Definition: geometryhelper.hh:160