1#ifndef DUNE_FEM_SPACE_TEST_CHECKBASISFUNCTIONSET_HH
2#define DUNE_FEM_SPACE_TEST_CHECKBASISFUNCTIONSET_HH
10#include <dune/fem/common/coordinate.hh>
21 template<
class BasisFunctionSet,
class Quadrature >
23 checkQuadratureConsistency (
const BasisFunctionSet &basisFunctionSet,
const Quadrature &quadrature,
24 bool supportsHessians =
false )
33 typedef typename BasisFunctionSet::Geometry Geometry;
39 const EntityType &entity = basisFunctionSet.entity();
40 const Geometry& geometry = basisFunctionSet.geometry();
41 const auto refElement = basisFunctionSet.referenceElement();
43 if( entity.type() != refElement.type() || geometry.type() != refElement.type() )
46 int order = basisFunctionSet.order();
50 const std::size_t
size = basisFunctionSet.size();
53 std::uniform_real_distribution< RangeFieldType > distribution( 1e-3, 1e3 );
54 std::default_random_engine randomEngine;
55 auto random = std::bind( distribution, randomEngine );
57 std::vector< RangeFieldType > dofs(
size );
58 for( RangeFieldType &dof : dofs )
61 RangeType valueFactor;
62 for( RangeFieldType &v : valueFactor )
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();
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();
79 const std::size_t nop = quadrature.nop();
81 std::vector< RangeType > aVec( nop );
82 std::vector< JacobianRangeType > aJac( nop );
84 basisFunctionSet.evaluateAll( quadrature, dofs, aVec );
85 basisFunctionSet.jacobianAll( quadrature, dofs, aJac );
87 for( std::size_t qp = 0; qp < nop; ++qp )
92 RangeType& a = aVec[ qp ];
94 std::vector< RangeType > values( basisFunctionSet.size() );
95 basisFunctionSet.evaluateAll( quadrature[ qp ], values );
96 for( std::size_t i = 0; i < values.size(); ++i )
98 b.axpy( dofs[ i ], values[ i ] );
101 ret[ 0 ] = std::max( ret[ 0 ], (a - b).two_norm() );
107 JacobianRangeType& a = aJac[ qp ];
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 ] );
115 ret[ 1 ] = std::max( ret[ 1 ], a.frobenius_norm() );
119 if( supportsHessians )
121 HessianRangeType a, b;
123 basisFunctionSet.hessianAll( quadrature[ qp ], dofs, a );
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 ] );
131 RangeFieldType error( 0 );
132 for(
int r = 0; r < dimRange; ++r )
134 for(
int r = 0; r < dimRange; ++r )
135 error += a[ r ].frobenius_norm2();
137 ret[ 2 ] = std::max( ret[ 2 ], std::sqrt( error ) );
142 std::vector< RangeFieldType > r1( dofs );
143 std::vector< RangeFieldType > r2( dofs );
145 basisFunctionSet.axpy( quadrature[ qp ], valueFactor, r1 );
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 ];
152 RangeFieldType error = 0;
153 for( std::size_t i = 0; i < values.size(); ++i )
154 error += std::abs( r2[ i ] - r1[ i ] );
156 ret[ 3 ] = std::max( ret[ 3 ], error );
161 DomainType x = coordinate( quadrature[ qp ] );
162 std::vector< RangeFieldType > r1( dofs );
163 std::vector< RangeFieldType > r2( dofs );
165 basisFunctionSet.axpy( x, jacobianFactor, r1 );
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 ];
173 RangeFieldType error = 0;
174 for( std::size_t i = 0; i < values.size(); ++i )
175 error += std::abs( r2[ i ] - r1[ i ] );
177 ret[ 4 ] = std::max( ret[ 4 ], error );
182 if( supportsHessians )
184 DomainType x = coordinate( quadrature[ qp ] );
185 std::vector< RangeFieldType > r1( dofs );
186 std::vector< RangeFieldType > r2( dofs );
188 basisFunctionSet.axpy( x, hessianFactor, r1 );
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 ];
198 RangeFieldType error = 0;
199 for( std::size_t i = 0; i < values.size(); ++i )
200 error += std::abs( r2[ i ] - r1[ i ] );
202 ret[ 5 ] = std::max( ret[ 5 ], error );
207 DomainType x = coordinate( quadrature[ qp ] );
208 std::vector< RangeFieldType > r1( dofs );
209 std::vector< RangeFieldType > r2( dofs );
211 basisFunctionSet.axpy( x, valueFactor, jacobianFactor, r1 );
213 std::vector< RangeType > values( basisFunctionSet.size() );
214 std::vector< JacobianRangeType > jacobians( basisFunctionSet.size() );
216 basisFunctionSet.evaluateAll( x, values );
217 basisFunctionSet.jacobianAll( x, jacobians );
219 for( std::size_t i = 0; i < values.size(); ++i )
221 r2[ i ] += valueFactor * values[ i ];
222 for(
int j = 0; j < JacobianRangeType::rows; ++j )
223 r2[ i ] += jacobianFactor[ j ] * jacobians[ i ][ j ];
226 RangeFieldType error = 0;
227 for( std::size_t i = 0; i < values.size(); ++i )
228 error += std::abs( r2[ i ] - r1[ i ] );
230 ret[ 6 ] = std::max( ret[ 6 ], error );
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