dune-geometry  2.2.1
cachedmapping.hh
Go to the documentation of this file.
1 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
2 #define DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
3 
4 #include <dune/geometry/type.hh>
9 //#include <dune/geometry/genericgeometry/traceprovider.hh>
10 
11 namespace Dune
12 {
13 
14  namespace GenericGeometry
15  {
16 
17  // Internal Forward Declarations
18  // -----------------------------
19 
20  template< unsigned int, class >
22 
23  template< unsigned int, class >
25 
26 
27 
28  // CachedStorage
29  // -------------
30 
31  template< unsigned int dim, class GeometryTraits >
33  {
34  friend class CachedJacobianTransposed< dim, GeometryTraits >;
35 
36  public:
37  static const unsigned int dimension = dim;
38  static const unsigned int dimWorld = GeometryTraits::dimWorld;
39 
41 
42  typedef typename GeometryTraits::Caching Caching;
43 
45  : affine( false ),
49  {}
50 
54 
55  bool affine : 1;
56 
57  bool jacobianTransposedComputed : 1; // = affine, if jacobian transposed was computed
58  bool jacobianInverseTransposedComputed : 1; // = affine, if jacobian inverse transposed was computed
59  bool integrationElementComputed : 1; // = affine, if integration element was computed
60  };
61 
62 
63 
64  // CachedJacobianTranposed
65  // -----------------------
66 
67  template< unsigned int dim, class GeometryTraits >
69  {
70  friend class CachedJacobianInverseTransposed< dim, GeometryTraits >;
71 
73  typedef typename Storage::Traits Traits;
74 
75  typedef typename Traits::MatrixHelper MatrixHelper;
76 
77  public:
78  typedef typename Traits::FieldType ctype;
79 
80  static const int rows = Traits::dimension;
81  static const int cols = Traits::dimWorld;
82 
84 
85  operator bool () const
86  {
87  return storage().jacobianTransposedComputed;
88  }
89 
90  operator const FieldMatrix & () const
91  {
92  return storage().jacobianTransposed;
93  }
94 
95  template< class X, class Y >
96  void mv ( const X &x, Y &y ) const
97  {
98  static_cast< const FieldMatrix & >( *this ).mv( x, y );
99  }
100 
101  template< class X, class Y >
102  void mtv ( const X &x, Y &y ) const
103  {
104  static_cast< const FieldMatrix & >( *this ).mtv( x, y );
105  }
106 
107  template< class X, class Y >
108  void umv ( const X &x, Y &y ) const
109  {
110  static_cast< const FieldMatrix & >( *this ).umv( x, y );
111  }
112 
113  template< class X, class Y >
114  void umtv ( const X &x, Y &y ) const
115  {
116  static_cast< const FieldMatrix & >( *this ).umtv( x, y );
117  }
118 
119  template< class X, class Y >
120  void mmv ( const X &x, Y &y ) const
121  {
122  static_cast< const FieldMatrix & >( *this ).mmv( x, y );
123  }
124 
125  template< class X, class Y >
126  void mmtv ( const X &x, Y &y ) const
127  {
128  static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
129  }
130 
131  ctype det () const
132  {
133  if( !storage().integrationElementComputed )
134  {
135  storage().integrationElement = MatrixHelper::template sqrtDetAAT< rows, cols >( storage().jacobianTransposed );
136  storage().integrationElementComputed = storage().affine;
137  }
138  return storage().integrationElement;
139  }
140 
141  private:
142  Storage &storage () const { return storage_; }
143 
144  mutable Storage storage_;
145  };
146 
147 
148 
149  // CachedJacobianInverseTransposed
150  // -------------------------------
151 
152  template< unsigned int dim, class GeometryTraits >
153  class CachedJacobianInverseTransposed
154  {
155  template< class, class > friend class CachedMapping;
156 
158  typedef typename JacobianTransposed::Storage Storage;
159  typedef typename JacobianTransposed::Traits Traits;
160 
161  typedef typename Traits::MatrixHelper MatrixHelper;
162 
163  public:
164  typedef typename Traits::FieldType ctype;
165 
166  static const int rows = Traits::dimWorld;
167  static const int cols = Traits::dimension;
168 
170 
171  operator bool () const
172  {
173  return storage().jacobianInverseTransposedComputed;
174  }
175 
176  operator const FieldMatrix & () const
177  {
178  return storage().jacobianInverseTransposed;
179  }
180 
181  template< class X, class Y >
182  void mv ( const X &x, Y &y ) const
183  {
184  static_cast< const FieldMatrix & >( *this ).mv( x, y );
185  }
186 
187  template< class X, class Y >
188  void mtv ( const X &x, Y &y ) const
189  {
190  static_cast< const FieldMatrix & >( *this ).mtv( x, y );
191  }
192 
193  template< class X, class Y >
194  void umv ( const X &x, Y &y ) const
195  {
196  static_cast< const FieldMatrix & >( *this ).umv( x, y );
197  }
198 
199  template< class X, class Y >
200  void umtv ( const X &x, Y &y ) const
201  {
202  static_cast< const FieldMatrix & >( *this ).umtv( x, y );
203  }
204 
205  template< class X, class Y >
206  void mmv ( const X &x, Y &y ) const
207  {
208  static_cast< const FieldMatrix & >( *this ).mmv( x, y );
209  }
210 
211  template< class X, class Y >
212  void mmtv ( const X &x, Y &y ) const
213  {
214  static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
215  }
216 
217  ctype det () const
218  {
219  // integrationElement is always computed with jacobianInverseTransposed
220  return ctype( 1 ) / storage().integrationElement;
221  }
222 
223  private:
224  JacobianTransposed &jacobianTransposed () { return jacobianTransposed_; }
225  const JacobianTransposed &jacobianTransposed () const { return jacobianTransposed_; }
226 
227  Storage &storage () const { return jacobianTransposed().storage(); }
228 
229  JacobianTransposed jacobianTransposed_;
230  };
231 
232 
233 
234  // CachedMapping
235  // -------------
236 
237  template< class Topology, class GeometryTraits >
239  {
241 
242  typedef typename GeometryTraits::template Mapping< Topology >::type MappingImpl;
243 
244  public:
246 
248 
249  static const unsigned int dimension = Traits::dimension;
250  static const unsigned int dimWorld = Traits::dimWorld;
251 
252  typedef typename Traits::FieldType FieldType;
255 
259 
261 
262  // can we safely assume that this mapping is always affine?
263  static const bool alwaysAffine = Mapping::alwaysAffine;
264 
265 #if 0
266  template< unsigned int codim >
267  struct Codim
268  {
270  };
271 #endif
272 
273  typedef typename GeometryTraits::Caching Caching;
274 
275  private:
276  typedef typename Traits::MatrixHelper MatrixHelper;
277 
278  public:
279  template< class CoordVector >
280  explicit CachedMapping ( const CoordVector &coords )
281  : mapping_( coords )
282  {
283  if( alwaysAffine )
284  storage().affine = true;
285  else
286  computeJacobianTransposed( baryCenter() );
287  preCompute();
288  }
289 
290  template< class CoordVector >
291  explicit CachedMapping ( const std::pair< const CoordVector &, bool > &coords )
292  : mapping_( coords.first )
293  {
294  storage().affine = coords.second;
295  preCompute();
296  }
297 
298  bool affine () const { return (alwaysAffine || storage().affine); }
300 
301 #if 0
302  unsigned int topologyId () const DUNE_DEPRECATED { return type().id(); }
303 #endif
304 
305  int numCorners () const { return ReferenceElement::numCorners; }
306  GlobalCoordinate corner ( int i ) const { return mapping().corner( i ); }
307  GlobalCoordinate center () const { return global( ReferenceElement::template baryCenter< 0 >( 0 ) ); }
308 
309  static bool checkInside ( const LocalCoordinate &x ) { return ReferenceElement::checkInside( x ); }
310 
312  {
314  if( jacobianTransposed() )
315  {
316  y = corner( 0 );
317  jacobianTransposed().umtv( x, y );
318  //MatrixHelper::template ATx< dimension, dimWorld >( jacobianTransposed_, x, y );
319  //y += corner( 0 );
320  }
321  else
322  mapping().global( x, y );
323  return y;
324  }
325 
327  {
328  LocalCoordinate x;
330  {
331  GlobalCoordinate z = y - corner( 0 );
332  jacobianInverseTransposed().mtv( z, x );
333  // MatrixHelper::template ATx< dimWorld, dimension >( jacobianInverseTransposed(), z, x );
334  }
335  else if( affine() )
336  {
337  const JacobianTransposed &JT = jacobianTransposed( baryCenter() );
338  GlobalCoordinate z = y - corner( 0 );
339  MatrixHelper::template xTRightInvA< dimension, dimWorld >( JT, z, x );
340  }
341  else
342  mapping().local( y, x );
343  return x;
344  }
345 
347  {
348  const EvaluationType evaluateI = Caching::evaluateIntegrationElement;
349  const EvaluationType evaluateJ = Caching::evaluateJacobianInverseTransposed;
350  if( ((evaluateI == PreCompute) || (evaluateJ == PreCompute)) && alwaysAffine )
351  return storage().integrationElement;
352  else
353  return jacobianTransposed( x ).det();
354  }
355 
356  FieldType volume () const
357  {
358  // do we need a quadrature of higher order, here?
359  const FieldType refVolume = ReferenceElement::volume();
360  return refVolume * integrationElement( baryCenter() );
361  }
362 
364  {
365  const EvaluationType evaluate = Caching::evaluateJacobianTransposed;
366  if( (evaluate == PreCompute) && alwaysAffine )
367  return jacobianTransposed();
368 
369  if( !jacobianTransposed() )
370  computeJacobianTransposed( x );
371  return jacobianTransposed();
372  }
373 
376  {
377  const EvaluationType evaluate = Caching::evaluateJacobianInverseTransposed;
378  if( (evaluate == PreCompute) && alwaysAffine )
379  return jacobianInverseTransposed();
380 
382  computeJacobianInverseTransposed( x );
383  return jacobianInverseTransposed();
384  }
385 
386 #if 0
387  This *clone () const
388  {
389  return new This( *this );
390  }
391 
392  This* clone ( char *mappingStorage ) const
393  {
394  return new( mappingStorage ) This( *this );
395  }
396 
397  template< unsigned int codim, bool hybrid >
398  typename TraceProvider< Topology, GeometryTraits, codim, hybrid >::Trace*
399  trace ( unsigned int i, char *mappingStorage ) const
400  {
402  }
403 #endif
404 
405  const Mapping &mapping () const { return mapping_; }
406 
407  private:
408  static const LocalCoordinate &baryCenter ()
409  {
410  return ReferenceElement::template baryCenter< 0 >( 0 );
411  }
412 
413  Storage &storage () const
414  {
415  return jacobianInverseTransposed().storage();
416  }
417 
418  const JacobianTransposed &jacobianTransposed () const
419  {
420  return jacobianInverseTransposed().jacobianTransposed();
421  }
422 
424  {
425  return jacobianInverseTransposed_;
426  }
427 
428  void preCompute ()
429  {
430  assert( affine() == mapping().jacobianTransposed( baryCenter(), storage().jacobianTransposed ) );
431  if( !affine() )
432  return;
433 
434  if( (Caching::evaluateJacobianTransposed == PreCompute) && !jacobianTransposed() )
435  computeJacobianTransposed( baryCenter() );
436 
437  if( Caching::evaluateJacobianInverseTransposed == PreCompute )
438  computeJacobianInverseTransposed( baryCenter() );
439  else if( Caching::evaluateIntegrationElement == PreCompute )
441  }
442 
443  void computeJacobianTransposed ( const LocalCoordinate &x ) const
444  {
445  storage().affine = mapping().jacobianTransposed( x, storage().jacobianTransposed );
446  storage().jacobianTransposedComputed = affine();
447  }
448 
449  void computeJacobianInverseTransposed ( const LocalCoordinate &x ) const
450  {
451  storage().integrationElement
452  = MatrixHelper::template rightInvA< dimension, dimWorld >( jacobianTransposed( x ), storage().jacobianInverseTransposed );
453  storage().integrationElementComputed = affine();
455  }
456 
457  private:
458  Mapping mapping_;
459  JacobianInverseTransposed jacobianInverseTransposed_;
460  };
461 
462  } // namespace GenericGeometry
463 
464 } // namespace Dune
465 
466 #endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH