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

generic operators updated

parent 5b7e5ed5
...@@ -111,7 +111,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm ...@@ -111,7 +111,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad, NULL);
} }
...@@ -120,7 +120,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm ...@@ -120,7 +120,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad, NULL);
} }
...@@ -160,7 +160,7 @@ struct GenericFirstOrderTerm_1 : public FirstOrderTerm ...@@ -160,7 +160,7 @@ struct GenericFirstOrderTerm_1 : public FirstOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad);
} }
...@@ -169,7 +169,7 @@ struct GenericFirstOrderTerm_1 : public FirstOrderTerm ...@@ -169,7 +169,7 @@ struct GenericFirstOrderTerm_1 : public FirstOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
} }
...@@ -209,7 +209,7 @@ struct GenericFirstOrderTerm_i : public FirstOrderTerm ...@@ -209,7 +209,7 @@ struct GenericFirstOrderTerm_i : public FirstOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad);
} }
...@@ -218,7 +218,7 @@ struct GenericFirstOrderTerm_i : public FirstOrderTerm ...@@ -218,7 +218,7 @@ struct GenericFirstOrderTerm_i : public FirstOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
} }
...@@ -250,7 +250,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm ...@@ -250,7 +250,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad);
} }
...@@ -259,7 +259,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm ...@@ -259,7 +259,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
} }
...@@ -294,7 +294,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm ...@@ -294,7 +294,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad);
} }
...@@ -303,7 +303,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm ...@@ -303,7 +303,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
} }
/// Implements SecondOrderTerm::getLALt(). /// Implements SecondOrderTerm::getLALt().
...@@ -365,7 +365,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm ...@@ -365,7 +365,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad);
} }
...@@ -374,7 +374,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm ...@@ -374,7 +374,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
} }
void getLALt(const ElInfo *elInfo, void getLALt(const ElInfo *elInfo,
...@@ -453,7 +453,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm ...@@ -453,7 +453,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, elInfo, subAssembler, quad); term.initElement(this, elInfo, subAssembler, quad);
} }
...@@ -462,7 +462,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm ...@@ -462,7 +462,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
{ {
term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad); term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
} }
void getLALt(const ElInfo *elInfo, void getLALt(const ElInfo *elInfo,
...@@ -767,6 +767,21 @@ inline void addSOT(Operator& op, double term, int I, int J) ...@@ -767,6 +767,21 @@ inline void addSOT(Operator& op, double term, int I, int J)
addSOT(op, constant(term), I, J); addSOT(op, constant(term), I, J);
} }
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
transformDOF(Term term, DOFVector<typename Term::value_type>* result);
template<typename Term>
typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
DOFVector<typename Term::value_type>& >::type
operator<<(DOFVector<typename Term::value_type>& result, const Term& term);
} }
#include "GenericOperatorTerm.hh"
#endif // AMDIS_GENERIC_OPERATOR_TERM_H #endif // AMDIS_GENERIC_OPERATOR_TERM_H
/******************************************************************************
*
* AMDiS - Adaptive multidimensional simulations
*
* Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
* Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
*
* Authors:
* Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
*
* This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
* WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
*
* This file is part of AMDiS
*
* See also license.opensource.txt in the distribution.
*
******************************************************************************/
/** \file GenericOperatorTerm.hh */
namespace AMDiS {
namespace detail {
template<typename Term, typename T>
struct GenericOperatorTerm { };
template<typename Term>
struct GenericOperatorTerm<Term, double> {
typedef GenericZeroOrderTerm<Term> type;
};
template<typename Term>
struct GenericOperatorTerm<Term, WorldVector<double> > {
typedef GenericFirstOrderTerm_b<Term> type;
};
template<typename Term >
struct GenericOperatorTerm<Term, WorldMatrix<double> > {
typedef GenericSecondOrderTerm_A<Term, true> type;
};
}
template<typename Term>
struct GenericOperatorTerm {
typedef typename detail::GenericOperatorTerm<Term, typename Term::value_type>::type type;
};
template<typename Term>
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, typename Term::value_type>::type
integrate(Term term)
{
typedef typename Term::value_type TOut;
typename GenericOperatorTerm<Term>::type ot(term);
std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
Mesh* mesh = (*feSpaces.begin())->getMesh();
int deg = term.getDegree();
int dim = mesh->getDim();
Quadrature* quad = Quadrature::provideQuadrature(dim, deg);
TOut value; nullify(value);
Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_DET;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo) {
term.initElement(&ot, elInfo, NULL, quad);
TOut tmp; nullify(tmp);
for (int iq = 0; iq < quad->getNumPoints(); iq++) {
tmp += quad->getWeight(iq) * term(iq);
}
value += tmp * elInfo->getDet();
elInfo = stack.traverseNext(elInfo);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::mpi::globalAdd(value);
#endif
return value;
}
// 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
transformDOF(Term term, DOFVector<typename Term::value_type>* result)
{
typedef typename Term::value_type TOut;
TOut tmp; nullify(tmp);
typename GenericOperatorTerm<Term>::type ot(term);
std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
Mesh* mesh = (*feSpaces.begin())->getMesh();
const FiniteElemSpace* resultFeSpace = result->getFeSpace();
const BasisFunction *basisFcts = resultFeSpace->getBasisFcts();
int nBasisFcts = basisFcts->getNumber();
DOFVector<short int> assigned(resultFeSpace, "assigned");
assigned.set(0);
result->set(tmp);
std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS | Mesh::FILL_GRD_LAMBDA);
while (elInfo) {
term.initElement(&ot, elInfo, NULL, NULL, basisFcts);
basisFcts->getLocalIndices(elInfo->getElement(), resultFeSpace->getAdmin(), localIndices);
for (int i = 0; i < nBasisFcts; i++) {
(*result)[localIndices[i]] += term(i);
assigned[localIndices[i]]++;
}
elInfo = stack.traverseNext(elInfo);
}
DOFIterator<typename Term::value_type> 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);
}
}
template<typename Term>
typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
DOFVector<typename Term::value_type>& >::type
operator<<(DOFVector<typename Term::value_type>& result, const Term& term)
{
transformDOF(term, &result);
return result;
}
}
\ No newline at end of file
...@@ -55,18 +55,53 @@ namespace result_of { ...@@ -55,18 +55,53 @@ namespace result_of {
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* elInfo, void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
subAssembler->getCoordsAtQPs(elInfo, quad, x); if (subAssembler) {
subAssembler->getCoordsAtQPs(elInfo, quad, x);
}
else if (quad) {
const int nPoints = quad->getNumPoints();
x.change_dim(nPoints);
for (int i = 0; i < nPoints; i++)
elInfo->coordToWorld(quad->getLambda(i), x[i]);
}
else if (basisFct) {
const int nBasisFct = basisFct->getNumber();
x.change_dim(nBasisFct);
for (int i = 0; i < nBasisFct; i++)
elInfo->coordToWorld(*basisFct->getCoords(i), x[i]);
MSG("elInfo->coordToWorld(basisFct)\n");
}
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo, void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
subAssembler->getCoordsAtQPs(smallElInfo, quad, x); if (subAssembler) {
subAssembler->getCoordsAtQPs(smallElInfo, quad, x);
}
else if (quad) {
const int nPoints = quad->getNumPoints();
x.change_dim(nPoints);
for (int i = 0; i < nPoints; i++)
smallElInfo->coordToWorld(quad->getLambda(i), x[i]);
}
else if (basisFct) {
const int nBasisFct = basisFct->getNumber();
x.change_dim(nBasisFct);
for (int i = 0; i < nBasisFct; i++)
smallElInfo->coordToWorld(*basisFct->getCoords(i), x[i]);
}
} }
inline value_type operator()(const int& iq) const { return x[iq]; } inline value_type operator()(const int& iq) const { return x[iq]; }
...@@ -91,17 +126,51 @@ namespace result_of { ...@@ -91,17 +126,51 @@ namespace result_of {
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* elInfo, void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
subAssembler->getCoordsAtQPs(elInfo, quad, x); if (subAssembler) {
subAssembler->getCoordsAtQPs(elInfo, quad, x);
}
else if (quad) {
const int nPoints = quad->getNumPoints();
x.change_dim(nPoints);
for (int i = 0; i < nPoints; i++)
elInfo->coordToWorld(quad->getLambda(i), x[i]);
}
else if (basisFct) {
const int nBasisFct = basisFct->getNumber();
x.change_dim(nBasisFct);
for (int i = 0; i < nBasisFct; i++)
elInfo->coordToWorld(*basisFct->getCoords(i), x[i]);
}
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo, void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
subAssembler->getCoordsAtQPs(smallElInfo, quad, x); if (subAssembler) {
subAssembler->getCoordsAtQPs(smallElInfo, quad, x);
}
else if (quad) {
const int nPoints = quad->getNumPoints();
x.change_dim(nPoints);
for (int i = 0; i < nPoints; i++)
smallElInfo->coordToWorld(quad->getLambda(i), x[i]);
}
else if (basisFct) {
const int nBasisFct = basisFct->getNumber();
x.change_dim(nBasisFct);
for (int i = 0; i < nBasisFct; i++)
smallElInfo->coordToWorld(*basisFct->getCoords(i), x[i]);
}
} }
inline double operator()(const int& iq) const { return x[iq][I]; } inline double operator()(const int& iq) const { return x[iq][I]; }
...@@ -127,16 +196,18 @@ namespace result_of { ...@@ -127,16 +196,18 @@ namespace result_of {
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* elInfo, void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
elInfo->getNormal(side, normal); elInfo->getNormal(side, normal);
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo, void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
smallElInfo->getNormal(side, normal); smallElInfo->getNormal(side, normal);
} }
...@@ -164,22 +235,24 @@ namespace result_of { ...@@ -164,22 +235,24 @@ namespace result_of {
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* elInfo, void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
elInfo->getNormal(side, normal); elInfo->getNormal(side, normal);
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo, void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
smallElInfo->getNormal(side, normal); smallElInfo->getNormal(side, normal);
} }
inline value_type operator()(const int& iq) const { return normal[I]; } inline value_type operator()(const int& iq) const { return normal[I]; }
inline value_type derivative(const int& iq, int identifier) const { 0.0; } inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
}; };
...@@ -201,16 +274,18 @@ namespace result_of { ...@@ -201,16 +274,18 @@ namespace result_of {
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* elInfo, void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
elInfo->getElementNormal(elementNormal); elInfo->getElementNormal(elementNormal);
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo, void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {
smallElInfo->getElementNormal(elementNormal); smallElInfo->getElementNormal(elementNormal);
} }
...@@ -237,22 +312,24 @@ namespace result_of { ...@@ -237,22 +312,24 @@ namespace result_of {
} }
template<typename OT> template<typename OT>
void initElement(OT& ot, const ElInfo* elInfo, void initElement(OT* ot, const ElInfo* elInfo,
SubAssembler* subAssembler, Quadrature *quad) SubAssembler* subAssembler, Quadrature *quad,
const BasisFunction *basisFct = NULL)
{ {