9#ifndef DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
10#define DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
34 public Dune::Function<const typename FE::Traits::Basis::Traits::DomainLocal&, typename FE::Traits::Basis::Traits::Range>
36 typedef typename FE::Traits::Basis::Traits::DomainLocal DomainLocal;
37 typedef typename FE::Traits::Basis::Traits::Range Range;
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 void evaluate (
const DomainLocal& x, Range& y)
const {
61 std::vector<Range> yy;
62 fe.basis().evaluateFunction(x, yy);
65 for (std::size_t i=0; i<yy.size(); ++i)
66 y.axpy(coeff[i], yy[i]);
84bool testInterpolation(
const FE& fe,
double eps,
int n=5)
89 std::vector<typename FEFunction<FE>::CT> coeff;
90 for(
int i=0; i<n && success; ++i) {
95 fe.interpolation().interpolate(f, coeff);
98 if (coeff.size() != fe.basis().size()) {
99 std::cout <<
"Bug in LocalInterpolation for finite element type "
100 << Dune::className<FE>() <<
":" << std::endl;
101 std::cout <<
" Interpolation vector has size " << coeff.size()
103 std::cout <<
" Basis has size " << fe.basis().size() << std::endl;
104 std::cout << std::endl;
112 for(std::size_t j=0; j<coeff.size() && success; ++j) {
113 if ( std::abs(coeff[j]-f.coeff[j]) >
114 eps*(
std::max(std::abs(f.coeff[j]), 1.0)) )
116 std::cout << std::setprecision(16);
117 std::cout <<
"Bug in LocalInterpolation for finite element type "
118 << Dune::className<FE>() <<
":" << std::endl;
119 std::cout <<
" Interpolation weight " << j <<
" differs by "
120 << std::abs(coeff[j]-f.coeff[j]) <<
" from coefficient of "
121 <<
"linear combination." << std::endl;
122 std::cout << std::endl;
142template<
class Geo,
class FE>
143bool testJacobian(
const Geo &geo,
const FE& fe,
double eps,
double delta,
144 std::size_t order = 2)
146 typedef typename FE::Traits::Basis Basis;
148 typedef typename Basis::Traits::DomainField DF;
149 static const std::size_t dimDLocal = Basis::Traits::dimDomainLocal;
150 typedef typename Basis::Traits::DomainLocal DomainLocal;
151 static const std::size_t dimDGlobal = Basis::Traits::dimDomainGlobal;
153 static const std::size_t dimR = Basis::Traits::dimRange;
154 typedef typename Basis::Traits::Range Range;
156 typedef typename Basis::Traits::Jacobian Jacobian;
170 for (std::size_t i=0; i < quad.size(); i++) {
173 const DomainLocal& testPoint = quad[i].position();
176 std::vector<Jacobian> jacobians;
177 fe.basis().evaluateJacobian(testPoint, jacobians);
178 if(jacobians.size() != fe.basis().size()) {
179 std::cout <<
"Bug in evaluateJacobianGlobal() for finite element type "
180 << Dune::className<FE>() <<
":" << std::endl;
181 std::cout <<
" Jacobian vector has size " << jacobians.size()
183 std::cout <<
" Basis has size " << fe.basis().size() << std::endl;
184 std::cout << std::endl;
189 geo.jacobianTransposed(testPoint);
192 for (std::size_t j=0; j<fe.basis().size(); ++j) {
198 for(std::size_t k = 0; k < dimR; ++k)
199 for(std::size_t l = 0; l < dimDGlobal; ++l)
200 for(std::size_t m = 0; m < dimDLocal; ++m)
201 localJacobian[k][m] += jacobians[j][k][l] * geoJT[m][l];
204 for (std::size_t m = 0; m < dimDLocal; ++m) {
207 DomainLocal upPos = testPoint;
208 DomainLocal downPos = testPoint;
213 std::vector<Range> upValues, downValues;
215 fe.basis().evaluateFunction(upPos, upValues);
216 fe.basis().evaluateFunction(downPos, downValues);
219 for(std::size_t k = 0; k < dimR; ++k) {
222 double derivative = localJacobian[k][m];
224 double finiteDiff = (upValues[j][k] - downValues[j][k]) / (2*delta);
227 if ( std::abs(derivative-finiteDiff) >
228 eps/delta*(
std::max(std::abs(finiteDiff), 1.0)) )
230 std::cout << std::setprecision(16);
231 std::cout <<
"Bug in evaluateJacobian() for finite element type "
232 << Dune::className<FE>() <<
":" << std::endl;
233 std::cout <<
" Shape function derivative does not agree with "
234 <<
"FD approximation" << std::endl;
235 std::cout <<
" Shape function " << j <<
" component " << k
236 <<
" at position " << testPoint <<
": derivative in "
237 <<
"local direction " << m <<
" is "
238 << derivative <<
", but " << finiteDiff <<
" is "
239 <<
"expected." << std::endl;
240 std::cout << std::endl;
252template<
class Geo,
class FE>
253bool testFE(
const Geo &geo,
const FE& fe,
double eps,
double delta,
258 success = testInterpolation(fe, eps) and success;
259 success = testJacobian(geo, fe, eps, delta, order) and success;
A dense n x m matrix.
Definition: fmatrix.hh:69
Base class template for function classes.
Definition: function.hh:29
void evaluate(const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Function evaluation.
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:126
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:254
A free function to provide the demangled class name of a given object or type as a string.
Simple base class templates for functions.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79