Commit 4a9832fa authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Implementation of dirichlet BC with expressions, some errors fixed and the ellipt-demo modified

parent b4adaf9e
......@@ -113,34 +113,4 @@ namespace AMDiS {
}
}
void DirichletExpressionBC::fillBoundaryCondition(
DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts)
{
bool isBoundaryElement = false;
for (int i = 0; i < nBasFcts; i++) {
if (localBound[i] == boundaryType) {
isBoundaryElement = true;
break;
}
}
if (!isBoundaryElement)
return;
WorldVector<double> worldCoords;
const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
term.initElement(&ot, elInfo, NULL, NULL, basFcts);
for (int i = 0; i < nBasFcts; i++) {
if (localBound[i] == boundaryType) {
typename Term::value_type value = term(i);
vector->setDirichletDofValue(dofIndices[i], value);
(*vector)[dofIndices[i]] = value;
}
}
}
}
} // end namespace AMDiS
......@@ -29,15 +29,14 @@
#include "BoundaryCondition.h"
#include "AbstractFunction.h"
#include "FixVec.h"
#include "Expressions.h"
namespace AMDiS
{
// some tags used for tag-dispatching
struct _value_by_abstractfunction {};
struct _value_by_dofvector {};
struct _value_by_function {};
struct _value_by_expression {};
namespace detail
{
......@@ -94,132 +93,130 @@ namespace AMDiS
} // end namespace detail
/**
* \ingroup Assembler
*
* \brief
* Sub class of BoundaryCondition. Implements Dirichlet boundary conditions.
* A DOFVectors is set to a given value at a Dirichlet dof and in a DOFMatrix
* the row corresponding to a Dirichlet dof is replaced by a row containing
* only a 1.0 in the diagonal.
*/
template <class ValueTag>
class DirichletBC {};
// specialization for AbstractFunctions as value container
template <>
class DirichletBC<_value_by_abstractfunction> : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true)
: super(type, rowFeSpace, colFeSpace, apply),
container(*fct)
{ }
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
AbstractFunction<double, WorldVector<double> > &container;
};
// specialization for DOFVectors as value container
template <>
class DirichletBC<_value_by_dofvector> : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletBC(BoundaryType type,
DOFVectorBase<double> *vec,
bool apply = true);
/**
* \ingroup Assembler
*
* \brief
* Sub class of BoundaryCondition. Implements Dirichlet boundary conditions.
* A DOFVectors is set to a given value at a Dirichlet dof and in a DOFMatrix
* the row corresponding to a Dirichlet dof is replaced by a row containing
* only a 1.0 in the diagonal.
*/
template <class ValueTag>
class DirichletBC {};
// specialization for AbstractFunctions as value container
template <>
class DirichletBC<_value_by_abstractfunction> : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true)
: super(type, rowFeSpace, colFeSpace, apply),
container(*fct)
{ }
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
AbstractFunction<double, WorldVector<double> > &container;
};
// specialization for DOFVectors as value container
template <>
class DirichletBC<_value_by_dofvector> : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletBC(BoundaryType type,
DOFVectorBase<double> *vec,
bool apply = true);
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
DOFVectorBase<double> *container;
};
protected:
DOFVectorBase<double> *container;
};
#if HAS_CXX11
// specialization for std::function or lambdas as value container
template <>
class DirichletBC<_value_by_function> : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletBC(BoundaryType type,
std::function<double(WorldVector<double>)> fct,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true)
: super(type, rowFeSpace, colFeSpace, apply),
container(fct)
{ }
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
std::function<double(WorldVector<double>)> container;
};
// specialization for std::function or lambdas as value container
template <>
class DirichletBC<_value_by_function> : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletBC(BoundaryType type,
std::function<double(WorldVector<double>)> fct,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true)
: super(type, rowFeSpace, colFeSpace, apply),
container(fct)
{ }
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
std::function<double(WorldVector<double>)> container;
};
#endif
// specialization for Expressions as value container
template <class Term>
class DirichletExpressionBC : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletExpressionBC(BoundaryType type,
Term const& term,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true)
: super(type, rowFeSpace, colFeSpace, apply),
term(term),
ot(term)
{ }
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
Term term;
GenericOperatorTerm<Term> ot;
};
// specialization for Expressions as value container
template <class Term>
class DirichletExpressionBC : public detail::DirichletBC
{
typedef detail::DirichletBC super;
public:
/// Constructor.
DirichletExpressionBC(BoundaryType type,
Term const& term,
const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL,
bool apply = true)
: super(type, rowFeSpace, colFeSpace, apply),
term(term),
ot(term)
{ }
/// Implementation of BoundaryCondition::fillBoundaryCondition().
virtual void fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
const BoundaryType* localBound,
int nBasFcts) override;
protected:
Term term;
GenericOperatorTerm<Term> ot;
};
} // end namespace AMDiS
......
......@@ -22,6 +22,7 @@
#define AMDIS_FUNCTORS_H
#include "AbstractFunction.h"
#include "Traits.h"
#include <boost/math/special_functions/cbrt.hpp>
#include <boost/math/special_functions/pow.hpp>
......
......@@ -34,11 +34,12 @@
#include "Boundary.h"
#include "MatrixVector.h"
#include "StandardProblemIteration.h"
#include "io/ElementFileWriter.h"
#include "ComponentTraverseInfo.h"
#include "AbstractFunction.h"
#include "solver/SolverMatrix.h"
#include "SystemVector.h"
#include "expressions/expr_traits.hpp"
#include "io/ElementFileWriter.h"
#include "solver/SolverMatrix.h"
namespace AMDiS {
......@@ -205,9 +206,9 @@ namespace AMDiS {
#endif
template <class Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type>::type
void ProblemStatSeq::addDirichletBC(BoundaryType type, int row, int col,
DOFVector<double> *vec);
typename boost::enable_if<typename traits::is_expr<Term>::type>::type
addDirichletBC(BoundaryType type, int row, int col,
Term const& term);
/// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
/// vector.
......
......@@ -18,12 +18,14 @@
*
******************************************************************************/
#include "DirichletBC.h"
namespace AMDiS
{
template <class Term>
typename boost::enable_if<typename traits::is_expr<Term>::type>::type
void ProblemStatSeq::addDirichletBC(BoundaryType type, int row, int col,
Term const& term)
ProblemStatSeq::addDirichletBC(BoundaryType type, int row, int col,
Term const& term)
{
FUNCNAME("ProblemStat::addDirichletBC()");
......@@ -32,8 +34,8 @@ namespace AMDiS
boundaryConditionSet = true;
DirichletExpressionBC<Term> *dirichletApply = new DirichletExpressionBC<Term>(type, term, true);
DirichletExpressionBC<Term> *dirichletNotApply = new DirichletExpressionBC<Term>(type, term, false);
DirichletExpressionBC<Term> *dirichletApply = new DirichletExpressionBC<Term>(type, term, componentSpaces[row], componentSpaces[col], true);
DirichletExpressionBC<Term> *dirichletNotApply = new DirichletExpressionBC<Term>(type, term, componentSpaces[row], componentSpaces[col], false);
for (int i = 0; i < nComponents; i++)
if (systemMatrix && (*systemMatrix)[row][i]) {
......
......@@ -26,7 +26,9 @@
#define AMDIS_COORDS_EXPRESSION_HPP
#include "AMDiS_fwd.h"
#include "ElInfo.h"
#include "LazyOperatorTerm.h"
#include "Mesh.h"
namespace AMDiS
{
......
......@@ -27,6 +27,8 @@
#include "LazyOperatorTerm.h"
#include <boost/lexical_cast.hpp>
namespace AMDiS
{
namespace expressions
......
......@@ -190,10 +190,10 @@ namespace AMDiS {
newTimestep = adaptInfo->getTimestep();
} else {
if (firstTimestep || succRejection) {
fac = pow((timeTol / errorEst), 1.0 / order);
fac = std::pow((timeTol / errorEst), 1.0 / order);
} else {
fac = adaptInfo->getTimestep() / tauAcc *
pow((timeTol * estAcc / (errorEst * errorEst)), 1.0 / order);
std::pow((timeTol * estAcc / (errorEst * errorEst)), 1.0 / order);
}
fac = std::max(0.1, std::min(fac, 3.0));
newTimestep = fac * adaptInfo->getTimestep();
......@@ -210,13 +210,13 @@ namespace AMDiS {
succRejection = true;
double reducedP =
log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej);
std::log(errorEst / estRej) / std::log(adaptInfo->getTimestep() / tauRej);
if (reducedP < order && reducedP > 0.0)
order = reducedP;
}
newTimestep = pow((timeTol / errorEst), 1.0 / order) * adaptInfo->getTimestep();
newTimestep = std::pow((timeTol / errorEst), 1.0 / order) * adaptInfo->getTimestep();
newTimestep *= 0.95;
tauRej = adaptInfo->getTimestep();
......
#include "AMDiS.h"
#include "boost/date_time/posix_time/posix_time.hpp"
using namespace AMDiS;
using namespace boost::posix_time;
// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return exp(-10.0 * (x * x));
}
};
/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >
{
......@@ -60,18 +45,18 @@ int main(int argc, char* argv[])
// ===== create matrix operator =====
Operator matrixOperator(ellipt.getFeSpace());
matrixOperator.addTerm(new Simple_SOT);
matrixOperator.addSot(1.0);
ellipt.addMatrixOperator(matrixOperator, 0, 0);
int degree = ellipt.getFeSpace()->getBasisFcts()->getDegree();
// ===== create rhs operator =====
Operator rhsOperator(ellipt.getFeSpace());
rhsOperator.addTerm(new CoordsAtQP_ZOT(new F(degree)));
rhsOperator.addZot(eval(new F(degree)));
ellipt.addVectorOperator(rhsOperator, 0);
// ===== add boundary conditions =====
ellipt.addDirichletBC(1, 0, 0, new G);
ellipt.addDirichletBC(1, 0, 0, exp(-10.0 * (X() * X())) );
// ===== start adaption loop =====
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment