3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_TEST_BASISTEST_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_TEST_BASISTEST_HH
13#include <dune/common/test/testsuite.hh>
16#include <dune/common/hybridutilities.hh>
20#include <dune/functions/functionspacebases/concepts.hh>
22struct CheckBasisFlag {};
23struct AllowZeroBasisFunctions {};
25template<
class T,
class... S>
26struct IsContained :
public std::disjunction<std::is_same<T,S>...>
34template<
class Element,
class Gr
idView>
35std::string elementStr(
const Element& element,
const GridView& gridView)
38 s << element.type() <<
"#" << gridView.indexSet().index(element);
46template<
class MultiIndex>
47bool multiIndicesConsecutive(
const MultiIndex& a,
const MultiIndex& b)
52 for (; (i<a.size()) and (i<b.size()) and (a[i] == b[i]); ++i)
56 if ((i<a.size()) and (i==b.size()))
60 if ((i<a.size()) and (i<b.size()))
70 for (; i<b.size(); ++i)
86template<
class MultiIndexSet>
87Dune::TestSuite checkBasisIndexTreeConsistency(
const MultiIndexSet& multiIndexSet)
93 auto it = multiIndexSet.begin();
94 auto end = multiIndexSet.end();
97 auto lastMultiIndex = *it;
100 test.require(lastMultiIndex.size()>0,
"multi-index size check")
101 <<
"empty multi-index found";
104 for (
decltype(lastMultiIndex.size()) i = 0; i<lastMultiIndex.size(); ++i)
106 test.require(lastMultiIndex[i] == 0,
"smallest index check")
107 <<
"smallest index contains non-zero entry " << lastMultiIndex[i] <<
" in position " << i;
111 for(; it != end; ++it)
113 auto multiIndex = *it;
116 test.require(multiIndex.size()>0,
"multi-index size check")
117 <<
"empty multi-index found";
120 test.check(multiIndicesConsecutive(lastMultiIndex, multiIndex),
"consecutive index check")
121 <<
"multi-indices " << lastMultiIndex <<
" and " << multiIndex <<
" are subsequent but not consecutive";
123 lastMultiIndex = multiIndex;
134template<
class Basis,
class MultiIndexSet>
135Dune::TestSuite checkBasisSizeConsistency(
const Basis& basis,
const MultiIndexSet& multiIndexSet)
141 using Prefix =
typename Basis::SizePrefix;
142 auto prefixSet = std::map<Prefix, std::size_t>();
143 for(
const auto& index : multiIndexSet)
145 auto prefix = Prefix();
146 for (
const auto& i: index)
148 prefixSet[prefix] =
std::max(prefixSet[prefix], i+1);
151 prefixSet[prefix] = 0;
156 for(
const auto& [prefix, size] : prefixSet)
158 auto prefixSize = basis.size(prefix);
159 test.check(prefixSize == size,
"basis.size(prefix) check")
160 <<
"basis.size(" << prefix <<
")=" << prefixSize <<
", but should be " << size;
179 using MultiIndex =
typename Basis::MultiIndex;
183 auto compare = [](
const auto& a,
const auto& b) {
184 return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end());
187 auto multiIndexSet = std::set<MultiIndex,
decltype(compare)>{compare};
189 auto localView = basis.localView();
190 for (
const auto& e :
elements(basis.gridView()))
194 for (
decltype(localView.size()) i=0; i< localView.size(); ++i)
196 auto multiIndex = localView.index(i);
197 for(
auto mi: multiIndex)
199 <<
"Global multi-index contains negative entry for shape function " << i
200 <<
" in element " << elementStr(localView.element(), basis.gridView());
201 multiIndexSet.insert(multiIndex);
205 test.subTest(checkBasisIndexTreeConsistency(multiIndexSet));
206 test.subTest(checkBasisSizeConsistency(basis, multiIndexSet));
207 test.check(basis.dimension() == multiIndexSet.size())
208 <<
"basis.dimension() does not equal the total number of basis functions.";
219template<
class LocalFiniteElement>
220Dune::TestSuite checkNonZeroShapeFunctions(
const LocalFiniteElement& fe, std::size_t order = 5,
double tol = 1e-10)
223 static const int dimension = LocalFiniteElement::Traits::LocalBasisType::Traits::dimDomain;
227 std::vector<typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeType> values;
228 std::vector<bool> isNonZero;
229 isNonZero.resize(fe.size(),
false);
230 for (
const auto& qp : quadRule)
232 fe.localBasis().evaluateFunction(qp.position(), values);
233 for(std::size_t i=0; i<fe.size(); ++i)
234 isNonZero[i] = (isNonZero[i] or (values[i].infinity_norm() > tol));
236 for(std::size_t i=0; i<fe.size(); ++i)
237 test.check(isNonZero[i])
238 <<
"Found a constant zero basis function";
248template<
class Basis,
class LocalView,
class... Flags>
249Dune::TestSuite checkLocalView(
const Basis& basis,
const LocalView& localView, Flags... flags)
251 Dune::TestSuite test(std::string(
"LocalView on ") + elementStr(localView.element(), basis.gridView()));
253 test.check(localView.size() <= localView.maxSize(),
"localView.size() check")
254 <<
"localView.size() is " << localView.size() <<
" but localView.maxSize() is " << localView.maxSize();
257 std::vector<std::size_t> localIndices;
258 localIndices.resize(localView.size(), 0);
260 test.check(node.size() == node.finiteElement().size())
261 <<
"Size of leaf node and finite element are different.";
262 for(std::size_t i=0; i<node.size(); ++i)
264 test.check(node.localIndex(i) < localView.size())
265 <<
"Local index exceeds localView.size().";
266 if (node.localIndex(i) < localView.size())
267 ++(localIndices[node.localIndex(i)]);
272 for(std::size_t i=0; i<localView.size(); ++i)
275 test.check(localIndices[i]>=1)
276 <<
"Local index " << i <<
" did not appear";
277 test.check(localIndices[i]<=1)
278 <<
"Local index " << i <<
" appears multiple times";
282 if (not IsContained<AllowZeroBasisFunctions, Flags...>::value)
285 test.subTest(checkNonZeroShapeFunctions(node.finiteElement()));
299struct EnableContinuityCheck
301 std::size_t order_ = 5;
304 template<
class JumpEvaluator>
305 auto localJumpContinuityCheck(
const JumpEvaluator& jumpEvaluator, std::size_t order,
double tol)
const
307 return [=](
const auto& intersection,
const auto&
treePath,
const auto& insideNode,
const auto& outsideNode,
const auto& insideToOutside) {
308 using Intersection = std::decay_t<
decltype(intersection)>;
309 using Node = std::decay_t<
decltype(insideNode)>;
311 std::vector<int> isContinuous(insideNode.size(),
true);
314 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
315 std::vector<std::vector<Range>> values;
316 std::vector<std::vector<Range>> neighborValues;
319 values.resize(quadRule.size());
320 neighborValues.resize(quadRule.size());
321 for(std::size_t k=0; k<quadRule.size(); ++k)
323 auto pointInElement = intersection.geometryInInside().global(quadRule[k].position());
324 auto pointInNeighbor = intersection.geometryInOutside().global(quadRule[k].position());
325 insideNode.finiteElement().localBasis().evaluateFunction(pointInElement, values[k]);
326 outsideNode.finiteElement().localBasis().evaluateFunction(pointInNeighbor, neighborValues[k]);
330 for(std::size_t i=0; i<insideNode.size(); ++i)
332 for(std::size_t k=0; k<quadRule.size(); ++k)
334 auto jump = values[k][i];
335 if (insideToOutside[i].has_value())
336 jump -= neighborValues[k][insideToOutside[i].value()];
337 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection, quadRule[k].position()) < tol);
344 auto localContinuityCheck()
const {
345 auto jumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
346 return jump.infinity_norm();
348 return localJumpContinuityCheck(jumpNorm, order_, tol_);
358struct EnableNormalContinuityCheck :
public EnableContinuityCheck
360 auto localContinuityCheck()
const {
361 auto normalJump = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
362 return jump * intersection.unitOuterNormal(x);
364 return localJumpContinuityCheck(normalJump, order_, tol_);
375struct EnableTangentialContinuityCheck :
public EnableContinuityCheck
377 auto localContinuityCheck()
const {
378 auto tangentialJumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
379 auto tangentialJump = jump - (jump * intersection.unitOuterNormal(x)) * intersection.unitOuterNormal(x);
380 return tangentialJump.two_norm();
382 return localJumpContinuityCheck(tangentialJumpNorm, order_, tol_);
393struct EnableCenterContinuityCheck :
public EnableContinuityCheck
395 template<
class JumpEvaluator>
396 auto localJumpCenterContinuityCheck(
const JumpEvaluator& jumpEvaluator,
double tol)
const
398 return [=](
const auto& intersection,
const auto&
treePath,
const auto& insideNode,
const auto& outsideNode,
const auto& insideToOutside) {
399 using Node = std::decay_t<
decltype(insideNode)>;
400 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
402 std::vector<int> isContinuous(insideNode.size(),
true);
403 std::vector<Range> insideValues;
404 std::vector<Range> outsideValues;
406 insideNode.finiteElement().localBasis().evaluateFunction(intersection.geometryInInside().center(), insideValues);
407 outsideNode.finiteElement().localBasis().evaluateFunction(intersection.geometryInOutside().center(), outsideValues);
409 auto centerLocal = intersection.geometry().local(intersection.geometry().center());
412 for(std::size_t i=0; i<insideNode.size(); ++i)
414 auto jump = insideValues[i];
415 if (insideToOutside[i].has_value())
416 jump -= outsideValues[insideToOutside[i].value()];
417 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection, centerLocal) < tol);
423 auto localContinuityCheck()
const {
424 auto jumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
425 return jump.infinity_norm();
427 return localJumpCenterContinuityCheck(jumpNorm, tol_);
443template<
class Basis,
class LocalCheck>
444Dune::TestSuite checkBasisContinuity(
const Basis& basis,
const LocalCheck& localCheck)
449 auto localView = basis.localView();
450 auto neighborLocalView = basis.localView();
452 for (
const auto& e :
elements(basis.gridView()))
455 for(
const auto& intersection :
intersections(basis.gridView(), e))
457 if (intersection.neighbor())
459 neighborLocalView.bind(intersection.outside());
462 const auto& outsideNode = Dune::TypeTree::child(neighborLocalView.tree(), treePath);
464 std::vector<std::optional<int>> insideToOutside;
465 insideToOutside.resize(insideNode.size());
468 for(std::size_t i=0; i<insideNode.size(); ++i)
470 for(std::size_t j=0; j<outsideNode.size(); ++j)
472 if (localView.index(insideNode.localIndex(i)) == neighborLocalView.index(outsideNode.localIndex(j)))
475 test.check(not insideToOutside[i].has_value())
476 <<
"Basis function " << localView.index(insideNode.localIndex(i))
477 <<
" appears twice in element " << elementStr(neighborLocalView.element(), basis.gridView());
478 insideToOutside[i] = j;
484 auto isContinuous = localCheck(intersection,
treePath, insideNode, outsideNode, insideToOutside);
486 for(std::size_t i=0; i<insideNode.size(); ++i)
488 test.check(isContinuous[i])
489 <<
"Basis function " << localView.index(insideNode.localIndex(i))
490 <<
" is discontinuous across intersection of elements "
491 << elementStr(localView.element(), basis.gridView())
492 <<
" and " << elementStr(neighborLocalView.element(), basis.gridView());
501template<
class Basis,
class... Flags>
506 using GridView =
typename Basis::GridView;
509 test.check(
Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, Basis>(),
"global basis concept check")
510 <<
"type passed to checkBasis() does not model the GlobalBasis concept";
513 auto localView = basis.localView();
514 for (
const auto& e :
elements(basis.gridView()))
517 test.subTest(checkLocalView(basis, localView, flags...));
521 test.subTest(checkBasisIndices(basis));
525 auto flagTuple = std::tie(flags...);
527 using Flag = std::decay_t<
decltype(flag)>;
528 if constexpr (std::is_base_of_v<EnableContinuityCheck, Flag>)
529 test.subTest(checkBasisContinuity(basis, flag.localContinuityCheck()));
535template<
class Basis,
class... Flags>
541 test.subTest(checkConstBasis(basis,flags...));
546 test.subTest(checkConstBasis(copy,flags...));
548 test.subTest(checkConstBasis(copy,flags...));
552 auto gridView = basis.gridView();
553 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 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