7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_TEST_BASISTEST_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_TEST_BASISTEST_HH
17#include <dune/common/test/testsuite.hh>
18#include <dune/common/concept.hh>
19#include <dune/common/typetraits.hh>
20#include <dune/common/hybridutilities.hh>
22#include <dune/geometry/quadraturerules.hh>
24#include <dune/functions/functionspacebases/concepts.hh>
26struct CheckBasisFlag {};
27struct AllowZeroBasisFunctions {};
29template<
class T,
class... S>
30struct IsContained :
public std::disjunction<std::is_same<T,S>...>
38template<
class Element,
class Gr
idView>
39std::string elementStr(
const Element& element,
const GridView& gridView)
42 s << element.type() <<
"#" << gridView.indexSet().index(element);
50template<
class MultiIndex>
51bool multiIndicesConsecutive(
const MultiIndex& a,
const MultiIndex& b)
56 for (; (i<a.size()) and (i<b.size()) and (a[i] == b[i]); ++i)
60 if ((i<a.size()) and (i==b.size()))
64 if ((i<a.size()) and (i<b.size()))
74 for (; i<b.size(); ++i)
90template<
class MultiIndexSet>
91Dune::TestSuite checkBasisIndexTreeConsistency(
const MultiIndexSet& multiIndexSet)
93 Dune::TestSuite test(
"index tree consistency check");
97 auto it = multiIndexSet.begin();
98 auto end = multiIndexSet.end();
101 auto lastMultiIndex = *it;
104 test.require(lastMultiIndex.size()>0,
"multi-index size check")
105 <<
"empty multi-index found";
108 for (
decltype(lastMultiIndex.size()) i = 0; i<lastMultiIndex.size(); ++i)
110 test.require(lastMultiIndex[i] == 0,
"smallest index check")
111 <<
"smallest index contains non-zero entry " << lastMultiIndex[i] <<
" in position " << i;
115 for(; it != end; ++it)
117 auto multiIndex = *it;
120 test.require(multiIndex.size()>0,
"multi-index size check")
121 <<
"empty multi-index found";
124 test.check(multiIndicesConsecutive(lastMultiIndex, multiIndex),
"consecutive index check")
125 <<
"multi-indices " << lastMultiIndex <<
" and " << multiIndex <<
" are subsequent but not consecutive";
127 lastMultiIndex = multiIndex;
138template<
class Basis,
class MultiIndexSet>
139Dune::TestSuite checkBasisSizeConsistency(
const Basis& basis,
const MultiIndexSet& multiIndexSet)
141 Dune::TestSuite test(
"index size consistency check");
145 using Prefix =
typename Basis::SizePrefix;
146 auto prefixSet = std::map<Prefix, std::size_t>();
147 for(
const auto& index : multiIndexSet)
149 auto prefix = Prefix();
150 for (
const auto& i: index)
152 prefixSet[prefix] = std::max(prefixSet[prefix], i+1);
155 prefixSet[prefix] = 0;
160 for(
const auto& [prefix, size] : prefixSet)
162 auto prefixSize = basis.size(prefix);
163 test.check(prefixSize == size,
"basis.size(prefix) check")
164 <<
"basis.size(" << prefix <<
")=" << prefixSize <<
", but should be " << size;
179Dune::TestSuite checkBasisIndices(
const Basis& basis)
181 Dune::TestSuite test(
"basis index check");
183 using MultiIndex =
typename Basis::MultiIndex;
185 static_assert(Dune::IsIndexable<MultiIndex>(),
"MultiIndex must support operator[]");
187 auto compare = [](
const auto& a,
const auto& b) {
188 return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end());
191 auto multiIndexSet = std::set<MultiIndex,
decltype(compare)>{compare};
193 auto localView = basis.localView();
194 for (
const auto& e : elements(basis.gridView()))
198 for (
decltype(localView.size()) i=0; i< localView.size(); ++i)
200 auto multiIndex = localView.index(i);
201 for(
auto mi: multiIndex)
203 <<
"Global multi-index contains negative entry for shape function " << i
204 <<
" in element " << elementStr(localView.element(), basis.gridView());
205 multiIndexSet.insert(multiIndex);
209 test.subTest(checkBasisIndexTreeConsistency(multiIndexSet));
210 test.subTest(checkBasisSizeConsistency(basis, multiIndexSet));
211 test.check(basis.dimension() == multiIndexSet.size())
212 <<
"basis.dimension() does not equal the total number of basis functions.";
223template<
class LocalFiniteElement>
224Dune::TestSuite checkNonZeroShapeFunctions(
const LocalFiniteElement& fe, std::size_t order = 5,
double tol = 1e-10)
226 Dune::TestSuite test;
227 static const int dimension = LocalFiniteElement::Traits::LocalBasisType::Traits::dimDomain;
229 auto quadRule = Dune::QuadratureRules<double, dimension>::rule(fe.type(), order);
231 std::vector<typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeType> values;
232 std::vector<bool> isNonZero;
233 isNonZero.resize(fe.size(),
false);
234 for (
const auto& qp : quadRule)
236 fe.localBasis().evaluateFunction(qp.position(), values);
237 for(std::size_t i=0; i<fe.size(); ++i)
238 isNonZero[i] = (isNonZero[i] or (values[i].infinity_norm() > tol));
240 for(std::size_t i=0; i<fe.size(); ++i)
241 test.check(isNonZero[i])
242 <<
"Found a constant zero basis function";
252template<
class Basis,
class LocalView,
class... Flags>
253Dune::TestSuite checkLocalView(
const Basis& basis,
const LocalView& localView, Flags... flags)
255 Dune::TestSuite test(std::string(
"LocalView on ") + elementStr(localView.element(), basis.gridView()));
257 test.check(localView.size() <= localView.maxSize(),
"localView.size() check")
258 <<
"localView.size() is " << localView.size() <<
" but localView.maxSize() is " << localView.maxSize();
261 std::vector<std::size_t> localIndices;
262 localIndices.resize(localView.size(), 0);
263 Dune::TypeTree::forEachLeafNode(localView.tree(), [&](
const auto& node,
auto&& treePath) {
264 test.check(node.size() == node.finiteElement().size())
265 <<
"Size of leaf node and finite element are different.";
266 for(std::size_t i=0; i<node.size(); ++i)
268 test.check(node.localIndex(i) < localView.size())
269 <<
"Local index exceeds localView.size().";
270 if (node.localIndex(i) < localView.size())
271 ++(localIndices[node.localIndex(i)]);
276 for(std::size_t i=0; i<localView.size(); ++i)
279 test.check(localIndices[i]>=1)
280 <<
"Local index " << i <<
" did not appear";
281 test.check(localIndices[i]<=1)
282 <<
"Local index " << i <<
" appears multiple times";
286 if (not IsContained<AllowZeroBasisFunctions, Flags...>::value)
288 Dune::TypeTree::forEachLeafNode(localView.tree(), [&](
const auto& node,
auto&& treePath) {
289 test.subTest(checkNonZeroShapeFunctions(node.finiteElement()));
303struct EnableContinuityCheck
305 std::size_t order_ = 5;
308 template<
class JumpEvaluator>
309 auto localJumpContinuityCheck(
const JumpEvaluator& jumpEvaluator, std::size_t order,
double tol)
const
311 return [=](
const auto& intersection,
const auto& treePath,
const auto& insideNode,
const auto& outsideNode,
const auto& insideToOutside) {
312 using Intersection = std::decay_t<
decltype(intersection)>;
313 using Node = std::decay_t<
decltype(insideNode)>;
315 std::vector<int> isContinuous(insideNode.size(),
true);
316 const auto& quadRule = Dune::QuadratureRules<double, Intersection::mydimension>::rule(intersection.type(), order);
318 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
319 std::vector<std::vector<Range>> values;
320 std::vector<std::vector<Range>> neighborValues;
323 values.resize(quadRule.size());
324 neighborValues.resize(quadRule.size());
325 for(std::size_t k=0; k<quadRule.size(); ++k)
327 auto pointInElement = intersection.geometryInInside().global(quadRule[k].position());
328 auto pointInNeighbor = intersection.geometryInOutside().global(quadRule[k].position());
329 insideNode.finiteElement().localBasis().evaluateFunction(pointInElement, values[k]);
330 outsideNode.finiteElement().localBasis().evaluateFunction(pointInNeighbor, neighborValues[k]);
334 for(std::size_t i=0; i<insideNode.size(); ++i)
336 for(std::size_t k=0; k<quadRule.size(); ++k)
338 auto jump = values[k][i];
339 if (insideToOutside[i].has_value())
340 jump -= neighborValues[k][insideToOutside[i].value()];
341 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection, quadRule[k].position()) < tol);
348 auto localContinuityCheck()
const {
349 auto jumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
350 return jump.infinity_norm();
352 return localJumpContinuityCheck(jumpNorm, order_, tol_);
362struct EnableNormalContinuityCheck :
public EnableContinuityCheck
364 auto localContinuityCheck()
const {
365 auto normalJump = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
366 return jump * intersection.unitOuterNormal(x);
368 return localJumpContinuityCheck(normalJump, order_, tol_);
379struct EnableTangentialContinuityCheck :
public EnableContinuityCheck
381 auto localContinuityCheck()
const {
382 auto tangentialJumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
383 auto tangentialJump = jump - (jump * intersection.unitOuterNormal(x)) * intersection.unitOuterNormal(x);
384 return tangentialJump.two_norm();
386 return localJumpContinuityCheck(tangentialJumpNorm, order_, tol_);
397struct EnableCenterContinuityCheck :
public EnableContinuityCheck
399 template<
class JumpEvaluator>
400 auto localJumpCenterContinuityCheck(
const JumpEvaluator& jumpEvaluator,
double tol)
const
402 return [=](
const auto& intersection,
const auto& treePath,
const auto& insideNode,
const auto& outsideNode,
const auto& insideToOutside) {
403 using Node = std::decay_t<
decltype(insideNode)>;
404 using Range =
typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
406 std::vector<int> isContinuous(insideNode.size(),
true);
407 std::vector<Range> insideValues;
408 std::vector<Range> outsideValues;
410 insideNode.finiteElement().localBasis().evaluateFunction(intersection.geometryInInside().center(), insideValues);
411 outsideNode.finiteElement().localBasis().evaluateFunction(intersection.geometryInOutside().center(), outsideValues);
413 auto centerLocal = intersection.geometry().local(intersection.geometry().center());
416 for(std::size_t i=0; i<insideNode.size(); ++i)
418 auto jump = insideValues[i];
419 if (insideToOutside[i].has_value())
420 jump -= outsideValues[insideToOutside[i].value()];
421 isContinuous[i] = isContinuous[i] and (jumpEvaluator(jump, intersection, centerLocal) < tol);
427 auto localContinuityCheck()
const {
428 auto jumpNorm = [](
auto&&jump,
auto&& intersection,
auto&& x) ->
double {
429 return jump.infinity_norm();
431 return localJumpCenterContinuityCheck(jumpNorm, tol_);
447template<
class Basis,
class LocalCheck>
448Dune::TestSuite checkBasisContinuity(
const Basis& basis,
const LocalCheck& localCheck)
450 Dune::TestSuite test(
"Global continuity check of basis functions");
453 auto localView = basis.localView();
454 auto neighborLocalView = basis.localView();
456 for (
const auto& e : elements(basis.gridView()))
459 for(
const auto& intersection : intersections(basis.gridView(), e))
461 if (intersection.neighbor())
463 neighborLocalView.bind(intersection.outside());
465 Dune::TypeTree::forEachLeafNode(localView.tree(), [&](
const auto& insideNode,
auto&& treePath) {
466 const auto& outsideNode = Dune::TypeTree::child(neighborLocalView.tree(), treePath);
468 std::vector<std::optional<int>> insideToOutside;
469 insideToOutside.resize(insideNode.size());
472 for(std::size_t i=0; i<insideNode.size(); ++i)
474 for(std::size_t j=0; j<outsideNode.size(); ++j)
476 if (localView.index(insideNode.localIndex(i)) == neighborLocalView.index(outsideNode.localIndex(j)))
479 test.check(not insideToOutside[i].has_value())
480 <<
"Basis function " << localView.index(insideNode.localIndex(i))
481 <<
" appears twice in element " << elementStr(neighborLocalView.element(), basis.gridView());
482 insideToOutside[i] = j;
488 auto isContinuous = localCheck(intersection, treePath, insideNode, outsideNode, insideToOutside);
490 for(std::size_t i=0; i<insideNode.size(); ++i)
492 test.check(isContinuous[i])
493 <<
"Basis function " << localView.index(insideNode.localIndex(i))
494 <<
" is discontinuous across intersection of elements "
495 << elementStr(localView.element(), basis.gridView())
496 <<
" and " << elementStr(neighborLocalView.element(), basis.gridView());
505template<
class Basis,
class... Flags>
506Dune::TestSuite checkConstBasis(
const Basis& basis, Flags... flags)
508 Dune::TestSuite test(
"const basis check");
510 using GridView =
typename Basis::GridView;
513 test.check(Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, Basis>(),
"global basis concept check")
514 <<
"type passed to checkBasis() does not model the GlobalBasis concept";
517 auto localView = basis.localView();
518 for (
const auto& e : elements(basis.gridView()))
521 test.subTest(checkLocalView(basis, localView, flags...));
525 test.subTest(checkBasisIndices(basis));
529 auto flagTuple = std::tie(flags...);
530 Dune::Hybrid::forEach(flagTuple, [&](
auto&& flag) {
531 using Flag = std::decay_t<
decltype(flag)>;
532 if constexpr (std::is_base_of_v<EnableContinuityCheck, Flag>)
533 test.subTest(checkBasisContinuity(basis, flag.localContinuityCheck()));
539template<
class Basis,
class... Flags>
540Dune::TestSuite checkBasis(Basis& basis, Flags... flags)
542 Dune::TestSuite test(
"basis check");
545 test.subTest(checkConstBasis(basis,flags...));
550 test.subTest(checkConstBasis(copy,flags...));
552 test.subTest(checkConstBasis(copy,flags...));
556 auto gridView = basis.gridView();
557 basis.update(gridView);
Definition: polynomial.hh:17