5#ifndef DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
6#define DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
22#include <dune/geometry/referenceelements.hh>
24#include <dune/localfunctions/common/virtualinterface.hh>
25#include <dune/localfunctions/common/virtualwrappers.hh>
30double jacobianTOL = 1e-5;
35class ShapeFunctionAsCallable
38 typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType;
39 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
42 typedef typename FE::Traits::LocalBasisType::Traits::RangeFieldType CT;
44 ShapeFunctionAsCallable(
const FE& fe,
int shapeFunction) :
46 shapeFunction_(shapeFunction)
49 RangeType operator() (DomainType x)
const
51 std::vector<RangeType> yy;
52 fe_.localBasis().evaluateFunction(x, yy);
53 return yy[shapeFunction_];
66bool testLocalInterpolation(
const FE& fe)
68 std::vector<typename ShapeFunctionAsCallable<FE>::CT> coeff;
69 for(
size_t i=0; i<fe.size(); ++i)
76 ShapeFunctionAsCallable<FE> sfAsCallable(fe, i);
80 fe.localInterpolation().interpolate(sfAsCallable, coeff);
83 if (coeff.size() != fe.localBasis().size())
85 std::cout <<
"Bug in LocalInterpolation for finite element type "
87 std::cout <<
" Interpolation produces " << coeff.size() <<
" degrees of freedom" << std::endl;
88 std::cout <<
" Basis has size " << fe.localBasis().size() << std::endl;
89 std::cout << std::endl;
94 for(std::size_t j=0; j<coeff.size(); ++j)
96 if ( std::abs(coeff[j] - (i==j)) > TOL)
98 std::cout << std::setprecision(16);
99 std::cout <<
"Bug in LocalInterpolation for finite element type "
101 std::cout <<
" Degree of freedom " << j <<
" applied to shape function " << i
102 <<
" yields value " << coeff[j] <<
", not the expected value " << (i==j) << std::endl;
103 std::cout << std::endl;
115bool testCanRepresentConstants(
const FE& fe,
118 typedef typename FE::Traits::LocalBasisType LB;
119 using RangeType =
typename LB::Traits::RangeType;
124 auto constantOne = [](
const typename LB::Traits::DomainType& xi) {
return RangeType(1.0); };
127 std::vector<double> coefficients;
128 fe.localInterpolation().interpolate(constantOne, coefficients);
134 for (
size_t i=0; i<quad.size(); i++) {
137 const auto& testPoint = quad[i].position();
140 std::vector<RangeType> values;
141 fe.localBasis().evaluateFunction(testPoint, values);
144 for (
size_t j=0; j<values.size(); j++)
145 sum += coefficients[j] * values[j];
147 if ((RangeType(1.0)-sum).two_norm() > TOL)
150 <<
" cannot represent constant functions!" << std::endl;
151 std::cout <<
" At position: " << testPoint <<
","
153 std::cout <<
" discrete approximation of the '1' function has value " << sum
155 std::cout << std::endl;
166bool testJacobian(
const FE& fe,
168 const std::function<
bool(
const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip =
nullptr)
170 typedef typename FE::Traits::LocalBasisType LB;
184 for (
size_t i=0; i<quad.size(); i++) {
191 std::vector<typename LB::Traits::JacobianType> jacobians;
192 fe.localBasis().evaluateJacobian(testPoint, jacobians);
193 if(jacobians.size() != fe.localBasis().size()) {
194 std::cout <<
"Bug in evaluateJacobian() for finite element type "
196 std::cout <<
" Jacobian vector has size " << jacobians.size()
198 std::cout <<
" Basis has size " << fe.localBasis().size()
200 std::cout << std::endl;
205 if (derivativePointSkip && derivativePointSkip(quad[i].position()))
209 for (
int k=0; k<LB::Traits::dimDomain; k++) {
215 upPos[k] += jacobianTOL;
216 downPos[k] -= jacobianTOL;
218 std::vector<typename LB::Traits::RangeType> upValues, downValues;
220 fe.localBasis().evaluateFunction(upPos, upValues);
221 fe.localBasis().evaluateFunction(downPos, downValues);
224 for (
unsigned int j=0; j<fe.localBasis().size(); ++j) {
226 for(
int l=0; l < LB::Traits::dimRange; ++l) {
229 double derivative = jacobians[j][l][k];
231 double finiteDiff = (upValues[j][l] - downValues[j][l])
235 if ( std::abs(derivative-finiteDiff) >
236 TOL/jacobianTOL*((std::abs(finiteDiff)>1) ? std::abs(finiteDiff) : 1.) )
238 std::cout << std::setprecision(16);
239 std::cout <<
"Bug in evaluateJacobian() for finite element type "
241 std::cout <<
" Shape function derivative does not agree with "
242 <<
"FD approximation" << std::endl;
243 std::cout <<
" Shape function " << j <<
" component " << l
244 <<
" at position " << testPoint <<
": derivative in "
245 <<
"direction " << k <<
" is " << derivative <<
", but "
246 << finiteDiff <<
" is expected." << std::endl;
247 std::cout << std::endl;
265 static bool test(
const FE& fe,
266 double eps,
double delta,
unsigned int diffOrder, std::size_t order = 2,
267 const std::function<
bool(
const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip =
nullptr)
272 std::cout <<
"No test for differentiability orders larger than 2!" << std::endl;
275 success = success and
testOrder2(fe, eps, delta, order, derivativePointSkip);
278 success = success and
testOrder1(fe, eps, delta, order, derivativePointSkip);
280 success = success and
testOrder0(fe, eps, delta, order);
290 std::size_t order = 2)
292 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
293 constexpr auto dimDomain = FE::Traits::LocalBasisType::Traits::dimDomain;
307 for (
size_t i = 0; i < quad.size(); i++)
313 std::vector<RangeType> partialValues;
314 std::array<unsigned int, dimDomain> multiIndex;
315 std::fill(multiIndex.begin(), multiIndex.end(), 0);
317 fe.localBasis().partial(multiIndex, testPoint, partialValues);
319 if (partialValues.size() != fe.localBasis().size())
321 std::cout <<
"Bug in partial() for finite element type "
323 std::cout <<
" values vector has size "
324 << partialValues.size() << std::endl;
325 std::cout <<
" Basis has size " << fe.localBasis().size()
327 std::cout << std::endl;
332 std::vector<RangeType> referenceValues;
333 fe.localBasis().evaluateFunction(testPoint, referenceValues);
336 for (
unsigned int j = 0; j < fe.localBasis().size(); ++j)
339 for (
int l = 0; l < FE::Traits::LocalBasisType::Traits::dimRange; ++l)
342 if (std::abs(partialValues[j][l] - referenceValues[j][l])
344 * ((std::abs(referenceValues[j][l]) > 1) ? std::abs(referenceValues[j][l]) : 1.))
346 std::cout << std::setprecision(16);
347 std::cout <<
"Bug in partial() for finite element type "
349 std::cout <<
" Shape function value does not agree with "
350 <<
"output of method evaluateFunction." << std::endl;
351 std::cout <<
" Shape function " << j <<
" component " << l
352 <<
" at position " << testPoint <<
": value is " << partialValues[j][l]
353 <<
", but " << referenceValues[j][l] <<
" is expected." << std::endl;
354 std::cout << std::endl;
370 std::size_t order = 2,
371 const std::function<
bool(
const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip =
nullptr)
373 typedef typename FE::Traits::LocalBasisType LB;
374 typedef typename LB::Traits::RangeFieldType RangeField;
389 for (
size_t i = 0; i < quad.size(); i++)
395 if (derivativePointSkip && derivativePointSkip(testPoint))
399 for (
int k = 0; k < LB::Traits::dimDomain; k++)
402 std::vector<typename LB::Traits::RangeType> firstPartialDerivatives;
403 std::array<unsigned int, LB::Traits::dimDomain> multiIndex;
404 std::fill(multiIndex.begin(), multiIndex.end(), 0);
406 fe.localBasis().partial(multiIndex, testPoint, firstPartialDerivatives);
407 if (firstPartialDerivatives.size() != fe.localBasis().size())
409 std::cout <<
"Bug in partial() for finite element type "
411 std::cout <<
" firstPartialDerivatives vector has size "
412 << firstPartialDerivatives.size() << std::endl;
413 std::cout <<
" Basis has size " << fe.localBasis().size()
415 std::cout << std::endl;
423 upPos[k] += jacobianTOL;
424 downPos[k] -= jacobianTOL;
426 std::vector<typename LB::Traits::RangeType> upValues, downValues;
428 fe.localBasis().evaluateFunction(upPos, upValues);
429 fe.localBasis().evaluateFunction(downPos, downValues);
432 for (
unsigned int j = 0; j < fe.localBasis().size(); ++j)
435 for (
int l = 0; l < LB::Traits::dimRange; ++l)
437 RangeField finiteDiff = (upValues[j][l] - downValues[j][l])
441 RangeField partialDerivative = firstPartialDerivatives[j][l];
442 if (std::abs(partialDerivative - finiteDiff)
444 * ((std::abs(finiteDiff) > 1) ? std::abs(finiteDiff) : 1.))
446 std::cout << std::setprecision(16);
447 std::cout <<
"Bug in partial() for finite element type "
449 std::cout <<
" Shape function derivative does not agree with "
450 <<
"FD approximation" << std::endl;
451 std::cout <<
" Shape function " << j <<
" component " << l
452 <<
" at position " << testPoint <<
": derivative in "
453 <<
"direction " << k <<
" is " << partialDerivative <<
", but "
454 << finiteDiff <<
" is expected." << std::endl;
455 std::cout << std::endl;
472 std::size_t order = 2,
473 const std::function<
bool(
const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip =
nullptr)
475 typedef typename FE::Traits::LocalBasisType LocalBasis;
476 typedef typename LocalBasis::Traits::DomainFieldType DF;
477 typedef typename LocalBasis::Traits::DomainType Domain;
478 static const int dimDomain = LocalBasis::Traits::dimDomain;
480 static const std::size_t dimR = LocalBasis::Traits::dimRange;
481 typedef typename LocalBasis::Traits::RangeType Range;
482 typedef typename LocalBasis::Traits::RangeFieldType RangeField;
496 for (std::size_t i = 0; i < quad.size(); i++)
499 const Domain& testPoint = quad[i].position();
502 if (derivativePointSkip && derivativePointSkip(testPoint))
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,
615 DisableRepresentConstants = 16
629bool testFE(
const FE& fe,
630 char disabledTests = DisableNone,
631 unsigned int diffOrder = 0,
632 const std::function<
bool(
const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip =
nullptr)
635 unsigned int quadOrder = 2;
639 if (FE::Traits::LocalBasisType::Traits::dimDomain != fe.type().dim())
641 std::cout <<
"Bug in type() for finite element type "
643 std::cout <<
" Coordinate dimension is " << FE::Traits::LocalBasisType::Traits::dimDomain << std::endl;
644 std::cout <<
" but GeometryType is " << fe.type() <<
" with dimension " << fe.type().dim() << std::endl;
648 if (fe.size() != fe.localBasis().size())
650 std::cout <<
"Bug in finite element type "
652 std::cout <<
" Size reported by LocalFiniteElement is " << fe.size() << std::endl;
653 std::cout <<
" but size reported by LocalBasis is " << fe.localBasis().size() << std::endl;
658 std::vector<typename FE::Traits::LocalBasisType::Traits::RangeType> values;
661 if (values.size() != fe.size())
663 std::cout <<
"Bug in finite element type "
665 std::cout <<
" LocalFiniteElement.size() returns " << fe.size() <<
"," << std::endl;
666 std::cout <<
" but LocalBasis::evaluateFunction returns " << values.size() <<
" values!" << std::endl;
670 if (fe.size() != fe.localCoefficients().size())
672 std::cout <<
"Bug in finite element type "
674 std::cout <<
" Size reported by LocalFiniteElement is " << fe.size() << std::endl;
675 std::cout <<
" but size reported by LocalCoefficients is " << fe.localCoefficients().size() << std::endl;
679 const auto& lc = fe.localCoefficients();
680 for(
size_t i=0; i<lc.size(); i++)
682 const auto& lk = lc.localKey(i);
683 if (lk.codim() > fe.type().dim())
685 std::cout <<
"Bug in finite element type "
687 std::cout <<
" Codimension reported by localKey(" << i <<
") is " << lk.codim() << std::endl;
688 std::cout <<
" but geometry is " << fe.type() <<
" with dimension " << fe.type().dim() << std::endl;
693 if (not (disabledTests & DisableLocalInterpolation))
695 success = testLocalInterpolation<FE>(fe) and success;
698 if (not (disabledTests & DisableRepresentConstants))
700 success = testCanRepresentConstants<FE>(fe) and success;
703 if (not (disabledTests & DisableJacobian))
705 success = testJacobian<FE>(fe, quadOrder, derivativePointSkip) and success;
710 success = (diffOrder == 0) and success;
713 if (not (disabledTests & DisableEvaluate))
715 success = TestPartial::test(fe, TOL, jacobianTOL, diffOrder, quadOrder, derivativePointSkip) and success;
718 if (not (disabledTests & DisableVirtualInterface))
720 typedef typename FE::Traits::LocalBasisType::Traits ImplementationLBTraits;
724 const VirtualFEImp virtualFE(fe);
725 if (not (disabledTests & DisableLocalInterpolation))
726 success = testLocalInterpolation<VirtualFEInterface>(virtualFE) and success;
727 if (not (disabledTests & DisableJacobian))
729 success = testJacobian<VirtualFEInterface>(virtualFE, quadOrder, derivativePointSkip) and success;
734 success = (diffOrder == 0) and success;
741#define TEST_FE(A) { bool b = testFE(A); std::cout << "testFE(" #A ") " << (b?"succeeded\n":"failed\n"); success &= b; }
742#define TEST_FE2(A,B) { bool b = testFE(A, B); std::cout << "testFE(" #A ", " #B ") " << (b?"succeeded\n":"failed\n"); success &= b; }
743#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; }
744#define TEST_FE4(A,B,C,D) { bool b = testFE(A, B, C, D); std::cout << "testFE(" #A ", " #B ", " #C ", " #D ") " << (b?"succeeded\n":"failed\n"); success &= b; }
vector space out of a tensor product of fields.
Definition: fvector.hh:91
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:225
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.
std::string className()
Provide the demangled class name of a type T as a string.
Definition: classname.hh:47
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:128
Helper class to test the 'partial' method.
Definition: test-localfe.hh:263
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:287
static bool testOrder1(const FE &fe, double eps, double delta, std::size_t order=2, const std::function< bool(const typename FE::Traits::LocalBasisType::Traits::DomainType &)> derivativePointSkip=nullptr)
Test the 'partial' method for first-order partial derivatives.
Definition: test-localfe.hh:367
static bool testOrder2(const FE &fe, double eps, double delta, std::size_t order=2, const std::function< bool(const typename FE::Traits::LocalBasisType::Traits::DomainType &)> derivativePointSkip=nullptr)
Test second-order partial derivatives.
Definition: test-localfe.hh:469