dune-geometry  2.3.1
cachedmapping.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 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
4 #define DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
5 
6 #include <dune/geometry/type.hh>
11 
12 namespace Dune
13 {
14 
15  namespace GenericGeometry
16  {
17 
18  // Internal Forward Declarations
19  // -----------------------------
20 
21  template< unsigned int, class >
23 
24  template< unsigned int, class >
26 
27 
28 
29  // CachedStorage
30  // -------------
31 
32  template< unsigned int dim, class GeometryTraits >
34  {
35  friend class CachedJacobianTransposed< dim, GeometryTraits >;
36 
37  public:
38  static const unsigned int dimension = dim;
39  static const unsigned int dimWorld = GeometryTraits::dimWorld;
40 
42 
43  typedef typename GeometryTraits::Caching Caching;
44 
46  : affine( false ),
50  {}
51 
55 
56  bool affine : 1;
57 
58  bool jacobianTransposedComputed : 1; // = affine, if jacobian transposed was computed
59  bool jacobianInverseTransposedComputed : 1; // = affine, if jacobian inverse transposed was computed
60  bool integrationElementComputed : 1; // = affine, if integration element was computed
61  };
62 
63 
64 
65  // CachedJacobianTranposed
66  // -----------------------
67 
68  template< unsigned int dim, class GeometryTraits >
70  {
71  friend class CachedJacobianInverseTransposed< dim, GeometryTraits >;
72 
74  typedef typename Storage::Traits Traits;
75 
76  typedef typename Traits::MatrixHelper MatrixHelper;
77 
78  public:
79  typedef typename Traits::FieldType ctype;
80 
81  static const int rows = Traits::dimension;
82  static const int cols = Traits::dimWorld;
83 
85 
86  operator bool () const
87  {
88  return storage().jacobianTransposedComputed;
89  }
90 
91  operator const FieldMatrix & () const
92  {
93  return storage().jacobianTransposed;
94  }
95 
96  template< class X, class Y >
97  void mv ( const X &x, Y &y ) const
98  {
99  static_cast< const FieldMatrix & >( *this ).mv( x, y );
100  }
101 
102  template< class X, class Y >
103  void mtv ( const X &x, Y &y ) const
104  {
105  static_cast< const FieldMatrix & >( *this ).mtv( x, y );
106  }
107 
108  template< class X, class Y >
109  void umv ( const X &x, Y &y ) const
110  {
111  static_cast< const FieldMatrix & >( *this ).umv( x, y );
112  }
113 
114  template< class X, class Y >
115  void umtv ( const X &x, Y &y ) const
116  {
117  static_cast< const FieldMatrix & >( *this ).umtv( x, y );
118  }
119 
120  template< class X, class Y >
121  void mmv ( const X &x, Y &y ) const
122  {
123  static_cast< const FieldMatrix & >( *this ).mmv( x, y );
124  }
125 
126  template< class X, class Y >
127  void mmtv ( const X &x, Y &y ) const
128  {
129  static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
130  }
131 
132  ctype det () const
133  {
134  if( !storage().integrationElementComputed )
135  {
136  storage().integrationElement = MatrixHelper::template sqrtDetAAT< rows, cols >( storage().jacobianTransposed );
137  storage().integrationElementComputed = storage().affine;
138  }
139  return storage().integrationElement;
140  }
141 
142  private:
143  Storage &storage () const { return storage_; }
144 
145  mutable Storage storage_;
146  };
147 
148 
149 
150  // CachedJacobianInverseTransposed
151  // -------------------------------
152 
153  template< unsigned int dim, class GeometryTraits >
154  class CachedJacobianInverseTransposed
155  {
156  template< class, class > friend class CachedMapping;
157 
159  typedef typename JacobianTransposed::Storage Storage;
160  typedef typename JacobianTransposed::Traits Traits;
161 
162  typedef typename Traits::MatrixHelper MatrixHelper;
163 
164  public:
165  typedef typename Traits::FieldType ctype;
166 
167  static const int rows = Traits::dimWorld;
168  static const int cols = Traits::dimension;
169 
171 
172  operator bool () const
173  {
174  return storage().jacobianInverseTransposedComputed;
175  }
176 
177  operator const FieldMatrix & () const
178  {
179  return storage().jacobianInverseTransposed;
180  }
181 
182  template< class X, class Y >
183  void mv ( const X &x, Y &y ) const
184  {
185  static_cast< const FieldMatrix & >( *this ).mv( x, y );
186  }
187 
188  template< class X, class Y >
189  void mtv ( const X &x, Y &y ) const
190  {
191  static_cast< const FieldMatrix & >( *this ).mtv( x, y );
192  }
193 
194  template< class X, class Y >
195  void umv ( const X &x, Y &y ) const
196  {
197  static_cast< const FieldMatrix & >( *this ).umv( x, y );
198  }
199 
200  template< class X, class Y >
201  void umtv ( const X &x, Y &y ) const
202  {
203  static_cast< const FieldMatrix & >( *this ).umtv( x, y );
204  }
205 
206  template< class X, class Y >
207  void mmv ( const X &x, Y &y ) const
208  {
209  static_cast< const FieldMatrix & >( *this ).mmv( x, y );
210  }
211 
212  template< class X, class Y >
213  void mmtv ( const X &x, Y &y ) const
214  {
215  static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
216  }
217 
218  ctype det () const
219  {
220  // integrationElement is always computed with jacobianInverseTransposed
221  return ctype( 1 ) / storage().integrationElement;
222  }
223 
224  private:
225  JacobianTransposed &jacobianTransposed () { return jacobianTransposed_; }
226  const JacobianTransposed &jacobianTransposed () const { return jacobianTransposed_; }
227 
228  Storage &storage () const { return jacobianTransposed().storage(); }
229 
230  JacobianTransposed jacobianTransposed_;
231  };
232 
233 
234 
235  // CachedMapping
236  // -------------
237 
238  template< class Topology, class GeometryTraits >
240  {
242 
243  typedef typename GeometryTraits::template Mapping< Topology >::type MappingImpl;
244 
245  public:
247 
249 
250  static const unsigned int dimension = Traits::dimension;
251  static const unsigned int dimWorld = Traits::dimWorld;
252 
253  typedef typename Traits::FieldType FieldType;
256 
260 
262 
263  // can we safely assume that this mapping is always affine?
264  static const bool alwaysAffine = Mapping::alwaysAffine;
265 
266  typedef typename GeometryTraits::Caching Caching;
267 
268  private:
269  typedef typename Traits::MatrixHelper MatrixHelper;
270 
271  public:
272  template< class CoordVector >
273  explicit CachedMapping ( const CoordVector &coords )
274  : mapping_( coords )
275  {
276  if( alwaysAffine )
277  storage().affine = true;
278  else
279  computeJacobianTransposed( baryCenter() );
280  preCompute();
281  }
282 
283  template< class CoordVector >
284  explicit CachedMapping ( const std::pair< const CoordVector &, bool > &coords )
285  : mapping_( coords.first )
286  {
287  storage().affine = coords.second;
288  preCompute();
289  }
290 
291  bool affine () const { return (alwaysAffine || storage().affine); }
293 
294  int numCorners () const { return ReferenceElement::numCorners; }
295  GlobalCoordinate corner ( int i ) const { return mapping().corner( i ); }
297 
298  static bool checkInside ( const LocalCoordinate &x ) { return ReferenceElement::checkInside( x ); }
299 
301  {
303  if( jacobianTransposed() )
304  {
305  y = corner( 0 );
306  jacobianTransposed().umtv( x, y );
307  //MatrixHelper::template ATx< dimension, dimWorld >( jacobianTransposed_, x, y );
308  //y += corner( 0 );
309  }
310  else
311  mapping().global( x, y );
312  return y;
313  }
314 
316  {
317  LocalCoordinate x;
319  {
320  GlobalCoordinate z = y - corner( 0 );
321  jacobianInverseTransposed().mtv( z, x );
322  // MatrixHelper::template ATx< dimWorld, dimension >( jacobianInverseTransposed(), z, x );
323  }
324  else if( affine() )
325  {
326  const JacobianTransposed &JT = jacobianTransposed( baryCenter() );
327  GlobalCoordinate z = y - corner( 0 );
328  MatrixHelper::template xTRightInvA< dimension, dimWorld >( JT, z, x );
329  }
330  else
331  mapping().local( y, x );
332  return x;
333  }
334 
336  {
337  const EvaluationType evaluateI = Caching::evaluateIntegrationElement;
338  const EvaluationType evaluateJ = Caching::evaluateJacobianInverseTransposed;
339  if( ((evaluateI == PreCompute) || (evaluateJ == PreCompute)) && alwaysAffine )
340  return storage().integrationElement;
341  else
342  return jacobianTransposed( x ).det();
343  }
344 
345  FieldType volume () const
346  {
347  // do we need a quadrature of higher order, here?
348  const FieldType refVolume = ReferenceElement::volume();
349  return refVolume * integrationElement( baryCenter() );
350  }
351 
353  {
354  const EvaluationType evaluate = Caching::evaluateJacobianTransposed;
355  if( (evaluate == PreCompute) && alwaysAffine )
356  return jacobianTransposed();
357 
358  if( !jacobianTransposed() )
359  computeJacobianTransposed( x );
360  return jacobianTransposed();
361  }
362 
365  {
366  const EvaluationType evaluate = Caching::evaluateJacobianInverseTransposed;
367  if( (evaluate == PreCompute) && alwaysAffine )
368  return jacobianInverseTransposed();
369 
371  computeJacobianInverseTransposed( x );
372  return jacobianInverseTransposed();
373  }
374 
375  const Mapping &mapping () const { return mapping_; }
376 
377  private:
378  static const LocalCoordinate &baryCenter ()
379  {
381  }
382 
383  Storage &storage () const
384  {
385  return jacobianInverseTransposed().storage();
386  }
387 
388  const JacobianTransposed &jacobianTransposed () const
389  {
390  return jacobianInverseTransposed().jacobianTransposed();
391  }
392 
394  {
395  return jacobianInverseTransposed_;
396  }
397 
398  void preCompute ()
399  {
400  assert( affine() == mapping().jacobianTransposed( baryCenter(), storage().jacobianTransposed ) );
401  if( !affine() )
402  return;
403 
404  if( (Caching::evaluateJacobianTransposed == PreCompute) && !jacobianTransposed() )
405  computeJacobianTransposed( baryCenter() );
406 
407  if( Caching::evaluateJacobianInverseTransposed == PreCompute )
408  computeJacobianInverseTransposed( baryCenter() );
409  else if( Caching::evaluateIntegrationElement == PreCompute )
411  }
412 
413  void computeJacobianTransposed ( const LocalCoordinate &x ) const
414  {
415  storage().affine = mapping().jacobianTransposed( x, storage().jacobianTransposed );
416  storage().jacobianTransposedComputed = affine();
417  }
418 
419  void computeJacobianInverseTransposed ( const LocalCoordinate &x ) const
420  {
421  storage().integrationElement
422  = MatrixHelper::template rightInvA< dimension, dimWorld >( jacobianTransposed( x ), storage().jacobianInverseTransposed );
423  storage().integrationElementComputed = affine();
425  }
426 
427  private:
428  Mapping mapping_;
429  JacobianInverseTransposed jacobianInverseTransposed_;
430  };
431 
432  } // namespace GenericGeometry
433 
434 } // namespace Dune
435 
436 #endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
static const unsigned int dimWorld
Definition: geometrytraits.hh:57
Traits::JacobianType FieldMatrix
Definition: cachedmapping.hh:170
static const unsigned int numCorners
Definition: genericgeometry/referenceelements.hh:33
Traits::GlobalCoordinate GlobalCoordinate
Definition: cachedmapping.hh:255
static bool checkInside(const CoordinateType &x)
Definition: genericgeometry/referenceelements.hh:70
Definition: genericgeometry/referenceelements.hh:28
Definition: matrixhelper.hh:33
GenericGeometry::Mapping< typename GeometryTraits::CoordTraits, Topology, GeometryTraits::dimWorld, MappingImpl > Mapping
Definition: cachedmapping.hh:248
int numCorners() const
Definition: cachedmapping.hh:294
Traits::FieldType ctype
Definition: cachedmapping.hh:79
CoordTraits::template Matrix< dimension, dimWorld >::type JacobianTransposedType
Definition: geometrytraits.hh:66
GeometryTraits::Caching Caching
Definition: cachedmapping.hh:266
Traits::JacobianType jacobianInverseTransposed
Definition: cachedmapping.hh:53
EvaluationType
If not affine only volume is cached (based on intElCompute) otherwise all quantities can be cached...
Definition: geometrytraits.hh:76
CachedMapping(const std::pair< const CoordVector &, bool > &coords)
Definition: cachedmapping.hh:284
static const unsigned int dimension
Definition: cachedmapping.hh:250
A unique label for each type of element that can occur in a grid.
static ctype volume()
Definition: genericgeometry/referenceelements.hh:82
GlobalCoordinate center() const
Definition: cachedmapping.hh:296
static const unsigned int dimension
Definition: geometrytraits.hh:56
void umv(const X &x, Y &y) const
Definition: cachedmapping.hh:109
MappingTraits< typename GeometryTraits::CoordTraits, dimension, dimWorld > Traits
Definition: cachedmapping.hh:41
ctype det() const
Definition: cachedmapping.hh:132
static const bool alwaysAffine
Definition: cachedmapping.hh:264
Traits::JacobianTransposedType jacobianTransposed
Definition: cachedmapping.hh:52
CachedStorage()
Definition: cachedmapping.hh:45
void umtv(const X &x, Y &y) const
Definition: cachedmapping.hh:201
GlobalCoordinate global(const LocalCoordinate &x) const
Definition: cachedmapping.hh:300
const Mapping & mapping() const
Definition: cachedmapping.hh:375
bool affine() const
Definition: cachedmapping.hh:291
Dune::GeometryType type() const
Definition: cachedmapping.hh:292
static const FieldVector< ctype, dimension > & baryCenter()
Return the element barycenter.
Definition: genericgeometry/referenceelements.hh:59
GeometryTraits::Caching Caching
Definition: cachedmapping.hh:43
CoordTraits::template Vector< dimWorld >::type GlobalCoordinate
Definition: geometrytraits.hh:61
static const unsigned int dimWorld
Definition: cachedmapping.hh:39
CachedJacobianInverseTransposed< dimension, GeometryTraits > JacobianInverseTransposed
Definition: cachedmapping.hh:259
Definition: cachedmapping.hh:239
static const unsigned int dimWorld
Definition: cachedmapping.hh:251
void mtv(const X &x, Y &y) const
Definition: cachedmapping.hh:189
void umtv(const X &x, Y &y) const
Definition: cachedmapping.hh:115
Traits::FieldType integrationElement
Definition: cachedmapping.hh:54
assign in constructor using barycenter
Definition: geometrytraits.hh:81
static bool checkInside(const LocalCoordinate &x)
Definition: cachedmapping.hh:298
void global(const LocalCoordinate &x, GlobalCoordinate &y) const
Definition: mapping.hh:81
FieldType volume() const
Definition: cachedmapping.hh:345
CachedStorage< dimension, GeometryTraits > Storage
Definition: cachedmapping.hh:257
void mv(const X &x, Y &y) const
Definition: cachedmapping.hh:97
Definition: cachedmapping.hh:33
static const int cols
Definition: cachedmapping.hh:168
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &x) const
Definition: cachedmapping.hh:364
bool jacobianTransposedComputed
Definition: cachedmapping.hh:58
Traits::FieldType FieldType
Definition: cachedmapping.hh:253
bool integrationElementComputed
Definition: cachedmapping.hh:60
static const int rows
Definition: cachedmapping.hh:81
Implements some reference element functionality needed by the generic geometries. ...
static const int cols
Definition: cachedmapping.hh:82
static const unsigned int dimension
Definition: cachedmapping.hh:38
LocalCoordinate local(const GlobalCoordinate &y) const
Definition: cachedmapping.hh:315
Default mapping traits using Dune::FieldVector and Dune::FieldMatrix.
Definition: cornermapping.hh:21
static const int rows
Definition: cachedmapping.hh:167
void umv(const X &x, Y &y) const
Definition: cachedmapping.hh:195
CoordTraits::template Vector< dimension >::type LocalCoordinate
Definition: geometrytraits.hh:60
Traits::JacobianTransposedType FieldMatrix
Definition: cachedmapping.hh:84
const GlobalCoordinate & corner(int i) const
Definition: mapping.hh:76
Traits::LocalCoordinate LocalCoordinate
Definition: cachedmapping.hh:254
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &x) const
Definition: cachedmapping.hh:352
void mmv(const X &x, Y &y) const
Definition: cachedmapping.hh:121
ctype det() const
Definition: cachedmapping.hh:218
bool affine
Definition: cachedmapping.hh:56
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
GenericGeometry::ReferenceElement< Topology, FieldType > ReferenceElement
Definition: cachedmapping.hh:261
bool jacobianTransposed(const LocalCoordinate &x, JacobianTransposedType &JT) const
Definition: mapping.hh:104
void mv(const X &x, Y &y) const
Definition: cachedmapping.hh:183
Traits::FieldType ctype
Definition: cachedmapping.hh:165
FieldType integrationElement(const LocalCoordinate &x) const
Definition: cachedmapping.hh:335
CoordTraits::ctype FieldType
Definition: geometrytraits.hh:59
void mmtv(const X &x, Y &y) const
Definition: cachedmapping.hh:213
MappingTraits< typename GeometryTraits::CoordTraits, Topology::dimension, GeometryTraits::dimWorld > Traits
Definition: cachedmapping.hh:246
bool jacobianInverseTransposedComputed
Definition: cachedmapping.hh:59
CachedMapping(const CoordVector &coords)
Definition: cachedmapping.hh:273
void mmtv(const X &x, Y &y) const
Definition: cachedmapping.hh:127
void mtv(const X &x, Y &y) const
Definition: cachedmapping.hh:103
GlobalCoordinate corner(int i) const
Definition: cachedmapping.hh:295
void local(const GlobalCoordinate &y, LocalCoordinate &x) const
Definition: mapping.hh:86
CoordTraits::template Matrix< dimWorld, dimension >::type JacobianType
Definition: geometrytraits.hh:64
Definition: topologytypes.hh:271
CachedJacobianTransposed< dimension, GeometryTraits > JacobianTransposed
Definition: cachedmapping.hh:258
void mmv(const X &x, Y &y) const
Definition: cachedmapping.hh:207