7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_TEST_BASISTEST_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_TEST_BASISTEST_HH
17#include <dune/common/test/testsuite.hh>
20#include <dune/common/hybridutilities.hh>
24#include <dune/functions/functionspacebases/concepts.hh>
25#include <dune/functions/functionspacebases/test/enabledifferentiabilitycheck.hh>
27struct CheckBasisFlag {};
28struct AllowZeroBasisFunctions {};
30template<
class T,
class... S>
31struct IsContained :
public std::disjunction<std::is_same<T,S>...>
37template<
class Element,
class Gr
idView>
38std::string elementStr(
const Element& element,
const GridView& gridView)
41 s << element.type() <<
"#" << gridView.indexSet().index(element);
49template<
class MultiIndex>
50bool multiIndicesConsecutive(
const MultiIndex& a,
const MultiIndex& b)
55 for (; (i<a.size()) and (i<b.size()) and (a[i] == b[i]); ++i)
59 if ((i<a.size()) and (i==b.size()))
63 if ((i<a.size()) and (i<b.size()))
73 for (; i<b.size(); ++i)
89template<
class MultiIndexSet>
90Dune::TestSuite checkBasisIndexTreeConsistency(
const MultiIndexSet& multiIndexSet)
96 auto it = multiIndexSet.begin();
97 auto end = multiIndexSet.end();
100 auto lastMultiIndex = *it;
103 test.require(lastMultiIndex.size()>0,
"multi-index size check")
104 <<
"empty multi-index found";
107 for (
decltype(lastMultiIndex.size()) i = 0; i<lastMultiIndex.size(); ++i)
109 test.require(lastMultiIndex[i] == 0,
"smallest index check")
110 <<
"smallest index contains non-zero entry " << lastMultiIndex[i] <<
" in position " << i;
114 for(; it != end; ++it)
116 auto multiIndex = *it;
119 test.require(multiIndex.size()>0,
"multi-index size check")
120 <<
"empty multi-index found";
123 test.check(multiIndicesConsecutive(lastMultiIndex, multiIndex),
"consecutive index check")
124 <<
"multi-indices " << lastMultiIndex <<
" and " << multiIndex <<
" are subsequent but not consecutive";
126 lastMultiIndex = multiIndex;
137template<
class Basis,
class MultiIndexSet>
138Dune::TestSuite checkBasisSizeConsistency(
const Basis& basis,
const MultiIndexSet& multiIndexSet)
144 using Prefix =
typename Basis::SizePrefix;
145 auto prefixSet = std::map<Prefix, std::size_t>();
146 for(
const auto& index : multiIndexSet)
148 auto prefix = Prefix();
149 for (
const auto& i: index)
151 prefixSet[prefix] =
std::max(prefixSet[prefix], i+1);
154 prefixSet[prefix] = 0;
159 for(
const auto& [prefix, size] : prefixSet)
161 auto prefixSize = basis.size(prefix);
162 test.check(prefixSize == size,
"basis.size(prefix) check")
163 <<
"basis.size(" << prefix <<
")=" << prefixSize <<
", but should be " << size;
182 using MultiIndex =
typename Basis::MultiIndex;
186 auto compare = [](
const auto& a,
const auto& b) {
187 return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end());
190 auto multiIndexSet = std::set<MultiIndex,
decltype(compare)>{compare};
192 auto localView = basis.localView();
193 for (
const auto& e :
elements(basis.gridView()))
197 for (
decltype(localView.size()) i=0; i< localView.size(); ++i)
199 auto multiIndex = localView.index(i);
200 for(
auto mi: multiIndex)
202 <<
"Global multi-index contains negative entry for shape function " << i
203 <<
" in element " << elementStr(localView.element(), basis.gridView());
204 multiIndexSet.insert(multiIndex);
208 test.subTest(checkBasisIndexTreeConsistency(multiIndexSet));
209 test.subTest(checkBasisSizeConsistency(basis, multiIndexSet));
210 test.check(basis.dimension() == multiIndexSet.size())
211 <<
"basis.dimension() does not equal the total number of basis functions.";
222template<
class LocalFiniteElement>
223Dune::TestSuite checkNonZeroShapeFunctions(
const LocalFiniteElement& fe, std::size_t order = 5,
double tol = 1e-10)
226 static const int dimension = LocalFiniteElement::Traits::LocalBasisType::Traits::dimDomain;
230 std::vector<typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeType> values;
231 std::vector<bool> isNonZero;
232 isNonZero.resize(fe.size(),
false);
233 for (
const auto& qp : quadRule)
235 fe.localBasis().evaluateFunction(qp.position(), values);
236 for(std::size_t i=0; i<fe.size(); ++i)
237 isNonZero[i] = (isNonZero[i] or (values[i].infinity_norm() > tol));
239 for(std::size_t i=0; i<fe.size(); ++i)
240 test.check(isNonZero[i])
241 <<
"Found a constant zero basis function";
251template<
class Basis,
class LocalView,
class... Flags>
252Dune::TestSuite checkLocalView(
const Basis& basis,
const LocalView& localView, Flags... flags)
254 Dune::TestSuite test(std::string(
"LocalView on ") + elementStr(localView.element(), basis.gridView()));
256 test.check(localView.size() <= localView.maxSize(),
"localView.size() check")
257 <<
"localView.size() is " << localView.size() <<
" but localView.maxSize() is " << localView.maxSize();
260 std::vector<std::size_t> localIndices;
261 localIndices.resize(localView.size(), 0);
263 test.check(node.size() == node.finiteElement().size())
264 <<
"Size of leaf node and finite element are different.";
265 for(std::size_t i=0; i<node.size(); ++i)
267 test.check(node.localIndex(i) < localView.size())
268 <<
"Local index exceeds localView.size().";
269 if (node.localIndex(i) < localView.size())
270 ++(localIndices[node.localIndex(i)]);
275 for(std::size_t i=0; i<localView.size(); ++i)
278 test.check(localIndices[i]>=1)
279 <<
"Local index " << i <<
" did not appear";
280 test.check(localIndices[i]<=1)
281 <<
"Local index " << i <<
" appears multiple times";
285 if (not IsContained<AllowZeroBasisFunctions, Flags...>::value)
288 test.subTest(checkNonZeroShapeFunctions(node.finiteElement()));
302struct EnableContinuityCheck
304 std::size_t order_ = 5;
307 template<
class JumpEvaluator>
308 auto localJumpContinuityCheck(
const JumpEvaluator& jumpEvaluator, std::size_t order,
double tol)
const
310 return [=](
const auto& intersection,
const auto&
treePath,
const auto& insideNode,
const auto& outsideNode,
const auto& insideToOutside) {
311 using Intersection = std::decay_t<
decltype(intersection)>;
312 using Node = std::decay_t<
decltype(insideNode)>;
314 std::vector<int> isContinuous(insideNode.size(),
true);
317 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
318 std::vector<std::vector<Range>> values;
319 std::vector<std::vector<Range>> neighborValues;
322 values.resize(quadRule.size());
323 neighborValues.resize(quadRule.size());
324 for(std::size_t k=0; k<quadRule.size(); ++k)
326 auto pointInElement = intersection.geometryInInside().global(quadRule[k].position());
327 auto pointInNeighbor = intersection.geometryInOutside().global(quadRule[k].position());
328 insideNode.finiteElement().localBasis().evaluateFunction(pointInElement, values[k]);
329 outsideNode.finiteElement().localBasis().evaluateFunction(pointInNeighbor, neighborValues[k]);
333 for(std::size_t i=0; i<insideNode.size(); ++i)
335 for(std::size_t k=0; k<quadRule.size(); ++k)
337 auto jump = values[k][i];
338 if (insideToOutside[i].has_value())
339 jump -= neighborValues[k][insideToOutside[i].value()];
340 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection, quadRule[k].position()) < tol);
347 auto localContinuityCheck()
const {
348 auto jumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
349 return jump.infinity_norm();
351 return localJumpContinuityCheck(jumpNorm, order_, tol_);
361struct EnableNormalContinuityCheck :
public EnableContinuityCheck
363 auto localContinuityCheck()
const {
364 auto normalJump = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
365 return jump * intersection.unitOuterNormal(x);
367 return localJumpContinuityCheck(normalJump, order_, tol_);
378struct EnableTangentialContinuityCheck :
public EnableContinuityCheck
380 auto localContinuityCheck()
const {
381 auto tangentialJumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
382 auto tangentialJump = jump - (jump * intersection.unitOuterNormal(x)) * intersection.unitOuterNormal(x);
383 return tangentialJump.two_norm();
385 return localJumpContinuityCheck(tangentialJumpNorm, order_, tol_);
396struct EnableCenterContinuityCheck :
public EnableContinuityCheck
398 template<
class JumpEvaluator>
399 auto localJumpCenterContinuityCheck(
const JumpEvaluator& jumpEvaluator,
double tol)
const
401 return [=](
const auto& intersection,
const auto&
treePath,
const auto& insideNode,
const auto& outsideNode,
const auto& insideToOutside) {
402 using Node = std::decay_t<
decltype(insideNode)>;
403 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
405 std::vector<int> isContinuous(insideNode.size(),
true);
406 std::vector<Range> insideValues;
407 std::vector<Range> outsideValues;
409 insideNode.finiteElement().localBasis().evaluateFunction(intersection.geometryInInside().center(), insideValues);
410 outsideNode.finiteElement().localBasis().evaluateFunction(intersection.geometryInOutside().center(), outsideValues);
412 auto centerLocal = intersection.geometry().local(intersection.geometry().center());
415 for(std::size_t i=0; i<insideNode.size(); ++i)
417 auto jump = insideValues[i];
418 if (insideToOutside[i].has_value())
419 jump -= outsideValues[insideToOutside[i].value()];
420 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection, centerLocal) < tol);
426 auto localContinuityCheck()
const {
427 auto jumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
428 return jump.infinity_norm();
430 return localJumpCenterContinuityCheck(jumpNorm, tol_);
442 struct EnableVertexContinuityCheck:
public EnableContinuityCheck
444 template <
class JumpEvaluator>
445 auto localJumpVertexContinuityCheck(
const JumpEvaluator &jumpEvaluator,
double tol)
const
447 return [=](
const auto &intersection,
const auto &
treePath,
const auto &insideNode,
448 const auto &outsideNode,
const auto &insideToOutside)
450 using Node = std::decay_t<
decltype(insideNode)>;
451 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
453 std::vector<int> isContinuous(insideNode.size(),
true);
454 std::vector<Range> insideValues;
455 std::vector<Range> outsideValues;
457 std::vector<
typename std::decay_t<
decltype(intersection)>::LocalCoordinate> vertices;
459 for (
int i = 0; i < intersection.geometry().corners(); ++i)
460 vertices.push_back(intersection.geometry().local(intersection.geometry().corner(i)));
462 for (
auto const &
vertex : vertices)
465 insideNode.finiteElement().localBasis().evaluateFunction(
466 intersection.geometryInInside().global(
vertex), insideValues);
467 outsideNode.finiteElement().localBasis().evaluateFunction(
468 intersection.geometryInOutside().global(
vertex), outsideValues);
470 for (std::size_t i = 0; i < insideNode.size(); ++i)
472 auto jump = insideValues[i];
473 if (insideToOutside[i].has_value())
474 jump -= outsideValues[insideToOutside[i].value()];
475 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection,
vertex) < tol);
482 auto localContinuityCheck()
const
484 auto jumpNorm = [](
auto &&jump,
auto &&intersection,
auto &&x) ->
double
485 {
return jump.infinity_norm(); };
486 return localJumpVertexContinuityCheck(jumpNorm, tol_);
501template<
class Basis,
class LocalCheck>
502Dune::TestSuite checkBasisContinuity(
const Basis& basis,
const LocalCheck& localCheck)
507 auto localView = basis.localView();
508 auto neighborLocalView = basis.localView();
510 for (
const auto& e :
elements(basis.gridView()))
513 for(
const auto& intersection :
intersections(basis.gridView(), e))
515 if (intersection.neighbor())
517 neighborLocalView.bind(intersection.outside());
520 const auto& outsideNode = Dune::TypeTree::child(neighborLocalView.tree(), treePath);
522 std::vector<std::optional<int>> insideToOutside;
523 insideToOutside.resize(insideNode.size());
526 for(std::size_t i=0; i<insideNode.size(); ++i)
528 for(std::size_t j=0; j<outsideNode.size(); ++j)
530 if (localView.index(insideNode.localIndex(i)) == neighborLocalView.index(outsideNode.localIndex(j)))
533 test.check(not insideToOutside[i].has_value())
534 <<
"Basis function " << localView.index(insideNode.localIndex(i))
535 <<
" appears twice in element " << elementStr(neighborLocalView.element(), basis.gridView());
536 insideToOutside[i] = j;
542 auto isContinuous = localCheck(intersection,
treePath, insideNode, outsideNode, insideToOutside);
544 for(std::size_t i=0; i<insideNode.size(); ++i)
546 test.check(isContinuous[i])
547 <<
"Basis function " << localView.index(insideNode.localIndex(i))
548 <<
" is discontinuous across intersection of elements "
549 << elementStr(localView.element(), basis.gridView())
550 <<
" and " << elementStr(neighborLocalView.element(), basis.gridView());
559template<
class Basis,
class... Flags>
564 using GridView =
typename Basis::GridView;
567 test.check(
Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, Basis>(),
"global basis concept check")
568 <<
"type passed to checkBasis() does not model the GlobalBasis concept";
571 auto localView = basis.localView();
572 for (
const auto& e :
elements(basis.gridView()))
575 test.subTest(checkLocalView(basis, localView, flags...));
579 test.subTest(checkBasisIndices(basis));
583 auto flagTuple = std::tie(flags...);
585 using Flag = std::decay_t<
decltype(flag)>;
586 if constexpr (std::is_base_of_v<EnableContinuityCheck, Flag>)
587 test.subTest(checkBasisContinuity(basis, flag.localContinuityCheck()));
588 else if constexpr (std::is_base_of_v<EnableDifferentiabilityCheck, Flag>)
589 test.subTest(checkBasisDifferentiability(basis, flag));
595template<
class Basis,
class... Flags>
601 test.subTest(checkConstBasis(basis,flags...));
606 test.subTest(checkConstBasis(copy,flags...));
608 test.subTest(checkConstBasis(copy,flags...));
612 auto gridView = basis.gridView();
613 basis.update(gridView);
Grid view abstract base class.
Definition: gridview.hh:66
Intersection of a mesh entity of codimension 0 ("element") with a "neighboring" element or with the d...
Definition: intersection.hh:164
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 Simple helper class to organize your test suite.
Definition: testsuite.hh:31
Infrastructure for concepts.
Traits for type conversions and type information.
constexpr auto models()
Check if concept is modeled by given types.
Definition: concept.hh:184
IteratorRange<... > intersections(const GV &gv, const Entity &e)
Iterates over all Intersections of an Entity with respect to the given GridView.
IteratorRange<... > elements(const GV &gv)
Iterates over all elements / cells (entities with codimension 0) of a GridView.
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:326
void forEachLeafNode(Tree &&tree, LeafFunc &&leafFunc)
Traverse tree and visit each leaf node.
Definition: traversal.hh:269
Dune namespace.
Definition: alignedallocator.hh:13
Type trait to determine whether an instance of T has an operator[](I), i.e. whether it can be indexed...
Definition: typetraits.hh:250