...
 
Commits (37)
......@@ -32,7 +32,7 @@ am
test-driver
dependencies.m4
dune.css
build-cmake/
build*/
*~
compile
CMakeFiles
......
This diff is collapsed.
......@@ -4,10 +4,10 @@
#Name of the module
Module: dune-functions
Version: 3.0-dev
Version: 2.4
Maintainer: dune-devel@dune-project.org
#depending on
Depends: dune-localfunctions (>= 3.0) dune-grid (>= 3.0) dune-typetree
Depends: dune-localfunctions (>= 2.4) dune-grid (>= 2.4.1) dune-typetree
# We suggest dune-istl, because it is used by the example programs.
Suggests: dune-istl
Whitespace-Hook: Yes
add_subdirectory(common)
add_subdirectory(functions)
install(FILES
concept.hh
typelist.hh
typeutilities.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/common)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_COMMON_CONCEPT_HH
#define DUNE_COMMON_CONCEPT_HH
#include <type_traits>
#include <utility>
#include <tuple>
#include <dune/common/typeutilities.hh>
#include <dune/common/typelist.hh>
#include <dune/common/tupleutility.hh>
namespace Dune {
// Forward declaration
template<class C, class... T>
constexpr bool models();
/**
* \brief Namespace for concepts
*
* This namespace contains helper functions for
* concept definitions and the concept definitions
* themselves.
*/
namespace Concept {
/**
* \brief Base class for refined concepts.
*
* If a new concept should refine one or more existing concepts,
* this can be achieved by deriving the new concept from
* Refines<C1,...,CN> where C1, ..., CN are the concepts
* to be refined. If you want to refine several concepts
* they should all be put in a single Refines<...> base
* class.
*
* \tparam BaseConcepts The list of concepts to be refined.
*/
template<class... BaseConcepts>
struct Refines
{
typedef TypeList<BaseConcepts...> BaseConceptList;
};
namespace Imp {
// #############################################################################
// # All functions following here are implementation details
// # for the models() function below.
// #############################################################################
// Here is the implementation of the concept checking.
// The first two overloads do the magic for checking
// if the requirements of a concept are satisfied.
// The rest is just for checking base concepts in case
// of refinement.
// This overload is present if type substitution for
// C::require(T...) is successful, i.e., if the T...
// matches the requirement of C. In this case this
// overload is selected because PriorityTag<1>
// is a better match for PrioriryTag<42> than
// PriorityTag<0> in the default overload.
template<class C, class... T,
decltype(std::declval<C>().require(std::declval<T>()...), 0) =0>
constexpr std::true_type matchesRequirement(PriorityTag<1>)
{ return {}; }
// If the above overload is ruled out by SFINAE because
// the T... does not match the requirements of C, then
// this default overload drops in.
template<class C, class... T>
constexpr std::false_type matchesRequirement(PriorityTag<0>)
{ return {}; }
// An empty list C of concepts is always matched by T...
template<class...T>
constexpr bool modelsConceptList(TypeList<>)
{ return true; }
// A nonempty list C0,..,CN of concepts is modeled
// by T... if it models the concept C0
// and all concepts in the list C1,..,CN.
template<class...T, class C0, class... CC>
constexpr bool modelsConceptList(TypeList<C0, CC...>)
{ return models<C0, T...>() and modelsConceptList<T...>(TypeList<CC...>()); }
// If C is an unrefined concept, then T... models C
// if it matches the requirement of C.
template<class C, class... T>
constexpr bool modelsConcept(PriorityTag<0>)
{ return matchesRequirement<C, T...>(PriorityTag<42>()); }
// If C is a refined concept, then T... models C
// if it matches the requirement of C and of
// all base concepts.
//
// This overload is used if C::BaseConceptList exists
// due to its higher priority.
template<class C, class... T,
decltype(typename C::BaseConceptList(), 0) = 0>
constexpr bool modelsConcept(PriorityTag<1>)
{ return matchesRequirement<C, T...>(PriorityTag<42>()) and modelsConceptList<T...>(typename C::BaseConceptList()); }
// #############################################################################
// # All functions following here are implementation details for the
// # for the tupleEntriesModel() function below.
// #############################################################################
template<class C, class Tuple>
struct TupleEntriesModelHelper
{
template<class Accumulated, class T>
struct AccumulateFunctor
{
using type = typename std::integral_constant<bool, Accumulated::value and models<C, T>()>;
};
using Result = typename ReduceTuple<AccumulateFunctor, Tuple, std::true_type>::type;
};
} // namespace Dune::Concept::Imp
// #############################################################################
// # The method tupleEntriesModel() does the actual check if the types in a tuple
// # model a concept using the implementation details above.
// #############################################################################
template<class C, class Tuple>
constexpr auto tupleEntriesModel()
-> typename Imp::TupleEntriesModelHelper<C, Tuple>::Result
{
return {};
}
// #############################################################################
// # The following require*() functions are just helpers that allow to
// # propagate a failed check as substitution failure. This is usefull
// # inside of a concept definition.
// #############################################################################
// Helper function for use in concept definitions.
// If the passed value b is not true, the concept will to be satisfied.
template<bool b, typename std::enable_if<b, int>::type = 0>
constexpr bool requireTrue()
{
return true;
}
// Helper function for use in concept definitions.
template<class C, class... T, typename std::enable_if<models<C, T...>(), int>::type = 0>
constexpr bool requireConcept()
{
return true;
}
// Helper function for use in concept definitions.
// This allows to avoid using decltype
template<class C, class... T, typename std::enable_if<models<C, T...>(), int>::type = 0>
constexpr bool requireConcept(T&&... t)
{
return true;
}
// Helper function for use in concept definitions.
// This checks if the concept given as first type is modelled by all types in the tuple passed as argument
template<class C, class Tuple, typename std::enable_if<tupleEntriesModel<C, Tuple>(), int>::type = 0>
constexpr bool requireConceptForTupleEntries()
{
return true;
}
// Helper function for use in concept definitions.
// If the first passed type is not convertible to the second, the concept will not be satisfied.
template<class From, class To,
typename std::enable_if< std::is_convertible<From, To>::value, int>::type = 0>
constexpr bool requireConvertible()
{
return true;
}
// Helper function for use in concept definitions.
// If passed argument is not convertible to the first passed type, the concept will not be satisfied.
template<class To, class From,
typename std::enable_if< std::is_convertible<From, To>::value, int>::type = 0>
constexpr bool requireConvertible(const From&)
{
return true;
}
// Helper function for use in concept definitions.
// This will always evaluate to true. If just allow
// to turn a type into an expression. The failure happens
// already during substitution for the type argument.
template<typename T>
constexpr bool requireType()
{
return true;
}
// Helper function for use in concept definitions.
// If first passed type is not a base class of second type, the concept will not be satisfied.
template<class Base, class Derived,
typename std::enable_if< std::is_base_of<Base, Derived>::value, int>::type = 0>
constexpr bool requireBaseOf()
{
return true;
}
// Helper function for use in concept definitions.
// If first passed type is not a base class of first arguments type, the concept will not be satisfied.
template<class Base, class Derived,
typename std::enable_if< std::is_base_of<Base, Derived>::value, int>::type = 0>
constexpr bool requireBaseOf(const Derived&)
{
return true;
}
// Helper function for use in concept definitions.
// If the passed types are not the same, the concept will not be satisfied.
template<class A, class B,
typename std::enable_if< std::is_same<A, B>::value, int>::type = 0>
constexpr bool requireSameType()
{
return true;
}
} // namespace Dune::Concept
/**
* \brief Check if concept is modeled by given types
*
* This will check if the given concept is modeled by the given
* list of types. This is true if the list of types models all
* the base concepts that are refined by the given concept
* and if it satisfies all additional requirements of the latter.
*
* Notice that a concept may be defined for a list of interacting types.
* The function will check if the given list of types matches the requirements
* on the whole list. It does not check if each individual type in the list
* satisfies the concept.
*
* This concept check mechanism is inspired by the concept checking
* facility in Eric Nieblers range-v3. For more information please
* refer to the libraries project page https://github.com/ericniebler/range-v3
* or this blog entry: http://ericniebler.com/2013/11/23/concept-checking-in-c11.
* In fact the interface provided here is almost exactly the same as in range-v3.
* However the implementation differs, because range-v3 uses its own meta-programming
* library whereas our implementation is more straight forward.
*
* \tparam C The concept to check
* \tparam T The list of type to check against the concept
*
*/
template<class C, class... T>
constexpr bool models()
{
return Concept::Imp::modelsConcept<C, T...>(PriorityTag<42>());
}
} // namespace Dune
#endif // DUNE_COMMON_CONCEPT_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_COMMON_TEST_COLLECTORSTREAM_HH
#define DUNE_COMMON_TEST_COLLECTORSTREAM_HH
#include <sstream>
#include <string>
#include <functional>
#include <dune/common/typeutilities.hh>
namespace Dune {
/**
* \brief Data collector stream
*
* A class derived from std::ostringstream that allows to
* collect data via a temporary returned object. To facilitate
* this it stores a callback that is used to pass the collected
* data to its creator on destruction.
*
* In order to avoid passing the same data twice, copy construction
* is forbidden and only move construction is allowed.
*/
class CollectorStream : public std::ostringstream
{
public:
/**
* \brief Create from callback
*
* \tparam CallBack Type of callback. Must be convertible to std::function<void(std::string)>
* \param callBack A copy of this function will be stored and called on destruction.
*/
template<class CallBack,
Dune::disableCopyMove<CollectorStream, CallBack> = 0>
CollectorStream(CallBack&& callBack) :
callBack_(callBack)
{}
CollectorStream(const CollectorStream& other) = delete;
/**
* \brief Move constructor
*
* This will take over the data and callback from the
* moved from CollectorStream and disable the callback
* in the latter.
*/
CollectorStream(CollectorStream&& other) :
callBack_(other.callBack_)
{
(*this) << other.str();
other.callBack_ = [](std::string){};
}
/**
* \brief Destructor
*
* This calls the callback function given on creation
* passing all collected data as a single string argument.
*/
~CollectorStream()
{
callBack_(this->str());
}
private:
std::function<void(std::string)> callBack_;
};
} // namespace Dune
#endif // DUNE_COMMON_TEST_COLLECTORSTREAM_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_COMMON_TEST_TESTSUITE_HH
#define DUNE_COMMON_TEST_TESTSUITE_HH
#include <iostream>
#include <sstream>
#include <string>
#include <dune/common/exceptions.hh>
#include <dune/common/test/collectorstream.hh>
namespace Dune {
/**
* \brief A Simple helper class to organize your test suite
*
* Usage: Construct a TestSuite and call check() or require()
* with the condition to check and probably a name for this check.
* These methods return a stream such that you can pipe in an
* explanantion accomponied by respective data to give a reason
* for a test failure.
*/
class TestSuite
{
public:
enum ThrowPolicy
{
AlwaysThrow,
ThrowOnRequired
};
/**
* \brief Create TestSuite
*
* \param name A name to identify this TestSuite. Defaults to "".
* \param policy If AlwaysThrow any failing check will throw, otherwise only required checks will do.
*/
TestSuite(ThrowPolicy policy, std::string name="") :
name_(name),
checks_(0),
failedChecks_(0),
throwPolicy_(policy==AlwaysThrow)
{}
/**
* \brief Create TestSuite
*
* \param name A name to identify this TestSuite. Defaults to "".
* \param policy If AlwaysThrow any failing check will throw, otherwise only required checks will do. Defaults to ThrowOnRequired
*/
TestSuite(std::string name="", ThrowPolicy policy=ThrowOnRequired) :
name_(name),
checks_(0),
failedChecks_(0),
throwPolicy_(policy==AlwaysThrow)
{}
/**
* \brief Check condition
*
* This will throw an exception if the check fails and if the AlwaysThrow policy was used on creation.
*
* \param conditon Checks if this is true and increases the failure counter if not.
* \param name A name to identify this check. Defaults to ""
* \returns A CollectorStream that can be used to create a diagnostic message to be printed on failure.
*/
CollectorStream check(bool condition, std::string name="")
{
++checks_;
if (not condition)
++failedChecks_;
return CollectorStream([=](std::string reason) {
if (not condition)
this->announceCheckResult(throwPolicy_, "CHECK ", name, reason);
});
}
/**
* \brief Check a required condition condition
*
* This will always throw an exception if the check fails.
*
* \param conditon Checks if this is true and increases the failure counter if not.
* \param name A name to identify this check. Defaults to ""
* \returns A CollectorStream that can be used to create a diagnostic message to be printed on failure.
*/
CollectorStream require(bool condition, std::string name="")
{
++checks_;
if (not condition)
++failedChecks_;
return CollectorStream([=](std::string reason) {
if (not condition)
this->announceCheckResult(true, "REQUIRED CHECK", name, reason);
});
}
/**
* \brief Collect data from a sub-TestSuite
*
* This will incorporate the accumulated results of the sub-TestSuite
* into this one. If the sub-TestSuite failed, i.e., contained failed
* checks, a summary will be printed.
*/
void subTest(const TestSuite& subTest)
{
checks_ += subTest.checks_;
failedChecks_ += subTest.failedChecks_;
if (not subTest)
announceCheckResult(throwPolicy_, "SUBTEST", subTest.name(), std::to_string(subTest.failedChecks_)+"/"+std::to_string(subTest.checks_) + " checks failed in this subtest.");
}
/**
* \brief Check if this TestSuite failed
*
* \returns False if any of the executed tests failed, otherwise true.
*/
operator const bool () const
{
return (failedChecks_==0);
}
/**
* \brief Query name
*
* \returns Name of this TestSuite
*/
std::string name() const
{
return name_;
}
/**
* \brief Print a summary of this TestSuite
*
* \returns False if any of the executed tests failed, otherwise true.
*/
bool report() const
{
if (failedChecks_>0)
std::cout << composeMessage("TEST ", name(), std::to_string(failedChecks_)+"/"+std::to_string(checks_) + " checks failed in this test.") << std::endl;
return (failedChecks_==0);
}
/**
* \brief Exit the test.
*
* This wil print a summary of the test and return an integer
* to be used on program exit.
*
* \returns 1 if any of the executed tests failed, otherwise 0.
*/
int exit() const
{
return (report() ? 0: 1);
}
protected:
// Compose a diagnostic message
static std::string composeMessage(std::string type, std::string name, std::string reason)
{
std::ostringstream s;
s << type << " FAILED";
if (name!="")
s << "(" << name << ")";
s << ": ";
if (reason!="")
s << reason;
return s.str();
}
// Announce check results. To be called on failed checks
static void announceCheckResult(bool throwException, std::string type, std::string name, std::string reason)
{
std::string message = composeMessage(type, name, reason);
std::cout << message << std::endl;
if (throwException)
{
Dune::Exception ex;
ex.message(message);
throw ex;
}
}
std::string name_;
std::size_t checks_;
std::size_t failedChecks_;
bool throwPolicy_;
};
} // namespace Dune
#endif // DUNE_COMMON_TEST_TESTSUITE_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_COMMON_TYPELIST_HH
#define DUNE_COMMON_TYPELIST_HH
#include <type_traits>
#include <tuple>
namespace Dune {
/**
* \brief A simple type list
*
* \ingroup TypeUtilities
*
* The purpose of this is to encapsulate a list of types.
* This allows, e.g., to pack an argument-pack into one type.
* In contrast to a std::tuple a TypeList can be created
* without creating any object of the stored types.
*
* This can, e.g., be used for overload resolution
* with tag-dispatch where TypeList is used as tag.
* In combination with PriorityTag this allows to emulate
* partial specialization of function templates in
* a sane way, i.e., without the hassle of classic
* specialization of function templates
*/
template<class... T>
struct TypeList
{};
/**
* \brief Check if given type is a TypeList
*
* \ingroup TypeUtilities
*
* The result of the check is encoded in the
* base class of type std::integral_constant<bool, result>.
*/
template<class T>
struct IsTypeList : std::false_type {};
/**
* \copydoc IsTypeList
*
* \ingroup TypeUtilities
*/
template<class... T>
struct IsTypeList<TypeList<T...> > : std::true_type {};
/**
* \brief Check if given type is an empty TypeList
*
* \ingroup TypeUtilities
*
* The result of the check is encoded in the
* base class of type std::integral_constant<bool, result>.
*/
template<class T>
struct IsEmptyTypeList : std::integral_constant<bool, IsTypeList<T>() and std::is_same<T, TypeList<> >() > {};
template<class T>
struct TypeListSize {};
/**
* \brief Get size of TypeList
*
* \ingroup TypeUtilities
*
* The result of is encoded in the base class of
* type std::integral_constant<std::size_t, result>.
*/
template<class... T>
struct TypeListSize<TypeList<T...>> : std::integral_constant<std::size_t, sizeof...(T)> {};
template<std::size_t i, class T>
struct TypeListElement {};
/**
* \brief Get element of TypeList
*
* \ingroup TypeUtilities
*/
template<std::size_t i, class... T>
struct TypeListElement<i, TypeList<T...>>
{
/**
* \brief Export type of i-th element in TypeList
*
* \todo Implement without using std::tuple.
*/
using type = typename std::tuple_element<i, std::tuple<T...>>::type;
/**
* \brief Export type of i-th element in TypeList
*
* \todo Implement without using std::tuple.
*/
using Type = type;
};
/**
* \brief Shortcut for TypeListElement<i, T>::type;
*/
template<std::size_t i, class T>
using TypeListEntry_t = typename TypeListElement<i, T>::type;
} // namespace Dune
#endif // DUNE_COMMON_TYPELIST_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_COMMON_TYPEUTILITIES_HH
#define DUNE_COMMON_TYPEUTILITIES_HH
#include <cstddef>
#include <type_traits>
#include <tuple>
namespace Dune {
/**
* \file
* \brief Utilities for type computations, constarining overloads, ...
* \author Carsten Gräser
*/
/**
* \brief Helper to disable constructor as copy and move constructor
*
* \ingroup TypeUtilities
*
* Helper typedef to remove constructor with forwarding reference from
* overload set for copy and move constructor or assignment.
*/
template<class This, class... T>
using disableCopyMove = typename std::enable_if<
(not(std::is_same<This, typename std::tuple_element<0, std::tuple<typename std::decay<T>::type...> >::type >::value)
and not(std::is_base_of<This, typename std::tuple_element<0, std::tuple<typename std::decay<T>::type...> >::type >::value)), int>::type;
/**
* \brief Helper class for tagging priorities.
*
* \ingroup TypeUtilities
*
* When using multiple overloads of a function
* where some are removed from the overload set
* via SFINAE, the remaining overloads may be ambiguous.
* A prototypic example would be a default overload
* that should be used if the others do not apply.
*
* By adding additional arguments of type PriorityTag<k>
* with increasing priority k to all overloads and calling
* the method with PriorityTag<m> where m is larger or equal
* to the maximal used priority, those can be made unambiguous.
*
* In this case the matching overload with highest priority
* will be used. This is achieved by the fact that PriorityTag<k>
* derives from all types PriorityTag<i> with i les than k.
*
* \tparam priority The priority of this tag.
*/
template<std::size_t priority>
struct PriorityTag : public PriorityTag<priority-1>
{};
/**
* \brief Helper class for tagging priorities.
*
* \ingroup TypeUtilities
*
* PriorityTag<0> does not derive from any
* other PriorityTag.
*/
template<>
struct PriorityTag<0>
{};
} // namespace Dune
#endif // DUNE_COMMON_TYPEUTILITIES_HH
# tests that should build and run successfully
set(TESTS differentiablefunctiontest)
dune_add_test(SOURCES differentiablefunctiontest.cc)
set(TESTPROGS ${TESTS})
# We do not want want to build the tests during make all,
# but just build them on demand
add_directory_test_target(_test_target)
add_dependencies(${_test_target} ${TESTPROGS})
add_executable("differentiablefunctiontest"
differentiablefunctiontest.cc)
target_link_libraries("differentiablefunctiontest" "dunecommon")
# Add the tests to be executed
foreach(_TEST ${TESTPROGS})
add_test(${_TEST} ${_TEST})
endforeach(_TEST)
......@@ -7,6 +7,9 @@
#include <dune/common/shared_ptr.hh>
#include <dune/typetree/visitor.hh>
#include <dune/typetree/traversal.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
......
......@@ -8,6 +8,10 @@
#include <dune/common/typeutilities.hh>
namespace Dune {
// forward declaration
template< class K, int SIZE > class FieldVector;
namespace Functions {
......@@ -53,6 +57,14 @@ namespace Imp {
return {};
}
// Specialization for FieldVector
template<class T, int N>
constexpr auto staticSize(const FieldVector<T,N>*, const PriorityTag<10>&)
-> decltype(std::integral_constant<std::size_t,N>())
{
return {};
}
template<class T>
constexpr std::false_type hasStaticSize(const T* t, const PriorityTag<0>& p)
{
......
......@@ -21,4 +21,5 @@ install(FILES
pqknodalbasis.hh
nodes.hh
sizeinfo.hh
taylorhoodbasis.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/functions/functionspacebases)
......@@ -5,6 +5,8 @@
* \brief The B-spline global function space basis
*/
#include <numeric>
/** \todo Don't use this matrix */
#include <dune/common/dynmatrix.hh>
......
......@@ -14,7 +14,7 @@
#include <dune/functions/common/utility.hh>
#include <dune/functions/common/staticforloop.hh>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/sizeinfo.hh>
namespace Dune {
......@@ -104,7 +104,7 @@ public:
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
using SizePrefix = Dune::ReservedVector<size_type, MultiIndex::max_size()+1>;
using SizePrefix = Dune::ReservedVector<size_type, PrefixSizeTraits<MultiIndex>::max_size()+1>;
/** \brief Constructor for a given grid view object */
......@@ -119,7 +119,7 @@ public:
void initializeIndices()
{
staticForLoop<0, sizeof...(SF)>([&](const auto& i) {
staticForLoop<0, sizeof...(SF)>([&](auto i) {
std::get<i.value>(subFactories_).initializeIndices();
});
}
......@@ -135,7 +135,7 @@ public:
Node<TP> node(const TP& tp) const
{
auto node = Node<TP>(tp);
staticForLoop<0, sizeof...(SF)>([&](const auto& i){
staticForLoop<0, sizeof...(SF)>([&](auto i){
node.setChild(std::get<i.value>(subFactories_).node(TypeTree::push_back(tp, i)), i);
});
return node;
......
......@@ -188,8 +188,8 @@ protected:
* \param nodeToRangeEntry Polymorphic functor mapping local ansatz nodes to range-indices of given function
* \param bitVector A vector with flags marking all DOFs that should be interpolated
*/
template <class B, class TP, class NTRE, class C, class F, class BV>
void interpolateTreeSubset(const B& basis, const TP& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry, const BV& bv)
template <class B, class... TreeIndices, class NTRE, class C, class F, class BV>
void interpolateTreeSubset(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry, const BV& bv)
{
using GridView = typename B::GridView;
using Element = typename GridView::template Codim<0>::Entity;
......@@ -233,22 +233,22 @@ void interpolateTreeSubset(const B& basis, const TP& treePath, C&& coeff, const
}
template <class B, class TP, class C, class F, class BV>
void interpolateTreeSubset(const B& basis, const TP& treePath, C&& coeff, const F& f, const BV& bitVector)
template <class B, class... TreeIndices, class C, class F, class BV>
void interpolateTreeSubset(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const BV& bitVector)
{
interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), bitVector);
}
template <class B, class TP, class NTRE, class C, class F>
void interpolateTree(const B& basis, const TP& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry)
template <class B, class... TreeIndices, class NTRE, class C, class F>
void interpolateTree(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry)
{
interpolateTreeSubset(basis, treePath, coeff, f, nodeToRangeEntry, Imp::AllTrueBitSetVector());
}
template <class B, class TP, class C, class F>
void interpolateTree(const B& basis, const TP& treePath, C&& coeff, const F& f)
template <class B, class... TreeIndices, class C, class F>
void interpolateTree(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f)
{
interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), Imp::AllTrueBitSetVector());
}
......@@ -271,13 +271,51 @@ void interpolateTree(const B& basis, const TP& treePath, C&& coeff, const F& f)
* \param f Function to interpolate
* \param bitVector A vector with flags marking ald DOFs that should be interpolated
*/
template <class B, class TP, class C, class F, class BV>
void interpolate(const B& basis, const TP& treePath, C&& coeff, const F& f, const BV& bitVector)
template <class B, class... TreeIndices, class C, class F, class BV>
void interpolate(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const BV& bitVector)
{
interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), bitVector);
}
namespace Imp {
template<class... T>
constexpr std::true_type isHybridTreePath(const TypeTree::HybridTreePath<T...>&)
{ return {}; }
template<class T>
constexpr std::false_type isHybridTreePath(const T&)
{ return {}; }
template<class T>
constexpr auto isHybridTreePath() -> decltype(isHybridTreePath(std::declval<std::decay_t<T>>()))
{ return {}; }
}
/**
* \brief Interpolate given function in discrete function space
*
* Interpolation is done wrt the leaf node of the ansatz tree
* corresponding to the given tree path. Only vector coefficients marked as 'true' in the
* bitVector argument are interpolated. Use this, e.g., to interpolate Dirichlet boundary values.
*
* Notice that this will only work if the range type of f and
* the block type of coeff are compatible and supported by
* FlatVectorBackend.
*
* \param basis Global function space basis of discrete function space
* \param coeff Coefficient vector to represent the interpolation
* \param f Function to interpolate
* \param bitVector A vector with flags marking all DOFs that should be interpolated
*/
template <class B, class C, class F, class BV,
std::enable_if_t<not Imp::isHybridTreePath<C>(), int> = 0>
void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
{
auto root = Dune::TypeTree::hybridTreePath();
interpolateTreeSubset(basis, root, coeff, f, makeDefaultNodeToRangeMap(basis, root), bitVector);
}
......@@ -316,8 +354,8 @@ void interpolate(const B& basis, C&& coeff, const F& f)
* \param coeff Coefficient vector to represent the interpolation
* \param f Function to interpolate
*/
template <class B, class TreePath, class C, class F>
void interpolate(const B& basis, const TreePath& treePath, C&& coeff, const F& f)
template <class B, class... TreeIndices, class C, class F>
void interpolate(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f)
{
interpolate (basis, treePath, coeff, f, Imp::AllTrueBitSetVector());
}
......
......@@ -12,7 +12,7 @@
#include <dune/functions/common/utility.hh>
#include <dune/functions/common/staticforloop.hh>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/sizeinfo.hh>
namespace Dune {
......@@ -78,7 +78,7 @@ public:
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
using SizePrefix = Dune::ReservedVector<size_type, SubFactory::SizePrefix::max_size()+1>;
using SizePrefix = Dune::ReservedVector<size_type, PrefixSizeTraits<typename SubFactory::SizePrefix>::max_size()+1>;
private:
......@@ -112,7 +112,7 @@ public:
Node<TP> node(const TP& tp) const
{
auto node = Node<TP>(tp);
for(int i=0; i<children; ++i)
for (std::size_t i=0; i<children; ++i)
node.setChild(i, subFactory_.node(TypeTree::push_back(tp, i)));
return node;
}
......
......@@ -12,6 +12,7 @@
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
namespace Dune {
......@@ -237,10 +238,10 @@ protected:
/** \brief Nodal basis of a scalar first-order Lagrangian finite element space
*
* \tparam GV The GridView that the space is defined on
* \tparam ST The type used for local indices; global indices are std::array<ST,1>
* \tparam ST The type used for local indices; global indices are FlatMultiIndex<ST>
*/
template<typename GV, class ST = std::size_t>
using PQ1NodalBasis = DefaultGlobalBasis<PQ1NodeFactory<GV, std::array<ST, 1>, ST> >;
using PQ1NodalBasis = DefaultGlobalBasis<PQ1NodeFactory<GV, FlatMultiIndex<ST>, ST> >;
} // end namespace Functions
} // end namespace Dune
......
......@@ -56,14 +56,22 @@ public:
// Precompute the number of dofs per entity type
const static int dofsPerEdge = k-1;
const static int dofsPerTriangle = (k-1)*(k-2)/2;
const static int dofsPerQuad = (k-1)*(k-1);
const static int dofsPerTetrahedron = ((k-3)*(k-2)*(k-1)/6 > 0) ? (k-3)*(k-2)*(k-1)/6 : 0;
const static int dofsPerPrism = (k-1)*(k-1)*(k-2)/2;
const static int dofsPerHexahedron = (k-1)*(k-1)*(k-1);
const static int dofsPerPyramid = ((k-2)*(k-1)*(2*k-3))/6;
const static int dofsPerVertex =
k == 0 ? (dim == 0 ? 1 : 0) : 1;
const static int dofsPerEdge =
k == 0 ? (dim == 1 ? 1 : 0) : k-1;
const static int dofsPerTriangle =
k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-2)/2;
const static int dofsPerQuad =
k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-1);
const static int dofsPerTetrahedron =
k == 0 ? (dim == 3 ? 1 : 0) : (k-3)*(k-2)*(k-1)/6;
const static int dofsPerPrism =
k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-2)/2;
const static int dofsPerHexahedron =
k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-1);
const static int dofsPerPyramid =
k == 0 ? (dim == 3 ? 1 : 0) : (k-2)*(k-1)*(2*k-3)/6;
template<class TP>
using Node = PQkNode<GV, k, size_type, TP>;
......@@ -85,7 +93,7 @@ public:
void initializeIndices()
{
vertexOffset_ = 0;
edgeOffset_ = vertexOffset_ + gridView_.size(dim);
edgeOffset_ = vertexOffset_ + dofsPerVertex * gridView_.size(dim);
triangleOffset_ = edgeOffset_ + dofsPerEdge * gridView_.size(dim-1);
GeometryType triangle;
......@@ -135,14 +143,18 @@ public:
switch (dim)
{
case 1:
return gridView_.size(1) + dofsPerEdge*gridView_.size(0);
return dofsPerVertex * gridView_.size(dim)
+ dofsPerEdge*gridView_.size(dim-1);
case 2:
{
GeometryType triangle, quad;
triangle.makeTriangle();
quad.makeQuadrilateral();
return gridView_.size(dim) + dofsPerEdge*gridView_.size(1)
+ dofsPerTriangle*gridView_.size(triangle) + dofsPerQuad*gridView_.size(quad);
return dofsPerVertex * gridView_.size(dim)
+ dofsPerEdge * gridView_.size(dim-1)
+ dofsPerTriangle * gridView_.size(triangle)
+ dofsPerQuad * gridView_.size(quad);
}
case 3:
{
......@@ -153,10 +165,14 @@ public:
pyramid.makePyramid();
prism.makePrism();
hexahedron.makeCube(3);
return gridView_.size(dim) + dofsPerEdge*gridView_.size(2)
+ dofsPerTriangle*gridView_.size(triangle) + dofsPerQuad*gridView_.size(quad)
+ dofsPerTetrahedron*gridView_.size(tetrahedron) + dofsPerPyramid*gridView_.size(pyramid)
+ dofsPerPrism*gridView_.size(prism) + dofsPerHexahedron*gridView_.size(hexahedron);
return dofsPerVertex * gridView_.size(dim)
+ dofsPerEdge * gridView_.size(dim-1)
+ dofsPerTriangle * gridView_.size(triangle)
+ dofsPerQuad * gridView_.size(quad)
+ dofsPerTetrahedron * gridView_.size(tetrahedron)
+ dofsPerPyramid * gridView_.size(pyramid)
+ dofsPerPrism * gridView_.size(prism)
+ dofsPerHexahedron * gridView_.size(hexahedron);
}
}
DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
......@@ -306,13 +322,15 @@ public:
size_t dofDim = dim - localKey.codim();
if (dofDim==0) { // vertex dof
return { gridIndexSet.subIndex(element,localKey.subEntity(),dim) };
return {{ gridIndexSet.subIndex(element,localKey.subEntity(),dim) }};
}
if (dofDim==1)
{ // edge dof
if (dim==1) // element dof -- any local numbering is fine
return { nodeFactory_->edgeOffset_ + (k-1)*gridIndexSet.subIndex(element,0,0) + localKey.index() };
if (dim==1) // element dof -- any local numbering is fine
return {{ nodeFactory_->edgeOffset_
+ nodeFactory_->dofsPerEdge * gridIndexSet.subIndex(element,0,0)
+ localKey.index() }};
else
{
const Dune::ReferenceElement<double,dim>& refElement
......@@ -323,9 +341,13 @@ public:
size_t v0 = gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
size_t v1 = gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
bool flip = (v0 > v1);
return { (flip)
? nodeFactory_->edgeOffset_ + (k-1)*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) + (k-2)-localKey.index()
: nodeFactory_->edgeOffset_ + (k-1)*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) + localKey.index() };
return {{ (flip)
? nodeFactory_->edgeOffset_
+ nodeFactory_->dofsPerEdge*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())
+ (nodeFactory_->dofsPerEdge-1)-localKey.index()
: nodeFactory_->edgeOffset_
+ nodeFactory_->dofsPerEdge*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())
+ localKey.index() }};
}
}
......@@ -336,12 +358,12 @@ public:
if (element.type().isTriangle())
{
const int interiorLagrangeNodesPerTriangle = (k-1)*(k-2)/2;
return { nodeFactory_->triangleOffset_ + interiorLagrangeNodesPerTriangle*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->triangleOffset_ + interiorLagrangeNodesPerTriangle*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
}
else if (element.type().isQuadrilateral())
{
const int interiorLagrangeNodesPerQuadrilateral = (k-1)*(k-1);
return { nodeFactory_->quadrilateralOffset_ + interiorLagrangeNodesPerQuadrilateral*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->quadrilateralOffset_ + interiorLagrangeNodesPerQuadrilateral*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
}
else
DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
......@@ -356,7 +378,7 @@ public:
if (k==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
DUNE_THROW(Dune::NotImplemented, "PQkNodalBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
return { nodeFactory_->triangleOffset_ + gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) };
return {{ nodeFactory_->triangleOffset_ + gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) }};
}
}
......@@ -365,13 +387,13 @@ public:
if (dim==3) // element dof -- any local numbering is fine
{
if (element.type().isTetrahedron())
return { nodeFactory_->tetrahedronOffset_ + NodeFactory::dofsPerTetrahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->tetrahedronOffset_ + NodeFactory::dofsPerTetrahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else if (element.type().isHexahedron())
return { nodeFactory_->hexahedronOffset_ + NodeFactory::dofsPerHexahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->hexahedronOffset_ + NodeFactory::dofsPerHexahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else if (element.type().isPrism())
return { nodeFactory_->prismOffset_ + NodeFactory::dofsPerPrism*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->prismOffset_ + NodeFactory::dofsPerPrism*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else if (element.type().isPyramid())
return { nodeFactory_->pyramidOffset_ + NodeFactory::dofsPerPyramid*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->pyramidOffset_ + NodeFactory::dofsPerPyramid*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else
DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
} else
......
......@@ -3,12 +3,38 @@
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_SIZEINFO_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_SIZEINFO_HH
#include <dune/common/reservedvector.hh>
#include <array>
namespace Dune {
namespace Functions {
/**
* \brief Traits class encapsulating size information of the prefix
*/
template<typename P>
struct PrefixSizeTraits
{
using size_type = typename P::size_type;
//! Returns the maximum length of the prefix.
static constexpr size_type max_size()
{
return P::max_size();
}
};
template<class T, int n>
struct PrefixSizeTraits<Dune::ReservedVector<T,n>>
{
using size_type = typename Dune::ReservedVector<T,n>::size_type;
//! Returns the maximum length of the prefix.
static constexpr size_type max_size()
{
return n;
}
};
/**
* \brief A class encapsulating size information
......
# tests that should build and run successfully
# Path to the example grid files in dune-grid
add_definitions(-DDUNE_GRID_EXAMPLE_GRIDS_PATH=\"${DUNE_GRID_EXAMPLE_GRIDS_PATH}\")
dune_add_test(SOURCES gridviewfunctionspacebasistest.cc)
# tests that should build and run successfully
set(TESTS compositebasistest gridviewfunctionspacebasistest taylorhoodbasistest hierarchicvectorwrappertest)
set(TESTPROGS ${TESTS})
# We do not want want to build the tests during make all,
# but just build them on demand
add_directory_test_target(_test_target)
add_dependencies(${_test_target} ${TESTPROGS})
add_executable("compositebasistest"
compositebasistest.cc)
add_executable("gridviewfunctionspacebasistest"
gridviewfunctionspacebasistest.cc)
add_executable("taylorhoodbasistest"
taylorhoodbasistest.cc)
add_executable("hierarchicvectorwrappertest"
hierarchicvectorwrappertest.cc)
target_link_libraries("compositebasistest" "dunecommon" "dunegeometry")
target_link_libraries("gridviewfunctionspacebasistest" "dunecommon" "dunegeometry")
target_link_libraries("taylorhoodbasistest" "dunecommon" "dunegeometry")
target_link_libraries("hierarchicvectorwrappertest" "dunecommon" "dunegeometry")
dune_add_test(SOURCES taylorhoodbasistest.cc)
add_dune_ug_flags("gridviewfunctionspacebasistest")
dune_add_test(SOURCES hierarchicvectorwrappertest.cc)
# Add the tests to be executed
foreach(_TEST ${TESTPROGS})
add_test(${_TEST} ${_TEST})
endforeach(_TEST)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <dune/common/function.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
using namespace Dune;
int main (int argc, char *argv[]) try
{
// Set up MPI, if available
MPIHelper::instance(argc, argv);
///////////////////////////////////
// Generate the grid
///////////////////////////////////
const int dim = 2;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> l(1);
std::array<int,dim> elements = {4, 4};
GridType grid(l,elements);
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
/////////////////////////////////////////////////////////
// Choose a finite element space
/////////////////////////////////////////////////////////
using namespace Functions::BasisTags;
using namespace Functions::BasisBuilder;
auto basis = makeBasis(
gridView,
composite(
pq<1>(),
pq<1>(),
pq<1>()
));
return 0;
}
// Error handling
catch (Exception e)
{
std::cout << e << std::endl;
}
......@@ -3,6 +3,7 @@
#include <config.h>
#include <iostream>
#include <numeric>
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
......@@ -269,6 +270,10 @@ void testOnStructuredGrid()
if (dim<3) // Currently not implemented for dim >= 3
testScalarBasis(pq4Basis, true);
// Test PQkNodalBasis for the corner case k == 0
PQkNodalBasis<GridView, 0> pq0Basis(gridView);
testScalarBasis(pq0Basis, true);
// Test LagrangeDGBasis for k==1
LagrangeDGBasis<GridView, 1> lagrangeDG1Basis(gridView);
testScalarBasis(lagrangeDG1Basis, true);
......