DUNE-FEM (unstable)

test-localfe.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
6#define DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
7
16#include <iomanip>
17#include <iostream>
18#include <typeinfo>
19
22#include <dune/geometry/referenceelements.hh>
23
24#include <dune/localfunctions/common/virtualinterface.hh>
25#include <dune/localfunctions/common/virtualwrappers.hh>
26
27double TOL = 1e-9;
28// The FD approximation used for checking the Jacobian uses half of the
29// precision -- so we have to be a little bit more tolerant here.
30double jacobianTOL = 1e-5; // sqrt(TOL)
31
32// This class wraps one shape function of a local finite element as a callable
33// that can be fed to the LocalInterpolation::interpolate method.
34template<class FE>
35class ShapeFunctionAsCallable
36{
37 // These types are deliberately private: They are not part of a Callable API
38 typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType;
39 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
40public:
41
42 typedef typename FE::Traits::LocalBasisType::Traits::RangeFieldType CT;
43
44 ShapeFunctionAsCallable(const FE& fe, int shapeFunction) :
45 fe_(fe),
46 shapeFunction_(shapeFunction)
47 {}
48
49 RangeType operator() (DomainType x) const
50 {
51 std::vector<RangeType> yy;
52 fe_.localBasis().evaluateFunction(x, yy);
53 return yy[shapeFunction_];
54 }
55
56private:
57 const FE& fe_;
58 int shapeFunction_;
59};
60
61
62// Check whether the degrees of freedom computed by LocalInterpolation
63// are dual to the shape functions. See Ciarlet, "The Finite Element Method
64// for Elliptic Problems", 1978, for details.
65template<class FE>
66bool testLocalInterpolation(const FE& fe)
67{
68 std::vector<typename ShapeFunctionAsCallable<FE>::CT> coeff;
69 for(size_t i=0; i<fe.size(); ++i)
70 {
72 // Feed the shape functions to the 'interpolate' method in form of a callable.
74
75 // The i-th shape function as a function that 'interpolate' can deal with
76 ShapeFunctionAsCallable<FE> sfAsCallable(fe, i);
77
78 // Compute degrees of freedom for that shape function
79 // We expect the result to be the i-th unit vector
80 fe.localInterpolation().interpolate(sfAsCallable, coeff);
81
82 // Check size of weight vector
83 if (coeff.size() != fe.localBasis().size())
84 {
85 std::cout << "Bug in LocalInterpolation for finite element type "
86 << Dune::className(fe) << std::endl;
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;
90 return false;
91 }
92
93 // Check if interpolation weights are equal to coefficients
94 for(std::size_t j=0; j<coeff.size(); ++j)
95 {
96 if ( std::abs(coeff[j] - (i==j)) > TOL)
97 {
98 std::cout << std::setprecision(16);
99 std::cout << "Bug in LocalInterpolation for finite element type "
100 << Dune::className(fe) << std::endl;
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;
104 return false;
105 }
106 }
107 }
108 return true;
109}
110
111
112// Check whether the space spanned by the shape functions
113// contains the constant functions
114template<class FE>
115bool testCanRepresentConstants(const FE& fe,
116 unsigned order = 5)
117{
118 typedef typename FE::Traits::LocalBasisType LB;
119 using RangeType = typename LB::Traits::RangeType;
120
121 bool success = true;
122
123 // Construct the constant '1' function
124 auto constantOne = [](const typename LB::Traits::DomainType& xi) { return RangeType(1.0); };
125
126 // Project the constant function onto the FE space
127 std::vector<double> coefficients;
128 fe.localInterpolation().interpolate(constantOne, coefficients);
129
130 // A set of test points
131 const auto& quad = Dune::QuadratureRules<double,LB::Traits::dimDomain>::rule(fe.type(),order);
132
133 // Loop over all quadrature points
134 for (size_t i=0; i<quad.size(); i++) {
135
136 // Get a test point
137 const auto& testPoint = quad[i].position();
138
139 // Compute value of the representation of constantOne at the test point
140 std::vector<RangeType> values;
141 fe.localBasis().evaluateFunction(testPoint, values);
142
143 RangeType sum(0);
144 for (size_t j=0; j<values.size(); j++)
145 sum += coefficients[j] * values[j];
146
147 if ((RangeType(1.0)-sum).two_norm() > TOL)
148 {
149 std::cout << "Finite element type " << Dune::className(fe)
150 << " cannot represent constant functions!" << std::endl;
151 std::cout << " At position: " << testPoint << ","
152 << std::endl;
153 std::cout << " discrete approximation of the '1' function has value " << sum
154 << std::endl;
155 std::cout << std::endl;
156 success = false;
157 }
158
159 } // Loop over all quadrature points
160
161 return success;
162}
163
164// check whether Jacobian agrees with FD approximation
165template<class FE>
166bool testJacobian(const FE& fe,
167 unsigned order = 2,
168 const std::function<bool(const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip = nullptr)
169{
170 typedef typename FE::Traits::LocalBasisType LB;
171
172 bool success = true;
173
174 // ////////////////////////////////////////////////////////////
175 // Check the partial derivatives by comparing them
176 // to finite difference approximations
177 // ////////////////////////////////////////////////////////////
178
179 // A set of test points
182
183 // Loop over all quadrature points
184 for (size_t i=0; i<quad.size(); i++) {
185
186 // Get a test point
188 quad[i].position();
189
190 // Get the shape function derivatives there
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 "
195 << Dune::className(fe) << std::endl;
196 std::cout << " Jacobian vector has size " << jacobians.size()
197 << std::endl;
198 std::cout << " Basis has size " << fe.localBasis().size()
199 << std::endl;
200 std::cout << std::endl;
201 return false;
202 }
203
204 // Skip this test point if we are supposed to
205 if (derivativePointSkip && derivativePointSkip(quad[i].position()))
206 continue;
207
208 // Loop over all directions
209 for (int k=0; k<LB::Traits::dimDomain; k++) {
210
211 // Compute an approximation to the derivative by finite differences
214
215 upPos[k] += jacobianTOL;
216 downPos[k] -= jacobianTOL;
217
218 std::vector<typename LB::Traits::RangeType> upValues, downValues;
219
220 fe.localBasis().evaluateFunction(upPos, upValues);
221 fe.localBasis().evaluateFunction(downPos, downValues);
222
223 // Loop over all shape functions in this set
224 for (unsigned int j=0; j<fe.localBasis().size(); ++j) {
225 //Loop over all components
226 for(int l=0; l < LB::Traits::dimRange; ++l) {
227
228 // The current partial derivative, just for ease of notation
229 double derivative = jacobians[j][l][k];
230
231 double finiteDiff = (upValues[j][l] - downValues[j][l])
232 / (2*jacobianTOL);
233
234 // Check
235 if ( std::abs(derivative-finiteDiff) >
236 TOL/jacobianTOL*((std::abs(finiteDiff)>1) ? std::abs(finiteDiff) : 1.) )
237 {
238 std::cout << std::setprecision(16);
239 std::cout << "Bug in evaluateJacobian() for finite element type "
240 << Dune::className(fe) << std::endl;
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;
248 success = false;
249 }
250 } //Loop over all components
251 } // Loop over all shape functions in this set
252 } // Loop over all directions
253 } // Loop over all quadrature points
254
255 return success;
256}
257
263{
264 template <class FE>
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)
268 {
269 bool success = true;
270
271 if (diffOrder > 2)
272 std::cout << "No test for differentiability orders larger than 2!" << std::endl;
273
274 if (diffOrder >= 2)
275 success = success and testOrder2(fe, eps, delta, order, derivativePointSkip);
276
277 if (diffOrder >= 1)
278 success = success and testOrder1(fe, eps, delta, order, derivativePointSkip);
279
280 success = success and testOrder0(fe, eps, delta, order);
281
282 return success;
283 }
284
286 template <class FE>
287 static bool testOrder0(const FE& fe,
288 double eps,
289 double delta,
290 std::size_t order = 2)
291 {
292 typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType;
293 constexpr auto dimDomain = FE::Traits::LocalBasisType::Traits::dimDomain;
294
295 bool success = true;
296
298 // Check the partial derivatives by comparing them
299 // to finite difference approximations
301
302 // A set of test points
303 const auto& quad = Dune::QuadratureRules<double, dimDomain>::rule(fe.type(),
304 order);
305
306 // Loop over all quadrature points
307 for (size_t i = 0; i < quad.size(); i++)
308 {
309 // Get a test point
310 const Dune::FieldVector<double, dimDomain>& testPoint = quad[i].position();
311
312 // Get the shape function values there using the 'partial' method
313 std::vector<RangeType> partialValues;
314 std::array<unsigned int, dimDomain> multiIndex;
315 std::fill(multiIndex.begin(), multiIndex.end(), 0);
316
317 fe.localBasis().partial(multiIndex, testPoint, partialValues);
318
319 if (partialValues.size() != fe.localBasis().size())
320 {
321 std::cout << "Bug in partial() for finite element type "
322 << Dune::className(fe) << std::endl;
323 std::cout << " values vector has size "
324 << partialValues.size() << std::endl;
325 std::cout << " Basis has size " << fe.localBasis().size()
326 << std::endl;
327 std::cout << std::endl;
328 return false;
329 }
330
331 // Get reference values
332 std::vector<RangeType> referenceValues;
333 fe.localBasis().evaluateFunction(testPoint, referenceValues);
334
335 // Loop over all shape functions in this set
336 for (unsigned int j = 0; j < fe.localBasis().size(); ++j)
337 {
338 // Loop over all components
339 for (int l = 0; l < FE::Traits::LocalBasisType::Traits::dimRange; ++l)
340 {
341 // Check the 'partial' method
342 if (std::abs(partialValues[j][l] - referenceValues[j][l])
343 > TOL / jacobianTOL
344 * ((std::abs(referenceValues[j][l]) > 1) ? std::abs(referenceValues[j][l]) : 1.))
345 {
346 std::cout << std::setprecision(16);
347 std::cout << "Bug in partial() for finite element type "
348 << Dune::className(fe) << std::endl;
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;
355 success = false;
356 }
357
358 } //Loop over all components
359 } // Loop over all shape functions in this set
360 } // Loop over all quadrature points
361
362 return success;
363 }
364
366 template <class FE>
367 static bool testOrder1(const FE& fe,
368 double eps,
369 double delta,
370 std::size_t order = 2,
371 const std::function<bool(const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip = nullptr)
372 {
373 typedef typename FE::Traits::LocalBasisType LB;
374 typedef typename LB::Traits::RangeFieldType RangeField;
375
376 bool success = true;
377
379 // Check the partial derivatives by comparing them
380 // to finite difference approximations
382
383 // A set of test points
386 order);
387
388 // Loop over all quadrature points
389 for (size_t i = 0; i < quad.size(); i++)
390 {
391 // Get a test point
392 const Dune::FieldVector<double, LB::Traits::dimDomain>& testPoint = quad[i].position();
393
394 // Skip the test points we are supposed to skip
395 if (derivativePointSkip && derivativePointSkip(testPoint))
396 continue;
397
398 // Loop over all directions
399 for (int k = 0; k < LB::Traits::dimDomain; k++)
400 {
401 // Get the shape function derivatives there using the 'partial' method
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);
405 multiIndex[k]++;
406 fe.localBasis().partial(multiIndex, testPoint, firstPartialDerivatives);
407 if (firstPartialDerivatives.size() != fe.localBasis().size())
408 {
409 std::cout << "Bug in partial() for finite element type "
410 << Dune::className(fe) << std::endl;
411 std::cout << " firstPartialDerivatives vector has size "
412 << firstPartialDerivatives.size() << std::endl;
413 std::cout << " Basis has size " << fe.localBasis().size()
414 << std::endl;
415 std::cout << std::endl;
416 return false;
417 }
418
419 // Compute an approximation to the derivative by finite differences
422
423 upPos[k] += jacobianTOL;
424 downPos[k] -= jacobianTOL;
425
426 std::vector<typename LB::Traits::RangeType> upValues, downValues;
427
428 fe.localBasis().evaluateFunction(upPos, upValues);
429 fe.localBasis().evaluateFunction(downPos, downValues);
430
431 // Loop over all shape functions in this set
432 for (unsigned int j = 0; j < fe.localBasis().size(); ++j)
433 {
434 // Loop over all components
435 for (int l = 0; l < LB::Traits::dimRange; ++l)
436 {
437 RangeField finiteDiff = (upValues[j][l] - downValues[j][l])
438 / (2 * jacobianTOL);
439
440 // Check the 'partial' method
441 RangeField partialDerivative = firstPartialDerivatives[j][l];
442 if (std::abs(partialDerivative - finiteDiff)
443 > TOL / jacobianTOL
444 * ((std::abs(finiteDiff) > 1) ? std::abs(finiteDiff) : 1.))
445 {
446 std::cout << std::setprecision(16);
447 std::cout << "Bug in partial() for finite element type "
448 << Dune::className(fe) << std::endl;
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;
456 success = false;
457 }
458
459 } // Loop over all directions
460 } // Loop over all shape functions in this set
461 } //Loop over all components
462 } // Loop over all quadrature points
463
464 return success;
465 }
466
468 template <class FE>
469 static bool testOrder2(const FE& fe,
470 double eps,
471 double delta,
472 std::size_t order = 2,
473 const std::function<bool(const typename FE::Traits::LocalBasisType::Traits::DomainType&)> derivativePointSkip = nullptr)
474 {
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;
479
480 static const std::size_t dimR = LocalBasis::Traits::dimRange;
481 typedef typename LocalBasis::Traits::RangeType Range;
482 typedef typename LocalBasis::Traits::RangeFieldType RangeField;
483
484 bool success = true;
485
487 // Check the partial derivatives by comparing them
488 // to finite difference approximations
490
491 // A set of test points
494
495 // Loop over all quadrature points
496 for (std::size_t i = 0; i < quad.size(); i++)
497 {
498 // Get a test point
499 const Domain& testPoint = quad[i].position();
500
501 // Skip the test points we are supposed to skip
502 if (derivativePointSkip && derivativePointSkip(testPoint))
503 continue;
504
505 // For testing the 'partial' method
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());
509
510 //loop over all local directions
511 for (int dir0 = 0; dir0 < dimDomain; dir0++)
512 {
513 for (int dir1 = 0; dir1 < dimDomain; dir1++)
514 {
515 // Get the shape function derivatives there using the 'partial' method
516 std::vector<Range> secondPartialDerivative;
517 std::array<unsigned int,dimDomain> multiIndex;
518 std::fill(multiIndex.begin(), multiIndex.end(), 0);
519 multiIndex[dir0]++;
520 multiIndex[dir1]++;
521 fe.localBasis().partial(multiIndex, testPoint, secondPartialDerivative);
522
523 if (secondPartialDerivative.size() != fe.localBasis().size())
524 {
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()
528 << std::endl;
529 std::cout << " Basis has size " << fe.localBasis().size()
530 << std::endl;
531 std::cout << std::endl;
532 return false;
533 }
534
535 //combine to Hesse matrices
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];
539
540 }
541 } //loop over all directions
542
543 // Loop over all local directions
544 for (std::size_t dir0 = 0; dir0 < dimDomain; ++dir0)
545 {
546 for (unsigned int dir1 = 0; dir1 < dimDomain; dir1++)
547 {
548 // Compute an approximation to the derivative by finite differences
549 std::array<Domain,4> neighbourPos;
550 std::fill(neighbourPos.begin(), neighbourPos.end(), testPoint);
551
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;
560
561 std::array<std::vector<Range>, 4> neighbourValues;
562 for (int k = 0; k < 4; k++)
563 fe.localBasis().evaluateFunction(neighbourPos[k],
564 neighbourValues[k]);
565
566 // Loop over all shape functions in this set
567 for (std::size_t j = 0; j < fe.localBasis().size(); ++j)
568 {
569 //Loop over all components
570 for (std::size_t k = 0; k < dimR; ++k)
571 {
572 RangeField finiteDiff = (neighbourValues[0][j][k]
573 - neighbourValues[1][j][k] - neighbourValues[2][j][k]
574 + neighbourValues[3][j][k]) / (4 * delta * delta);
575
576 // The current partial derivative, just for ease of notation, evaluated by the 'partial' method
577 RangeField partialDerivative = partialHessians[k][j][dir0][dir1];
578
579 // Check
580 if (std::abs(partialDerivative - finiteDiff)
581 > eps / delta * (std::max(std::abs(finiteDiff), 1.0)))
582 {
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;
594 success = false;
595 }
596
597 } //Loop over all components
598 }
599 } // Loop over all local directions
600 } // Loop over all shape functions in this set
601 } // Loop over all quadrature points
602
603 return success;
604 }
605
606};
607
608// Flags for disabling parts of testFE
609enum {
610 DisableNone = 0,
611 DisableLocalInterpolation = 1,
612 DisableVirtualInterface = 2,
613 DisableJacobian = 4,
614 DisableEvaluate = 8,
615 DisableRepresentConstants = 16
616};
617
628template<class FE>
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)
633{
634 // Order of the quadrature rule used to generate test points
635 unsigned int quadOrder = 2;
636
637 bool success = true;
638
639 if (FE::Traits::LocalBasisType::Traits::dimDomain != fe.type().dim())
640 {
641 std::cout << "Bug in type() for finite element type "
642 << Dune::className(fe) << std::endl;
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;
645 success = false;
646 }
647
648 if (fe.size() != fe.localBasis().size())
649 {
650 std::cout << "Bug in finite element type "
651 << Dune::className(fe) << std::endl;
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;
654 success = false;
655 }
656
657 // Make sure evaluateFunction returns the correct number of values
658 std::vector<typename FE::Traits::LocalBasisType::Traits::RangeType> values;
659 fe.localBasis().evaluateFunction(Dune::ReferenceElements<double,FE::Traits::LocalBasisType::Traits::dimDomain>::general(fe.type()).position(0,0), values);
660
661 if (values.size() != fe.size())
662 {
663 std::cout << "Bug in finite element type "
664 << Dune::className(fe) << std::endl;
665 std::cout << " LocalFiniteElement.size() returns " << fe.size() << "," << std::endl;
666 std::cout << " but LocalBasis::evaluateFunction returns " << values.size() << " values!" << std::endl;
667 success = false;
668 }
669
670 if (fe.size() != fe.localCoefficients().size())
671 {
672 std::cout << "Bug in finite element type "
673 << Dune::className(fe) << std::endl;
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;
676 success = false;
677 }
678
679 const auto& lc = fe.localCoefficients();
680 for(size_t i=0; i<lc.size(); i++)
681 {
682 const auto& lk = lc.localKey(i);
683 if (lk.codim() > fe.type().dim())
684 {
685 std::cout << "Bug in finite element type "
686 << Dune::className(fe) << std::endl;
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;
689 success = false;
690 }
691 }
692
693 if (not (disabledTests & DisableLocalInterpolation))
694 {
695 success = testLocalInterpolation<FE>(fe) and success;
696 }
697
698 if (not (disabledTests & DisableRepresentConstants))
699 {
700 success = testCanRepresentConstants<FE>(fe) and success;
701 }
702
703 if (not (disabledTests & DisableJacobian))
704 {
705 success = testJacobian<FE>(fe, quadOrder, derivativePointSkip) and success;
706 }
707 else
708 {
709 // make sure diffOrder is 0
710 success = (diffOrder == 0) and success;
711 }
712
713 if (not (disabledTests & DisableEvaluate))
714 {
715 success = TestPartial::test(fe, TOL, jacobianTOL, diffOrder, quadOrder, derivativePointSkip) and success;
716 }
717
718 if (not (disabledTests & DisableVirtualInterface))
719 {
720 typedef typename FE::Traits::LocalBasisType::Traits ImplementationLBTraits;
722 typedef typename Dune::LocalFiniteElementVirtualImp<FE> VirtualFEImp;
723
724 const VirtualFEImp virtualFE(fe);
725 if (not (disabledTests & DisableLocalInterpolation))
726 success = testLocalInterpolation<VirtualFEInterface>(virtualFE) and success;
727 if (not (disabledTests & DisableJacobian))
728 {
729 success = testJacobian<VirtualFEInterface>(virtualFE, quadOrder, derivativePointSkip) and success;
730 }
731 else
732 {
733 // make sure diffOrder is 0
734 success = (diffOrder == 0) and success;
735 }
736 }
737
738 return success;
739}
740
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; }
745
746#endif // DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)