dune-geometry  2.8.0
axisalignedcubegeometry.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
5 #define DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
6 
11 #include <bitset>
12 
13 #include <dune/common/fvector.hh>
14 #include <dune/common/fmatrix.hh>
15 #include <dune/common/diagonalmatrix.hh>
16 
18 #include <dune/geometry/type.hh>
19 
20 
21 namespace Dune {
22 
46  template <class CoordType, unsigned int dim, unsigned int coorddim>
48  {
49 
50 
51  public:
52 
54  enum {mydimension = dim};
55 
57  enum {coorddimension = coorddim};
58 
60  typedef CoordType ctype;
61 
63  typedef FieldVector<ctype,dim> LocalCoordinate;
64 
66  typedef FieldVector<ctype,coorddim> GlobalCoordinate;
67 
69  typedef ctype Volume;
70 
77  typedef typename std::conditional<dim==coorddim,
78  DiagonalMatrix<ctype,dim>,
79  FieldMatrix<ctype,dim,coorddim> >::type JacobianTransposed;
80 
87  typedef typename std::conditional<dim==coorddim,
88  DiagonalMatrix<ctype,dim>,
89  FieldMatrix<ctype,coorddim,dim> >::type JacobianInverseTransposed;
90 
95  AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower,
96  const Dune::FieldVector<ctype,coorddim> upper)
97  : lower_(lower),
98  upper_(upper),
99  axes_()
100  {
101  static_assert(dim==coorddim, "Use this constructor only if dim==coorddim!");
102  // all 'true', but is never actually used
103  axes_ = (1<<coorddim)-1;
104  }
105 
113  AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower,
114  const Dune::FieldVector<ctype,coorddim> upper,
115  const std::bitset<coorddim>& axes)
116  : lower_(lower),
117  upper_(upper),
118  axes_(axes)
119  {
120  assert(axes.count()==dim);
121  for (size_t i=0; i<coorddim; i++)
122  if (not axes_[i])
123  upper_[i] = lower_[i];
124  }
125 
130  AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower)
131  : lower_(lower)
132  {}
133 
136  {
137  return GeometryTypes::cube(dim);
138  }
139 
142  {
143  GlobalCoordinate result;
144  if (dim == coorddim) { // fast case
145  for (size_t i=0; i<coorddim; i++)
146  result[i] = lower_[i] + local[i]*(upper_[i] - lower_[i]);
147  } else if (dim == 0) { // a vertex -- the other fast case
148  result = lower_; // hope for named-return-type-optimization
149  } else { // slow case
150  size_t lc=0;
151  for (size_t i=0; i<coorddim; i++)
152  result[i] = (axes_[i])
153  ? lower_[i] + local[lc++]*(upper_[i] - lower_[i])
154  : lower_[i];
155  }
156  return result;
157  }
158 
161  {
162  LocalCoordinate result;
163  if (dim == coorddim) { // fast case
164  for (size_t i=0; i<dim; i++)
165  result[i] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
166  } else if (dim != 0) { // slow case
167  size_t lc=0;
168  for (size_t i=0; i<coorddim; i++)
169  if (axes_[i])
170  result[lc++] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
171  }
172  return result;
173  }
174 
177  {
178  JacobianTransposed result;
179 
180  // Actually compute the result. Uses different methods depending
181  // on what kind of matrix JacobianTransposed is.
182  jacobianTransposed(result);
183 
184  return result;
185  }
186 
189  {
191 
192  // Actually compute the result. Uses different methods depending
193  // on what kind of matrix JacobianTransposed is.
195 
196  return result;
197  }
198 
202  ctype integrationElement([[maybe_unused]] const LocalCoordinate& local) const
203  {
204  return volume();
205  }
206 
209  {
210  GlobalCoordinate result;
211  if (dim==0)
212  result = lower_;
213  else {
214  // Since lower_==upper_ for unused coordinates, this always does the right thing
215  for (size_t i=0; i<coorddim; i++)
216  result[i] = CoordType(0.5) * (lower_[i] + upper_[i]);
217  }
218  return result;
219  }
220 
222  int corners() const
223  {
224  return 1<<dim;
225  }
226 
229  {
230  GlobalCoordinate result;
231  if (dim == coorddim) { // fast case
232  for (size_t i=0; i<coorddim; i++)
233  result[i] = (k & (1<<i)) ? upper_[i] : lower_[i];
234  } else if (dim == 0) { // vertex
235  result = lower_; // rely on named return-type optimization
236  } else { // slow case
237  unsigned int mask = 1;
238 
239  for (size_t i=0; i<coorddim; i++) {
240  if (not axes_[i])
241  result[i] = lower_[i];
242  else {
243  result[i] = (k & mask) ? upper_[i] : lower_[i];
244  mask = (mask<<1);
245  }
246  }
247  }
248 
249 
250  return result;
251  }
252 
254  Volume volume() const
255  {
256  ctype vol = 1;
257  if (dim == coorddim) { // fast case
258  for (size_t i=0; i<dim; i++)
259  vol *= upper_[i] - lower_[i];
260  // do nothing if dim == 0
261  } else if (dim != 0) { // slow case
262  for (size_t i=0; i<coorddim; i++)
263  if (axes_[i])
264  vol *= upper_[i] - lower_[i];
265  }
266  return vol;
267  }
268 
270  bool affine() const
271  {
272  return true;
273  }
274 
276  {
278  }
279 
280  private:
281  // jacobianTransposed: fast case --> diagonal matrix
282  void jacobianTransposed ( DiagonalMatrix<ctype,dim> &jacobianTransposed ) const
283  {
284  for (size_t i=0; i<dim; i++)
285  jacobianTransposed.diagonal()[i] = upper_[i] - lower_[i];
286  }
287 
288  // jacobianTransposed: slow case --> dense matrix
289  void jacobianTransposed ( FieldMatrix<ctype,dim,coorddim> &jacobianTransposed ) const
290  {
291  if (dim==0)
292  return;
293 
294  size_t lc = 0;
295  for (size_t i=0; i<coorddim; i++)
296  if (axes_[i])
297  jacobianTransposed[lc++][i] = upper_[i] - lower_[i];
298  }
299 
300  // jacobianInverseTransposed: fast case --> diagonal matrix
301  void jacobianInverseTransposed ( DiagonalMatrix<ctype,dim> &jacobianInverseTransposed ) const
302  {
303  for (size_t i=0; i<dim; i++)
304  jacobianInverseTransposed.diagonal()[i] = CoordType(1.0) / (upper_[i] - lower_[i]);
305  }
306 
307  // jacobianInverseTransposed: slow case --> dense matrix
308  void jacobianInverseTransposed ( FieldMatrix<ctype,coorddim,dim> &jacobianInverseTransposed ) const
309  {
310  if (dim==0)
311  return;
312 
313  size_t lc = 0;
314  for (size_t i=0; i<coorddim; i++)
315  if (axes_[i])
316  jacobianInverseTransposed[i][lc++] = CoordType(1.0) / (upper_[i] - lower_[i]);
317  }
318 
319  Dune::FieldVector<ctype,coorddim> lower_;
320 
321  Dune::FieldVector<ctype,coorddim> upper_;
322 
323  std::bitset<coorddim> axes_;
324  };
325 
326 } // namespace Dune
327 #endif
A unique label for each type of element that can occur in a grid.
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:470
unspecified-type ReferenceElement
Returns the type of reference element for the argument types T...
Definition: referenceelements.hh:415
Definition: affinegeometry.hh:19
static const ReferenceElement & cube()
get hypercube reference elements
Definition: referenceelements.hh:208
A geometry implementation for axis-aligned hypercubes.
Definition: axisalignedcubegeometry.hh:48
Volume volume() const
Return the element volume.
Definition: axisalignedcubegeometry.hh:254
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower, const Dune::FieldVector< ctype, coorddim > upper, const std::bitset< coorddim > &axes)
Constructor from a lower left and an upper right corner.
Definition: axisalignedcubegeometry.hh:113
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower, const Dune::FieldVector< ctype, coorddim > upper)
Constructor from a lower left and an upper right corner.
Definition: axisalignedcubegeometry.hh:95
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower)
Constructor from a single point only.
Definition: axisalignedcubegeometry.hh:130
GlobalCoordinate corner(int k) const
Return world coordinates of the k-th corner of the element.
Definition: axisalignedcubegeometry.hh:228
ctype Volume
Type used for volume.
Definition: axisalignedcubegeometry.hh:69
@ mydimension
Definition: axisalignedcubegeometry.hh:54
FieldVector< ctype, dim > LocalCoordinate
Type used for a vector of element coordinates.
Definition: axisalignedcubegeometry.hh:63
JacobianInverseTransposed jacobianInverseTransposed([[maybe_unused]] const LocalCoordinate &local) const
Jacobian transposed of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:188
friend Dune::Transitional::ReferenceElement< ctype, Dim< dim > > referenceElement(const AxisAlignedCubeGeometry &)
Definition: axisalignedcubegeometry.hh:275
JacobianTransposed jacobianTransposed([[maybe_unused]] const LocalCoordinate &local) const
Jacobian transposed of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:176
FieldVector< ctype, coorddim > GlobalCoordinate
Type used for a vector of world coordinates.
Definition: axisalignedcubegeometry.hh:66
LocalCoordinate local(const GlobalCoordinate &global) const
Map a point in global (world) coordinates to element coordinates.
Definition: axisalignedcubegeometry.hh:160
CoordType ctype
Type used for single coordinate coefficients.
Definition: axisalignedcubegeometry.hh:60
@ coorddimension
Definition: axisalignedcubegeometry.hh:57
std::conditional< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, dim, coorddim > >::type JacobianTransposed
Return type of jacobianTransposed.
Definition: axisalignedcubegeometry.hh:79
GeometryType type() const
Type of the cube. Here: a hypercube of the correct dimension.
Definition: axisalignedcubegeometry.hh:135
int corners() const
Return the number of corners of the element.
Definition: axisalignedcubegeometry.hh:222
ctype integrationElement([[maybe_unused]] const LocalCoordinate &local) const
Return the integration element, i.e., the determinant term in the integral transformation formula.
Definition: axisalignedcubegeometry.hh:202
std::conditional< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, coorddim, dim > >::type JacobianInverseTransposed
Return type of jacobianInverseTransposed.
Definition: axisalignedcubegeometry.hh:89
GlobalCoordinate center() const
Return center of mass of the element.
Definition: axisalignedcubegeometry.hh:208
bool affine() const
Return if the element is affine. Here: yes.
Definition: axisalignedcubegeometry.hh:270
GlobalCoordinate global(const LocalCoordinate &local) const
Map a point in local (element) coordinates to world coordinates.
Definition: axisalignedcubegeometry.hh:141
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123