Skip to content

GitLab

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

new Reader for XYZ files and a lot of corrections in GenericOperatorTerms

parent 449a699f
Expand all Hide whitespace changes
Inline Side-by-side
AMDiS/src/GenericOperatorTerm.h
View file @ baa63e82
......@@ -36,6 +36,7 @@
#include "expressions/LazyOperatorTerm.h"
#include "expressions/AtomicExpression.h"
#include "expressions/LazyExpression.h"
#include "expressions/cmath_expr.h"
/** generic operator-terms provide an easy way of automated generation of
* 'arbitrary' operator-terms out of some elementary operations, by using a
......@@ -767,13 +768,28 @@ inline void addSOT(Operator& op, double term, int I, int J)
addSOT(op, constant(term), I, J);
}
// _____________________________________________________________________________
template<typename TOut>
inline result_of::Function1<result_of::Wrapper<TOut,WorldVector<double> >, result_of::Coords>
eval(AbstractFunction<TOut, WorldVector<double> >* fct) { return function_(wrap(fct), X()); }
// C++11:
// template<typename TOut>
// inline auto eval(AbstractFunction<TOut, WorldVector<double> >* fct) -> decltype( function_(wrap(fct), X()) )
// { return function_(wrap(fct), X()); }
template<typename Term>
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, typename Term::value_type>::type
integrate(Term term);
template<typename Term>
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, typename Term::value_type>::type
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, void>::type
transformDOF(Term term, DOFVector<typename Term::value_type>* result);
template<typename Term>
......
AMDiS/src/GenericOperatorTerm.hh
View file @ baa63e82
......@@ -53,13 +53,15 @@ namespace AMDiS {
template<typename Term>
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, typename Term::value_type>::type
integrate(Term term)
integrate(Term term, Mesh* mesh_ = NULL)
{
typedef typename Term::value_type TOut;
typename GenericOperatorTerm<Term>::type ot(term);
std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
Mesh* mesh = (*feSpaces.begin())->getMesh();
TEST_EXIT(mesh_ || !feSpaces.empty())("The Expression must contain a DOFVector or FeSpace depended value!\n");
Mesh* mesh = mesh_ ? mesh_ : (*feSpaces.begin())->getMesh();
int deg = term.getDegree();
int dim = mesh->getDim();
......@@ -91,23 +93,24 @@ integrate(Term term)
// works only for nodal basis functions!
template<typename Term>
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, typename Term::value_type>::type
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, void>::type
transformDOF(Term term, DOFVector<typename Term::value_type>* result)
{
typedef typename Term::value_type TOut;
TOut tmp; nullify(tmp);
DOFVector<TOut> temp(result->getFeSpace(), "temp");
DOFVector<short int> assigned(result->getFeSpace(), "assigned");
typename GenericOperatorTerm<Term>::type ot(term);
std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
Mesh* mesh = (*feSpaces.begin())->getMesh();
Mesh* mesh = result->getFeSpace()->getMesh();
const FiniteElemSpace* resultFeSpace = result->getFeSpace();
const FiniteElemSpace* resultFeSpace = temp.getFeSpace();
const BasisFunction *basisFcts = resultFeSpace->getBasisFcts();
int nBasisFcts = basisFcts->getNumber();
DOFVector<short int> assigned(resultFeSpace, "assigned");
assigned.set(0);
result->set(tmp);
temp.set(tmp);
std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
TraverseStack stack;
......@@ -120,16 +123,17 @@ transformDOF(Term term, DOFVector<typename Term::value_type>* result)
basisFcts->getLocalIndices(elInfo->getElement(), resultFeSpace->getAdmin(), localIndices);
for (int i = 0; i < nBasisFcts; i++) {
(*result)[localIndices[i]] += term(i);
temp[localIndices[i]] += term(i);
assigned[localIndices[i]]++;
}
elInfo = stack.traverseNext(elInfo);
}
}
DOFIterator<typename Term::value_type> resultIter(result, USED_DOFS);
DOFIterator<TOut> tempIter(&temp, USED_DOFS);
DOFIterator<TOut> resultIter(result, USED_DOFS);
DOFIterator<short int> assignedIter(&assigned, USED_DOFS);
for (resultIter.reset(), assignedIter.reset(); !resultIter.end(); ++resultIter, ++assignedIter) {
*resultIter *= 1.0/static_cast<double>(*assignedIter);
for (tempIter.reset(), resultIter.reset(), assignedIter.reset(); !resultIter.end(); ++tempIter, ++resultIter, ++assignedIter) {
*resultIter = (*tempIter)/static_cast<double>(*assignedIter);
}
}
......
AMDiS/src/expressions/AtomicExpression.h
View file @ baa63e82
......@@ -36,6 +36,34 @@
namespace AMDiS {
namespace traits {
template<typename T>
struct ValueType {
typedef T type;
static type eval(const T& t) { return t; }
};
template<typename T>
struct ValueType<T&> {
typedef typename ValueType<T>::type type;
static type eval(T& t) { return t; }
};
template<typename T>
struct ValueType<T*> {
typedef typename ValueType<T>::type type;
static type eval(T* t) { return ValueType<T>::eval(*t); }
};
template<typename T>
struct ValueType<const T> {
typedef typename ValueType<T>::type type;
static type eval(const T& t) { return ValueType<T>::eval(t); }
};
}
namespace result_of {
/// initialize a coordinate vector at quadrature points
......@@ -75,7 +103,6 @@ namespace result_of {
x.change_dim(nBasisFct);
for (int i = 0; i < nBasisFct; i++)
elInfo->coordToWorld(*basisFct->getCoords(i), x[i]);
MSG("elInfo->coordToWorld(basisFct)\n");
}
}
......@@ -105,7 +132,9 @@ namespace result_of {
}
inline value_type operator()(const int& iq) const { return x[iq]; }
inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
};
......@@ -174,7 +203,9 @@ namespace result_of {
}
inline double operator()(const int& iq) const { return x[iq][I]; }
inline double derivative(const int& iq, int identifier) const { return 0.0; }
template<typename Identifier>
inline double derivative(const int& iq) const { return 0.0; }
};
......@@ -213,7 +244,9 @@ namespace result_of {
}
inline value_type operator()(const int& iq) const { return normal; }
inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
};
......@@ -252,7 +285,9 @@ namespace result_of {
}
inline value_type operator()(const int& iq) const { return normal[I]; }
inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { return 0.0; }
};
......@@ -291,7 +326,9 @@ namespace result_of {
}
inline value_type operator()(const int& iq) const { return elementNormal; }
inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
};
......@@ -329,7 +366,9 @@ namespace result_of {
}
inline value_type operator()(const int& iq) const { return elementNormal[I]; }
inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { return 0.0; }
};
......@@ -359,22 +398,54 @@ namespace result_of {
const BasisFunction *basisFct = NULL) {}
inline value_type operator()(const int& iq) const { return value; }
inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { return 0.0; }
};
template<typename T>
struct Reference : public LazyOperatorTermBase
{
typedef typename ValueType<T>::type value_type;
const value_type& value;
Reference(const T& value_) : value(value_) {}
Reference(const T* value_) : value(*value_) {}
template<typename List>
void insertFeSpaces(List& feSpaces) const {}
int getDegree() const
{
return 0;
}
template<typename OT>
void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL) {}
template<typename OT>
void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL) {}
inline value_type operator()(const int& iq) const { return value; }
template<typename Identifier>
inline value_type derivative(const int& iq) const { return 0.0; }
};
///
template<typename Term>
template<typename Term, typename Identifier>
struct Jacobian : public LazyOperatorTerm1<Term>
{
typedef LazyOperatorTerm1<Term> super;
typedef typename Term::value_type value_type;
const int I;
Jacobian(const Term& term_, int I_) : super(term_), I(I_) {}
Jacobian(const Term& term_) : super(term_) {}
int getDegree() const
{
......@@ -382,7 +453,7 @@ namespace result_of {
return super::term.getDegree();
}
inline value_type operator()(const int& iq) const { return super::term.derivative(iq, I); }
inline value_type operator()(const int& iq) const { return super::term.template derivative<Identifier>(iq); }
};
}
......@@ -400,13 +471,18 @@ inline result_of::ElementNormal M(int I) { return result_of::ElementNormal(I); }
template<typename T>
inline result_of::Const<T> constant(const T& value) { return result_of::Const<T>(value); }
template<typename T>
inline result_of::Reference<T> ref_(T& value) { return result_of::Reference<T>(value); }
template<typename T>
inline result_of::Reference<T> ref_(T* value) { return result_of::Reference<T>(value); }
namespace Private {
template<typename Term>
template<typename Identifier, typename Term>
inline typename boost::enable_if<
typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
result_of::Jacobian<Term> >::type
jacobian(const Term& t, int I) { return result_of::Jacobian<Term>(t, I); }
result_of::Jacobian<Term, Identifier> >::type
jacobian(const Term& t) { return result_of::Jacobian<Term, Identifier>(t); }
}
......
AMDiS/src/expressions/LazyExpression.h
View file @ baa63e82
......@@ -52,118 +52,6 @@ protected:
namespace result_of {
///
template<int I, typename Term>
struct Pow : public LazyOperatorTerm1<Term>
{
typedef LazyOperatorTerm1<Term> super;
typedef typename AMDiS::detail::Pow<I,typename Term::value_type>::result_type value_type;
Pow(const Term& term_) : super(term_) {}
int getDegree() const
{
return I * (super::term.getDegree());
}
inline value_type operator()(const int& iq) const { return AMDiS::detail::Pow<I,typename Term::value_type>::eval(super::term(iq)); }
inline value_type derivative(const int& iq, int identifier) const
{
return super::term.derivative(iq, identifier) * I*AMDiS::detail::Pow<I-1,typename Term::value_type>::eval(super::term(iq));
}
};
///
template<typename Term>
struct Sqrt : public LazyOperatorTerm1<Term>
{
typedef LazyOperatorTerm1<Term> super;
typedef typename Term::value_type value_type;
Sqrt(const Term& term_) : super(term_) {}
int getDegree() const
{
return 2 * (super::term.getDegree()); // stimmt nicht ganz
}
inline value_type operator()(const int& iq) const { return sqrt(super::term(iq)); }
inline value_type derivative(const int& iq, int identifier) const
{
return super::term.derivative(iq, identifier) / (2.0 * sqrt(super::term(iq)));
}
};
///
template<typename Term>
struct Exp : public LazyOperatorTerm1<Term>
{
typedef LazyOperatorTerm1<Term> super;
typedef typename Term::value_type value_type;
int degree;
Exp(const Term& term_, int degree_ = 1) : super(term_), degree(degree_) {}
int getDegree() const
{
return degree * (super::term.getDegree());
}
inline value_type operator()(const int& iq) const { return exp(super::term(iq)); }
inline value_type derivative(const int& iq, int identifier) const
{
return super::term.derivative(iq, identifier) * exp(super::term(iq));
}
};
///
template<typename Term>
struct Abs : public LazyOperatorTerm1<Term>
{
typedef LazyOperatorTerm1<Term> super;
typedef typename Term::value_type value_type;
Abs(const Term& term_) : super(term_) {}
int getDegree() const
{
return super::term.getDegree();
}
inline value_type operator()(const int& iq) const { return std::abs(super::term(iq)); }
inline value_type derivative(const int& iq, int identifier) const
{
return super::term(iq) > 0.0 ? super::term.derivative(iq, identifier) : -super::term.derivative(iq, identifier);
}
};
///
template<typename Term>
struct Signum : public LazyOperatorTerm1<Term>
{
typedef LazyOperatorTerm1<Term> super;
typedef typename Term::value_type value_type;
Signum(const Term& term_) : super(term_) {}
int getDegree() const
{
return 0;
}
inline value_type operator()(const int& iq) const { return (super::term(iq) > 0.0 ? 1.0 : -1.0); }
inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
};
///
template<typename Term1, typename Term2>
......@@ -182,9 +70,24 @@ namespace result_of {
inline value_type operator()(const int& iq) const { return super::term1(iq) + super::term2(iq); }
inline value_type derivative(const int& iq, int identifier) const
template<typename Identifier>
inline value_type derivative(const int& iq) const
{
return derivative<Identifier>(iq, typename traits::DependsOn<Identifier, super>::type());
}
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::true_) const
{
return super::term1.derivative(iq, identifier) + super::term2.derivative(iq, identifier);
return super::term1.derivative<Identifier>(iq) + super::term2.derivative<Identifier>(iq);
}
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::false_) const
{
value_type result;
nullify(result);
return result;
}
};
......@@ -207,9 +110,24 @@ namespace result_of {
inline value_type operator()(const int& iq) const { return super::term1(iq) * super::term2(iq); }
inline value_type derivative(const int& iq, int identifier) const
template<typename Identifier>
inline value_type derivative(const int& iq) const
{
return derivative<Identifier>(iq, typename traits::DependsOn<Identifier, super>::type());
}
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::true_) const
{
return super::term1.derivative(iq, identifier) * super::term2(iq) + super::term1(iq) * super::term2.derivative(iq, identifier);
return super::term1.derivative<Identifier>(iq) * super::term2(iq) + super::term1(iq) * super::term2.derivative<Identifier>(iq);
}
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::false_) const
{
value_type result;
nullify(result);
return result;
}
};
......@@ -232,64 +150,25 @@ namespace result_of {
inline value_type operator()(const int& iq) const { return super::term1(iq) / super::term2(iq); }
inline value_type derivative(const int& iq, int identifier) const
template<typename Identifier>
inline value_type derivative(const int& iq) const
{
typename Term2::value_type term2_eval = super::term2(iq);
return (super::term1.derivative(iq, identifier) * term2_eval - super::term1(iq) * super::term2.derivative(iq, identifier)) / (sqr(term2_eval));
}
};
///
template<typename Term1, typename Term2>
struct Max : public LazyOperatorTerm2<Term1, Term2>
{
typedef LazyOperatorTerm2<Term1, Term2> super;
typedef typename Term1::value_type value_type;
Max(const Term1& term1_, const Term2& term2_)
: super(term1_, term2_) {}
int getDegree() const
{
return std::max(super::term1.getDegree(), super::term2.getDegree());
}
inline value_type operator()(const int& iq) const
{
return std::max(super::term1(iq), super::term2(iq));
return derivative<Identifier>(iq, typename traits::DependsOn<Identifier, super>::type());
}
inline value_type derivative(const int& iq, int identifier) const
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::true_) const
{
return super::term1(iq) > super::term2(iq) ? super::term1.derivative(iq, identifier) : super::term2.derivative(iq, identifier);
}
};
///
template<typename Term1, typename Term2>
struct Min : public LazyOperatorTerm2<Term1, Term2>
{
typedef LazyOperatorTerm2<Term1, Term2> super;
typedef typename Term1::value_type value_type;
Min(const Term1& term1_, const Term2& term2_)
: super(term1_, term2_) {}
int getDegree() const
{
return std::max(super::term1.getDegree(), super::term2.getDegree());
}
inline value_type operator()(const int& iq) const
{
return std::min(super::term1(iq), super::term2(iq));
typename Term2::value_type term2_eval = super::term2(iq);
return (super::term1.derivative<Identifier>(iq) * term2_eval - super::term1(iq) * super::term2.derivative<Identifier>(iq)) / (sqr(term2_eval));
}
inline value_type derivative(const int& iq, int identifier) const
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::false_) const
{
return super::term1(iq) < super::term2(iq) ? super::term1.derivative(iq, identifier) : super::term2.derivative(iq, identifier);
value_type result;
nullify(result);
return result;
}
};
......@@ -314,7 +193,25 @@ namespace result_of {
inline value_type operator()(const int& iq) const { return f(super::term(iq)); }
inline value_type derivative(const int& iq, int identifier) const { return super::term.derivative(iq, identifier) * f.derivative<0>(super::term(iq)); }
template<typename Identifier>
inline value_type derivative(const int& iq) const
{
return derivative<Identifier>(iq, typename traits::DependsOn<Identifier, super>::type());
}
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::true_) const
{
return super::term.derivative<Identifier>(iq) * f.derivative<0>(super::term(iq));
}
template<typename Identifier>
inline value_type derivative(const int& iq, boost::mpl::false_) const
{
value_type result;
nullify(result);
return result;
}
};
......@@ -338,11 +235,26 @@ namespace result_of {
}
inline value_type operator()(const int& iq) const { return f(super::term1(iq), super::term2(iq)); }
template<typename Identifier>
inline value_type derivative(const int& iq) const
{
return derivative<Identifier>(iq, typename traits::DependsOn<Identifier, super>::type());
}