dune-geometry  2.2.1
cornermapping.hh
Go to the documentation of this file.
1 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CORNERMAPPING_HH
2 #define DUNE_GEOMETRY_GENERICGEOMETRY_CORNERMAPPING_HH
3 
7 
8 namespace Dune
9 {
10 
11  namespace GenericGeometry
12  {
13 
14  // External Forward Declarations
15  // -----------------------------
16 
17  template< class CT, unsigned int dim, unsigned int dimW >
18  struct MappingTraits;
19 
20 
21 
22  // GenericCornerMapping
23  // --------------------
24 
25  template< class Topology, class Traits, bool affine, unsigned int offset = 0 >
26  class GenericCornerMapping;
27 
28  template< class Traits, bool affine, unsigned int offset >
29  class GenericCornerMapping < Point, Traits, affine, offset >
30  {
31  typedef Point Topology;
32 
33  public:
34  static const unsigned int dim = Topology::dimension;
35  static const unsigned int dimW = Traits::dimWorld;
36 
37  typedef typename Traits::FieldType FieldType;
38  typedef typename Traits::LocalCoordinate LocalCoordinate;
39  typedef typename Traits::GlobalCoordinate GlobalCoordinate;
40  typedef typename Traits::JacobianTransposedType JacobianTransposedType;
41 
42  static const bool alwaysAffine = true;
43 
44  template< class CoordStorage >
45  static const GlobalCoordinate &origin ( const CoordStorage &coords )
46  {
47  dune_static_assert( CoordStorage::size, "Invalid offset." );
48  return coords[ offset ];
49  }
50 
51  template< class CoordStorage >
52  static void phi_set ( const CoordStorage &coords,
53  const LocalCoordinate &x,
54  const FieldType &factor,
55  GlobalCoordinate &p )
56  {
57  const GlobalCoordinate &y = origin( coords );
58  for( unsigned int i = 0; i < dimW; ++i )
59  p[ i ] = factor * y[ i ];
60  }
61 
62  template< class CoordStorage >
63  static void phi_add ( const CoordStorage &coords,
64  const LocalCoordinate &x,
65  const FieldType &factor,
66  GlobalCoordinate &p )
67  {
68  const GlobalCoordinate &y = origin( coords );
69  for( unsigned int i = 0; i < dimW; ++i )
70  p[ i ] += factor * y[ i ];
71  }
72 
73  template< class CoordStorage >
74  static bool Dphi_set ( const CoordStorage &coords,
75  const LocalCoordinate &x,
76  const FieldType &factor,
78  {
79  return true;
80  }
81 
82  template< class CoordStorage >
83  static bool Dphi_add ( const CoordStorage &coords,
84  const LocalCoordinate &x,
85  const FieldType &factor,
87  {
88  return true;
89  }
90  };
91 
92 
93  template< class BaseTopology, class Traits, bool affine, unsigned int offset >
94  class GenericCornerMapping< Prism< BaseTopology >, Traits, affine, offset >
95  {
97 
98  typedef GenericCornerMapping< BaseTopology, Traits, affine, offset >
99  BottomMapping;
100  typedef GenericCornerMapping
101  < BaseTopology, Traits, affine, offset + BaseTopology::numCorners >
102  TopMapping;
103 
104  public:
105  static const unsigned int dim = Topology::dimension;
106  static const unsigned int dimW = Traits::dimWorld;
107 
108  typedef typename Traits::FieldType FieldType;
109  typedef typename Traits::LocalCoordinate LocalCoordinate;
110  typedef typename Traits::GlobalCoordinate GlobalCoordinate;
111  typedef typename Traits::JacobianTransposedType JacobianTransposedType;
112 
113  static const bool alwaysAffine = ((dim < 2) || affine);
114 
115  template< class CoordStorage >
116  static const GlobalCoordinate &origin ( const CoordStorage &coords )
117  {
118  return BottomMapping::origin( coords );
119  }
120 
121  template< class CoordStorage >
122  static void phi_set ( const CoordStorage &coords,
123  const LocalCoordinate &x,
124  const FieldType &factor,
125  GlobalCoordinate &p )
126  {
127  const FieldType xn = x[ dim-1 ];
128  const FieldType cxn = FieldType( 1 ) - xn;
129  BottomMapping::phi_set( coords, x, factor * cxn, p );
130  TopMapping::phi_add( coords, x, factor * xn, p );
131  }
132 
133  template< class CoordStorage >
134  static void phi_add ( const CoordStorage &coords,
135  const LocalCoordinate &x,
136  const FieldType &factor,
137  GlobalCoordinate &p )
138  {
139  const FieldType xn = x[ dim-1 ];
140  const FieldType cxn = FieldType( 1 ) - xn;
141  BottomMapping::phi_add( coords, x, factor * cxn, p );
142  TopMapping::phi_add( coords, x, factor * xn, p );
143  }
144 
145  template< class CoordStorage >
146  static bool Dphi_set ( const CoordStorage &coords,
147  const LocalCoordinate &x,
148  const FieldType &factor,
150  {
151  const FieldType xn = x[ dim-1 ];
152  bool isAffine = true;
153  if( alwaysAffine )
154  {
155  const FieldType cxn = FieldType( 1 ) - xn;
156  BottomMapping::Dphi_set( coords, x, factor * cxn, J );
157  TopMapping::Dphi_add( coords, x, factor * xn, J );
158  }
159  else
160  {
162  isAffine &= BottomMapping::Dphi_set( coords, x, factor, J );
163  isAffine &= TopMapping::Dphi_set( coords, x, factor, Jtop );
164 
165  FieldType norm = FieldType( 0 );
166  for( unsigned int i = 0; i < dim-1; ++i )
167  {
168  Jtop[ i ] -= J[ i ];
169  norm += Jtop[ i ].two_norm2();
170  J[ i ].axpy( xn, Jtop[ i ] );
171  }
172  isAffine &= (norm < 1e-12);
173  }
174  BottomMapping::phi_set( coords, x, -factor, J[ dim-1 ] );
175  TopMapping::phi_add( coords, x, factor, J[ dim-1 ] );
176  return isAffine;
177  }
178 
179  template< class CoordStorage >
180  static bool Dphi_add ( const CoordStorage &coords,
181  const LocalCoordinate &x,
182  const FieldType &factor,
184  {
185  const FieldType xn = x[ dim-1 ];
186  bool isAffine = true;
187  if( alwaysAffine )
188  {
189  const FieldType cxn = FieldType( 1 ) - xn;
190  BottomMapping::Dphi_add( coords, x, factor * cxn, J );
191  TopMapping::Dphi_add( coords, x, factor * xn, J );
192  }
193  else
194  {
195  JacobianTransposedType Jbottom, Jtop;
196  isAffine &= BottomMapping::Dphi_set( coords, x, FieldType( 1 ), Jbottom );
197  isAffine &= TopMapping::Dphi_set( coords, x, FieldType( 1 ), Jtop );
198 
199  FieldType norm = FieldType( 0 );
200  for( unsigned int i = 0; i < dim-1; ++i )
201  {
202  Jtop[ i ] -= Jbottom[ i ];
203  norm += Jtop[ i ].two_norm2();
204  J[ i ].axpy( factor, Jbottom[ i ] );
205  J[ i ].axpy( factor*xn, Jtop[ i ] );
206  }
207  isAffine &= (norm < 1e-12);
208  }
209  BottomMapping::phi_add( coords, x, -factor, J[ dim-1 ] );
210  TopMapping::phi_add( coords, x, factor, J[ dim-1 ] );
211  return isAffine;
212  }
213  };
214 
215 
216  template< class BaseTopology, class Traits, bool affine, unsigned int offset >
217  class GenericCornerMapping < Pyramid< BaseTopology >, Traits, affine, offset >
218  {
220 
221  typedef GenericCornerMapping< BaseTopology, Traits, affine, offset >
222  BottomMapping;
223  typedef GenericCornerMapping
224  < Point, Traits, affine, offset + BaseTopology::numCorners >
225  TopMapping;
226 
227  public:
228  static const unsigned int dim = Topology::dimension;
229  static const unsigned int dimW = Traits::dimWorld;
230 
231  typedef typename Traits::FieldType FieldType;
232  typedef typename Traits::LocalCoordinate LocalCoordinate;
233  typedef typename Traits::GlobalCoordinate GlobalCoordinate;
234  typedef typename Traits::JacobianTransposedType JacobianTransposedType;
235 
236  static const bool alwaysAffine = (BottomMapping::alwaysAffine || affine);
237 
238  template< class CoordStorage >
239  static const GlobalCoordinate &origin ( const CoordStorage &coords )
240  {
241  return BottomMapping::origin( coords );
242  }
243 
244  template< class CoordStorage >
245  static void phi_set ( const CoordStorage &coords,
246  const LocalCoordinate &x,
247  const FieldType &factor,
248  GlobalCoordinate &p )
249  {
250  const FieldType xn = x[ dim-1 ];
251  if( alwaysAffine )
252  {
253  const GlobalCoordinate &top = TopMapping::origin( coords );
254  const GlobalCoordinate &bottom = BottomMapping::origin( coords );
255 
256  BottomMapping::phi_set( coords, x, factor, p );
257  for( unsigned int i = 0; i < dimW; ++i )
258  p[ i ] += (factor * xn) * (top[ i ] - bottom[ i ]);
259  }
260  else
261  {
262  TopMapping::phi_set( coords, x, factor * xn, p );
263  const FieldType cxn = FieldType( 1 ) - xn;
264  if( cxn > 1e-12 )
265  {
266  const FieldType icxn = FieldType( 1 ) / cxn;
267  LocalCoordinate xb;
268  for( unsigned int i = 0; i < dim-1; ++i )
269  xb[ i ] = icxn * x[ i ];
270 
271  BottomMapping::phi_add( coords, xb, factor * cxn, p );
272  }
273  }
274  }
275 
276  template< class CoordStorage >
277  static void phi_add ( const CoordStorage &coords,
278  const LocalCoordinate &x,
279  const FieldType &factor,
280  GlobalCoordinate &p )
281  {
282  const FieldType xn = x[ dim-1 ];
283  if( alwaysAffine )
284  {
285  const GlobalCoordinate &top = TopMapping::origin( coords );
286  const GlobalCoordinate &bottom = BottomMapping::origin( coords );
287 
288  BottomMapping::phi_add( coords, x, factor, p );
289  for( unsigned int i = 0; i < dimW; ++i )
290  p[ i ] += (factor * xn) * (top[ i ] - bottom[ i ]);
291  }
292  else
293  {
294  TopMapping::phi_add( coords, x, factor * xn, p );
295  const FieldType cxn = FieldType( 1 ) - xn;
296  if( cxn > 1e-12 )
297  {
298  const FieldType icxn = FieldType( 1 ) / cxn;
299  LocalCoordinate xb;
300  for( unsigned int i = 0; i < dim-1; ++i )
301  xb[ i ] = icxn * x[ i ];
302 
303  BottomMapping::phi_add( coords, xb, factor * cxn, p );
304  }
305  }
306  }
307 
308  template< class CoordStorage >
309  static bool Dphi_set ( const CoordStorage &coords,
310  const LocalCoordinate &x,
311  const FieldType &factor,
313  {
314  GlobalCoordinate &q = J[ dim-1 ];
315  bool isAffine;
316  if( alwaysAffine )
317  {
318  const GlobalCoordinate &top = TopMapping::origin( coords );
319  const GlobalCoordinate &bottom = BottomMapping::origin( coords );
320 
321  isAffine = BottomMapping::Dphi_set( coords, x, factor, J );
322  for( unsigned int i = 0; i < dimW; ++i )
323  q[ i ] = factor * (top[ i ] - bottom[ i ]);
324  }
325  else
326  {
327  const FieldType xn = x[ dim-1 ];
328  const FieldType cxn = FieldType( 1 ) - xn;
329  const FieldType icxn = FieldType( 1 ) / cxn;
330  LocalCoordinate xb;
331  for( unsigned int i = 0; i < dim-1; ++i )
332  xb[ i ] = icxn * x[ i ];
333  isAffine = BottomMapping::Dphi_set( coords, xb, factor, J );
334 
335  TopMapping::phi_set( coords, x, factor, q );
336  BottomMapping::phi_add( coords, xb, -factor, q );
337  xb *= factor;
338  for( unsigned int j = 0; j < dim-1; ++j )
339  {
340  for( unsigned int i = 0; i < dimW; ++i )
341  q[ i ] += J[ j ][ i ] * xb[ j ];
342  }
343  }
344  return isAffine;
345  }
346 
347  template< class CoordStorage >
348  static bool Dphi_add ( const CoordStorage &coords,
349  const LocalCoordinate &x,
350  const FieldType &factor,
352  {
353  GlobalCoordinate &q = J[ dim-1 ];
354  bool isAffine;
355  if( alwaysAffine )
356  {
357  const GlobalCoordinate &top = TopMapping::origin( coords );
358  const GlobalCoordinate &bottom = BottomMapping::origin( coords );
359 
360  isAffine = BottomMapping::Dphi_add( coords, x, factor, J );
361  for( unsigned int i = 0; i < dimW; ++i )
362  q[ i ] += factor * (top[ i ] - bottom[ i ]);
363  }
364  else
365  {
366  const FieldType xn = x[ dim-1 ];
367  const FieldType cxn = FieldType( 1 ) - xn;
368  const FieldType icxn = FieldType( 1 ) / cxn;
369  LocalCoordinate xb;
370  for( unsigned int i = 0; i < dim-1; ++i )
371  xb[ i ] = icxn * x[ i ];
372  isAffine = BottomMapping::Dphi_add( coords, xb, factor, J );
373 
374  TopMapping::phi_add( coords, x, factor, q );
375  BottomMapping::phi_add( coords, xb, -factor, q );
376  xb *= factor;
377  for( unsigned int j = 0; j < dim-1; ++j )
378  {
379  for( unsigned int i = 0; i < dimW; ++i )
380  q[ i ] += J[ j ][ i ] * xb[ j ];
381  }
382  }
383  return isAffine;
384  }
385  };
386 
387 
388 
389  // SubMappingCoords
390  // ----------------
391 
392  template< class Mapping, unsigned int codim >
394  {
395  typedef typename Mapping::GlobalCoordinate GlobalCoordinate;
397 
398  static const unsigned int dimension = ReferenceElement::dimension;
399 
400  const Mapping &mapping_;
401  const unsigned int i_;
402 
403  public:
404  SubMappingCoords ( const Mapping &mapping, unsigned int i )
405  : mapping_( mapping ), i_( i )
406  {}
407 
408  const GlobalCoordinate &operator[] ( unsigned int j ) const
409  {
410  const unsigned int k
411  = ReferenceElement::template subNumbering< codim, dimension - codim >( i_, j );
412  return mapping_.corner( k );
413  }
414  };
415 
416 
417 
418  // CoordStorage
419  // ------------
420 
425  template< class CoordTraits, class Topology, unsigned int dimW >
427  {
429 
430  public:
431  static const unsigned int size = Topology::numCorners;
432 
433  static const unsigned int dimWorld = dimW;
434 
435  typedef typename CoordTraits::template Vector< dimWorld >::type GlobalCoordinate;
436 
437  template< class SubTopology >
438  struct SubStorage
439  {
441  };
442 
443  template< class CoordVector >
444  explicit CoordStorage ( const CoordVector &coords )
445  {
446  for( unsigned int i = 0; i < size; ++i )
447  coords_[ i ] = coords[ i ];
448  }
449 
450  const GlobalCoordinate &operator[] ( unsigned int i ) const
451  {
452  return coords_[ i ];
453  }
454 
455  private:
456  GlobalCoordinate coords_[ size ];
457  };
458 
459 
460 
461  // CoordPointerStorage
462  // -------------------
463 
468  template< class CoordTraits, class Topology, unsigned int dimW >
470  {
472 
473  public:
474  static const unsigned int size = Topology::numCorners;
475 
476  static const unsigned int dimWorld = dimW;
477 
478  typedef typename CoordTraits::template Vector< dimWorld >::type GlobalCoordinate;
479 
480  template< class SubTopology >
481  struct SubStorage
482  {
484  };
485 
486  template< class CoordVector >
487  explicit CoordPointerStorage ( const CoordVector &coords )
488  {
489  for( unsigned int i = 0; i < size; ++i )
490  coords_[ i ] = &(coords[ i ]);
491  }
492 
493  const GlobalCoordinate &operator[] ( unsigned int i ) const
494  {
495  return *(coords_[ i ]);
496  }
497 
498  private:
499  const GlobalCoordinate *coords_[ size ];
500  };
501 
502 
503 
504  // CornerMapping
505  // -------------
506 
512  template< class CoordTraits, class Topo, unsigned int dimW,
513  class CStorage = CoordStorage< CoordTraits, Topo, dimW >,
514  bool affine = false >
516  {
518 
519  public:
520  typedef Topo Topology;
521  typedef CStorage CornerStorage;
523 
524  static const unsigned int dimension = Traits::dimension;
525  static const unsigned int dimWorld = Traits::dimWorld;
526 
527  typedef typename Traits::FieldType FieldType;
532 
534 
535  template< unsigned int codim, unsigned int i >
536  struct SubTopology
537  {
539  typedef typename CStorage::template SubStorage< Topology >::type CornerStorage;
541  };
542 
543  private:
544  typedef GenericGeometry::GenericCornerMapping< Topology, Traits, affine > GenericMapping;
545 
546  public:
548 
549  template< class CoordVector >
550  explicit CornerMapping ( const CoordVector &coords )
551  : coords_( coords )
552  {}
553 
554  const GlobalCoordinate &corner ( int i ) const
555  {
556  return coords_[ i ];
557  }
558 
559  void global ( const LocalCoordinate &x, GlobalCoordinate &y ) const
560  {
561  GenericMapping::phi_set( coords_, x, FieldType( 1 ), y );
562  }
563 
565  JacobianTransposedType &JT ) const
566  {
567  return GenericMapping::Dphi_set( coords_, x, FieldType( 1 ), JT );
568  }
569 
570  template< unsigned int codim, unsigned int i >
572  {
573  typedef typename SubTopology< codim, i >::Trace Trace;
574  typedef SubMappingCoords< This, codim > CoordVector;
575  return Trace( CoordVector( *this, i ) );
576  }
577 
578  protected:
580  };
581 
582  } // namespace GenericGeometry
583 
584 } // namespace Dune
585 
586 #endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CORNERMAPPING_HH