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 
34 namespace Ewoms {
38 template <class Scalar, unsigned dim>
40 {
41 public:
42  enum { numCorners = (1 << dim) };
43 
44  typedef Dune::FieldVector<Scalar, dim> LocalPosition;
45  typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
46 
47  Dune::GeometryType type() const
48  { return Dune::GeometryType(Dune::GeometryType::cube, dim); }
49 
50  template <class CornerContainer>
51  void setCorners(const CornerContainer& corners, unsigned numCorners)
52  {
53  unsigned cornerIdx;
54  for (cornerIdx = 0; cornerIdx < numCorners; ++cornerIdx) {
55  for (unsigned j = 0; j < dim; ++j)
56  corners_[cornerIdx][j] = corners[cornerIdx][j];
57  }
58  assert(cornerIdx == numCorners);
59 
60  center_ = 0;
61  for (cornerIdx = 0; cornerIdx < numCorners; ++cornerIdx)
62  center_ += corners_[cornerIdx];
63  center_ /= numCorners;
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 
147 private:
148  GlobalPosition corners_[numCorners];
149  GlobalPosition center_;
150 };
151 
152 } // namespace Ewoms
153 
154 #endif // EWOMS_QUADRATURE_GEOMETRY_HH
Definition: baseauxiliarymodule.hh:37
const GlobalPosition & corner(unsigned cornerIdx) const
Return the position of the corner with a given index.
Definition: quadraturegeometries.hh:127
Quadrature geometry for quadrilaterals.
Definition: quadraturegeometries.hh:39
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
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
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition: quadraturegeometries.hh:69
GlobalPosition global(const LocalPosition &localPos) const
Convert a local coordinate into a global one.
Definition: quadraturegeometries.hh:75
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