3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
geometry/refinementquadraturerule.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 *****************************************************************************/
26#ifndef DUMUX_GEOMETRY_REFINEMENT_QUADRATURERULE_HH
27#define DUMUX_GEOMETRY_REFINEMENT_QUADRATURERULE_HH
28
29#include <dune/geometry/quadraturerules.hh>
30#include <dune/geometry/virtualrefinement.hh>
31#include <dumux/common/math.hh>
32
33namespace Dumux {
34
39template<typename ct, int mydim>
40class RefinementQuadratureRule final : public Dune::QuadratureRule<ct, mydim>
41{
42public:
44 static constexpr int dim = mydim;
45
47 static constexpr int highest_order = 1;
48
50
51 RefinementQuadratureRule(Dune::GeometryType type, int levels)
52 : Dune::QuadratureRule<ct, dim>(type, highest_order)
53 {
54 auto& refinement = Dune::buildRefinement<dim, ct>(type, type);
55 const auto tag = Dune::refinementLevels(levels);
56
57 // get the vertices
58 const auto numVertices = refinement.nVertices(tag);
59 std::vector<Dune::FieldVector<ct, dim>> points; points.reserve(numVertices);
60 auto vIt = refinement.vBegin(tag);
61 const auto vItEnd = refinement.vEnd(tag);
62 for (; vIt != vItEnd; ++vIt)
63 points.emplace_back(vIt.coords());
64
65 // go over all elements and get center and volume
66 const auto numElements = refinement.nElements(tag);
67 this->reserve(numElements);
68 auto eIt = refinement.eBegin(tag);
69 const auto eItEnd = refinement.eEnd(tag);
70 for (; eIt != eItEnd; ++eIt)
71 {
72 // quadrature point in the centroid of the sub-element and weight is its volume
73 const auto weight = computeVolume_(type, points, eIt.vertexIndices());
74 this->emplace_back(eIt.coords(), weight);
75 }
76 }
77
78private:
79 template<class VertexIndices>
80 ct computeVolume_(Dune::GeometryType t, const std::vector<Dune::FieldVector<ct, dim>>& points, const VertexIndices& indices) const
81 {
82 if (t != Dune::GeometryTypes::simplex(2))
83 DUNE_THROW(Dune::NotImplemented, "Only implemented for 2d simplices");
84
85 const auto ab = points[indices[1]] - points[indices[0]];
86 const auto ac = points[indices[2]] - points[indices[0]];
87 using std::abs;
88 return abs(crossProduct(ab, ac))*0.5;
89 }
90};
91
92} // end namespace Dumux
93
94#endif
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:640
Definition: adapt.hh:29
Definition: common/pdesolver.hh:35
A "quadrature" based on virtual refinement.
Definition: geometry/refinementquadraturerule.hh:41
~RefinementQuadratureRule() final
Definition: geometry/refinementquadraturerule.hh:49
RefinementQuadratureRule(Dune::GeometryType type, int levels)
Definition: geometry/refinementquadraturerule.hh:51
static constexpr int dim
The space dimension.
Definition: geometry/refinementquadraturerule.hh:44
static constexpr int highest_order
The highest quadrature order available.
Definition: geometry/refinementquadraturerule.hh:47