Dune Core Modules (2.3.1)

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