3 #ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH 4 #define DUNE_GEOMETRY_AFFINEGEOMETRY_HH 13 #include <dune/common/fmatrix.hh> 14 #include <dune/common/fvector.hh> 24 template<
class ctype,
int dim >
27 template<
class ctype,
int dim >
43 template<
int m,
int n >
44 static void Ax (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, n > &x, FieldVector< ctype, m > &ret )
46 for(
int i = 0; i < m; ++i )
48 ret[ i ] =
ctype( 0 );
49 for(
int j = 0; j < n; ++j )
50 ret[ i ] += A[ i ][ j ] * x[ j ];
54 template<
int m,
int n >
55 static void ATx (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, m > &x, FieldVector< ctype, n > &ret )
57 for(
int i = 0; i < n; ++i )
59 ret[ i ] =
ctype( 0 );
60 for(
int j = 0; j < m; ++j )
61 ret[ i ] += A[ j ][ i ] * x[ j ];
65 template<
int m,
int n,
int p >
66 static void AB (
const FieldMatrix< ctype, m, n > &A,
const FieldMatrix< ctype, n, p > &B, FieldMatrix< ctype, m, p > &ret )
68 for(
int i = 0; i < m; ++i )
70 for(
int j = 0; j < p; ++j )
72 ret[ i ][ j ] =
ctype( 0 );
73 for(
int k = 0; k < n; ++k )
74 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
79 template<
int m,
int n,
int p >
80 static void ATBT (
const FieldMatrix< ctype, m, n > &A,
const FieldMatrix< ctype, p, m > &B, FieldMatrix< ctype, n, p > &ret )
82 for(
int i = 0; i < n; ++i )
84 for(
int j = 0; j < p; ++j )
86 ret[ i ][ j ] =
ctype( 0 );
87 for(
int k = 0; k < m; ++k )
88 ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
93 template<
int m,
int n >
94 static void ATA_L (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
96 for(
int i = 0; i < n; ++i )
98 for(
int j = 0; j <= i; ++j )
100 ret[ i ][ j ] =
ctype( 0 );
101 for(
int k = 0; k < m; ++k )
102 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
107 template<
int m,
int n >
108 static void ATA (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
110 for(
int i = 0; i < n; ++i )
112 for(
int j = 0; j <= i; ++j )
114 ret[ i ][ j ] =
ctype( 0 );
115 for(
int k = 0; k < m; ++k )
116 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
117 ret[ j ][ i ] = ret[ i ][ j ];
120 ret[ i ][ i ] =
ctype( 0 );
121 for(
int k = 0; k < m; ++k )
122 ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
126 template<
int m,
int n >
127 static void AAT_L (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
137 for(
int i = 0; i < m; ++i )
139 for(
int j = 0; j <= i; ++j )
141 ctype &retij = ret[ i ][ j ];
142 retij = A[ i ][ 0 ] * A[ j ][ 0 ];
143 for(
int k = 1; k < n; ++k )
144 retij += A[ i ][ k ] * A[ j ][ k ];
149 template<
int m,
int n >
150 static void AAT (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
152 for(
int i = 0; i < m; ++i )
154 for(
int j = 0; j < i; ++j )
156 ret[ i ][ j ] =
ctype( 0 );
157 for(
int k = 0; k < n; ++k )
158 ret[ i ][ j ] += A[ i ][ k ] * A[ j ][ k ];
159 ret[ j ][ i ] = ret[ i ][ j ];
161 ret[ i ][ i ] =
ctype( 0 );
162 for(
int k = 0; k < n; ++k )
163 ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
168 static void Lx (
const FieldMatrix< ctype, n, n > &L,
const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
170 for(
int i = 0; i < n; ++i )
172 ret[ i ] =
ctype( 0 );
173 for(
int j = 0; j <= i; ++j )
174 ret[ i ] += L[ i ][ j ] * x[ j ];
179 static void LTx (
const FieldMatrix< ctype, n, n > &L,
const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
181 for(
int i = 0; i < n; ++i )
183 ret[ i ] =
ctype( 0 );
184 for(
int j = i; j < n; ++j )
185 ret[ i ] += L[ j ][ i ] * x[ j ];
190 static void LTL (
const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
192 for(
int i = 0; i < n; ++i )
194 for(
int j = 0; j < i; ++j )
196 ret[ i ][ j ] =
ctype( 0 );
197 for(
int k = i; k < n; ++k )
198 ret[ i ][ j ] += L[ k ][ i ] * L[ k ][ j ];
199 ret[ j ][ i ] = ret[ i ][ j ];
201 ret[ i ][ i ] =
ctype( 0 );
202 for(
int k = i; k < n; ++k )
203 ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
208 static void LLT (
const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
210 for(
int i = 0; i < n; ++i )
212 for(
int j = 0; j < i; ++j )
214 ret[ i ][ j ] =
ctype( 0 );
215 for(
int k = 0; k <= j; ++k )
216 ret[ i ][ j ] += L[ i ][ k ] * L[ j ][ k ];
217 ret[ j ][ i ] = ret[ i ][ j ];
219 ret[ i ][ i ] =
ctype( 0 );
220 for(
int k = 0; k <= i; ++k )
221 ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
226 static void cholesky_L (
const FieldMatrix< ctype, n, n > &A, FieldMatrix< ctype, n, n > &ret )
228 for(
int i = 0; i < n; ++i )
230 ctype &rii = ret[ i ][ i ];
232 ctype xDiag = A[ i ][ i ];
233 for(
int j = 0; j < i; ++j )
234 xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
235 assert( xDiag >
ctype( 0 ) );
238 ctype invrii =
ctype( 1 ) / rii;
239 for(
int k = i+1; k < n; ++k )
241 ctype x = A[ k ][ i ];
242 for(
int j = 0; j < i; ++j )
243 x -= ret[ i ][ j ] * ret[ k ][ j ];
244 ret[ k ][ i ] = invrii * x;
250 static ctype
detL (
const FieldMatrix< ctype, n, n > &L )
253 for(
int i = 0; i < n; ++i )
259 static ctype
invL ( FieldMatrix< ctype, n, n > &L )
262 for(
int i = 0; i < n; ++i )
264 ctype &lii = L[ i ][ i ];
266 lii =
ctype( 1 ) / lii;
267 for(
int j = 0; j < i; ++j )
269 ctype &lij = L[ i ][ j ];
270 ctype x = lij * L[ j ][ j ];
271 for(
int k = j+1; k < i; ++k )
272 x += L[ i ][ k ] * L[ k ][ j ];
281 static void invLx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
283 for(
int i = 0; i < n; ++i )
285 for(
int j = 0; j < i; ++j )
286 x[ i ] -= L[ i ][ j ] * x[ j ];
287 x[ i ] /= L[ i ][ i ];
293 static void invLTx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
295 for(
int i = n; i > 0; --i )
297 for(
int j = i; j < n; ++j )
298 x[ i-1 ] -= L[ j ][ i-1 ] * x[ j ];
299 x[ i-1 ] /= L[ i-1 ][ i-1 ];
304 static ctype
spdDetA (
const FieldMatrix< ctype, n, n > &A )
307 FieldMatrix< ctype, n, n > L;
313 static ctype
spdInvA ( FieldMatrix< ctype, n, n > &A )
315 FieldMatrix< ctype, n, n > L;
317 const ctype det =
invL( L );
324 static void spdInvAx ( FieldMatrix< ctype, n, n > &A, FieldVector< ctype, n > &x )
326 FieldMatrix< ctype, n, n > L;
332 template<
int m,
int n >
333 static ctype
detATA (
const FieldMatrix< ctype, m, n > &A )
337 FieldMatrix< ctype, n, n > ata;
350 template<
int m,
int n >
351 static ctype
sqrtDetAAT (
const FieldMatrix< ctype, m, n > &A )
358 if( (n == 2) && (m == 2) )
361 return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
363 else if( (n == 3) && (m == 3) )
366 const ctype v0 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 1 ][ 1 ] * A[ 0 ][ 2 ];
367 const ctype v1 = A[ 0 ][ 2 ] * A[ 1 ][ 0 ] - A[ 1 ][ 2 ] * A[ 0 ][ 0 ];
368 const ctype v2 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 1 ][ 0 ] * A[ 0 ][ 1 ];
369 return abs( v0 * A[ 2 ][ 0 ] + v1 * A[ 2 ][ 1 ] + v2 * A[ 2 ][ 2 ] );
371 else if ( (n == 3) && (m == 2) )
374 const ctype v0 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 0 ][ 1 ] * A[ 1 ][ 0 ];
375 const ctype v1 = A[ 0 ][ 0 ] * A[ 1 ][ 2 ] - A[ 1 ][ 0 ] * A[ 0 ][ 2 ];
376 const ctype v2 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 0 ][ 2 ] * A[ 1 ][ 1 ];
377 return sqrt( v0*v0 + v1*v1 + v2*v2);
382 FieldMatrix< ctype, m, m > aat;
392 template<
int m,
int n >
393 static ctype
leftInvA (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
395 static_assert((m >= n),
"Matrix has no left inverse.");
396 FieldMatrix< ctype, n, n > ata;
398 const ctype det =
spdInvA( ata );
403 template<
int m,
int n >
404 static void leftInvAx (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, m > &x, FieldVector< ctype, n > &y )
406 static_assert((m >= n),
"Matrix has no left inverse.");
407 FieldMatrix< ctype, n, n > ata;
414 template<
int m,
int n >
415 static ctype
rightInvA (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
417 static_assert((n >= m),
"Matrix has no right inverse.");
419 if( (n == 2) && (m == 2) )
421 const ctype det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
422 const ctype detInv =
ctype( 1 ) / det;
423 ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
424 ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
425 ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
426 ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
431 FieldMatrix< ctype, m , m > aat;
433 const ctype det =
spdInvA( aat );
434 ATBT( A , aat , ret );
439 template<
int m,
int n >
440 static void xTRightInvA (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, n > &x, FieldVector< ctype, m > &y )
442 static_assert((n >= m),
"Matrix has no right inverse.");
443 FieldMatrix< ctype, m, m > aat;
459 template<
class ct,
int mydim,
int cdim>
468 static const int mydimension= mydim;
471 static const int coorddimension = cdim;
496 AffineGeometry (
const ReferenceElement &refElement,
const GlobalCoordinate &origin,
497 const JacobianTransposed &jt )
498 : refElement_(&refElement), origin_(origin), jacobianTransposed_(jt)
500 integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
505 const JacobianTransposed &jt )
506 : refElement_( &ReferenceElements::general( gt ) ), origin_(origin), jacobianTransposed_( jt )
508 integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
512 template<
class CoordVector >
513 AffineGeometry (
const ReferenceElement &refElement,
const CoordVector &coordVector )
514 : refElement_(&refElement), origin_(coordVector[0])
516 for(
int i = 0; i < mydimension; ++i )
517 jacobianTransposed_[ i ] = coordVector[ i+1 ] - origin_;
518 integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
522 template<
class CoordVector >
524 : refElement_(&ReferenceElements::general( gt )), origin_(coordVector[0] )
526 for(
int i = 0; i < mydimension; ++i )
527 jacobianTransposed_[ i ] = coordVector[ i+1 ] - origin_;
528 integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
538 int corners ()
const {
return refElement_->size( mydimension ); }
543 return global( refElement_->position( i, mydimension ) );
547 GlobalCoordinate
center ()
const {
return global( refElement_->position( 0, 0 ) ); }
555 GlobalCoordinate
global (
const LocalCoordinate &local )
const 557 GlobalCoordinate global( origin_ );
558 jacobianTransposed_.umtv( local, global );
575 LocalCoordinate
local (
const GlobalCoordinate &global )
const 577 LocalCoordinate local;
578 jacobianInverseTransposed_.mtv( global - origin_, local );
594 DUNE_UNUSED_PARAMETER(local);
595 return integrationElement_;
601 return integrationElement_ * refElement_->volume();
612 DUNE_UNUSED_PARAMETER(local);
613 return jacobianTransposed_;
624 DUNE_UNUSED_PARAMETER(local);
625 return jacobianInverseTransposed_;
631 const ReferenceElement* refElement_;
632 GlobalCoordinate origin_;
633 JacobianTransposed jacobianTransposed_;
634 JacobianInverseTransposed jacobianInverseTransposed_;
635 ctype integrationElement_;
640 #endif // #ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH AffineGeometry(Dune::GeometryType gt, const CoordVector &coordVector)
Create affine geometry from GeometryType and a vector of vertex coordinates.
Definition: affinegeometry.hh:523
FieldVector< ctype, mydimension > LocalCoordinate
Type for local coordinate vector.
Definition: affinegeometry.hh:474
static void xTRightInvA(const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &y)
Definition: affinegeometry.hh:440
Class providing access to the singletons of the reference elements.
Definition: affinegeometry.hh:28
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:267
static void AAT_L(const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret)
Definition: affinegeometry.hh:127
friend const ReferenceElement & referenceElement(const AffineGeometry &geometry)
Definition: affinegeometry.hh:628
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian's inverse.
Definition: affinegeometry.hh:622
static ctype rightInvA(const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret)
Compute right pseudo-inverse of matrix A.
Definition: affinegeometry.hh:415
GlobalCoordinate corner(int i) const
Obtain coordinates of the i-th corner.
Definition: affinegeometry.hh:541
static void ATA_L(const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret)
Definition: affinegeometry.hh:94
static void ATx(const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &ret)
Definition: affinegeometry.hh:55
Implementation of the Geometry interface for affine geometries.
Definition: affinegeometry.hh:460
bool affine() const
Always true: this is an affine geometry.
Definition: affinegeometry.hh:532
This class provides access to geometric and topological properties of a reference element...
Definition: affinegeometry.hh:25
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian.
Definition: affinegeometry.hh:610
Definition: affinegeometry.hh:18
static ctype detATA(const FieldMatrix< ctype, m, n > &A)
Definition: affinegeometry.hh:333
A unique label for each type of element that can occur in a grid.
AffineGeometry(const ReferenceElement &refElement, const GlobalCoordinate &origin, const JacobianTransposed &jt)
Create affine geometry from reference element, one vertex, and the Jacobian matrix.
Definition: affinegeometry.hh:496
static void LLT(const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret)
Definition: affinegeometry.hh:208
static ctype spdDetA(const FieldMatrix< ctype, n, n > &A)
Definition: affinegeometry.hh:304
static ctype sqrtDetAAT(const FieldMatrix< ctype, m, n > &A)
Compute the square root of the determinant of A times A transposed.
Definition: affinegeometry.hh:351
ctype integrationElement(const LocalCoordinate &local) const
Obtain the integration element.
Definition: affinegeometry.hh:592
static ctype invL(FieldMatrix< ctype, n, n > &L)
Definition: affinegeometry.hh:259
static void leftInvAx(const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &y)
Definition: affinegeometry.hh:404
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the mapping.
Definition: affinegeometry.hh:555
static ctype leftInvA(const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret)
Definition: affinegeometry.hh:393
Dune::GeometryType type() const
Obtain the type of the reference element.
Definition: affinegeometry.hh:535
ct ctype
Type used for coordinates.
Definition: affinegeometry.hh:465
static void LTL(const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret)
Definition: affinegeometry.hh:190
static void cholesky_L(const FieldMatrix< ctype, n, n > &A, FieldMatrix< ctype, n, n > &ret)
Definition: affinegeometry.hh:226
ctype volume() const
Obtain the volume of the element.
Definition: affinegeometry.hh:599
static void AB(const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, n, p > &B, FieldMatrix< ctype, m, p > &ret)
Definition: affinegeometry.hh:66
static void ATBT(const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, p, m > &B, FieldMatrix< ctype, n, p > &ret)
Definition: affinegeometry.hh:80
static void invLTx(FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x)
Definition: affinegeometry.hh:293
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type for the transposed inverse Jacobian matrix.
Definition: affinegeometry.hh:483
static void LTx(const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret)
Definition: affinegeometry.hh:179
static void Ax(const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &ret)
Definition: affinegeometry.hh:44
LocalCoordinate local(const GlobalCoordinate &global) const
Evaluate the inverse mapping.
Definition: affinegeometry.hh:575
FieldVector< ctype, coorddimension > GlobalCoordinate
Type for coordinate vector in world space.
Definition: affinegeometry.hh:477
GlobalCoordinate center() const
Obtain the centroid of the mapping's image.
Definition: affinegeometry.hh:547
Definition: affinegeometry.hh:39
static void invLx(FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x)
Definition: affinegeometry.hh:281
int corners() const
Obtain number of corners of the corresponding reference element.
Definition: affinegeometry.hh:538
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type for the transposed Jacobian matrix.
Definition: affinegeometry.hh:480
AffineGeometry(Dune::GeometryType gt, const GlobalCoordinate &origin, const JacobianTransposed &jt)
Create affine geometry from GeometryType, one vertex, and the Jacobian matrix.
Definition: affinegeometry.hh:504
static void ATA(const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret)
Definition: affinegeometry.hh:108
static ctype detL(const FieldMatrix< ctype, n, n > &L)
Definition: affinegeometry.hh:250
static ctype spdInvA(FieldMatrix< ctype, n, n > &A)
Definition: affinegeometry.hh:313
static void spdInvAx(FieldMatrix< ctype, n, n > &A, FieldVector< ctype, n > &x)
Definition: affinegeometry.hh:324
static void Lx(const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret)
Definition: affinegeometry.hh:168
static void AAT(const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret)
Definition: affinegeometry.hh:150
AffineGeometry(const ReferenceElement &refElement, const CoordVector &coordVector)
Create affine geometry from reference element and a vector of vertex coordinates. ...
Definition: affinegeometry.hh:513
ct ctype
Definition: affinegeometry.hh:41