11#ifndef DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
12#define DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
39 using RangeType =
typename FE::Traits::Basis::Traits::Range;
40 using DomainType =
typename FE::Traits::Basis::Traits::DomainLocal;
42 typedef typename FE::Traits::Basis::Traits::RangeField CT;
44 std::vector<CT> coeff;
46 FEFunction(
const FE& fe_) : fe(fe_) { resetCoefficients(); }
48 void resetCoefficients() {
49 coeff.resize(fe.basis().size());
50 for(std::size_t i=0; i<coeff.size(); ++i)
54 void setRandom(
double max) {
55 coeff.resize(fe.basis().size());
56 for(std::size_t i=0; i<coeff.size(); ++i)
57 coeff[i] = ((1.0*std::rand()) / RAND_MAX - 0.5)*2.0*
max;
60 RangeType operator() (
const DomainType& x)
const {
62 std::vector<RangeType> yy;
63 fe.basis().evaluateFunction(x, yy);
65 for (std::size_t i=0; i<yy.size(); ++i)
66 y.axpy(coeff[i], yy[i]);
86bool testInterpolation(
const FE& fe,
double eps,
int n=5)
91 std::vector<typename FEFunction<FE>::CT> coeff;
92 for(
int i=0; i<n && success; ++i) {
101 fe.interpolation().interpolate(f, coeff);
104 if (coeff.size() != fe.basis().size()) {
105 std::cout <<
"Bug in LocalInterpolation for finite element type "
106 << Dune::className<FE>() <<
":" << std::endl;
107 std::cout <<
" Interpolation vector has size " << coeff.size()
109 std::cout <<
" Basis has size " << fe.basis().size() << std::endl;
110 std::cout << std::endl;
118 for(std::size_t j=0; j<coeff.size() && success; ++j) {
119 if ( std::abs(coeff[j]-f.coeff[j]) >
120 (2*coeff.size()*eps)*(std::max(std::abs(f.coeff[j]), 1.0)) )
122 std::cout << std::setprecision(16);
123 std::cout <<
"Bug in LocalInterpolation for finite element type "
124 << Dune::className<FE>() <<
":" << std::endl;
125 std::cout <<
" Interpolation weight " << j <<
" differs by "
126 << std::abs(coeff[j]-f.coeff[j]) <<
" from coefficient of "
127 <<
"linear combination." << std::endl;
128 std::cout << std::endl;
148template<
class Geo,
class FE>
149bool testJacobian(
const Geo &geo,
const FE& fe,
double eps,
double delta,
150 std::size_t order = 2)
152 typedef typename FE::Traits::Basis Basis;
154 typedef typename Basis::Traits::DomainField DF;
155 static const std::size_t dimDLocal = Basis::Traits::dimDomainLocal;
156 typedef typename Basis::Traits::DomainLocal DomainLocal;
157 static const std::size_t dimDGlobal = Basis::Traits::dimDomainGlobal;
159 static const std::size_t dimR = Basis::Traits::dimRange;
160 typedef typename Basis::Traits::Range Range;
162 typedef typename Basis::Traits::Jacobian Jacobian;
176 for (std::size_t i=0; i < quad.size(); i++) {
179 const DomainLocal& testPoint = quad[i].position();
182 std::vector<Jacobian> jacobians;
183 fe.basis().evaluateJacobian(testPoint, jacobians);
184 if(jacobians.size() != fe.basis().size()) {
185 std::cout <<
"Bug in evaluateJacobianGlobal() for finite element type "
186 << Dune::className<FE>() <<
":" << std::endl;
187 std::cout <<
" Jacobian vector has size " << jacobians.size()
189 std::cout <<
" Basis has size " << fe.basis().size() << std::endl;
190 std::cout << std::endl;
195 geo.jacobianTransposed(testPoint);
198 for (std::size_t j=0; j<fe.basis().size(); ++j) {
204 for(std::size_t k = 0; k < dimR; ++k)
205 for(std::size_t l = 0; l < dimDGlobal; ++l)
206 for(std::size_t m = 0; m < dimDLocal; ++m)
207 localJacobian[k][m] += jacobians[j][k][l] * geoJT[m][l];
210 for (std::size_t m = 0; m < dimDLocal; ++m) {
213 DomainLocal upPos = testPoint;
214 DomainLocal downPos = testPoint;
219 std::vector<Range> upValues, downValues;
221 fe.basis().evaluateFunction(upPos, upValues);
222 fe.basis().evaluateFunction(downPos, downValues);
225 for(std::size_t k = 0; k < dimR; ++k) {
228 double derivative = localJacobian[k][m];
230 double finiteDiff = (upValues[j][k] - downValues[j][k]) / (2*delta);
233 if ( std::abs(derivative-finiteDiff) >
234 eps/delta*(std::max(std::abs(finiteDiff), 1.0)) )
236 std::cout << std::setprecision(16);
237 std::cout <<
"Bug in evaluateJacobian() for finite element type "
238 << Dune::className<FE>() <<
":" << std::endl;
239 std::cout <<
" Shape function derivative does not agree with "
240 <<
"FD approximation" << std::endl;
241 std::cout <<
" Shape function " << j <<
" component " << k
242 <<
" at position " << testPoint <<
": derivative in "
243 <<
"local direction " << m <<
" is "
244 << derivative <<
", but " << finiteDiff <<
" is "
245 <<
"expected." << std::endl;
246 std::cout << std::endl;
258template<
class Geo,
class FE>
259bool testFE(
const Geo &geo,
const FE& fe,
double eps,
double delta,
264 success = testInterpolation(fe, eps) and success;
265 success = testJacobian(geo, fe, eps, delta, order) and success;
A dense n x m matrix.
Definition: fmatrix.hh:117
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
A free function to provide the demangled class name of a given object or type as a string.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484