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
11namespace 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.111.3 (Nov 12, 23:30, 2024)