DUNE-FEM (unstable)

checkbasisfunctionset.hh
1#ifndef DUNE_FEM_SPACE_TEST_CHECKBASISFUNCTIONSET_HH
2#define DUNE_FEM_SPACE_TEST_CHECKBASISFUNCTIONSET_HH
3
4#include <algorithm>
5#include <functional>
6#include <iostream>
7#include <random>
8
10#include <dune/fem/common/coordinate.hh>
11
12namespace Dune
13{
14
15 namespace Fem
16 {
17
18 // checkQuadratureConsistency
19 // --------------------------
20
21 template< class BasisFunctionSet, class Quadrature >
23 checkQuadratureConsistency ( const BasisFunctionSet &basisFunctionSet, const Quadrature &quadrature,
24 bool supportsHessians = false )
25 {
26 // get types
28 typedef typename BasisFunctionSet::DomainType DomainType;
29 typedef typename BasisFunctionSet::RangeType RangeType;
30 typedef typename BasisFunctionSet::JacobianRangeType JacobianRangeType;
31 typedef typename BasisFunctionSet::HessianRangeType HessianRangeType;
32 typedef typename BasisFunctionSet::EntityType EntityType;
33 typedef typename BasisFunctionSet::Geometry Geometry;
34
35 typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
36
37 static const int dimRange = FunctionSpaceType::dimRange;
38
39 const EntityType &entity = basisFunctionSet.entity();
40 const Geometry& geometry = basisFunctionSet.geometry();
41 const auto refElement = basisFunctionSet.referenceElement();
42
43 if( entity.type() != refElement.type() || geometry.type() != refElement.type() )
44 DUNE_THROW( Dune::InvalidStateException, "GeometryType of referenceElement and entity mismatch for this basisFunctionSet" );
45
46 int order = basisFunctionSet.order();
47 // prevent warning about unused params
48 (void) order;
49
50 const std::size_t size = basisFunctionSet.size();
51
52 // init random dof vectors
53 std::uniform_real_distribution< RangeFieldType > distribution( 1e-3, 1e3 );
54 std::default_random_engine randomEngine;
55 auto random = std::bind( distribution, randomEngine );
56
57 std::vector< RangeFieldType > dofs( size );
58 for( RangeFieldType &dof : dofs )
59 dof = random();
60
61 RangeType valueFactor;
62 for( RangeFieldType &v : valueFactor )
63 v = random();
64
65 JacobianRangeType jacobianFactor;
66 for( std::size_t j = 0; j < JacobianRangeType::rows; ++j )
67 for( std::size_t k = 0; k < JacobianRangeType::cols; ++k )
68 jacobianFactor[ j ][ k ] = random();
69
70 HessianRangeType hessianFactor;
71 for( std::size_t j = 0; j < JacobianRangeType::rows; ++j )
72 for( std::size_t k = 0; k < JacobianRangeType::cols; ++k )
73 for( std::size_t l = 0; l < JacobianRangeType::cols; ++l )
74 hessianFactor[ j ][ k ][ l ] = random();
75
76 // return value
78
79 const std::size_t nop = quadrature.nop();
80
81 std::vector< RangeType > aVec( nop );
82 std::vector< JacobianRangeType > aJac( nop );
83
84 basisFunctionSet.evaluateAll( quadrature, dofs, aVec );
85 basisFunctionSet.jacobianAll( quadrature, dofs, aJac );
86
87 for( std::size_t qp = 0; qp < nop; ++qp )
88 {
89 // check evaluate methods
90 {
91 RangeType b;
92 RangeType& a = aVec[ qp ];
93
94 std::vector< RangeType > values( basisFunctionSet.size() );
95 basisFunctionSet.evaluateAll( quadrature[ qp ], values );
96 for( std::size_t i = 0; i < values.size(); ++i )
97 {
98 b.axpy( dofs[ i ], values[ i ] );
99 }
100
101 ret[ 0 ] = std::max( ret[ 0 ], (a - b).two_norm() );
102 }
103
104 // check jacobian methods
105 {
106 JacobianRangeType b;
107 JacobianRangeType& a = aJac[ qp ];
108
109 std::vector< JacobianRangeType > values( basisFunctionSet.size() );
110 basisFunctionSet.jacobianAll( quadrature[ qp ], values );
111 for( std::size_t i = 0; i < values.size(); ++i )
112 b.axpy( dofs[ i ], values[ i ] );
113
114 a -= b;
115 ret[ 1 ] = std::max( ret[ 1 ], a.frobenius_norm() );
116 }
117
118 // check hessian methods
119 if( supportsHessians )
120 {
121 HessianRangeType a, b;
122
123 basisFunctionSet.hessianAll( quadrature[ qp ], dofs, a );
124
125 std::vector< HessianRangeType > values( basisFunctionSet.size() );
126 basisFunctionSet.hessianAll( quadrature[ qp ], values );
127 for( int r = 0; r < dimRange; ++r )
128 for( std::size_t i = 0; i < values.size(); ++i )
129 b[ r ].axpy( dofs[ i ], values[ i ][ r ] );
130
131 RangeFieldType error( 0 );
132 for( int r = 0; r < dimRange; ++r )
133 a[ r ] -= b[ r ];
134 for( int r = 0; r < dimRange; ++r )
135 error += a[ r ].frobenius_norm2();
136
137 ret[ 2 ] = std::max( ret[ 2 ], std::sqrt( error ) );
138 }
139
140 // check value axpy method
141 {
142 std::vector< RangeFieldType > r1( dofs );
143 std::vector< RangeFieldType > r2( dofs );
144
145 basisFunctionSet.axpy( quadrature[ qp ], valueFactor, r1 );
146
147 std::vector< RangeType > values( basisFunctionSet.size() );
148 basisFunctionSet.evaluateAll( quadrature[ qp ], values );
149 for( std::size_t i = 0; i < values.size(); ++i )
150 r2[ i ] += valueFactor * values[ i ];
151
152 RangeFieldType error = 0;
153 for( std::size_t i = 0; i < values.size(); ++i )
154 error += std::abs( r2[ i ] - r1[ i ] );
155
156 ret[ 3 ] = std::max( ret[ 3 ], error );
157 }
158
159 // check jacobian axpy method
160 {
161 DomainType x = coordinate( quadrature[ qp ] );
162 std::vector< RangeFieldType > r1( dofs );
163 std::vector< RangeFieldType > r2( dofs );
164
165 basisFunctionSet.axpy( x, jacobianFactor, r1 );
166
167 std::vector< JacobianRangeType > values( basisFunctionSet.size() );
168 basisFunctionSet.jacobianAll( x, values );
169 for( std::size_t i = 0; i < values.size(); ++i )
170 for( int j = 0; j < JacobianRangeType::rows; ++j )
171 r2[ i ] += jacobianFactor[ j ] * values[ i ][ j ];
172
173 RangeFieldType error = 0;
174 for( std::size_t i = 0; i < values.size(); ++i )
175 error += std::abs( r2[ i ] - r1[ i ] );
176
177 ret[ 4 ] = std::max( ret[ 4 ], error );
178 }
179
180
181 // check hessian axpy method
182 if( supportsHessians )
183 {
184 DomainType x = coordinate( quadrature[ qp ] );
185 std::vector< RangeFieldType > r1( dofs );
186 std::vector< RangeFieldType > r2( dofs );
187
188 basisFunctionSet.axpy( x, hessianFactor, r1 );
189
190 std::vector< HessianRangeType > values( basisFunctionSet.size() );
191 basisFunctionSet.hessianAll( x, values );
192 for( std::size_t i = 0; i < values.size(); ++i )
193 for( int j = 0; j < JacobianRangeType::rows; ++j )
194 for( std::size_t k = 0; k < JacobianRangeType::cols; ++k )
195 for( std::size_t l = 0; l < JacobianRangeType::cols; ++l )
196 r2[ i ] += hessianFactor[ j ][ k ][ l ] * values[ i ][ j ][ k ][ l ];
197
198 RangeFieldType error = 0;
199 for( std::size_t i = 0; i < values.size(); ++i )
200 error += std::abs( r2[ i ] - r1[ i ] );
201
202 ret[ 5 ] = std::max( ret[ 5 ], error );
203 }
204
205 // check value, jaacobian axpy method
206 {
207 DomainType x = coordinate( quadrature[ qp ] );
208 std::vector< RangeFieldType > r1( dofs );
209 std::vector< RangeFieldType > r2( dofs );
210
211 basisFunctionSet.axpy( x, valueFactor, jacobianFactor, r1 );
212
213 std::vector< RangeType > values( basisFunctionSet.size() );
214 std::vector< JacobianRangeType > jacobians( basisFunctionSet.size() );
215
216 basisFunctionSet.evaluateAll( x, values );
217 basisFunctionSet.jacobianAll( x, jacobians );
218
219 for( std::size_t i = 0; i < values.size(); ++i )
220 {
221 r2[ i ] += valueFactor * values[ i ];
222 for( int j = 0; j < JacobianRangeType::rows; ++j )
223 r2[ i ] += jacobianFactor[ j ] * jacobians[ i ][ j ];
224 }
225
226 RangeFieldType error = 0;
227 for( std::size_t i = 0; i < values.size(); ++i )
228 error += std::abs( r2[ i ] - r1[ i ] );
229
230 ret[ 6 ] = std::max( ret[ 6 ], error );
231 }
232 }
233
234 return ret;
235 }
236
237 } // namespace Fem
238
239} // namespace Dune
240
241#endif // #ifndef DUNE_FEM_SPACE_TEST_CHECKBASISFUNCTIONSET_HH
Entity EntityType
entity type
Definition: basisfunctionset.hh:35
FunctionSpace< typename Entity::Geometry::ctype, typename Range::value_type, Entity::Geometry::coorddimension, Range::dimension > FunctionSpaceType
function space type
Definition: basisfunctionset.hh:40
FunctionSpaceType::DomainType DomainType
range type
Definition: basisfunctionset.hh:43
FunctionSpaceType::JacobianRangeType JacobianRangeType
jacobian range type
Definition: basisfunctionset.hh:47
FunctionSpaceType::HessianRangeType HessianRangeType
hessian range type
Definition: basisfunctionset.hh:49
FunctionSpaceType::RangeType RangeType
range type
Definition: basisfunctionset.hh:45
@ dimRange
dimension of range vector space
Definition: functionspaceinterface.hh:48
FunctionSpaceTraits::RangeFieldType RangeFieldType
Intrinsic type used for values in the range field (usually a double)
Definition: functionspaceinterface.hh:63
A vector valued function space.
Definition: functionspace.hh:60
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
actual interface class for quadratures
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)