Dune Core Modules (2.9.0)

geometrycache.hh
1#ifndef DUNE_SPGRID_GEOMETRYCACHE_HH
2#define DUNE_SPGRID_GEOMETRYCACHE_HH
3
4#include <utility>
5
6#include <dune/common/exceptions.hh>
7#include <dune/common/fvector.hh>
8#include <dune/common/fmatrix.hh>
9
10#include <dune/grid/spgrid/direction.hh>
11#include <dune/grid/spgrid/normal.hh>
12
13namespace Dune
14{
15
16 // SPGeometryPattern
17 // -----------------
18
19 template< int dim, int codim >
20 struct SPGeometryPattern
21 {
22 typedef SPDirection< dim > Direction;
23
24 SPGeometryPattern ();
25 explicit SPGeometryPattern ( Direction dir );
26
27 int nonzero ( const int k ) const;
28 int zero ( const int k ) const;
29
30 private:
31 int nonzero_[ dim - codim ];
32 int zero_[ codim ];
33 };
34
35 template< int dim >
36 struct SPGeometryPattern< dim, 0 >
37 {
38 typedef SPDirection< dim > Direction;
39
40 SPGeometryPattern () = default;
41 explicit SPGeometryPattern ( Direction /*dir*/ ) {}
42
43 int nonzero ( const int k ) const;
44 int zero ( const int k ) const;
45 };
46
47 template< int dim >
48 struct SPGeometryPattern< dim, dim >
49 {
50 typedef SPDirection< dim > Direction;
51
52 SPGeometryPattern () = default;
53 explicit SPGeometryPattern ( Direction /*dir*/ ) {}
54
55 int nonzero ( const int k ) const;
56 int zero ( const int k ) const;
57 };
58
59 template<>
60 struct SPGeometryPattern< 0, 0 >
61 {
62 typedef SPDirection< 0 > Direction;
63
64 SPGeometryPattern () = default;
65 explicit SPGeometryPattern ( Direction /*dir*/ ) {}
66
67 int nonzero ( const int k ) const;
68 int zero ( const int k ) const;
69 };
70
71
72
73 // SPJacobianTransposed
74 // --------------------
75
76 template< class ct, int dim, int mydim >
77 class SPJacobianTransposed
78 : private SPGeometryPattern< dim, dim - mydim >
79 {
80 typedef SPGeometryPattern< dim, dim - mydim > Pattern;
81
82 public:
83 typedef ct field_type;
84 typedef ct value_type;
85
86 typedef std::size_t size_type;
87
88 static const int rows = mydim;
89 static const int cols = dim;
90
91 typedef FieldVector< field_type, dim > GlobalVector;
92 typedef FieldVector< field_type, mydim > LocalVector;
93
94 typedef Dune::FieldMatrix< field_type, rows, cols > FieldMatrix;
95 typedef typename Dune::FieldTraits< field_type >::real_type real_type;
96
97 SPJacobianTransposed () = default;
98
99 SPJacobianTransposed ( const GlobalVector &h, SPDirection< dim > dir )
100 : Pattern( std::move( dir ) )
101 {
102 for( int k = 0; k < rows; ++k )
103 h_[ k ] = h[ pattern().nonzero( k ) ];
104 }
105
106 operator FieldMatrix () const;
107
108 template< class X, class Y > void mv ( const X &x, Y &y ) const;
109 template< class X, class Y > void mtv ( const X &x, Y &y ) const;
110 template< class X, class Y > void mhv ( const X &x, Y &y ) const;
111
112 template< class X, class Y > void umv ( const X &x, Y &y ) const;
113 template< class X, class Y > void umtv ( const X &x, Y &y ) const;
114 template< class X, class Y > void umhv ( const X &x, Y &y ) const;
115
116 template< class X, class Y > void mmv ( const X &x, Y &y ) const;
117 template< class X, class Y > void mmtv ( const X &x, Y &y ) const;
118 template< class X, class Y > void mmhv ( const X &x, Y &y ) const;
119
120 template< class X, class Y > void usmv ( const field_type &alpha, const X &x, Y &y ) const;
121 template< class X, class Y > void usmtv ( const field_type &alpha, const X &x, Y &y ) const;
122 template< class X, class Y > void usmhv ( const field_type &alpha, const X &x, Y &y ) const;
123
124 field_type det () const;
125
126 field_type determinant () const
127 {
128 if( rows == cols )
129 return det();
130 else
131 DUNE_THROW( FMatrixError, "There is no determinant for a " << rows << "x" << cols << " matrix" );
132 }
133
134 real_type frobenius_norm () const { return h().two_norm(); }
135 real_type frobenius_norm2 () const { return h().two_norm2(); }
136 real_type infinity_norm () const { return h().infinity_norm(); }
137 real_type infinity_norm_real () const { return h().infinity_norm_real(); }
138
139 private:
140 const Pattern &pattern () const { return static_cast< const Pattern & >( *this ); }
141 const LocalVector &h () const { return h_; }
142
143 LocalVector h_;
144 };
145
146
147
148 // FieldTraits for SPJacobianTransposed
149 // ------------------------------------
150
151 template< class ct, int dim, int mydim >
152 struct FieldTraits< SPJacobianTransposed< ct, dim, mydim > >
153 {
154 typedef typename FieldTraits< ct >::field_type field_type;
155 typedef typename FieldTraits< ct >::real_type real_type;
156 };
157
158
159
160 // SPJacobianInverseTransposed
161 // ---------------------------
162
163 template< class ct, int dim, int mydim >
164 class SPJacobianInverseTransposed
165 : private SPGeometryPattern< dim, dim - mydim >
166 {
167 typedef SPGeometryPattern< dim, dim - mydim > Pattern;
168
169 public:
170 typedef ct field_type;
171 typedef ct value_type;
172
173 typedef std::size_t size_type;
174
175 static const int rows = dim;
176 static const int cols = mydim;
177
178 typedef FieldVector< field_type, dim > GlobalVector;
179 typedef FieldVector< field_type, mydim > LocalVector;
180
181 typedef Dune::FieldMatrix< field_type, rows, cols > FieldMatrix;
182 typedef typename Dune::FieldTraits< field_type >::real_type real_type;
183
184 SPJacobianInverseTransposed () = default;
185
186 SPJacobianInverseTransposed ( const GlobalVector &h, SPDirection< dim > dir )
187 : Pattern( dir )
188 {
189 for( int k = 0; k < cols; ++k )
190 hInv_[ k ] = field_type( 1 ) / h[ pattern().nonzero( k ) ];
191 }
192
193 operator FieldMatrix () const;
194
195 template< class X, class Y > void mv ( const X &x, Y &y ) const;
196 template< class X, class Y > void mtv ( const X &x, Y &y ) const;
197 template< class X, class Y > void mhv ( const X &x, Y &y ) const;
198
199 template< class X, class Y > void umv ( const X &x, Y &y ) const;
200 template< class X, class Y > void umtv ( const X &x, Y &y ) const;
201 template< class X, class Y > void umhv ( const X &x, Y &y ) const;
202
203 template< class X, class Y > void usmv ( const field_type &alpha, const X &x, Y &y ) const;
204 template< class X, class Y > void usmtv ( const field_type &alpha, const X &x, Y &y ) const;
205 template< class X, class Y > void usmhv ( const field_type &alpha, const X &x, Y &y ) const;
206
207 template< class X, class Y > void mmv ( const X &x, Y &y ) const;
208 template< class X, class Y > void mmtv ( const X &x, Y &y ) const;
209 template< class X, class Y > void mmhv ( const X &x, Y &y ) const;
210
211 field_type det () const;
212
213 field_type determinant () const
214 {
215 if( rows == cols )
216 return det();
217 else
218 DUNE_THROW( FMatrixError, "There is no determinant for a " << rows << "x" << cols << " matrix" );
219 }
220
221 real_type frobenius_norm () const { return hInv().two_norm(); }
222 real_type frobenius_norm2 () const { return hInv().two_norm2(); }
223 real_type infinity_norm () const { return hInv().infinity_norm(); }
224 real_type infinity_norm_real () const { return hInv().infinity_norm_real(); }
225
226 private:
227 const Pattern &pattern () const { return static_cast< const Pattern & >( *this ); }
228 const LocalVector &hInv () const { return hInv_; }
229
230 LocalVector hInv_;
231 };
232
233
234
235 // FieldTraits for SPJacobianInverseTransposed
236 // -------------------------------------------
237
238 template< class ct, int dim, int mydim >
239 struct FieldTraits< SPJacobianInverseTransposed< ct, dim, mydim > >
240 {
241 typedef typename FieldTraits< ct >::field_type field_type;
242 typedef typename FieldTraits< ct >::real_type real_type;
243 };
244
245
246
247 // SPGeometryCache
248 // ---------------
249
250 template< class ct, int dim, int codim >
251 class SPGeometryCache
252 {
253 typedef SPGeometryCache< ct, dim, codim > This;
254
255 public:
256 typedef ct ctype;
257
258 static const int dimension = dim;
259 static const int codimension = codim;
260 static const int mydimension = dimension - codimension;
261
262 typedef SPDirection< dimension > Direction;
263
264 typedef FieldVector< ctype, dimension > GlobalVector;
265 typedef FieldVector< ctype, mydimension > LocalVector;
266
267 typedef SPJacobianTransposed< ctype, dimension, mydimension > JacobianTransposed;
268 typedef SPJacobianInverseTransposed< ctype, dimension, mydimension > JacobianInverseTransposed;
269
270 SPGeometryCache ( const GlobalVector &h, Direction dir )
271 : jacobianTransposed_( h, dir ), jacobianInverseTransposed_( h, dir ), volume_( jacobianTransposed_.det() )
272 {}
273
274 const ctype &volume () const { return volume_; }
275 const JacobianTransposed &jacobianTransposed () const { return jacobianTransposed_; }
276 const JacobianInverseTransposed &jacobianInverseTransposed () const { return jacobianInverseTransposed_; }
277
278 JacobianTransposed jacobianTransposed_;
279 JacobianInverseTransposed jacobianInverseTransposed_;
280 ctype volume_;
281 };
282
283
284
285 // Implementation of SPGeometryPattern
286 // -----------------------------------
287
288 template< int dim, int codim >
289 inline SPGeometryPattern< dim, codim >::SPGeometryPattern ()
290 {
291 const int mydim = dim - codim;
292 for( int k = 0; k < mydim; ++k )
293 nonzero_[ k ] = k;
294 for( int k = 0; k < codim; ++k )
295 zero_[ k ] = mydim + k;
296 }
297
298 template< int dim, int codim >
299 inline SPGeometryPattern< dim, codim >::SPGeometryPattern ( Direction dir )
300 {
301 int k = 0;
302 for( int j = 0; j < dim; ++j )
303 {
304 if( dir[ j ] != 0 )
305 nonzero_[ k++ ] = j;
306 else
307 zero_[ j - k ] = j;
308 }
309 assert( k == dim - codim );
310 }
311
312
313 template< int dim, int codim >
314 inline int SPGeometryPattern< dim, codim >::nonzero ( const int k ) const
315 {
316 assert( (k >= 0) && (k < dim - codim) );
317 return nonzero_[ k ];
318 }
319
320
321 template< int dim >
322 inline int SPGeometryPattern< dim, 0 >::nonzero ( const int k ) const
323 {
324 assert( (k >= 0) && (k < dim) );
325 return k;
326 }
327
328
329 template< int dim >
330 inline int SPGeometryPattern< dim, dim >::nonzero ( const int k ) const
331 {
332 assert( false );
333 return k;
334 }
335
336
337 inline int SPGeometryPattern< 0, 0 >::nonzero ( const int k ) const
338 {
339 assert( false );
340 return k;
341 }
342
343
344 template< int dim, int codim >
345 inline int SPGeometryPattern< dim, codim >::zero ( const int k ) const
346 {
347 assert( (k >= 0) && (k < codim) );
348 return zero_[ k ];
349 }
350
351
352 template< int dim >
353 inline int SPGeometryPattern< dim, 0 >::zero ( const int k ) const
354 {
355 assert( false );
356 return k;
357 }
358
359
360 template< int dim >
361 inline int SPGeometryPattern< dim, dim >::zero ( const int k ) const
362 {
363 assert( (k >= 0) && (k < dim) );
364 return k;
365 }
366
367
368 inline int SPGeometryPattern< 0, 0 >::zero ( const int k ) const
369 {
370 assert( false );
371 return k;
372 }
373
374
375
376 // Implementation of SPJacobianTransposed
377 // --------------------------------------
378
379 template< class ct, int dim, int mydim >
380 inline SPJacobianTransposed< ct, dim, mydim >::operator FieldMatrix () const
381 {
382 FieldMatrix matrix( field_type( 0 ) );
383 for( int k = 0; k < rows; ++k )
384 matrix[ k ][ pattern().nonzero( k ) ] = h()[ k ];
385 return matrix;
386 }
387
388
389 template< class ct, int dim, int mydim >
390 template< class X, class Y >
391 inline void SPJacobianTransposed< ct, dim, mydim >::mv ( const X &x, Y &y ) const
392 {
393 for( int k = 0; k < rows; ++k )
394 y[ k ] = h()[ k ] * x[ pattern().nonzero( k ) ];
395 }
396
397
398 template< class ct, int dim, int mydim >
399 template< class X, class Y >
400 inline void SPJacobianTransposed< ct, dim, mydim >::mtv ( const X &x, Y &y ) const
401 {
402 for( int k = 0; k < rows; ++k )
403 y[ pattern().nonzero( k ) ] = h()[ k ] * x[ k ];
404 for( int k = 0; k < cols - rows; ++k )
405 y[ pattern().zero( k ) ] = field_type( 0 );
406 }
407
408
409 template< class ct, int dim, int mydim >
410 template< class X, class Y >
411 inline void SPJacobianTransposed< ct, dim, mydim >::mhv ( const X &x, Y &y ) const
412 {
413 for( int k = 0; k < rows; ++k )
414 y[ pattern().nonzero( k ) ] = conjugateComplex( h()[ k ] ) * x[ k ];
415 for( int k = 0; k < cols - rows; ++k )
416 y[ pattern().zero( k ) ] = field_type( 0 );
417 }
418
419
420 template< class ct, int dim, int mydim >
421 template< class X, class Y >
422 inline void SPJacobianTransposed< ct, dim, mydim >::umv ( const X &x, Y &y ) const
423 {
424 for( int k = 0; k < rows; ++k )
425 y[ k ] += h()[ k ] * x[ pattern().nonzero( k ) ];
426 }
427
428
429 template< class ct, int dim, int mydim >
430 template< class X, class Y >
431 inline void SPJacobianTransposed< ct, dim, mydim >::umtv ( const X &x, Y &y ) const
432 {
433 for( int k = 0; k < rows; ++k )
434 y[ pattern().nonzero( k ) ] += h()[ k ] * x[ k ];
435 }
436
437
438 template< class ct, int dim, int mydim >
439 template< class X, class Y >
440 inline void SPJacobianTransposed< ct, dim, mydim >::umhv ( const X &x, Y &y ) const
441 {
442 for( int k = 0; k < rows; ++k )
443 y[ pattern().nonzero( k ) ] += conjugateComplex( h()[ k ] ) * x[ k ];
444 }
445
446
447 template< class ct, int dim, int mydim >
448 template< class X, class Y >
449 inline void SPJacobianTransposed< ct, dim, mydim >::usmv ( const field_type &alpha, const X &x, Y &y ) const
450 {
451 for( int k = 0; k < rows; ++k )
452 y[ k ] += alpha * h()[ k ] * x[ pattern().nonzero( k ) ];
453 }
454
455
456 template< class ct, int dim, int mydim >
457 template< class X, class Y >
458 inline void SPJacobianTransposed< ct, dim, mydim >::usmtv ( const field_type &alpha, const X &x, Y &y ) const
459 {
460 for( int k = 0; k < rows; ++k )
461 y[ pattern().nonzero( k ) ] += alpha * h()[ k ] * x[ k ];
462 }
463
464
465 template< class ct, int dim, int mydim >
466 template< class X, class Y >
467 inline void SPJacobianTransposed< ct, dim, mydim >::usmhv ( const field_type &alpha, const X &x, Y &y ) const
468 {
469 for( int k = 0; k < rows; ++k )
470 y[ pattern().nonzero( k ) ] += alpha * conjugateComplex( h()[ k ] ) * x[ k ];
471 }
472
473
474 template< class ct, int dim, int mydim >
475 template< class X, class Y >
476 inline void SPJacobianTransposed< ct, dim, mydim >::mmv ( const X &x, Y &y ) const
477 {
478 for( int k = 0; k < rows; ++k )
479 y[ k ] -= h()[ k ] * x[ pattern().nonzero( k ) ];
480 }
481
482
483 template< class ct, int dim, int mydim >
484 template< class X, class Y >
485 inline void SPJacobianTransposed< ct, dim, mydim >::mmtv ( const X &x, Y &y ) const
486 {
487 for( int k = 0; k < rows; ++k )
488 y[ pattern().nonzero( k ) ] -= h()[ k ] * x[ k ];
489 }
490
491
492 template< class ct, int dim, int mydim >
493 template< class X, class Y >
494 inline void SPJacobianTransposed< ct, dim, mydim >::mmhv ( const X &x, Y &y ) const
495 {
496 for( int k = 0; k < rows; ++k )
497 y[ pattern().nonzero( k ) ] -= conjugateComplex( h()[ k ] ) * x[ k ];
498 }
499
500
501 template< class ct, int dim, int mydim >
502 inline typename SPJacobianTransposed< ct, dim, mydim >::field_type SPJacobianTransposed< ct, dim, mydim >::det () const
503 {
504 field_type det( 1 );
505 for( int k = 0; k < rows; ++k )
506 det *= h()[ k ];
507 return det;
508 }
509
510
511
512 // Implementation of SPJacobianInverseTransposed
513 // ---------------------------------------------
514
515 template< class ct, int dim, int mydim >
516 inline SPJacobianInverseTransposed< ct, dim, mydim >::operator FieldMatrix () const
517 {
518 FieldMatrix matrix( field_type( 0 ) );
519 for( int k = 0; k < cols; ++k )
520 matrix[ pattern().nonzero( k ) ][ k ] = hInv()[ k ];
521 return matrix;
522 }
523
524
525 template< class ct, int dim, int mydim >
526 template< class X, class Y >
527 inline void SPJacobianInverseTransposed< ct, dim, mydim >::mv ( const X &x, Y &y ) const
528 {
529 for( int k = 0; k < cols; ++k )
530 y[ pattern().nonzero( k ) ] = hInv()[ k ] * x[ k ];
531 for( int k = 0; k < rows - cols; ++k )
532 y[ pattern().zero( k ) ] = field_type( 0 );
533 }
534
535
536 template< class ct, int dim, int mydim >
537 template< class X, class Y >
538 inline void SPJacobianInverseTransposed< ct, dim, mydim >::mtv ( const X &x, Y &y ) const
539 {
540 for( int k = 0; k < cols; ++k )
541 y[ k ] = hInv()[ k ] * x[ pattern().nonzero( k ) ];
542 }
543
544
545 template< class ct, int dim, int mydim >
546 template< class X, class Y >
547 inline void SPJacobianInverseTransposed< ct, dim, mydim >::mhv ( const X &x, Y &y ) const
548 {
549 for( int k = 0; k < cols; ++k )
550 y[ k ] = conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
551 }
552
553
554 template< class ct, int dim, int mydim >
555 template< class X, class Y >
556 inline void SPJacobianInverseTransposed< ct, dim, mydim >::umv ( const X &x, Y &y ) const
557 {
558 for( int k = 0; k < cols; ++k )
559 y[ pattern().nonzero( k ) ] += hInv()[ k ] * x[ k ];
560 }
561
562
563 template< class ct, int dim, int mydim >
564 template< class X, class Y >
565 inline void SPJacobianInverseTransposed< ct, dim, mydim >::umtv ( const X &x, Y &y ) const
566 {
567 for( int k = 0; k < cols; ++k )
568 y[ k ] += hInv()[ k ] * x[ pattern().nonzero( k ) ];
569 }
570
571
572 template< class ct, int dim, int mydim >
573 template< class X, class Y >
574 inline void SPJacobianInverseTransposed< ct, dim, mydim >::umhv ( const X &x, Y &y ) const
575 {
576 for( int k = 0; k < cols; ++k )
577 y[ k ] += conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
578 }
579
580
581 template< class ct, int dim, int mydim >
582 template< class X, class Y >
583 inline void SPJacobianInverseTransposed< ct, dim, mydim >::usmv ( const field_type &alpha, const X &x, Y &y ) const
584 {
585 for( int k = 0; k < cols; ++k )
586 y[ pattern().nonzero( k ) ] += alpha * hInv()[ k ] * x[ k ];
587 }
588
589
590 template< class ct, int dim, int mydim >
591 template< class X, class Y >
592 inline void SPJacobianInverseTransposed< ct, dim, mydim >::usmtv ( const field_type &alpha, const X &x, Y &y ) const
593 {
594 for( int k = 0; k < cols; ++k )
595 y[ k ] += alpha * hInv()[ k ] * x[ pattern().nonzero( k ) ];
596 }
597
598
599 template< class ct, int dim, int mydim >
600 template< class X, class Y >
601 inline void SPJacobianInverseTransposed< ct, dim, mydim >::usmhv ( const field_type &alpha, const X &x, Y &y ) const
602 {
603 for( int k = 0; k < cols; ++k )
604 y[ k ] += alpha * conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
605 }
606
607
608 template< class ct, int dim, int mydim >
609 template< class X, class Y >
610 inline void SPJacobianInverseTransposed< ct, dim, mydim >::mmv ( const X &x, Y &y ) const
611 {
612 for( int k = 0; k < cols; ++k )
613 y[ pattern().nonzero( k ) ] -= hInv()[ k ] * x[ k ];
614 }
615
616
617 template< class ct, int dim, int mydim >
618 template< class X, class Y >
619 inline void SPJacobianInverseTransposed< ct, dim, mydim >::mmtv ( const X &x, Y &y ) const
620 {
621 for( int k = 0; k < cols; ++k )
622 y[ k ] -= hInv()[ k ] * x[ pattern().nonzero( k ) ];
623 }
624
625
626 template< class ct, int dim, int mydim >
627 template< class X, class Y >
628 inline void SPJacobianInverseTransposed< ct, dim, mydim >::mmhv ( const X &x, Y &y ) const
629 {
630 for( int k = 0; k < cols; ++k )
631 y[ k ] -= conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
632 }
633
634
635 template< class ct, int dim, int mydim >
636 inline typename SPJacobianInverseTransposed< ct, dim, mydim >::field_type SPJacobianInverseTransposed< ct, dim, mydim >::det () const
637 {
638 field_type det( 1 );
639 for( int k = 0; k < cols; ++k )
640 det *= hInv()[ k ];
641 return det;
642 }
643
644} // namespace Dune
645
646#endif // #ifndef DUNE_SPGRID_GEOMETRYCACHE_HH
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:117
exceptions.hh
A few common exception classes.
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_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::conjugateComplex
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:164
std
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)