1#ifndef DUNE_SPGRID_GEOMETRYCACHE_HH
2#define DUNE_SPGRID_GEOMETRYCACHE_HH
6#include <dune/common/exceptions.hh>
7#include <dune/common/fvector.hh>
8#include <dune/common/fmatrix.hh>
19 template<
int dim,
int codim >
28 int zero (
const int k )
const;
31 int nonzero_[ dim - codim ];
43 int nonzero (
const int k )
const;
44 int zero (
const int k )
const;
55 int nonzero (
const int k )
const;
56 int zero (
const int k )
const;
67 int nonzero (
const int k )
const;
68 int zero (
const int k )
const;
76 template<
class ct,
int dim,
int mydim >
88 static const int rows = mydim;
89 static const int cols = dim;
94 typedef Dune::FieldMatrix< field_type, rows, cols >
FieldMatrix;
95 typedef typename Dune::FieldTraits< field_type >::real_type
real_type;
102 for(
int k = 0; k <
rows; ++k )
103 h_[ k ] = h[ pattern().
nonzero( k ) ];
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;
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;
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;
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;
131 DUNE_THROW( FMatrixError,
"There is no determinant for a " <<
rows <<
"x" <<
cols <<
" matrix" );
140 const Pattern &pattern ()
const {
return static_cast< const Pattern &
>( *this ); }
151 template<
class ct,
int dim,
int mydim >
155 typedef typename FieldTraits< ct >::real_type
real_type;
163 template<
class ct,
int dim,
int mydim >
182 typedef typename Dune::FieldTraits< field_type >::real_type
real_type;
189 for(
int k = 0; k <
cols; ++k )
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;
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;
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;
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;
218 DUNE_THROW( FMatrixError,
"There is no determinant for a " <<
rows <<
"x" <<
cols <<
" matrix" );
227 const Pattern &pattern ()
const {
return static_cast< const Pattern &
>( *this ); }
238 template<
class ct,
int dim,
int mydim >
242 typedef typename FieldTraits< ct >::real_type
real_type;
250 template<
class ct,
int dim,
int codim >
288 template<
int dim,
int codim >
291 const int mydim = dim - codim;
292 for(
int k = 0; k < mydim; ++k )
294 for(
int k = 0; k < codim; ++k )
295 zero_[ k ] = mydim + k;
298 template<
int dim,
int codim >
302 for(
int j = 0; j < dim; ++j )
309 assert( k == dim - codim );
313 template<
int dim,
int codim >
316 assert( (k >= 0) && (k < dim - codim) );
317 return nonzero_[ k ];
324 assert( (k >= 0) && (k < dim) );
344 template<
int dim,
int codim >
347 assert( (k >= 0) && (k < codim) );
363 assert( (k >= 0) && (k < dim) );
379 template<
class ct,
int dim,
int mydim >
383 for(
int k = 0; k < rows; ++k )
384 matrix[ k ][ pattern().nonzero( k ) ] = h()[ k ];
389 template<
class ct,
int dim,
int mydim >
390 template<
class X,
class Y >
393 for(
int k = 0; k < rows; ++k )
394 y[ k ] = h()[ k ] * x[ pattern().nonzero( k ) ];
398 template<
class ct,
int dim,
int mydim >
399 template<
class X,
class Y >
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 )
409 template<
class ct,
int dim,
int mydim >
410 template<
class X,
class Y >
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 )
420 template<
class ct,
int dim,
int mydim >
421 template<
class X,
class Y >
424 for(
int k = 0; k < rows; ++k )
425 y[ k ] += h()[ k ] * x[ pattern().nonzero( k ) ];
429 template<
class ct,
int dim,
int mydim >
430 template<
class X,
class Y >
433 for(
int k = 0; k < rows; ++k )
434 y[ pattern().nonzero( k ) ] += h()[ k ] * x[ k ];
438 template<
class ct,
int dim,
int mydim >
439 template<
class X,
class Y >
442 for(
int k = 0; k < rows; ++k )
443 y[ pattern().nonzero( k ) ] += conjugateComplex( h()[ k ] ) * x[ k ];
447 template<
class ct,
int dim,
int mydim >
448 template<
class X,
class Y >
451 for(
int k = 0; k < rows; ++k )
452 y[ k ] += alpha * h()[ k ] * x[ pattern().nonzero( k ) ];
456 template<
class ct,
int dim,
int mydim >
457 template<
class X,
class Y >
460 for(
int k = 0; k < rows; ++k )
461 y[ pattern().nonzero( k ) ] += alpha * h()[ k ] * x[ k ];
465 template<
class ct,
int dim,
int mydim >
466 template<
class X,
class Y >
469 for(
int k = 0; k < rows; ++k )
470 y[ pattern().nonzero( k ) ] += alpha * conjugateComplex( h()[ k ] ) * x[ k ];
474 template<
class ct,
int dim,
int mydim >
475 template<
class X,
class Y >
478 for(
int k = 0; k < rows; ++k )
479 y[ k ] -= h()[ k ] * x[ pattern().nonzero( k ) ];
483 template<
class ct,
int dim,
int mydim >
484 template<
class X,
class Y >
487 for(
int k = 0; k < rows; ++k )
488 y[ pattern().nonzero( k ) ] -= h()[ k ] * x[ k ];
492 template<
class ct,
int dim,
int mydim >
493 template<
class X,
class Y >
496 for(
int k = 0; k < rows; ++k )
497 y[ pattern().nonzero( k ) ] -= conjugateComplex( h()[ k ] ) * x[ k ];
501 template<
class ct,
int dim,
int mydim >
505 for(
int k = 0; k < rows; ++k )
515 template<
class ct,
int dim,
int mydim >
519 for(
int k = 0; k < cols; ++k )
520 matrix[ pattern().nonzero( k ) ][ k ] = hInv()[ k ];
525 template<
class ct,
int dim,
int mydim >
526 template<
class X,
class Y >
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 )
536 template<
class ct,
int dim,
int mydim >
537 template<
class X,
class Y >
540 for(
int k = 0; k < cols; ++k )
541 y[ k ] = hInv()[ k ] * x[ pattern().nonzero( k ) ];
545 template<
class ct,
int dim,
int mydim >
546 template<
class X,
class Y >
549 for(
int k = 0; k < cols; ++k )
550 y[ k ] = conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
554 template<
class ct,
int dim,
int mydim >
555 template<
class X,
class Y >
558 for(
int k = 0; k < cols; ++k )
559 y[ pattern().nonzero( k ) ] += hInv()[ k ] * x[ k ];
563 template<
class ct,
int dim,
int mydim >
564 template<
class X,
class Y >
567 for(
int k = 0; k < cols; ++k )
568 y[ k ] += hInv()[ k ] * x[ pattern().nonzero( k ) ];
572 template<
class ct,
int dim,
int mydim >
573 template<
class X,
class Y >
576 for(
int k = 0; k < cols; ++k )
577 y[ k ] += conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
581 template<
class ct,
int dim,
int mydim >
582 template<
class X,
class Y >
585 for(
int k = 0; k < cols; ++k )
586 y[ pattern().nonzero( k ) ] += alpha * hInv()[ k ] * x[ k ];
590 template<
class ct,
int dim,
int mydim >
591 template<
class X,
class Y >
594 for(
int k = 0; k < cols; ++k )
595 y[ k ] += alpha * hInv()[ k ] * x[ pattern().nonzero( k ) ];
599 template<
class ct,
int dim,
int mydim >
600 template<
class X,
class Y >
603 for(
int k = 0; k < cols; ++k )
604 y[ k ] += alpha * conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
608 template<
class ct,
int dim,
int mydim >
609 template<
class X,
class Y >
612 for(
int k = 0; k < cols; ++k )
613 y[ pattern().nonzero( k ) ] -= hInv()[ k ] * x[ k ];
617 template<
class ct,
int dim,
int mydim >
618 template<
class X,
class Y >
621 for(
int k = 0; k < cols; ++k )
622 y[ k ] -= hInv()[ k ] * x[ pattern().nonzero( k ) ];
626 template<
class ct,
int dim,
int mydim >
627 template<
class X,
class Y >
630 for(
int k = 0; k < cols; ++k )
631 y[ k ] -= conjugateComplex( hInv()[ k ] ) * x[ pattern().nonzero( k ) ];
635 template<
class ct,
int dim,
int mydim >
639 for(
int k = 0; k < cols; ++k )
Definition: iostream.hh:7
Definition: direction.hh:18
Definition: geometrycache.hh:21
SPDirection< dim > Direction
Definition: geometrycache.hh:22
int zero(const int k) const
Definition: geometrycache.hh:345
SPGeometryPattern(Direction dir)
Definition: geometrycache.hh:299
SPGeometryPattern()
Definition: geometrycache.hh:289
int nonzero(const int k) const
Definition: geometrycache.hh:314
SPGeometryPattern()=default
SPDirection< dim > Direction
Definition: geometrycache.hh:38
SPGeometryPattern(Direction dir)
Definition: geometrycache.hh:41
SPDirection< dim > Direction
Definition: geometrycache.hh:50
SPGeometryPattern(Direction dir)
Definition: geometrycache.hh:53
SPGeometryPattern()=default
SPGeometryPattern()=default
SPDirection< 0 > Direction
Definition: geometrycache.hh:62
SPGeometryPattern(Direction dir)
Definition: geometrycache.hh:65
Definition: geometrycache.hh:79
void usmhv(const field_type &alpha, const X &x, Y &y) const
Definition: geometrycache.hh:467
static const int rows
Definition: geometrycache.hh:88
void mv(const X &x, Y &y) const
Definition: geometrycache.hh:391
ct field_type
Definition: geometrycache.hh:83
field_type det() const
Definition: geometrycache.hh:502
field_type determinant() const
Definition: geometrycache.hh:126
std::size_t size_type
Definition: geometrycache.hh:86
Dune::FieldTraits< field_type >::real_type real_type
Definition: geometrycache.hh:95
static const int cols
Definition: geometrycache.hh:89
real_type infinity_norm() const
Definition: geometrycache.hh:136
void mmtv(const X &x, Y &y) const
Definition: geometrycache.hh:485
SPJacobianTransposed()=default
FieldVector< field_type, dim > GlobalVector
Definition: geometrycache.hh:91
real_type infinity_norm_real() const
Definition: geometrycache.hh:137
Dune::FieldMatrix< field_type, rows, cols > FieldMatrix
Definition: geometrycache.hh:94
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: geometrycache.hh:449
void usmtv(const field_type &alpha, const X &x, Y &y) const
Definition: geometrycache.hh:458
real_type frobenius_norm2() const
Definition: geometrycache.hh:135
void mmv(const X &x, Y &y) const
Definition: geometrycache.hh:476
real_type frobenius_norm() const
Definition: geometrycache.hh:134
void umhv(const X &x, Y &y) const
Definition: geometrycache.hh:440
void mtv(const X &x, Y &y) const
Definition: geometrycache.hh:400
FieldVector< field_type, mydim > LocalVector
Definition: geometrycache.hh:92
void mmhv(const X &x, Y &y) const
Definition: geometrycache.hh:494
ct value_type
Definition: geometrycache.hh:84
SPJacobianTransposed(const GlobalVector &h, SPDirection< dim > dir)
Definition: geometrycache.hh:99
void umv(const X &x, Y &y) const
Definition: geometrycache.hh:422
void umtv(const X &x, Y &y) const
Definition: geometrycache.hh:431
void mhv(const X &x, Y &y) const
Definition: geometrycache.hh:411
FieldTraits< ct >::field_type field_type
Definition: geometrycache.hh:154
FieldTraits< ct >::real_type real_type
Definition: geometrycache.hh:155
Definition: geometrycache.hh:166
void umv(const X &x, Y &y) const
Definition: geometrycache.hh:556
std::size_t size_type
Definition: geometrycache.hh:173
void mmv(const X &x, Y &y) const
Definition: geometrycache.hh:610
void umhv(const X &x, Y &y) const
Definition: geometrycache.hh:574
Dune::FieldTraits< field_type >::real_type real_type
Definition: geometrycache.hh:182
void mmhv(const X &x, Y &y) const
Definition: geometrycache.hh:628
real_type infinity_norm_real() const
Definition: geometrycache.hh:224
real_type infinity_norm() const
Definition: geometrycache.hh:223
static const int cols
Definition: geometrycache.hh:176
void usmhv(const field_type &alpha, const X &x, Y &y) const
Definition: geometrycache.hh:601
real_type frobenius_norm() const
Definition: geometrycache.hh:221
FieldVector< field_type, mydim > LocalVector
Definition: geometrycache.hh:179
real_type frobenius_norm2() const
Definition: geometrycache.hh:222
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: geometrycache.hh:583
ct value_type
Definition: geometrycache.hh:171
Dune::FieldMatrix< field_type, rows, cols > FieldMatrix
Definition: geometrycache.hh:181
SPJacobianInverseTransposed()=default
void mtv(const X &x, Y &y) const
Definition: geometrycache.hh:538
void usmtv(const field_type &alpha, const X &x, Y &y) const
Definition: geometrycache.hh:592
ct field_type
Definition: geometrycache.hh:170
void umtv(const X &x, Y &y) const
Definition: geometrycache.hh:565
field_type determinant() const
Definition: geometrycache.hh:213
field_type det() const
Definition: geometrycache.hh:636
void mv(const X &x, Y &y) const
Definition: geometrycache.hh:527
FieldVector< field_type, dim > GlobalVector
Definition: geometrycache.hh:178
void mhv(const X &x, Y &y) const
Definition: geometrycache.hh:547
SPJacobianInverseTransposed(const GlobalVector &h, SPDirection< dim > dir)
Definition: geometrycache.hh:186
static const int rows
Definition: geometrycache.hh:175
void mmtv(const X &x, Y &y) const
Definition: geometrycache.hh:619
FieldTraits< ct >::field_type field_type
Definition: geometrycache.hh:241
FieldTraits< ct >::real_type real_type
Definition: geometrycache.hh:242
Definition: geometrycache.hh:252
static const int codimension
Definition: geometrycache.hh:259
static const int dimension
Definition: geometrycache.hh:258
const JacobianTransposed & jacobianTransposed() const
Definition: geometrycache.hh:275
static const int mydimension
Definition: geometrycache.hh:260
FieldVector< ctype, dimension > GlobalVector
Definition: geometrycache.hh:264
const JacobianInverseTransposed & jacobianInverseTransposed() const
Definition: geometrycache.hh:276
JacobianInverseTransposed jacobianInverseTransposed_
Definition: geometrycache.hh:279
SPGeometryCache(const GlobalVector &h, Direction dir)
Definition: geometrycache.hh:270
const ctype & volume() const
Definition: geometrycache.hh:274
SPDirection< dimension > Direction
Definition: geometrycache.hh:262
SPJacobianTransposed< ctype, dimension, mydimension > JacobianTransposed
Definition: geometrycache.hh:267
ctype volume_
Definition: geometrycache.hh:280
JacobianTransposed jacobianTransposed_
Definition: geometrycache.hh:278
ct ctype
Definition: geometrycache.hh:256
FieldVector< ctype, mydimension > LocalVector
Definition: geometrycache.hh:265
SPJacobianInverseTransposed< ctype, dimension, mydimension > JacobianInverseTransposed
Definition: geometrycache.hh:268