3#ifndef DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
4#define DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
20#include <dune/geometry/referenceelements.hh>
22#include <dune/localfunctions/common/virtualinterface.hh>
23#include <dune/localfunctions/common/virtualwrappers.hh>
28double jacobianTOL = 1e-5;
33class ShapeFunctionAsFunction :
39 typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType;
40 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
43 typedef typename FE::Traits::LocalBasisType::Traits::RangeFieldType CT;
45 ShapeFunctionAsFunction(
const FE& fe,
int shapeFunction) :
47 shapeFunction_(shapeFunction)
50 void evaluate (
const DomainType& x, RangeType& y)
const
52 std::vector<RangeType> yy;
53 fe_.localBasis().evaluateFunction(x, yy);
54 y = yy[shapeFunction_];
65class ShapeFunctionAsCallable
68 typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType;
69 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
72 ShapeFunctionAsCallable(
const FE& fe,
int shapeFunction) :
74 shapeFunction_(shapeFunction)
77 RangeType operator() (DomainType x)
const
79 std::vector<RangeType> yy;
80 fe_.localBasis().evaluateFunction(x, yy);
81 return yy[shapeFunction_];
94bool testLocalInterpolation(
const FE& fe)
96 std::vector<typename ShapeFunctionAsFunction<FE>::CT> coeff;
97 for(
size_t i=0; i<fe.size(); ++i)
106 ShapeFunctionAsFunction<FE> f(fe, i);
110 fe.localInterpolation().interpolate(f, coeff);
113 if (coeff.size() != fe.localBasis().size())
115 std::cout <<
"Bug in LocalInterpolation for finite element type "
117 std::cout <<
" Interpolation produces " << coeff.size() <<
" degrees of freedom" << std::endl;
118 std::cout <<
" Basis has size " << fe.localBasis().size() << std::endl;
119 std::cout << std::endl;
124 for(std::size_t j=0; j<coeff.size(); ++j)
126 if ( std::abs(coeff[j] - (i==j)) > TOL)
128 std::cout << std::setprecision(16);
129 std::cout <<
"Bug in LocalInterpolation for finite element type "
131 std::cout <<
" Degree of freedom " << j <<
" applied to shape function " << i
132 <<
" yields value " << coeff[j] <<
", not the expected value " << (i==j) << std::endl;
133 std::cout << std::endl;
145 ShapeFunctionAsCallable<FE> sfAsCallable(fe, i);
149 fe.localInterpolation().interpolate(sfAsCallable, coeff);
152 if (coeff.size() != fe.localBasis().size())
154 std::cout <<
"Bug in LocalInterpolation for finite element type "
156 std::cout <<
" Interpolation produces " << coeff.size() <<
" degrees of freedom" << std::endl;
157 std::cout <<
" Basis has size " << fe.localBasis().size() << std::endl;
158 std::cout << std::endl;
163 for(std::size_t j=0; j<coeff.size(); ++j)
165 if ( std::abs(coeff[j] - (i==j)) > TOL)
167 std::cout << std::setprecision(16);
168 std::cout <<
"Bug in LocalInterpolation for finite element type "
170 std::cout <<
" Degree of freedom " << j <<
" applied to shape function " << i
171 <<
" yields value " << coeff[j] <<
", not the expected value " << (i==j) << std::endl;
172 std::cout << std::endl;
183bool testJacobian(
const FE& fe,
unsigned order = 2)
185 typedef typename FE::Traits::LocalBasisType LB;
199 for (
size_t i=0; i<quad.size(); i++) {
206 std::vector<typename LB::Traits::JacobianType> jacobians;
207 fe.localBasis().evaluateJacobian(testPoint, jacobians);
208 if(jacobians.size() != fe.localBasis().size()) {
209 std::cout <<
"Bug in evaluateJacobian() for finite element type "
211 std::cout <<
" Jacobian vector has size " << jacobians.size()
213 std::cout <<
" Basis has size " << fe.localBasis().size()
215 std::cout << std::endl;
220 for (
int k=0; k<LB::Traits::dimDomain; k++) {
226 upPos[k] += jacobianTOL;
227 downPos[k] -= jacobianTOL;
229 std::vector<typename LB::Traits::RangeType> upValues, downValues;
231 fe.localBasis().evaluateFunction(upPos, upValues);
232 fe.localBasis().evaluateFunction(downPos, downValues);
235 for (
unsigned int j=0; j<fe.localBasis().size(); ++j) {
237 for(
int l=0; l < LB::Traits::dimRange; ++l) {
240 double derivative = jacobians[j][l][k];
242 double finiteDiff = (upValues[j][l] - downValues[j][l])
246 if ( std::abs(derivative-finiteDiff) >
247 TOL/jacobianTOL*((std::abs(finiteDiff)>1) ? std::abs(finiteDiff) : 1.) )
249 std::cout << std::setprecision(16);
250 std::cout <<
"Bug in evaluateJacobian() for finite element type "
252 std::cout <<
" Shape function derivative does not agree with "
253 <<
"FD approximation" << std::endl;
254 std::cout <<
" Shape function " << j <<
" component " << l
255 <<
" at position " << testPoint <<
": derivative in "
256 <<
"direction " << k <<
" is " << derivative <<
", but "
257 << finiteDiff <<
" is expected." << std::endl;
258 std::cout << std::endl;
276 static bool test(
const FE& fe,
277 double eps,
double delta,
unsigned int diffOrder, std::size_t order = 2)
282 std::cout <<
"No test for differentiability orders larger than 2!" << std::endl;
285 success = success and
testOrder2(fe, eps, delta, order);
288 success = success and
testOrder1(fe, eps, delta, order);
290 success = success and
testOrder0(fe, eps, delta, order);
300 std::size_t order = 2)
302 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
303 constexpr auto dimDomain = FE::Traits::LocalBasisType::Traits::dimDomain;
317 for (
size_t i = 0; i < quad.size(); i++)
323 std::vector<RangeType> partialValues;
324 std::array<unsigned int, dimDomain> multiIndex;
325 std::fill(multiIndex.begin(), multiIndex.end(), 0);
327 fe.localBasis().partial(multiIndex, testPoint, partialValues);
329 if (partialValues.size() != fe.localBasis().size())
331 std::cout <<
"Bug in partial() for finite element type "
333 std::cout <<
" values vector has size "
334 << partialValues.size() << std::endl;
335 std::cout <<
" Basis has size " << fe.localBasis().size()
337 std::cout << std::endl;
342 std::vector<RangeType> referenceValues;
343 fe.localBasis().evaluateFunction(testPoint, referenceValues);
346 for (
unsigned int j = 0; j < fe.localBasis().size(); ++j)
349 for (
int l = 0; l < FE::Traits::LocalBasisType::Traits::dimRange; ++l)
352 if (std::abs(partialValues[j][l] - referenceValues[j][l])
354 * ((std::abs(referenceValues[j][l]) > 1) ? std::abs(referenceValues[j][l]) : 1.))
356 std::cout << std::setprecision(16);
357 std::cout <<
"Bug in partial() for finite element type "
359 std::cout <<
" Shape function value does not agree with "
360 <<
"output of method evaluateFunction." << std::endl;
361 std::cout <<
" Shape function " << j <<
" component " << l
362 <<
" at position " << testPoint <<
": value is " << partialValues[j][l]
363 <<
", but " << referenceValues[j][l] <<
" is expected." << std::endl;
364 std::cout << std::endl;
380 std::size_t order = 2)
382 typedef typename FE::Traits::LocalBasisType LB;
383 typedef typename LB::Traits::RangeFieldType RangeField;
398 for (
size_t i = 0; i < quad.size(); i++)
404 for (
int k = 0; k < LB::Traits::dimDomain; k++)
407 std::vector<typename LB::Traits::RangeType> firstPartialDerivatives;
408 std::array<unsigned int, LB::Traits::dimDomain> multiIndex;
409 std::fill(multiIndex.begin(), multiIndex.end(), 0);
411 fe.localBasis().partial(multiIndex, testPoint, firstPartialDerivatives);
412 if (firstPartialDerivatives.size() != fe.localBasis().size())
414 std::cout <<
"Bug in partial() for finite element type "
416 std::cout <<
" firstPartialDerivatives vector has size "
417 << firstPartialDerivatives.size() << std::endl;
418 std::cout <<
" Basis has size " << fe.localBasis().size()
420 std::cout << std::endl;
428 upPos[k] += jacobianTOL;
429 downPos[k] -= jacobianTOL;
431 std::vector<typename LB::Traits::RangeType> upValues, downValues;
433 fe.localBasis().evaluateFunction(upPos, upValues);
434 fe.localBasis().evaluateFunction(downPos, downValues);
437 for (
unsigned int j = 0; j < fe.localBasis().size(); ++j)
440 for (
int l = 0; l < LB::Traits::dimRange; ++l)
442 RangeField finiteDiff = (upValues[j][l] - downValues[j][l])
446 RangeField partialDerivative = firstPartialDerivatives[j][l];
447 if (std::abs(partialDerivative - finiteDiff)
449 * ((std::abs(finiteDiff) > 1) ? std::abs(finiteDiff) : 1.))
451 std::cout << std::setprecision(16);
452 std::cout <<
"Bug in partial() for finite element type "
454 std::cout <<
" Shape function derivative does not agree with "
455 <<
"FD approximation" << std::endl;
456 std::cout <<
" Shape function " << j <<
" component " << l
457 <<
" at position " << testPoint <<
": derivative in "
458 <<
"direction " << k <<
" is " << partialDerivative <<
", but "
459 << finiteDiff <<
" is expected." << std::endl;
460 std::cout << std::endl;
477 std::size_t order = 2)
479 typedef typename FE::Traits::LocalBasisType LocalBasis;
480 typedef typename LocalBasis::Traits::DomainFieldType DF;
481 typedef typename LocalBasis::Traits::DomainType Domain;
482 static const int dimDomain = LocalBasis::Traits::dimDomain;
484 static const std::size_t dimR = LocalBasis::Traits::dimRange;
485 typedef typename LocalBasis::Traits::RangeType Range;
486 typedef typename LocalBasis::Traits::RangeFieldType RangeField;
500 for (std::size_t i = 0; i < quad.size(); i++)
503 const Domain& testPoint = quad[i].position();
506 std::array<std::vector<Dune::FieldMatrix<RangeField, dimDomain, dimDomain> >, dimR> partialHessians;
507 for (
size_t k = 0; k < dimR; k++)
508 partialHessians[k].resize(fe.size());
511 for (
int dir0 = 0; dir0 < dimDomain; dir0++)
513 for (
int dir1 = 0; dir1 < dimDomain; dir1++)
516 std::vector<Range> secondPartialDerivative;
517 std::array<unsigned int,dimDomain> multiIndex;
518 std::fill(multiIndex.begin(), multiIndex.end(), 0);
521 fe.localBasis().partial(multiIndex, testPoint, secondPartialDerivative);
523 if (secondPartialDerivative.size() != fe.localBasis().size())
525 std::cout <<
"Bug in partial() for finite element type "
526 << Dune::className<FE>() <<
":" << std::endl;
527 std::cout <<
" return vector has size " << secondPartialDerivative.size()
529 std::cout <<
" Basis has size " << fe.localBasis().size()
531 std::cout << std::endl;
536 for (
size_t k = 0; k < dimR; k++)
537 for (std::size_t j = 0; j < fe.localBasis().size(); ++j)
538 partialHessians[k][j][dir0][dir1] = secondPartialDerivative[j][k];
544 for (std::size_t dir0 = 0; dir0 < dimDomain; ++dir0)
546 for (
unsigned int dir1 = 0; dir1 < dimDomain; dir1++)
549 std::array<Domain,4> neighbourPos;
550 std::fill(neighbourPos.begin(), neighbourPos.end(), testPoint);
552 neighbourPos[0][dir0] += delta;
553 neighbourPos[0][dir1] += delta;
554 neighbourPos[1][dir0] -= delta;
555 neighbourPos[1][dir1] += delta;
556 neighbourPos[2][dir0] += delta;
557 neighbourPos[2][dir1] -= delta;
558 neighbourPos[3][dir0] -= delta;
559 neighbourPos[3][dir1] -= delta;
561 std::array<std::vector<Range>, 4> neighbourValues;
562 for (
int k = 0; k < 4; k++)
563 fe.localBasis().evaluateFunction(neighbourPos[k],
567 for (std::size_t j = 0; j < fe.localBasis().size(); ++j)
570 for (std::size_t k = 0; k < dimR; ++k)
572 RangeField finiteDiff = (neighbourValues[0][j][k]
573 - neighbourValues[1][j][k] - neighbourValues[2][j][k]
574 + neighbourValues[3][j][k]) / (4 * delta * delta);
577 RangeField partialDerivative = partialHessians[k][j][dir0][dir1];
580 if (std::abs(partialDerivative - finiteDiff)
581 > eps / delta * (
std::max(std::abs(finiteDiff), 1.0)))
583 std::cout << std::setprecision(16);
584 std::cout <<
"Bug in partial() for finite element type "
585 << Dune::className<FE>() <<
":" << std::endl;
586 std::cout <<
" Second shape function derivative does not agree with "
587 <<
"FD approximation" << std::endl;
588 std::cout <<
" Shape function " << j <<
" component " << k
589 <<
" at position " << testPoint <<
": derivative in "
590 <<
"local direction (" << dir0 <<
", " << dir1 <<
") is "
591 << partialDerivative <<
", but " << finiteDiff
592 <<
" is expected." << std::endl;
593 std::cout << std::endl;
611 DisableLocalInterpolation = 1,
612 DisableVirtualInterface = 2,
619bool testFE(
const FE& fe,
char disabledTests = DisableNone,
unsigned int diffOrder = 0)
622 unsigned int quadOrder = 2;
626 if (FE::Traits::LocalBasisType::Traits::dimDomain != fe.type().dim())
628 std::cout <<
"Bug in type() for finite element type "
630 std::cout <<
" Coordinate dimension is " << FE::Traits::LocalBasisType::Traits::dimDomain << std::endl;
631 std::cout <<
" but GeometryType is " << fe.type() <<
" with dimension " << fe.type().dim() << std::endl;
635 if (fe.size() != fe.localBasis().size())
637 std::cout <<
"Bug in finite element type "
639 std::cout <<
" Size reported by LocalFiniteElement is " << fe.size() << std::endl;
640 std::cout <<
" but size reported by LocalBasis is " << fe.localBasis().size() << std::endl;
645 std::vector<typename FE::Traits::LocalBasisType::Traits::RangeType> values;
648 if (values.size() != fe.size())
650 std::cout <<
"Bug in finite element type "
652 std::cout <<
" LocalFiniteElement.size() returns " << fe.size() <<
"," << std::endl;
653 std::cout <<
" but LocalBasis::evaluateFunction returns " << values.size() <<
" values!" << std::endl;
657 if (fe.size() != fe.localCoefficients().size())
659 std::cout <<
"Bug in finite element type "
661 std::cout <<
" Size reported by LocalFiniteElement is " << fe.size() << std::endl;
662 std::cout <<
" but size reported by LocalCoefficients is " << fe.localCoefficients().size() << std::endl;
666 const auto& lc = fe.localCoefficients();
667 for(
size_t i=0; i<lc.size(); i++)
669 const auto& lk = lc.localKey(i);
670 if (lk.codim() > fe.type().dim())
672 std::cout <<
"Bug in finite element type "
674 std::cout <<
" Codimension reported by localKey(" << i <<
") is " << lk.codim() << std::endl;
675 std::cout <<
" but geometry is " << fe.type() <<
" with dimension " << fe.type().dim() << std::endl;
680 if (not (disabledTests & DisableLocalInterpolation))
682 success = testLocalInterpolation<FE>(fe) and success;
684 if (not (disabledTests & DisableJacobian))
686 success = testJacobian<FE>(fe, quadOrder) and success;
691 success = (diffOrder == 0) and success;
694 if (not (disabledTests & DisableEvaluate))
696 success = TestPartial::test(fe, TOL, jacobianTOL, diffOrder, quadOrder) and success;
699 if (not (disabledTests & DisableVirtualInterface))
701 typedef typename FE::Traits::LocalBasisType::Traits ImplementationLBTraits;
705 const VirtualFEImp virtualFE(fe);
706 if (not (disabledTests & DisableLocalInterpolation))
707 success = testLocalInterpolation<VirtualFEInterface>(virtualFE) and success;
708 if (not (disabledTests & DisableJacobian))
710 success = testJacobian<VirtualFEInterface>(virtualFE) and success;
715 success = (diffOrder == 0) and success;
722#define TEST_FE(A) { bool b = testFE(A); std::cout << "testFE(" #A ") " << (b?"succeeded\n":"failed\n"); success &= b; }
723#define TEST_FE2(A,B) { bool b = testFE(A, B); std::cout << "testFE(" #A ", " #B ") " << (b?"succeeded\n":"failed\n"); success &= b; }
724#define TEST_FE3(A,B,C) { bool b = testFE(A, B, C); std::cout << "testFE(" #A ", " #B ", " #C ") " << (b?"succeeded\n":"failed\n"); success &= b; }
vector space out of a tensor product of fields.
Definition: fvector.hh:96
void evaluate(const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Function evaluation.
Return a proper base class for functions to use with LocalInterpolation.
Definition: virtualinterface.hh:37
class for wrapping a finite element using the virtual interface
Definition: virtualwrappers.hh:240
virtual base class for local finite elements with functions
Definition: virtualinterface.hh:260
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.
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
std::string className()
Provide the demangled class name of a type T as a string.
Definition: classname.hh:44
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:168
Helper class to test the 'partial' method.
Definition: test-localfe.hh:274
static bool testOrder0(const FE &fe, double eps, double delta, std::size_t order=2)
Test the 'partial' method for zero-order partial derivatives, i.e., values.
Definition: test-localfe.hh:297
static bool testOrder2(const FE &fe, double eps, double delta, std::size_t order=2)
Test second-order partial derivatives.
Definition: test-localfe.hh:474
static bool testOrder1(const FE &fe, double eps, double delta, std::size_t order=2)
Test the 'partial' method for first-order partial derivatives.
Definition: test-localfe.hh:377