dune-fem  2.8-git
simplegeometry.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_GRIDPART_COMMON_SIMPLEGEOMETRY_HH
2 #define DUNE_FEM_GRIDPART_COMMON_SIMPLEGEOMETRY_HH
3 
4 #include <type_traits>
5 #include <utility>
6 
7 #include <dune/common/fmatrix.hh>
8 #include <dune/common/fvector.hh>
9 
10 #include <dune/geometry/affinegeometry.hh>
11 #include <dune/geometry/referenceelements.hh>
12 
13 #include <dune/grid/common/geometry.hh>
14 
16 
17 namespace Dune
18 {
19 
20  // SimpleGeometry
21  // --------------
22 
23  template< class BasicGeometry >
25  : public BasicGeometry
26  {
27  typedef typename BasicGeometry::ctype ctype;
28 
29  using BasicGeometry::mydimension;
30  using BasicGeometry::coorddimension;
31 
32  using BasicGeometry::global;
33  using BasicGeometry::jacobianTransposed;
34  using BasicGeometry::quadrature;
35  using BasicGeometry::type;
36 
37  typedef FieldVector< ctype, mydimension > LocalCoordinate;
38  typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
39 
40  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
41  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
43 
44  // Helper class to compute a matrix pseudo inverse
45  typedef Impl::FieldMatrixHelper< ctype > MatrixHelper;
46 
47  template< class... Args, std::enable_if_t< std::is_constructible< BasicGeometry, Args &&... >::value, int > = 0 >
48  SimpleGeometry ( Args &&...args )
49  : BasicGeometry( std::forward< Args >( args )... )
50  {}
51 
52  int corners () const { return referenceElement().size( mydimension ); }
53  GlobalCoordinate corner ( int i ) const { return global( referenceElement().position( i, mydimension ) ); }
54 
55  bool affine () const { return false; }
56 
57  LocalCoordinate local ( const GlobalCoordinate &global ) const
58  {
59  const ctype tolerance = 1e-12; // use something better here e.g. Traits::tolerance();
60  LocalCoordinate local = referenceElement().position( 0, 0 );
61  LocalCoordinate dlocal;
62  do
63  {
64  // Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n
65  const GlobalCoordinate dglobal = this->global( local ) - global;
66  MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed( local ), dglobal, dlocal );
67  local -= dlocal;
68  assert( referenceElement().checkInside( local ) );
69  }
70  while( dlocal.two_norm2() > tolerance );
71  return local;
72  }
73 
74  ctype integrationElement ( const LocalCoordinate &local ) const
75  {
76  return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jacobianTransposed( local ) );
77  }
78 
80  {
81  JacobianInverseTransposed jacInverseTransposed( 0 );
82  MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed( local ), jacInverseTransposed );
83  return jacInverseTransposed;
84  }
85 
87  {
89  ctype volume( 0 );
90  for( const auto &qp : quadrature( 0 ) )
91  {
92  const ctype weight = qp.weight() * integrationElement( qp.position() );
93  center.axpy( weight, global( qp.position() ) );
94  volume += weight;
95  }
96  return center /= volume;
97  }
98 
99  ctype volume () const
100  {
101  ctype volume( 0 );
102  for( const auto &qp : quadrature( 0 ) )
103  volume += qp.weight() * integrationElement( qp.position() );
104  return volume;
105  }
106 
107  auto referenceElement () const { return Dune::referenceElement<double>( type(), Dune::Dim<mydimension>() ); }
108  };
109 
110 } // namespace Dune
111 
112 #endif // #ifndef DUNE_FEM_GRIDPART_COMMON_SIMPLEGEOMETRY_HH
Definition: bindguard.hh:11
Definition: simplegeometry.hh:26
FieldVector< ctype, mydimension > LocalCoordinate
Definition: simplegeometry.hh:37
SimpleGeometry(Args &&...args)
Definition: simplegeometry.hh:48
int corners() const
Definition: simplegeometry.hh:52
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Definition: simplegeometry.hh:41
FieldVector< ctype, coorddimension > GlobalCoordinate
Definition: simplegeometry.hh:38
ctype volume() const
Definition: simplegeometry.hh:99
JacobianInverseTransposed Jacobian
Definition: simplegeometry.hh:42
ctype integrationElement(const LocalCoordinate &local) const
Definition: simplegeometry.hh:74
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Definition: simplegeometry.hh:79
BasicGeometry::ctype ctype
Definition: simplegeometry.hh:27
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Definition: simplegeometry.hh:40
auto referenceElement() const
Definition: simplegeometry.hh:107
Impl::FieldMatrixHelper< ctype > MatrixHelper
Definition: simplegeometry.hh:45
GlobalCoordinate center() const
Definition: simplegeometry.hh:86
bool affine() const
Definition: simplegeometry.hh:55
GlobalCoordinate corner(int i) const
Definition: simplegeometry.hh:53