Commit d2aafda1 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

expression and DirichletBC with lambda functions

parent 7bd09e85
...@@ -97,7 +97,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc ...@@ -97,7 +97,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/DOFMatrix.cc ${SOURCE_DIR}/DOFMatrix.cc
${SOURCE_DIR}/DOFVector.cc ${SOURCE_DIR}/DOFVector.cc
${SOURCE_DIR}/Debug.cc ${SOURCE_DIR}/Debug.cc
${SOURCE_DIR}/DirichletBC.cc # ${SOURCE_DIR}/DirichletBC.cc
${SOURCE_DIR}/DualTraverse.cc ${SOURCE_DIR}/DualTraverse.cc
${SOURCE_DIR}/ElInfo.cc ${SOURCE_DIR}/ElInfo.cc
${SOURCE_DIR}/ElInfo1d.cc ${SOURCE_DIR}/ElInfo1d.cc
......
...@@ -28,6 +28,7 @@ ...@@ -28,6 +28,7 @@
#include <boost/numeric/mtl/mtl.hpp> #include <boost/numeric/mtl/mtl.hpp>
#include "OpenMP.h" #include "OpenMP.h"
#include "Config.h"
namespace AMDiS { namespace AMDiS {
...@@ -134,7 +135,11 @@ namespace AMDiS { ...@@ -134,7 +135,11 @@ namespace AMDiS {
struct BoundaryObject; struct BoundaryObject;
struct AtomicBoundary; struct AtomicBoundary;
#if HAS_VARIADIC_TEMPLATES
template<typename ReturnType, typename... Args > class AbstractFunction;
#else
template<typename ReturnType, typename ArgumentType> class AbstractFunction; template<typename ReturnType, typename ArgumentType> class AbstractFunction;
#endif
template<typename ReturnType, template<typename ReturnType,
typename ArgumentType1, typename ArgumentType1,
typename ArgumentType2> class BinaryAbstractFunction; typename ArgumentType2> class BinaryAbstractFunction;
...@@ -153,6 +158,7 @@ namespace AMDiS { ...@@ -153,6 +158,7 @@ namespace AMDiS {
template<typename T> class DOFVector; template<typename T> class DOFVector;
template<typename T> class DimVec; template<typename T> class DimVec;
template<typename T> class DimMat; template<typename T> class DimMat;
template<typename T> class DirichletBC;
// template<typename ITLSolver> class ITL_LinearSolverInterface; // template<typename ITLSolver> class ITL_LinearSolverInterface;
template<typename T, typename MatT, typename VecT > class ITL_Preconditioner; template<typename T, typename MatT, typename VecT > class ITL_Preconditioner;
template<typename T> class Matrix; template<typename T> class Matrix;
......
...@@ -33,43 +33,6 @@ ...@@ -33,43 +33,6 @@
namespace AMDiS { namespace AMDiS {
/**
* \ingroup Common
*
* \brief
* An AbstractFunction object represents a function
* f : ArgumentType -> ReturnType.
*
* AbstractFunction is a pure virtual class interface class.
* To create your own function you have to derive AbstractFunction and
* overload operator().
*/
// #ifndef HAS_VARIADIC_TEMPLATES
template<typename ReturnType, typename ArgumentType>
class AbstractFunction
{
public:
/// Constructor.
AbstractFunction(int degree = 0) :
degree_(degree)
{}
virtual ~AbstractFunction() {}
/// Returns \ref degree_.
inline int getDegree() const
{
return degree_;
}
/// Deligates the evaluation to overriden method f.
virtual ReturnType operator()(const ArgumentType& x) const = 0;
protected:
int degree_;
};
// #endif
/** /**
* \ingroup Common * \ingroup Common
* *
...@@ -175,7 +138,68 @@ namespace AMDiS { ...@@ -175,7 +138,68 @@ namespace AMDiS {
int degree_; int degree_;
}; };
// #ifndef HAS_VARIADIC_TEMPLATES #if HAS_VARIADIC_TEMPLATES
// C++11 implementation of AbstractFunction with arbitrary nr of arguments using variadic templates
/**
* \ingroup Common
*
* \brief
* An AbstractFunction object represents a function
* f : ArgumentTypes... -> ReturnType.
*
* AbstractFunction is a pure virtual class interface class.
* To create your own function you have to derive AbstractFunction and
* overload operator().
*/
template< typename ReturnType, typename... Args >
class AbstractFunction
{
public:
AbstractFunction(int degree = 0)
: degree_(degree)
{}
virtual ~AbstractFunction() {}
/// Returns \ref degree_.
inline int getDegree() const
{
return degree_;
}
/// function evaluation.
virtual ReturnType operator()(Args const&... args) const = 0;
protected:
int degree_;
};
#else
template<typename ReturnType, typename ArgumentType>
class AbstractFunction
{
public:
/// Constructor.
AbstractFunction(int degree = 0) :
degree_(degree)
{}
virtual ~AbstractFunction() {}
/// Returns \ref degree_.
inline int getDegree() const
{
return degree_;
}
/// Deligates the evaluation to overriden method f.
virtual ReturnType operator()(const ArgumentType& x) const = 0;
protected:
int degree_;
};
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
// test of AbstractFunction with arbitrary number of arguments // test of AbstractFunction with arbitrary number of arguments
...@@ -226,33 +250,7 @@ namespace AMDiS { ...@@ -226,33 +250,7 @@ namespace AMDiS {
BOOST_PP_REPEAT_FROM_TO(1, 11, ABSTRACT_FUNCTION_MACRO, nil) BOOST_PP_REPEAT_FROM_TO(1, 11, ABSTRACT_FUNCTION_MACRO, nil)
// #else #endif
//
// /// C++11 implementation of abstract functions with arbitrary nr of arguments using variadic templates
// template< typename ReturnType, typename... Ts >
// class AbstractFunction
// {
// public:
// AbstractFunction(int degree = 0)
// : degree_(degree)
// {}
//
// virtual ~AbstractFunction() {}
//
// /// Returns \ref degree_.
// inline int getDegree() const
// {
// return degree_;
// }
//
// /// function evaluation.
// virtual ReturnType operator()(Ts const&... args) const = 0;
//
// protected:
// int degree_;
// };
//
// #endif
} }
#endif // AMDIS_ABSTRACTFUNCTION_H #endif // AMDIS_ABSTRACTFUNCTION_H
...@@ -72,21 +72,23 @@ namespace AMDiS { ...@@ -72,21 +72,23 @@ namespace AMDiS {
const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
for (int i = 0; i < nBasFcts; i++) {
if (localBound[i] == boundaryType) {
double value = 0.0;
if (f) { if (f) {
for (int i = 0; i < nBasFcts; i++)
if (localBound[i] == boundaryType) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords); elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords);
value = (*f)(worldCoords); double value = (*f)(worldCoords);
} else { vector->setDirichletDofValue(dofIndices[i], value);
if (dofVec) (*vector)[dofIndices[i]] = value;
value = (*dofVec)[dofIndices[i]];
else
ERROR_EXIT("There is something wrong!\n");
} }
} else if (dofVec) {
for (int i = 0; i < nBasFcts; i++)
if (localBound[i] == boundaryType) {
double value = value = (*dofVec)[dofIndices[i]];
vector->setDirichletDofValue(dofIndices[i], value); vector->setDirichletDofValue(dofIndices[i], value);
(*vector)[dofIndices[i]] = value; (*vector)[dofIndices[i]] = value;
} }
} else {
ERROR_EXIT("No data provided to assemble DirichletBC!\n");
} }
} }
......
...@@ -30,7 +30,48 @@ ...@@ -30,7 +30,48 @@
#include "AbstractFunction.h" #include "AbstractFunction.h"
#include "FixVec.h" #include "FixVec.h"
namespace AMDiS { namespace AMDiS
{
// some tags used for tag-dispatching
struct _value_by_abstractfunction {};
struct _value_by_dofvector {};
struct _value_by_function {};
namespace detail
{
template <class Tag>
struct ValueContainer {};
template <>
struct ValueContainer<_value_by_abstractfunction>
{
ValueContainer(AbstractFunction<double, WorldVector<double> > *fct) : value(*fct) {};
ValueContainer(AbstractFunction<double, WorldVector<double> > &fct) : value(fct) {};
AbstractFunction<double, WorldVector<double> > &value;
};
template <>
struct ValueContainer<_value_by_dofvector>
{
ValueContainer(DOFVectorBase<double> *vec) : value(vec) {};
ValueContainer(DOFVectorBase<double> &vec) : value(&vec) {};
DOFVectorBase<double> *value;
};
#if __cplusplus > 199711L
template <>
struct ValueContainer<_value_by_function>
{
ValueContainer(std::function<double(WorldVector<double>)> fct) : value(fct) {};
std::function<double(WorldVector<double>)> value;
};
#endif
} // end namespace detail
/** /**
* \ingroup Assembler * \ingroup Assembler
...@@ -41,6 +82,7 @@ namespace AMDiS { ...@@ -41,6 +82,7 @@ namespace AMDiS {
* the row corresponding to a Dirichlet dof is replaced by a row containing * the row corresponding to a Dirichlet dof is replaced by a row containing
* only a 1.0 in the diagonal. * only a 1.0 in the diagonal.
*/ */
template <class ValueTag>
class DirichletBC : public BoundaryCondition class DirichletBC : public BoundaryCondition
{ {
public: public:
...@@ -51,24 +93,33 @@ namespace AMDiS { ...@@ -51,24 +93,33 @@ namespace AMDiS {
const FiniteElemSpace *colFeSpace = NULL, const FiniteElemSpace *colFeSpace = NULL,
bool apply = true); bool apply = true);
#if __cplusplus > 199711L
/// Constructor.
DirichletBC(BoundaryType type,
std::function<double(WorldVector<double>)> fct,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true);
#endif
/// Constructor. /// Constructor.
DirichletBC(BoundaryType type, DirichletBC(BoundaryType type,
DOFVectorBase<double> *vec, DOFVectorBase<double> *vec,
bool apply = true); bool apply = true);
/// Implementation of BoundaryCondition::fillBoundaryCondition(). /// Implementation of BoundaryCondition::fillBoundaryCondition().
void fillBoundaryCondition(DOFMatrix* matrix, virtual void fillBoundaryCondition(DOFMatrix* matrix,
ElInfo* elInfo, ElInfo* elInfo,
const DegreeOfFreedom* dofIndices, const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound, const BoundaryType* localBound,
int nBasFcts); int nBasFcts) override;
/// Implementation of BoundaryCondition::fillBoundaryCondition(). /// Implementation of BoundaryCondition::fillBoundaryCondition().
void fillBoundaryCondition(DOFVectorBase<double>* vector, virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo, ElInfo* elInfo,
const DegreeOfFreedom* dofIndices, const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound, const BoundaryType* localBound,
int nBasFcts); int nBasFcts) override;
/// ///
void initVector(DOFVectorBase<double>*); void initVector(DOFVectorBase<double>*);
...@@ -93,25 +144,30 @@ namespace AMDiS { ...@@ -93,25 +144,30 @@ namespace AMDiS {
return applyBC; return applyBC;
} }
inline AbstractFunction<double, WorldVector<double> > *getF()
{
return f;
}
inline DOFVectorBase<double> *getDOFVector()
{
return dofVec;
}
protected: protected:
/// Function which is evaluated at world coords of Dirichlet dofs. void fillBC(_value_by_abstractfunction,
AbstractFunction<double, WorldVector<double> > *f; DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts);
/// void fillBC(_value_by_function,
WorldVector<double> worldCoords; DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts);
void fillBC(_value_by_dofvector,
DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts);
/// DOFVector containing the boundary values protected:
DOFVectorBase<double> *dofVec; detail::ValueContainer<ValueTag> container;
/// Defines, if the boundary condition must be applied to the matrix. See /// Defines, if the boundary condition must be applied to the matrix. See
/// comment of \ref BoundaryCondition::applyBoundaryCondition. /// comment of \ref BoundaryCondition::applyBoundaryCondition.
...@@ -120,4 +176,6 @@ namespace AMDiS { ...@@ -120,4 +176,6 @@ namespace AMDiS {
} }
#include "DirichletBC.hh"
#endif #endif
...@@ -36,8 +36,6 @@ namespace AMDiS { ...@@ -36,8 +36,6 @@ namespace AMDiS {
#define MAX_DIM 3 #define MAX_DIM 3
#define MAX_DEGREE 4 #define MAX_DEGREE 4
template<typename ReturnType, typename ArgumentType> class AbstractFunction;
/** \ingroup FEMSpace /** \ingroup FEMSpace
* \brief * \brief
* Lagrange basis functions. Sub class of BasisFunction * Lagrange basis functions. Sub class of BasisFunction
......
...@@ -1383,10 +1383,10 @@ namespace AMDiS { ...@@ -1383,10 +1383,10 @@ namespace AMDiS {
boundaryConditionSet = true; boundaryConditionSet = true;
DirichletBC *dirichletApply = DirichletBC<_value_by_abstractfunction> *dirichletApply =
new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], true); new DirichletBC<_value_by_abstractfunction>(type, b, componentSpaces[row], componentSpaces[col], true);
DirichletBC *dirichletNotApply = DirichletBC<_value_by_abstractfunction> *dirichletNotApply =
new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], false); new DirichletBC<_value_by_abstractfunction>(type, b, componentSpaces[row], componentSpaces[col], false);
for (int i = 0; i < nComponents; i++) for (int i = 0; i < nComponents; i++)
if (systemMatrix && (*systemMatrix)[row][i]) { if (systemMatrix && (*systemMatrix)[row][i]) {
...@@ -1403,6 +1403,38 @@ namespace AMDiS { ...@@ -1403,6 +1403,38 @@ namespace AMDiS {
} }
#if __cplusplus > 199711L
void ProblemStatSeq::addDirichletBC(BoundaryType type, int row, int col,
std::function<double(WorldVector<double>)> b)
{
FUNCNAME("ProblemStat::addDirichletBC()");
TEST_EXIT(row >= 0 && row < nComponents)("Wrong row number: %d\n", row);
TEST_EXIT(col >= 0 && col < nComponents)("Wrong col number: %d\n", col);
boundaryConditionSet = true;
DirichletBC<_value_by_function> *dirichletApply =
new DirichletBC<_value_by_function>(type, b, componentSpaces[row], componentSpaces[col], true);
DirichletBC<_value_by_function> *dirichletNotApply =
new DirichletBC<_value_by_function>(type, b, componentSpaces[row], componentSpaces[col], false);
for (int i = 0; i < nComponents; i++)
if (systemMatrix && (*systemMatrix)[row][i]) {
if (i == col)
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletApply);
else
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletNotApply);
}
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
if (solution)
solution->getDOFVector(col)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
}
#endif
void ProblemStatSeq::addDirichletBC(BoundaryType type, int row, int col, void ProblemStatSeq::addDirichletBC(BoundaryType type, int row, int col,
DOFVector<double> *vec) DOFVector<double> *vec)
{ {
...@@ -1413,8 +1445,8 @@ namespace AMDiS { ...@@ -1413,8 +1445,8 @@ namespace AMDiS {
boundaryConditionSet = true; boundaryConditionSet = true;
DirichletBC *dirichletApply = new DirichletBC(type, vec, true); DirichletBC<_value_by_dofvector> *dirichletApply = new DirichletBC<_value_by_dofvector>(type, vec, true);
DirichletBC *dirichletNotApply = new DirichletBC(type, vec, false); DirichletBC<_value_by_dofvector> *dirichletNotApply = new DirichletBC<_value_by_dofvector>(type, vec, false);
for (int i = 0; i < nComponents; i++) for (int i = 0; i < nComponents; i++)
if (systemMatrix && (*systemMatrix)[row][i]) { if (systemMatrix && (*systemMatrix)[row][i]) {
......
...@@ -27,6 +27,9 @@ ...@@ -27,6 +27,9 @@
#include <vector> #include <vector>
#include <list> #include <list>
#if __cplusplus > 199711L
#include <functional>
#endif
#include "AMDiS_fwd.h" #include "AMDiS_fwd.h"
#include "ProblemStatBase.h" #include "ProblemStatBase.h"
#include "Initfile.h" #include "Initfile.h"
...@@ -189,6 +192,14 @@ namespace AMDiS { ...@@ -189,6 +192,14 @@ namespace AMDiS {
virtual void addDirichletBC(BoundaryType type, int row, int col, virtual void addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *b); AbstractFunction<double, WorldVector<double> > *b);
#if __cplusplus > 199711L
/// Adds a Dirichlet boundary condition, where the rhs is given by an
/// lambda function or a std::function object
virtual void addDirichletBC(BoundaryType type, int row, int col,
std::function<double(WorldVector<double>)> b);
#endif
/// Adds a Dirichlet boundary condition, where the rhs is given by a DOF /// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
/// vector. /// vector.
virtual void addDirichletBC(BoundaryType type, int row, int col, virtual void addDirichletBC(BoundaryType type, int row, int col,
......
...@@ -34,11 +34,7 @@ ...@@ -34,11 +34,7 @@
#include "add_expr.hpp" // add two expressions #include "add_expr.hpp" // add two expressions
#include "mult_expr.hpp" // multiply two expressions #include "mult_expr.hpp" // multiply two expressions
#if HAS_VARIADIC_TEMPLATES && HAS_ALIAS_TEMPLATES #include "functor_expr.hpp" // apply a functor with 1/2/3 arguments to expressions
#include "functorN_expr.hpp" // apply a functor with arbitrary nr. of arguments to expressions
#else
#include "functor_expr.hpp" // apply a functor with 1/2/3 arguments to expressions
#endif
#include "cmath_expr.hpp" // apply a cmath function to expressions #include "cmath_expr.hpp" // apply a cmath function to expressions
#include "vec_functors.hpp" // apply a vector function to expressions #include "vec_functors.hpp" // apply a vector function to expressions
......
...@@ -28,13 +28,14 @@ ...@@ -28,13 +28,14 @@
#include "AMDiS_fwd.h" #include "AMDiS_fwd.h"
#include "LazyOperatorTerm.h" #include "LazyOperatorTerm.h"
#include "Functors.h" #include "Functors.h"
#include "expressions/functor_expr.hpp" // #include "expressions/functor_expr.hpp"
#include "operations/functors.hpp" #include "operations/functors.hpp"
#include <boost/static_assert.hpp> #include <boost/static_assert.hpp>
#include <tuple> #include <tuple>
#include <utility> #include <utility>
#include <functional>
namespace AMDiS namespace AMDiS
{ {
...@@ -71,7 +72,66 @@ namespace AMDiS ...@@ -71,7 +72,66 @@ namespace AMDiS
template<typename... Int> template<typename... Int>
static int eval(F f, Int... d) { return f.getDegree(d...); } static int eval(F f, Int... d) { return f.getDegree(d...); }
}; };
}