dune-geometry  2.4
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 #include <dune/common/unused.hh>
17 
18 #include <dune/geometry/type.hh>
19 
20 
21 
22 
23 namespace Dune {
24 
48  template <class CoordType, unsigned int dim, unsigned int coorddim>
50  {
51 
52 
53  public:
54 
56  enum {mydimension = dim};
57 
59  enum {coorddimension = coorddim};
60 
62  typedef CoordType ctype;
63 
65  typedef FieldVector<ctype,dim> LocalCoordinate;
66 
68  typedef FieldVector<ctype,coorddim> GlobalCoordinate;
69 
76  typedef typename conditional<dim==coorddim,
77  DiagonalMatrix<ctype,dim>,
78  FieldMatrix<ctype,dim,coorddim> >::type JacobianTransposed;
79 
86  typedef typename conditional<dim==coorddim,
87  DiagonalMatrix<ctype,dim>,
88  FieldMatrix<ctype,coorddim,dim> >::type JacobianInverseTransposed;
89 
94  AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower,
95  const Dune::FieldVector<ctype,coorddim> upper)
96  : lower_(lower),
97  upper_(upper),
98  axes_(),
99  jacobianTransposed_(0),
100  jacobianInverseTransposed_(0)
101  {
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  jacobianTransposed_(0),
120  jacobianInverseTransposed_(0)
121  {
122  assert(axes.count()==dim);
123  for (size_t i=0; i<coorddim; i++)
124  if (not axes_[i])
125  upper_[i] = lower_[i];
126  }
127 
132  AxisAlignedCubeGeometry(const Dune::FieldVector<ctype,coorddim> lower)
133  : lower_(lower),
134  jacobianTransposed_(0),
135  jacobianInverseTransposed_(0)
136  {}
137 
140  {
141  lower_ = other.lower_;
142  upper_ = other.upper_;
143  axes_ = other.axes_;
144  jacobianTransposed_ = other.jacobianTransposed_;
145  jacobianInverseTransposed_ = other.jacobianInverseTransposed_;
146  return *this;
147  }
148 
151  {
152  return GeometryType(GeometryType::cube,dim);
153  }
154 
156  GlobalCoordinate global(const LocalCoordinate& local) const
157  {
158  GlobalCoordinate result;
159  if (dim == coorddim) { // fast case
160  for (size_t i=0; i<coorddim; i++)
161  result[i] = lower_[i] + local[i]*(upper_[i] - lower_[i]);
162  } if (dim == 0) { // a vertex -- the other fast case
163  result = lower_; // hope for named-return-type-optimization
164  } else { // slow case
165  size_t lc=0;
166  for (size_t i=0; i<coorddim; i++)
167  result[i] = (axes_[i])
168  ? lower_[i] + local[lc++]*(upper_[i] - lower_[i])
169  : lower_[i];
170  }
171  return result;
172  }
173 
175  LocalCoordinate local(const GlobalCoordinate& global) const
176  {
177  LocalCoordinate result;
178  if (dim == coorddim) { // fast case
179  for (size_t i=0; i<dim; i++)
180  result[i] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
181  } else if (dim != 0) { // slow case
182  size_t lc=0;
183  for (size_t i=0; i<coorddim; i++)
184  if (axes_[i])
185  result[lc++] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
186  }
187  return result;
188  }
189 
191  const JacobianTransposed& jacobianTransposed(DUNE_UNUSED const LocalCoordinate& local) const
192  {
193  jacobianTransposed( jacobianTransposed_ );
194  return jacobianTransposed_;
195  }
196 
198  const JacobianInverseTransposed& jacobianInverseTransposed(DUNE_UNUSED const LocalCoordinate& local) const
199  {
200  jacobianInverseTransposed( jacobianInverseTransposed_ );
201  return jacobianInverseTransposed_;
202  }
203 
207  ctype integrationElement(DUNE_UNUSED const LocalCoordinate& local) const
208  {
209  return volume();
210  }
211 
213  GlobalCoordinate center() const
214  {
215  GlobalCoordinate result;
216  if (dim==0)
217  result = lower_;
218  else {
219  // Since lower_==upper_ for unused coordinates, this always does the right thing
220  for (size_t i=0; i<coorddim; i++)
221  result[i] = 0.5 * (lower_[i] + upper_[i]);
222  }
223  return result;
224  }
225 
227  int corners() const
228  {
229  return 1<<dim;
230  }
231 
233  GlobalCoordinate corner(int k) const
234  {
235  GlobalCoordinate result;
236  if (dim == coorddim) { // fast case
237  for (size_t i=0; i<coorddim; i++)
238  result[i] = (k & (1<<i)) ? upper_[i] : lower_[i];
239  } if (dim == 0) { // vertex
240  result = lower_; // rely on named return-type optimization
241  } else { // slow case
242  unsigned int mask = 1;
243 
244  for (size_t i=0; i<coorddim; i++) {
245  if (not axes_[i])
246  result[i] = lower_[i];
247  else {
248  result[i] = (k & mask) ? upper_[i] : lower_[i];
249  mask = (mask<<1);
250  }
251  }
252  }
253 
254 
255  return result;
256  }
257 
259  ctype volume() const
260  {
261  ctype vol = 1;
262  if (dim == coorddim) { // fast case
263  for (size_t i=0; i<dim; i++)
264  vol *= upper_[i] - lower_[i];
265  // do nothing if dim == 0
266  } else if (dim != 0) { // slow case
267  for (size_t i=0; i<coorddim; i++)
268  if (axes_[i])
269  vol *= upper_[i] - lower_[i];
270  }
271  return vol;
272  }
273 
275  bool affine() const
276  {
277  return true;
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] = 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++] = 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  // Storage so method jacobianTransposed can return a const reference
326  mutable JacobianTransposed jacobianTransposed_;
327 
328  // Storage so method jacobianInverseTransposed can return a const reference
329  mutable JacobianInverseTransposed jacobianInverseTransposed_;
330 
331  };
332 
333 } // namespace Dune
334 #endif
const JacobianTransposed & jacobianTransposed(DUNE_UNUSED const LocalCoordinate &local) const
Jacobian transposed of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:191
Definition: affinegeometry.hh:18
conditional< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, coorddim, dim > >::type JacobianInverseTransposed
Return type of jacobianInverseTransposed.
Definition: axisalignedcubegeometry.hh:88
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
ctype volume() const
Return the element volume.
Definition: axisalignedcubegeometry.hh:259
FieldVector< ctype, coorddim > GlobalCoordinate
Type used for a vector of world coordinates.
Definition: axisalignedcubegeometry.hh:68
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower)
Constructor from a single point only.
Definition: axisalignedcubegeometry.hh:132
const JacobianInverseTransposed & jacobianInverseTransposed(DUNE_UNUSED const LocalCoordinate &local) const
Jacobian transposed of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:198
Definition: axisalignedcubegeometry.hh:59
GlobalCoordinate corner(int k) const
Return world coordinates of the k-th corner of the element.
Definition: axisalignedcubegeometry.hh:233
AxisAlignedCubeGeometry & operator=(const AxisAlignedCubeGeometry &other)
Assignment operator.
Definition: axisalignedcubegeometry.hh:139
conditional< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, dim, coorddim > >::type JacobianTransposed
Return type of jacobianTransposed.
Definition: axisalignedcubegeometry.hh:78
GlobalCoordinate center() const
Return center of mass of the element.
Definition: axisalignedcubegeometry.hh:213
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:94
FieldVector< ctype, dim > LocalCoordinate
Type used for a vector of element coordinates.
Definition: axisalignedcubegeometry.hh:65
GlobalCoordinate global(const LocalCoordinate &local) const
Map a point in local (element) coordinates to world coordinates.
Definition: axisalignedcubegeometry.hh:156
A unique label for each type of element that can occur in a grid.
int corners() const
Return the number of corners of the element.
Definition: axisalignedcubegeometry.hh:227
GeometryType type() const
Type of the cube. Here: a hypercube of the correct dimension.
Definition: axisalignedcubegeometry.hh:150
LocalCoordinate local(const GlobalCoordinate &global) const
Map a point in global (world) coordinates to element coordinates.
Definition: axisalignedcubegeometry.hh:175
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
bool affine() const
Return if the element is affine. Here: yes.
Definition: axisalignedcubegeometry.hh:275
Cube element in any nonnegative dimension.
Definition: type.hh:31
ctype integrationElement(DUNE_UNUSED const LocalCoordinate &local) const
Return the integration element, i.e., the determinant term in the integral transformation formula...
Definition: axisalignedcubegeometry.hh:207
Definition: axisalignedcubegeometry.hh:56
CoordType ctype
Type used for single coordinate coefficients.
Definition: axisalignedcubegeometry.hh:62
A geometry implementation for axis-aligned hypercubes.
Definition: axisalignedcubegeometry.hh:49