dune-fem 2.8.0
Loading...
Searching...
No Matches
agglomerationquadrature.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_AGGLOMERATIONQUADRATURE_HH
2#define DUNE_FEM_AGGLOMERATIONQUADRATURE_HH
3
4#include <stack>
5#include <dune/geometry/axisalignedcubegeometry.hh>
7
8#include "quadrature.hh"
10
11namespace Dune
12{
13
14 namespace Fem
15 {
16
17
20 template< typename GridPartImp, class IntegrationPointList >
22 {
23 public:
25 typedef GridPartImp GridPartType;
26
28 enum { codimension = 0 };
29
31 enum { dimension = GridPartType :: dimension };
32
34 enum { dimensionworld = GridPartType :: dimensionworld };
35
37 typedef typename GridPartType :: ctype RealType;
38
40
41 typedef typename IntegrationPointListType :: CoordinateType CoordinateType;
42
43 typedef int QuadratureKeyType;
44
45 // for compatibility
46 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
47
48 protected:
50 typedef typename IntegrationPointListType :: IntegrationPointListType IntegrationPointListImpl;
51
52 typedef std::stack< std::unique_ptr< PolyhedronQuadratureType > > PolyhedronQuadratureStorageType;
53
54 // PolyhedronQuadrature storage
56 {
58 return *storage;
59 }
60
61 // get object from stack or create new
62 static PolyhedronQuadratureType* getObject( const GeometryType& type )
63 {
65 if( storage.empty() )
66 {
67 return new PolyhedronQuadratureType( type, 0, IdProvider ::instance().newId() );
68 }
69 else
70 {
71 PolyhedronQuadratureType* quad = storage.top().release();
72 assert( quad );
73 storage.pop();
74 return quad;
75 }
76 }
77
78 // push object to stack or delete
79 static void pushObject( const PolyhedronQuadratureType* quad )
80 {
81 PolyhedronQuadratureType* polyQuad = const_cast< PolyhedronQuadratureType* > (quad);
83 if( storage.size() < 20 )
84 {
85 storage.emplace( polyQuad );
86 }
87 else
88 {
89 delete polyQuad;
90 }
91 }
92
93 // deleter object returning pointer to stack object
94 struct Deleter
95 {
97 {
98 const PolyhedronQuadratureType* polyQuad = static_cast< const PolyhedronQuadratureType* > ( quad );
99 pushObject( polyQuad );
100 }
101 };
102
103 public:
105 // this only works for 2d so far
106 template <int dimw = dimensionworld >
107 static std::enable_if_t<dimension == 2 && dimw == 2, IntegrationPointListType>
108 computeQuadrature( const EntityType &entity, const QuadratureKeyType& quadKey )
109 {
110 typedef ElementQuadrature< GridPartImp, 0 > QuadratureType;
111 Dune::GeometryType simplexType = Dune::GeometryTypes::simplex( dimension );
112
113 typedef AxisAlignedCubeGeometry< RealType, dimension, dimension> CubeGeometryType;
114
115 const auto& elemGeo = entity.geometry();
116
117 // compute bounding box for setting up a cube quadrature
118 CoordinateType lower = elemGeo.corner( 0 );
119 CoordinateType upper = lower;
120 const int corners = elemGeo.corners();
121 for( int c = 1; c < corners; ++c )
122 {
123 const auto& corner = elemGeo.corner( c );
124 for( int d=0; d< dimension; ++d )
125 {
126 lower[ d ] = std::min( lower[ d ], corner[ d ]);
127 upper[ d ] = std::max( upper[ d ], corner[ d ]);
128 }
129 }
130
131 //std::cout << "BBox = " << lower << " " << upper << std::endl;
132
133 CubeGeometryType cubeGeom( lower, upper );
134
135 QuadratureType quad( simplexType, quadKey );
136 // TODO needs generalization
137 const int subEntities = entity.subEntities( 1 );
138 const int quadNop = quad.nop();
139 int order = quad.order();
140
141 PolyhedronQuadratureType& quadImp = *(getObject( entity.type() ));
142
143 quadImp.reset( order, subEntities * quadNop );
144
145 Dune::FieldMatrix<double,dimension,dimension> A( 0 ) ;
146 const auto& center = entity.geometry().center();
147
148 for( int i = 0; i<subEntities; ++i )
149 {
150 const auto subEntity = entity.template subEntity< 1 >( i );
151 const auto& geom = subEntity.geometry();
152 assert( geom.corners() == dimension );
153 // setup transformation matrix, here setup transposed matrix
154 for( int c = 0; c<dimension; ++c )
155 {
156 A[ c ] = geom.corner( c );
157 A[ c ] -= center;
158 }
159
160 // compute simplex volume / ref volume (which removes the 0.5)
161 // abs is taken because determinant may be negative since we did not
162 // care about the orientation
163 double vol = std::abs(A.determinant()) / elemGeo.volume();
164
165 CoordinateType point;
166 for( int qp = 0; qp < quadNop; ++qp )
167 {
168 // p = A^T * xHat + p_0 (A is stored as transposed)
169 A.mtv( quad.point( qp ), point );
170 point += center;
171
172 point = cubeGeom.local( point );
173
174 // scale weights with number of sub-triangles
175 double weight = quad.weight( qp ) * vol;
176
177 quadImp.addQuadraturePoint( point, weight );
178 }
179 }
180 // return a shared pointer with the correct deleter
181 // removing the pointer to the stack
182 std::shared_ptr< const IntegrationPointListImpl > quadPtr( &quadImp, Deleter() );
183 return IntegrationPointListType( quadPtr );
184 }
185
186 template <int dimw = dimensionworld >
187 static std::enable_if_t<dimension != 2 || dimw != 2, IntegrationPointListType>
188 computeQuadrature( const EntityType &entity, const QuadratureKeyType& quadKey )
189 {
190 assert(false);
191 return IntegrationPointListType( {} );
192 }
193 };
194
195 } // namespace Fem
196
197} // namespace Dune
198
199#endif // #ifndef DUNE_FEM_ELEMENTQUADRATURE_HH
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
double min(const Dune::Fem::Double &v, const double p)
Definition: double.hh:953
Definition: bindguard.hh:11
ThreadSafeValue realizes thread safety for a given variable by creating an instance of this variable ...
Definition: threadsafevalue.hh:18
Agglomeration is a simple quadrature for polyhedral cells based on sub triangulation
Definition: agglomerationquadrature.hh:22
std::stack< std::unique_ptr< PolyhedronQuadratureType > > PolyhedronQuadratureStorageType
Definition: agglomerationquadrature.hh:52
PolyhedronQuadrature< RealType, dimension > PolyhedronQuadratureType
Definition: agglomerationquadrature.hh:49
static std::enable_if_t< dimension !=2||dimw !=2, IntegrationPointListType > computeQuadrature(const EntityType &entity, const QuadratureKeyType &quadKey)
Definition: agglomerationquadrature.hh:188
IntegrationPointListType::CoordinateType CoordinateType
Definition: agglomerationquadrature.hh:41
IntegrationPointListType::IntegrationPointListType IntegrationPointListImpl
Definition: agglomerationquadrature.hh:50
static void pushObject(const PolyhedronQuadratureType *quad)
Definition: agglomerationquadrature.hh:79
@ dimensionworld
Definition: agglomerationquadrature.hh:34
GridPartType::ctype RealType
type for reals (usually double)
Definition: agglomerationquadrature.hh:37
GridPartImp GridPartType
type of the grid partition
Definition: agglomerationquadrature.hh:25
@ codimension
Definition: agglomerationquadrature.hh:28
static std::enable_if_t< dimension==2 &&dimw==2, IntegrationPointListType > computeQuadrature(const EntityType &entity, const QuadratureKeyType &quadKey)
returns quadrature points for polyhedral cells
Definition: agglomerationquadrature.hh:108
static PolyhedronQuadratureStorageType & quadStorage()
Definition: agglomerationquadrature.hh:55
static PolyhedronQuadratureType * getObject(const GeometryType &type)
Definition: agglomerationquadrature.hh:62
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: agglomerationquadrature.hh:46
@ dimension
Definition: agglomerationquadrature.hh:31
int QuadratureKeyType
Definition: agglomerationquadrature.hh:43
IntegrationPointList IntegrationPointListType
Definition: agglomerationquadrature.hh:39
Definition: agglomerationquadrature.hh:95
void operator()(const IntegrationPointListImpl *quad)
Definition: agglomerationquadrature.hh:96
quadrature on the codim-0 reference element
Definition: elementquadrature.hh:58
Definition: femquadratures.hh:287
void reset(const int order, const int nop)
Definition: femquadratures.hh:338
void addQuadraturePoint(const CoordinateType &point, const FieldType weight)
Adds a point-weight pair to the quadrature.
Definition: quadratureimp.hh:270
actual interface class for integration point lists
Definition: quadrature.hh:161