quadraturegeometries.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 This file is part of the Open Porous Media project (OPM).
5
6 OPM 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 2 of the License, or
9 (at your option) any later version.
10
11 OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef EWOMS_QUADRATURE_GEOMETRIES_HH
28#define EWOMS_QUADRATURE_GEOMETRIES_HH
29
30#include <dune/common/fvector.hh>
31#include <dune/common/fmatrix.hh>
32#include <dune/geometry/type.hh>
33
34namespace Opm {
38template <class Scalar, unsigned dim>
40{
41public:
42 enum { numCorners = (1 << dim) };
43
44 using LocalPosition = Dune::FieldVector<Scalar, dim>;
45 using GlobalPosition = Dune::FieldVector<Scalar, dim>;
46
47 Dune::GeometryType type() const
48 { return Dune::GeometryType(/*topologyId=*/(1 << dim) - 1, dim); }
49
50 template <class CornerContainer>
51 void setCorners(const CornerContainer& corners, unsigned nCorners)
52 {
53 unsigned cornerIdx;
54 for (cornerIdx = 0; cornerIdx < nCorners; ++cornerIdx) {
55 for (unsigned j = 0; j < dim; ++j)
56 corners_[cornerIdx][j] = corners[cornerIdx][j];
57 }
58 assert(cornerIdx == nCorners);
59
60 center_ = 0;
61 for (cornerIdx = 0; cornerIdx < nCorners; ++cornerIdx)
62 center_ += corners_[cornerIdx];
63 center_ /= nCorners;
64 }
65
69 const GlobalPosition& center() const
70 { return center_; }
71
75 GlobalPosition global(const LocalPosition& localPos) const
76 {
77 GlobalPosition globalPos(0.0);
78
79 for (unsigned cornerIdx = 0; cornerIdx < numCorners; ++cornerIdx)
80 globalPos.axpy(cornerWeight(localPos, cornerIdx),
81 corners_[cornerIdx]);
82
83 return globalPos;
84 }
85
90 void jacobian(Dune::FieldMatrix<Scalar, dim, dim>& jac,
91 const LocalPosition& localPos) const
92 {
93 jac = 0.0;
94 for (unsigned cornerIdx = 0; cornerIdx < numCorners; ++cornerIdx) {
95 for (unsigned k = 0; k < dim; ++k) {
96 Scalar dWeight_dk = (cornerIdx& (1 << k)) ? 1 : -1;
97 for (unsigned j = 0; j < dim; ++j) {
98 if (k != j) {
99 if (cornerIdx& (1 << j))
100 dWeight_dk *= localPos[j];
101 else
102 dWeight_dk *= 1 - localPos[j];
103 ;
104 }
105 }
106
107 jac[k].axpy(dWeight_dk, corners_[cornerIdx]);
108 }
109 }
110 }
111
117 Scalar integrationElement(const LocalPosition& localPos) const
118 {
119 Dune::FieldMatrix<Scalar, dim, dim> jac;
120 jacobian(jac, localPos);
121 return jac.determinant();
122 }
123
127 const GlobalPosition& corner(unsigned cornerIdx) const
128 { return corners_[cornerIdx]; }
129
134 Scalar cornerWeight(const LocalPosition& localPos, unsigned cornerIdx) const
135 {
136 GlobalPosition globalPos(0.0);
137
138 // this code is based on the Q1 finite element code from
139 // dune-localfunctions
140 Scalar weight = 1.0;
141 for (unsigned j = 0; j < dim; ++j)
142 weight *= (cornerIdx& (1 << j)) ? localPos[j] : (1 - localPos[j]);
143
144 return weight;
145 }
146
147private:
148 GlobalPosition corners_[numCorners];
149 GlobalPosition center_;
150};
151
152} // namespace Opm
153
154#endif // EWOMS_QUADRATURE_GEOMETRY_HH
Quadrature geometry for quadrilaterals.
Definition: quadraturegeometries.hh:40
Scalar integrationElement(const LocalPosition &localPos) const
Return the determinant of the Jacobian of the mapping from local to global coordinates at a given loc...
Definition: quadraturegeometries.hh:117
Dune::GeometryType type() const
Definition: quadraturegeometries.hh:47
Scalar cornerWeight(const LocalPosition &localPos, unsigned cornerIdx) const
Return the weight of an individual corner for the local to global mapping.
Definition: quadraturegeometries.hh:134
const GlobalPosition & corner(unsigned cornerIdx) const
Return the position of the corner with a given index.
Definition: quadraturegeometries.hh:127
GlobalPosition global(const LocalPosition &localPos) const
Convert a local coordinate into a global one.
Definition: quadraturegeometries.hh:75
Dune::FieldVector< Scalar, dim > LocalPosition
Definition: quadraturegeometries.hh:44
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition: quadraturegeometries.hh:69
void setCorners(const CornerContainer &corners, unsigned nCorners)
Definition: quadraturegeometries.hh:51
void jacobian(Dune::FieldMatrix< Scalar, dim, dim > &jac, const LocalPosition &localPos) const
Returns the Jacobian matrix of the local to global mapping at a given local position.
Definition: quadraturegeometries.hh:90
@ numCorners
Definition: quadraturegeometries.hh:42
Dune::FieldVector< Scalar, dim > GlobalPosition
Definition: quadraturegeometries.hh:45
static const int dim
Definition: structuredgridvanguard.hh:67
Definition: blackoilboundaryratevector.hh:37