3#ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH
4#define DUNE_GEOMETRY_AFFINEGEOMETRY_HH
27 template<
typename Implementation >
30 template<
class ctype,
int dim >
31 class ReferenceElementImplementation;
33 template<
class ctype,
int dim >
34 struct ReferenceElements;
46 struct FieldMatrixHelper
50 template<
int m,
int n >
51 static void Ax (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, n > &x, FieldVector< ctype, m > &ret )
53 for(
int i = 0; i < m; ++i )
55 ret[ i ] = ctype( 0 );
56 for(
int j = 0; j < n; ++j )
57 ret[ i ] += A[ i ][ j ] * x[ j ];
61 template<
int m,
int n >
62 static void ATx (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, m > &x, FieldVector< ctype, n > &ret )
64 for(
int i = 0; i < n; ++i )
66 ret[ i ] = ctype( 0 );
67 for(
int j = 0; j < m; ++j )
68 ret[ i ] += A[ j ][ i ] * x[ j ];
72 template<
int m,
int n,
int p >
73 static void AB (
const FieldMatrix< ctype, m, n > &A,
const FieldMatrix< ctype, n, p > &B, FieldMatrix< ctype, m, p > &ret )
75 for(
int i = 0; i < m; ++i )
77 for(
int j = 0; j < p; ++j )
79 ret[ i ][ j ] = ctype( 0 );
80 for(
int k = 0; k < n; ++k )
81 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
86 template<
int m,
int n,
int p >
87 static void ATBT (
const FieldMatrix< ctype, m, n > &A,
const FieldMatrix< ctype, p, m > &B, FieldMatrix< ctype, n, p > &ret )
89 for(
int i = 0; i < n; ++i )
91 for(
int j = 0; j < p; ++j )
93 ret[ i ][ j ] = ctype( 0 );
94 for(
int k = 0; k < m; ++k )
95 ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
100 template<
int m,
int n >
101 static void ATA_L (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
103 for(
int i = 0; i < n; ++i )
105 for(
int j = 0; j <= i; ++j )
107 ret[ i ][ j ] = ctype( 0 );
108 for(
int k = 0; k < m; ++k )
109 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
114 template<
int m,
int n >
115 static void ATA (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
117 for(
int i = 0; i < n; ++i )
119 for(
int j = 0; j <= i; ++j )
121 ret[ i ][ j ] = ctype( 0 );
122 for(
int k = 0; k < m; ++k )
123 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
124 ret[ j ][ i ] = ret[ i ][ j ];
127 ret[ i ][ i ] = ctype( 0 );
128 for(
int k = 0; k < m; ++k )
129 ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
133 template<
int m,
int n >
134 static void AAT_L (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
144 for(
int i = 0; i < m; ++i )
146 for(
int j = 0; j <= i; ++j )
148 ctype &retij = ret[ i ][ j ];
149 retij = A[ i ][ 0 ] * A[ j ][ 0 ];
150 for(
int k = 1; k < n; ++k )
151 retij += A[ i ][ k ] * A[ j ][ k ];
156 template<
int m,
int n >
157 static void AAT (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
159 for(
int i = 0; i < m; ++i )
161 for(
int j = 0; j < i; ++j )
163 ret[ i ][ j ] = ctype( 0 );
164 for(
int k = 0; k < n; ++k )
165 ret[ i ][ j ] += A[ i ][ k ] * A[ j ][ k ];
166 ret[ j ][ i ] = ret[ i ][ j ];
168 ret[ i ][ i ] = ctype( 0 );
169 for(
int k = 0; k < n; ++k )
170 ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
175 static void Lx (
const FieldMatrix< ctype, n, n > &L,
const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
177 for(
int i = 0; i < n; ++i )
179 ret[ i ] = ctype( 0 );
180 for(
int j = 0; j <= i; ++j )
181 ret[ i ] += L[ i ][ j ] * x[ j ];
186 static void LTx (
const FieldMatrix< ctype, n, n > &L,
const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
188 for(
int i = 0; i < n; ++i )
190 ret[ i ] = ctype( 0 );
191 for(
int j = i; j < n; ++j )
192 ret[ i ] += L[ j ][ i ] * x[ j ];
197 static void LTL (
const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
199 for(
int i = 0; i < n; ++i )
201 for(
int j = 0; j < i; ++j )
203 ret[ i ][ j ] = ctype( 0 );
204 for(
int k = i; k < n; ++k )
205 ret[ i ][ j ] += L[ k ][ i ] * L[ k ][ j ];
206 ret[ j ][ i ] = ret[ i ][ j ];
208 ret[ i ][ i ] = ctype( 0 );
209 for(
int k = i; k < n; ++k )
210 ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
215 static void LLT (
const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
217 for(
int i = 0; i < n; ++i )
219 for(
int j = 0; j < i; ++j )
221 ret[ i ][ j ] = ctype( 0 );
222 for(
int k = 0; k <= j; ++k )
223 ret[ i ][ j ] += L[ i ][ k ] * L[ j ][ k ];
224 ret[ j ][ i ] = ret[ i ][ j ];
226 ret[ i ][ i ] = ctype( 0 );
227 for(
int k = 0; k <= i; ++k )
228 ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
233 static bool cholesky_L (
const FieldMatrix< ctype, n, n > &A, FieldMatrix< ctype, n, n > &ret,
const bool checkSingular =
false )
235 for(
int i = 0; i < n; ++i )
237 ctype &rii = ret[ i ][ i ];
239 ctype xDiag = A[ i ][ i ];
240 for(
int j = 0; j < i; ++j )
241 xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
245 if( checkSingular && ! ( xDiag > ctype( 0 )) )
249 assert( xDiag > ctype( 0 ) );
252 ctype invrii = ctype( 1 ) / rii;
253 for(
int k = i+1; k < n; ++k )
255 ctype x = A[ k ][ i ];
256 for(
int j = 0; j < i; ++j )
257 x -= ret[ i ][ j ] * ret[ k ][ j ];
258 ret[ k ][ i ] = invrii * x;
267 static ctype detL (
const FieldMatrix< ctype, n, n > &L )
270 for(
int i = 0; i < n; ++i )
276 static ctype invL ( FieldMatrix< ctype, n, n > &L )
279 for(
int i = 0; i < n; ++i )
281 ctype &lii = L[ i ][ i ];
283 lii = ctype( 1 ) / lii;
284 for(
int j = 0; j < i; ++j )
286 ctype &lij = L[ i ][ j ];
287 ctype x = lij * L[ j ][ j ];
288 for(
int k = j+1; k < i; ++k )
289 x += L[ i ][ k ] * L[ k ][ j ];
298 static void invLx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
300 for(
int i = 0; i < n; ++i )
302 for(
int j = 0; j < i; ++j )
303 x[ i ] -= L[ i ][ j ] * x[ j ];
304 x[ i ] /= L[ i ][ i ];
310 static void invLTx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
312 for(
int i = n; i > 0; --i )
314 for(
int j = i; j < n; ++j )
315 x[ i-1 ] -= L[ j ][ i-1 ] * x[ j ];
316 x[ i-1 ] /= L[ i-1 ][ i-1 ];
321 static ctype spdDetA (
const FieldMatrix< ctype, n, n > &A )
324 FieldMatrix< ctype, n, n > L;
330 static ctype spdInvA ( FieldMatrix< ctype, n, n > &A )
332 FieldMatrix< ctype, n, n > L;
334 const ctype det = invL( L );
341 static bool spdInvAx ( FieldMatrix< ctype, n, n > &A, FieldVector< ctype, n > &x,
const bool checkSingular =
false )
343 FieldMatrix< ctype, n, n > L;
344 const bool invertible = cholesky_L( A, L, checkSingular );
345 if( ! invertible )
return invertible ;
351 template<
int m,
int n >
352 static ctype detATA (
const FieldMatrix< ctype, m, n > &A )
356 FieldMatrix< ctype, n, n > ata;
358 return spdDetA( ata );
369 template<
int m,
int n >
370 static ctype sqrtDetAAT (
const FieldMatrix< ctype, m, n > &A )
377 if( (n == 2) && (m == 2) )
380 return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
382 else if( (n == 3) && (m == 3) )
385 const ctype v0 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 1 ][ 1 ] * A[ 0 ][ 2 ];
386 const ctype v1 = A[ 0 ][ 2 ] * A[ 1 ][ 0 ] - A[ 1 ][ 2 ] * A[ 0 ][ 0 ];
387 const ctype v2 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 1 ][ 0 ] * A[ 0 ][ 1 ];
388 return abs( v0 * A[ 2 ][ 0 ] + v1 * A[ 2 ][ 1 ] + v2 * A[ 2 ][ 2 ] );
390 else if ( (n == 3) && (m == 2) )
393 const ctype v0 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 0 ][ 1 ] * A[ 1 ][ 0 ];
394 const ctype v1 = A[ 0 ][ 0 ] * A[ 1 ][ 2 ] - A[ 1 ][ 0 ] * A[ 0 ][ 2 ];
395 const ctype v2 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 0 ][ 2 ] * A[ 1 ][ 1 ];
396 return sqrt( v0*v0 + v1*v1 + v2*v2);
401 FieldMatrix< ctype, m, m > aat;
403 return spdDetA( aat );
411 template<
int m,
int n >
412 static ctype leftInvA (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
414 static_assert((m >= n),
"Matrix has no left inverse.");
415 FieldMatrix< ctype, n, n > ata;
417 const ctype det = spdInvA( ata );
422 template<
int m,
int n >
423 static void leftInvAx (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, m > &x, FieldVector< ctype, n > &y )
425 static_assert((m >= n),
"Matrix has no left inverse.");
426 FieldMatrix< ctype, n, n > ata;
433 template<
int m,
int n >
434 static ctype rightInvA (
const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
436 static_assert((n >= m),
"Matrix has no right inverse.");
438 if( (n == 2) && (m == 2) )
440 const ctype det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
441 const ctype detInv = ctype( 1 ) / det;
442 ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
443 ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
444 ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
445 ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
450 FieldMatrix< ctype, m , m > aat;
452 const ctype det = spdInvA( aat );
453 ATBT( A , aat , ret );
458 template<
int m,
int n >
459 static bool xTRightInvA (
const FieldMatrix< ctype, m, n > &A,
const FieldVector< ctype, n > &x, FieldVector< ctype, m > &y )
461 static_assert((n >= m),
"Matrix has no right inverse.");
462 FieldMatrix< ctype, m, m > aat;
466 return spdInvAx( aat, y,
true );
479 template<
class ct,
int mydim,
int cdim>
515 typedef Impl::FieldMatrixHelper< ct > MatrixHelper;
521 : refElement_(refElement), origin_(origin), jacobianTransposed_(jt)
523 integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
533 template<
class CoordVector >
535 : refElement_(refElement), origin_(coordVector[0])
538 jacobianTransposed_[ i ] = coordVector[ i+1 ] - origin_;
539 integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
543 template<
class CoordVector >
595 jacobianInverseTransposed_.
mtv(
global - origin_, local );
612 return integrationElement_;
618 return integrationElement_ * refElement_.
volume();
630 return jacobianTransposed_;
642 return jacobianInverseTransposed_;
647 return geometry.refElement_;
651 ReferenceElement refElement_;
655 ctype integrationElement_;
Implementation of the Geometry interface for affine geometries.
Definition: affinegeometry.hh:481
AffineGeometry(const ReferenceElement &refElement, const CoordVector &coordVector)
Create affine geometry from reference element and a vector of vertex coordinates.
Definition: affinegeometry.hh:534
AffineGeometry(Dune::GeometryType gt, const GlobalCoordinate &origin, const JacobianTransposed &jt)
Create affine geometry from GeometryType, one vertex, and the Jacobian matrix.
Definition: affinegeometry.hh:527
FieldVector< ctype, mydimension > LocalCoordinate
Type for local coordinate vector.
Definition: affinegeometry.hh:494
Dune::GeometryType type() const
Obtain the type of the reference element.
Definition: affinegeometry.hh:552
static const int mydimension
Dimension of the geometry.
Definition: affinegeometry.hh:488
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:519
ctype Volume
Type used for volume.
Definition: affinegeometry.hh:500
AffineGeometry(Dune::GeometryType gt, const CoordVector &coordVector)
Create affine geometry from GeometryType and a vector of vertex coordinates.
Definition: affinegeometry.hh:544
ctype integrationElement(const LocalCoordinate &local) const
Obtain the integration element.
Definition: affinegeometry.hh:609
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian's inverse.
Definition: affinegeometry.hh:639
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type for the transposed Jacobian matrix.
Definition: affinegeometry.hh:503
GlobalCoordinate corner(int i) const
Obtain coordinates of the i-th corner.
Definition: affinegeometry.hh:558
int corners() const
Obtain number of corners of the corresponding reference element.
Definition: affinegeometry.hh:555
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type for the transposed inverse Jacobian matrix.
Definition: affinegeometry.hh:506
static const int coorddimension
Dimension of the world space.
Definition: affinegeometry.hh:491
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the mapping.
Definition: affinegeometry.hh:572
GlobalCoordinate center() const
Obtain the centroid of the mapping's image.
Definition: affinegeometry.hh:564
ct ctype
Type used for coordinates.
Definition: affinegeometry.hh:485
FieldVector< ctype, coorddimension > GlobalCoordinate
Type for coordinate vector in world space.
Definition: affinegeometry.hh:497
bool affine() const
Always true: this is an affine geometry.
Definition: affinegeometry.hh:549
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian.
Definition: affinegeometry.hh:627
Volume volume() const
Obtain the volume of the element.
Definition: affinegeometry.hh:616
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:414
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:446
vector space out of a tensor product of fields.
Definition: fvector.hh:96
This class provides access to geometric and topological properties of a reference element.
Definition: referenceelement.hh:56
CoordinateField volume() const
obtain the volume of the reference element
Definition: referenceelement.hh:245
decltype(auto) type(int i, int c) const
obtain the type of subentity (i,c)
Definition: referenceelement.hh:175
int size(int c) const
number of subentities of codimension c
Definition: referenceelement.hh:98
decltype(auto) position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelement.hh:207
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:279
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
unspecified-type ReferenceElement
Returns the type of reference element for the argument type T.
Definition: referenceelements.hh:495
Dune namespace.
Definition: alignedallocator.hh:14
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:168
A unique label for each type of element that can occur in a grid.