Dune Core Modules (2.3.1)

referencedomain.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_REFERENCEDOMAIN_HH
4#define DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
5
6#include <algorithm>
7
8#include <dune/common/array.hh>
9#include <dune/common/fmatrix.hh>
10#include <dune/common/fvector.hh>
11#include <dune/common/typetraits.hh>
12#include <dune/common/unused.hh>
13
14#include <dune/geometry/genericgeometry/topologytypes.hh>
15#include <dune/geometry/genericgeometry/subtopologies.hh>
16
17namespace Dune
18{
19
20 namespace GenericGeometry
21 {
22
23 // Internal Forward Declarations
24 // -----------------------------
25
26 template< class Topology >
27 struct ReferenceDomain;
28
29
30
31 // ReferenceDomain
32 // ---------------
33
34 template< class Topology >
35 class ReferenceDomainBase;
36
38 template<>
39 class ReferenceDomainBase< Point >
40 {
41 typedef Point Topology;
42
43 friend struct ReferenceDomain< Topology >;
44 friend class ReferenceDomainBase< Prism< Topology > >;
45 friend class ReferenceDomainBase< Pyramid< Topology > >;
46
47 static const unsigned int numNormals = 0;
48
49 template< class ctype, int dim >
50 static void corner ( unsigned int i, FieldVector< ctype, dim > &n )
51 {
52 DUNE_UNUSED_PARAMETER(i);
53 DUNE_UNUSED_PARAMETER(n);
54 assert( i < Topology::numCorners );
55 }
56
57 template< class ctype, int dim >
58 static bool
59 checkInside ( const FieldVector< ctype, dim > &x, ctype factor )
60 {
61 DUNE_UNUSED_PARAMETER(x);
62 DUNE_UNUSED_PARAMETER(factor);
63 return true;
64 }
65
66 template< class ctype, int dim >
67 static void
68 integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
69 {
70 DUNE_UNUSED_PARAMETER(i);
71 DUNE_UNUSED_PARAMETER(n);
72 assert( i < numNormals );
73 }
74
75 template< class ctype >
76 static ctype volume ()
77 {
78 return ctype( 1 );
79 }
80 };
81
82
83 template< class BaseTopology >
84 class ReferenceDomainBase< Prism< BaseTopology > >
85 {
86 typedef Prism< BaseTopology > Topology;
87
88 friend struct ReferenceDomain< Topology >;
89 friend class ReferenceDomainBase< Prism< Topology > >;
90 friend class ReferenceDomainBase< Pyramid< Topology > >;
91
92 static const unsigned int numNormals = Size< Topology, 1 >::value;
93
94 static const unsigned int dimension = Topology::dimension;
95 static const unsigned int myindex = dimension - 1;
96
97 template< class ctype, int dim >
98 static void corner ( unsigned int i, FieldVector< ctype, dim > &x )
99 {
100 assert( i < Topology::numCorners );
101 const unsigned int j = i % BaseTopology::numCorners;
102 ReferenceDomainBase< BaseTopology >::corner( j, x );
103 if( i >= BaseTopology::numCorners )
104 x[ myindex ] = ctype( 1 );
105 }
106
107 template< class ctype, int dim >
108 static bool
109 checkInside ( const FieldVector< ctype, dim > &x, ctype factor )
110 {
111 const ctype xn = x[ myindex ];
112 const ctype cxn = factor - xn;
113 return (xn > -1e-12) && (cxn > -1e-12)
114 && ReferenceDomainBase< BaseTopology >::checkInside( x, factor );
115 }
116
117 template< class ctype, int dim >
118 static void
119 integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
120 {
121 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
122
123 if( i >= BaseReferenceDomain::numNormals )
124 {
125 const unsigned int j = i - BaseReferenceDomain::numNormals;
126 n[ myindex ] = (j == 0 ? ctype( -1 ) : ctype( 1 ));
127 }
128 else
129 BaseReferenceDomain::integrationOuterNormal( i, n );
130 }
131
132 template< class ctype >
133 static ctype volume ()
134 {
135 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
136 return BaseReferenceDomain::template volume< ctype >();
137 }
138 };
139
140
141 template< class BaseTopology >
142 class ReferenceDomainBase< Pyramid< BaseTopology > >
143 {
144 typedef Pyramid< BaseTopology > Topology;
145
146 friend struct ReferenceDomain< Topology >;
147 friend class ReferenceDomainBase< Prism< Topology > >;
148 friend class ReferenceDomainBase< Pyramid< Topology > >;
149
150 static const unsigned int numNormals = Size< Topology, 1 >::value;
151
152 static const unsigned int dimension = Topology::dimension;
153 static const unsigned int myindex = dimension - 1;
154
155 template< bool >
156 struct MultiDimensional
157 {
158 template< class ctype, int dim >
159 static void
160 integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
161 {
162 multiDimensionalIntegrationOuterNormal( i, n );
163 }
164 };
165
166 template< bool >
167 struct OneDimensional
168 {
169 template< class ctype, int dim >
170 static void
171 integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
172 {
173 n[ myindex ] = (i > 0) ? ctype( 1 ) : ctype( -1 );
174 }
175 };
176
177 template< class ctype, int dim >
178 static void corner ( unsigned int i, FieldVector< ctype, dim > &x )
179 {
180 assert( i < Topology::numCorners );
181 if( i < BaseTopology::numCorners )
182 ReferenceDomainBase< BaseTopology >::corner( i, x );
183 else
184 x[ myindex ] = ctype( 1 );
185 }
186
187 template< class ctype, int dim >
188 static bool
189 checkInside ( const FieldVector< ctype, dim > &x, ctype factor )
190 {
191 const ctype xn = x[ myindex ];
192 const ctype cxn = factor - xn;
193 return (xn > -1e-12) && (cxn > -1e-12)
194 && ReferenceDomainBase< BaseTopology >::checkInside( x, cxn );
195 }
196
197 template< class ctype, int dim >
198 static void
199 multiDimensionalIntegrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
200 {
201 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
202 typedef SubTopologyNumbering< BaseTopology, 1, dimension-2 > Numbering;
203
204 if( i > 0 )
205 {
206 const unsigned int j = Numbering::number( i-1, 0 );
207 FieldVector< ctype, dim > x( ctype( 0 ) );
208 BaseReferenceDomain::corner( j, x );
209
210 BaseReferenceDomain::integrationOuterNormal ( i-1, n );
211 n[ myindex ] = (x * n);
212 }
213 else
214 n[ myindex ] = ctype( -1 );
215 }
216
217 template< class ctype, int dim >
218 static void
219 integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
220 {
221 conditional< (dimension > 1), MultiDimensional<true>, OneDimensional<false> > :: type
222 ::integrationOuterNormal( i, n );
223 }
224
225 template< class ctype >
226 static ctype volume ()
227 {
228 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
229 const ctype baseVolume = BaseReferenceDomain::template volume< ctype >();
230 return baseVolume / ctype( (unsigned int)(dimension) ); // linker problem when using dimension directly
231 }
232 };
237 // ReferenceDomain
238 // ---------------
239
240 template< class Topology >
241 struct ReferenceDomain
242 {
243 static const unsigned int numCorners = Topology::numCorners;
244 static const unsigned int dimension = Topology::dimension;
245
246 static const unsigned int numNormals
247 = ReferenceDomainBase< Topology >::numNormals;
248
249 template< class ctype >
250 static void corner ( unsigned int i, FieldVector< ctype, dimension > &x )
251 {
252 x = ctype( 0 );
253 ReferenceDomainBase< Topology >::corner( i, x );
254 }
255
256 template< class ctype >
257 static bool checkInside ( const FieldVector< ctype, dimension > &x )
258 {
259 return ReferenceDomainBase< Topology >::checkInside( x, ctype( 1 ) );
260 }
261
262 template< class ctype >
263 static void
264 integrationOuterNormal ( unsigned int i, FieldVector< ctype, dimension > &n )
265 {
266 n = ctype( 0 );
267 return ReferenceDomainBase< Topology >::integrationOuterNormal( i, n );
268 }
269
270 template< class ctype >
271 static ctype volume ()
272 {
273 return ReferenceDomainBase< Topology >::template volume< ctype >();
274 }
275 };
276
277
278
279 // checkInside
280 // -----------
281
282 template< class ct, int cdim >
283 inline bool
284 checkInside ( unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
285 {
286 assert( (dim >= 0) && (dim <= cdim) );
287 assert( topologyId < numTopologies( dim ) );
288
289 if( dim > 0 )
290 {
291 const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
292 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
293 return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
294 else
295 return false;
296 }
297 else
298 return true;
299 }
300
301
302
303 // referenceCorners
304 // ----------------
305
306 template< class ct, int cdim >
307 inline unsigned int
308 referenceCorners ( unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners )
309 {
310 assert( (dim >= 0) && (dim <= cdim) );
311 assert( topologyId < numTopologies( dim ) );
312
313 if( dim > 0 )
314 {
315 const unsigned int nBaseCorners
316 = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
317 assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
318 if( isPrism( topologyId, dim ) )
319 {
320 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
321 for( unsigned int i = 0; i < nBaseCorners; ++i )
322 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
323 return 2*nBaseCorners;
324 }
325 else
326 {
327 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
328 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
329 return nBaseCorners+1;
330 }
331 }
332 else
333 {
334 *corners = FieldVector< ct, cdim >( ct( 0 ) );
335 return 1;
336 }
337 }
338
339
340
341 // referenceVolume
342 // ---------------
343
344 unsigned long referenceVolumeInverse ( unsigned int topologyId, int dim );
345
346 template< class ct >
347 inline ct referenceVolume ( unsigned int topologyId, int dim )
348 {
349 return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
350 }
351
352
353
354 // referenceOrigins
355 // ----------------
356
357 template< class ct, int cdim >
358 inline unsigned int
359 referenceOrigins ( unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins )
360 {
361 assert( (dim >= 0) && (dim <= cdim) );
362 assert( topologyId < numTopologies( dim ) );
363 assert( (codim >= 0) && (codim <= dim) );
364
365 if( codim > 0 )
366 {
367 const unsigned int baseId = baseTopologyId( topologyId, dim );
368 if( isPrism( topologyId, dim ) )
369 {
370 const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
371 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
372 for( unsigned int i = 0; i < m; ++i )
373 {
374 origins[ n+m+i ] = origins[ n+i ];
375 origins[ n+m+i ][ dim-1 ] = ct( 1 );
376 }
377 return n+2*m;
378 }
379 else
380 {
381 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
382 if( codim == dim )
383 {
384 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
385 origins[ m ][ dim-1 ] = ct( 1 );
386 return m+1;
387 }
388 else
389 return m+referenceOrigins( baseId, dim-1, codim, origins+m );
390 }
391 }
392 else
393 {
394 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
395 return 1;
396 }
397 }
398
399
400
401 // referenceEmbeddings
402 // -------------------
403
404 template< class ct, int cdim, int mydim >
405 inline unsigned int
406 referenceEmbeddings ( unsigned int topologyId, int dim, int codim,
407 FieldVector< ct, cdim > *origins,
408 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
409 {
410 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
411 assert( (dim - codim <= mydim) && (mydim <= cdim) );
412 assert( topologyId < numTopologies( dim ) );
413
414 if( codim > 0 )
415 {
416 const unsigned int baseId = baseTopologyId( topologyId, dim );
417 if( isPrism( topologyId, dim ) )
418 {
419 const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
420 for( unsigned int i = 0; i < n; ++i )
421 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
422
423 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
424 std::copy( origins+n, origins+n+m, origins+n+m );
425 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
426 for( unsigned int i = 0; i < m; ++i )
427 origins[ n+m+i ][ dim-1 ] = ct( 1 );
428
429 return n+2*m;
430 }
431 else
432 {
433 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
434 if( codim == dim )
435 {
436 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
437 origins[ m ][ dim-1 ] = ct( 1 );
438 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
439 return m+1;
440 }
441 else
442 {
443 const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
444 for( unsigned int i = 0; i < n; ++i )
445 {
446 for( int k = 0; k < dim-1; ++k )
447 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
448 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
449 }
450 return m+n;
451 }
452 }
453 }
454 else
455 {
456 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
457 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
458 for( int k = 0; k < dim; ++k )
459 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
460 return 1;
461 }
462 }
463
464
465
466 // referenceIntegrationOuterNormals
467 // --------------------------------
468
469 template< class ct, int cdim >
470 inline unsigned int
471 referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
472 const FieldVector< ct, cdim > *origins,
473 FieldVector< ct, cdim > *normals )
474 {
475 assert( (dim > 0) && (dim <= cdim) );
476 assert( topologyId < numTopologies( dim ) );
477
478 if( dim > 1 )
479 {
480 const unsigned int baseId = baseTopologyId( topologyId, dim );
481 if( isPrism( topologyId, dim ) )
482 {
483 const unsigned int numBaseFaces
484 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
485
486 for( unsigned int i = 0; i < 2; ++i )
487 {
488 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
489 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*int( i )-1 );
490 }
491
492 return numBaseFaces+2;
493 }
494 else
495 {
496 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
497 normals[ 0 ][ dim-1 ] = ct( -1 );
498
499 const unsigned int numBaseFaces
500 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
501 for( unsigned int i = 1; i <= numBaseFaces; ++i )
502 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
503
504 return numBaseFaces+1;
505 }
506 }
507 else
508 {
509 for( unsigned int i = 0; i < 2; ++i )
510 {
511 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
512 normals[ i ][ 0 ] = ct( 2*int( i )-1 );
513 }
514
515 return 2;
516 }
517 }
518
519 template< class ct, int cdim >
520 inline unsigned int
521 referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
522 FieldVector< ct, cdim > *normals )
523 {
524 assert( (dim > 0) && (dim <= cdim) );
525
526 FieldVector< ct, cdim > *origins
527 = new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
528 referenceOrigins( topologyId, dim, 1, origins );
529
530 const unsigned int numFaces
531 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
532 assert( numFaces == size( topologyId, dim, 1 ) );
533
534 delete[] origins;
535
536 return numFaces;
537 }
538
539 } // namespace GenericGeometry
540
541} // namespace Dune
542
543#endif // DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
array.hh
Fallback implementation of the std::array class (a static array)
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune
Dune namespace.
Definition: alignment.hh:14
typetraits.hh
Traits for type conversions and type information.
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 (Jul 15, 22:36, 2024)