From ced7e4c94c8a6acb68daecf43095e3b5d0cdb642 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 14 Feb 2014 09:18:30 +0000
Subject: [PATCH] generic operators updated
---
AMDiS/src/GenericOperatorTerm.h | 43 +++--
AMDiS/src/GenericOperatorTerm.hh | 146 ++++++++++++++++
AMDiS/src/expressions/AtomicExpression.h | 147 ++++++++++++----
AMDiS/src/expressions/GradientOf.h | 128 ++++++++++++--
AMDiS/src/expressions/LazyExpression.h | 6 +
AMDiS/src/expressions/LazyOperatorTerm.h | 54 +++---
AMDiS/src/expressions/ValueOf.h | 204 +++++++++++++++--------
7 files changed, 568 insertions(+), 160 deletions(-)
create mode 100644 AMDiS/src/GenericOperatorTerm.hh
diff --git a/AMDiS/src/GenericOperatorTerm.h b/AMDiS/src/GenericOperatorTerm.h
index 198bcda0..e677e8b0 100644
--- a/AMDiS/src/GenericOperatorTerm.h
+++ b/AMDiS/src/GenericOperatorTerm.h
@@ -111,7 +111,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad, NULL);
}
@@ -120,7 +120,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm
SubAssembler* subAssembler,
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
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad);
}
@@ -169,7 +169,7 @@ struct GenericFirstOrderTerm_1 : public FirstOrderTerm
SubAssembler* subAssembler,
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
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad);
}
@@ -218,7 +218,7 @@ struct GenericFirstOrderTerm_i : public FirstOrderTerm
SubAssembler* subAssembler,
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
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad);
}
@@ -259,7 +259,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm
SubAssembler* subAssembler,
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
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad);
}
@@ -303,7 +303,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+ term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
/// Implements SecondOrderTerm::getLALt().
@@ -365,7 +365,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad);
}
@@ -374,7 +374,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+ term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
void getLALt(const ElInfo *elInfo,
@@ -453,7 +453,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, elInfo, subAssembler, quad);
+ term.initElement(this, elInfo, subAssembler, quad);
}
@@ -462,7 +462,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
SubAssembler* subAssembler,
Quadrature *quad)
{
- term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+ term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
void getLALt(const ElInfo *elInfo,
@@ -767,6 +767,21 @@ inline void addSOT(Operator& op, double term, int I, int 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
diff --git a/AMDiS/src/GenericOperatorTerm.hh b/AMDiS/src/GenericOperatorTerm.hh
new file mode 100644
index 00000000..4d33f0b9
--- /dev/null
+++ b/AMDiS/src/GenericOperatorTerm.hh
@@ -0,0 +1,146 @@
+/******************************************************************************
+ *
+ * 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
diff --git a/AMDiS/src/expressions/AtomicExpression.h b/AMDiS/src/expressions/AtomicExpression.h
index dbc7ad1d..93233b82 100644
--- a/AMDiS/src/expressions/AtomicExpression.h
+++ b/AMDiS/src/expressions/AtomicExpression.h
@@ -55,18 +55,53 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ 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>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ 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]; }
@@ -91,17 +126,51 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ 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>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ 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]; }
@@ -127,16 +196,18 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
elInfo->getNormal(side, normal);
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
smallElInfo->getNormal(side, normal);
}
@@ -164,22 +235,24 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
elInfo->getNormal(side, normal);
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
smallElInfo->getNormal(side, normal);
}
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 {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
elInfo->getElementNormal(elementNormal);
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
smallElInfo->getElementNormal(elementNormal);
}
@@ -237,22 +312,24 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
elInfo->getElementNormal(elementNormal);
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
smallElInfo->getElementNormal(elementNormal);
}
inline value_type operator()(const int& iq) const { return elementNormal[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; }
};
@@ -272,12 +349,14 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad) {}
+ 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) {}
+ 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; }
inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
diff --git a/AMDiS/src/expressions/GradientOf.h b/AMDiS/src/expressions/GradientOf.h
index d6aba4ff..29ca53b9 100644
--- a/AMDiS/src/expressions/GradientOf.h
+++ b/AMDiS/src/expressions/GradientOf.h
@@ -44,6 +44,7 @@ namespace result_of {
DOFVector<T>* vecDV;
mutable mtl::dense_vector<typename GradientType<T>::type> vec;
+ mutable mtl::dense_vector<T> coeff;
GradientOf(Vector& vector) : vecDV(&vector) {}
GradientOf(Vector* vector) : vecDV(vector) {}
@@ -60,18 +61,56 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
- {
- ot.getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
+ {
+ if (ot && subAssembler)
+ ot->getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ else if (quad)
+ vecDV->getGrdAtQPs(elInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ for (size_t i = 0; i < nBasisFct; i++)
+ localBasisFct->evalGrdUh(*basisFct->getCoords(i), grdLambda, coeff, vec[i]);
+ }
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ if (ot && subAssembler)
+ ot->getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ else if (quad)
+ vecDV->getGrdAtQPs(smallElInfo, largeElInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(smallElInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+
+ const DimVec<WorldVector<double> > &grdLambda = smallElInfo->getGrdLambda();
+ for (size_t i = 0; i < nBasisFct; i++)
+ localBasisFct->evalGrdUh(*basisFct->getCoords(i), grdLambda, coeff, vec[i]);
+ }
}
value_type operator()(const int& iq) const { return vec[iq]; }
@@ -85,7 +124,9 @@ namespace result_of {
typedef T value_type;
DOFVector<T>* vecDV;
- mutable mtl::dense_vector<typename GradientType<T>::type> vec;
+// mutable mtl::dense_vector<typename GradientType<T>::type> vec;
+ mutable mtl::dense_vector<T> vec;
+ mutable mtl::dense_vector<T> coeff;
int comp;
DerivativeOf(Vector& vector) : vecDV(&vector), comp(I) {}
@@ -114,21 +155,78 @@ namespace result_of {
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+// if (op && subAssembler)
+// ot->getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+// else
+// vecDV->getGrdAtQPs(elInfo, localQuad, NULL, vec);
+
+ if (ot && subAssembler)
+ subAssembler->getDerivativeAtQPs(vecDV, elInfo, quad, comp, vec);
+ else if (quad)
+ vecDV->getDerivativeAtQPs(elInfo, quad, NULL, comp, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ mtl::dense_vector<typename GradientType<T>::type> helper(nBasisFct);
+
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ for (size_t i = 0; i < nBasisFct; i++)
+ localBasisFct->evalGrdUh(*basisFct->getCoords(i), grdLambda, coeff, helper[i]);
+
+ vec.change_dim(nBasisFct);
+ for (size_t i = 0; i < num_rows(helper); i++)
+ vec[i] = helper[i][comp];
+ }
}
template<typename OT>
- void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+// if (op && subAssembler)
+// ot->getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+// else
+// vecDV->getGrdAtQPs(smallElInfo, largeElInfo, localQuad, NULL, vec);
+
+ if (ot && subAssembler)
+ subAssembler->getDerivativeAtQPs(vecDV, smallElInfo, largeElInfo, quad, comp, vec);
+ else if (quad)
+ vecDV->getDerivativeAtQPs(smallElInfo, largeElInfo, quad, NULL, comp, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(smallElInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ mtl::dense_vector<typename GradientType<T>::type> helper(nBasisFct);
+
+ const DimVec<WorldVector<double> > &grdLambda = smallElInfo->getGrdLambda();
+ for (size_t i = 0; i < nBasisFct; i++)
+ localBasisFct->evalGrdUh(*basisFct->getCoords(i), grdLambda, coeff, helper[i]);
+
+ vec.change_dim(nBasisFct);
+ for (size_t i = 0; i < num_rows(helper); i++)
+ vec[i] = helper[i][comp];
+ }
}
- value_type operator()(const int& iq) const { return vec[iq][comp]; }
+// value_type operator()(const int& iq) const { return vec[iq][comp]; }
+ value_type operator()(const int& iq) const { return vec[iq]; }
};
} // end namespace result_of
diff --git a/AMDiS/src/expressions/LazyExpression.h b/AMDiS/src/expressions/LazyExpression.h
index eb7a5b1b..f053bccb 100644
--- a/AMDiS/src/expressions/LazyExpression.h
+++ b/AMDiS/src/expressions/LazyExpression.h
@@ -40,6 +40,12 @@ struct NoType {};
struct FunctorBase
{
typedef NoType value_type;
+ int getDegree() const { return 0; }
+ int getDegree(int d0) const { return 0; }
+ int getDegree(int d0, int d1) const { return 0; }
+ int getDegree(int d0, int d1, int d2) const { return 0; }
+ int getDegree(int d0, int d1, int d2, int d3) const { return 0; }
+
protected:
virtual void helper() {};
};
diff --git a/AMDiS/src/expressions/LazyOperatorTerm.h b/AMDiS/src/expressions/LazyOperatorTerm.h
index 67b15eea..6dd75db8 100644
--- a/AMDiS/src/expressions/LazyOperatorTerm.h
+++ b/AMDiS/src/expressions/LazyOperatorTerm.h
@@ -49,17 +49,19 @@ struct LazyOperatorTerm1 : public LazyOperatorTermBase
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- term.initElement(ot, elInfo, subAssembler, quad);
+ term.initElement(ot, elInfo, subAssembler, quad, basisFct);
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- term.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ term.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad, basisFct);
}
inline double operator()(const int& iq) const;
@@ -80,19 +82,21 @@ struct LazyOperatorTerm2 : public LazyOperatorTermBase
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- term1.initElement(ot, elInfo, subAssembler, quad);
- term2.initElement(ot, elInfo, subAssembler, quad);
+ term1.initElement(ot, elInfo, subAssembler, quad, basisFct);
+ term2.initElement(ot, elInfo, subAssembler, quad, basisFct);
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- term1.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
- term2.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ term1.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad, basisFct);
+ term2.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad, basisFct);
}
double operator()(const int& iq) const;
@@ -116,21 +120,23 @@ struct LazyOperatorTerm3 : public LazyOperatorTermBase
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- term1.initElement(ot, elInfo, subAssembler, quad);
- term2.initElement(ot, elInfo, subAssembler, quad);
- term3.initElement(ot, elInfo, subAssembler, quad);
+ term1.initElement(ot, elInfo, subAssembler, quad, basisFct);
+ term2.initElement(ot, elInfo, subAssembler, quad, basisFct);
+ term3.initElement(ot, elInfo, subAssembler, quad, basisFct);
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- term1.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
- term2.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
- term3.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ term1.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad, basisFct);
+ term2.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad, basisFct);
+ term3.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad, basisFct);
}
inline double operator()(const int& iq) const;
diff --git a/AMDiS/src/expressions/ValueOf.h b/AMDiS/src/expressions/ValueOf.h
index a740d9b9..4af2ff3b 100644
--- a/AMDiS/src/expressions/ValueOf.h
+++ b/AMDiS/src/expressions/ValueOf.h
@@ -47,6 +47,7 @@ namespace result_of {
DOFVector<T>* vecDV;
mutable mtl::dense_vector<T> vec;
+ mutable mtl::dense_vector<T> coeff;
ValueOf(DOFVector<T>& vector) : vecDV(&vector) {}
ValueOf(DOFVector<T>* vector) : vecDV(vector) {}
@@ -63,18 +64,52 @@ namespace result_of {
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ if (ot && subAssembler)
+ ot->getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ else if (quad)
+ vecDV->getVecAtQPs(elInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+ for (size_t i = 0; i < nBasisFct; i++)
+ vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
+ }
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ if (ot && subAssembler)
+ ot->getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ else if (quad)
+ vecDV->getVecAtQPs(smallElInfo, largeElInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(smallElInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+ for (size_t i = 0; i < nBasisFct; i++)
+ vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
+ }
}
inline value_type operator()(const int& iq) const { return vec[iq]; }
@@ -89,8 +124,12 @@ namespace result_of {
Vector<DOFVector<T>*> vecDV;
mutable mtl::dense_vector<value_type> vec;
+ mutable Vector<mtl::dense_vector<T> > coeff;
- ValueOf(Vector<DOFVector<T>*>& vector) : vecDV(vector) {}
+ ValueOf(Vector<DOFVector<T>*>& vector) : vecDV(vector)
+ {
+ resize(coeff, num_rows(vecDV));
+ }
template<typename List>
inline void insertFeSpaces(List& feSpaces) const
@@ -105,12 +144,30 @@ namespace result_of {
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
Vector<mtl::dense_vector<T> > helper; resize(helper, num_rows(vecDV));
- for (size_t i = 0; i < num_rows(vecDV); i++)
- ot.getVectorAtQPs(vecDV[i], elInfo, subAssembler, quad, helper[i]);
+ for (size_t i = 0; i < num_rows(vecDV); i++) {
+ if (ot && subAssembler)
+ ot->getVectorAtQPs(vecDV[i], elInfo, subAssembler, quad, helper[i]);
+ else if (quad)
+ vecDV[i]->getVecAtQPs(elInfo, quad, NULL, helper[i]);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV[i]->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff[i].change_dim(localBasisFct->getNumber());
+ vecDV[i]->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ helper[i].change_dim(nBasisFct);
+ for (size_t j = 0; j < nBasisFct; j++)
+ helper[i][j] = localBasisFct->evalUh(*basisFct->getCoords(j), coeff[i]);
+ }
+ }
vec.change_dim(num_rows(helper[0]));
for (size_t iq = 0; iq < num_rows(helper[0]); iq++) {
value_type tmp; resize(tmp, num_rows(vecDV));
@@ -122,12 +179,30 @@ namespace result_of {
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
Vector<mtl::dense_vector<T> > helper; resize(helper, num_rows(vecDV));
- for (size_t i = 0; i < num_rows(vecDV); i++)
- ot.getVectorAtQPs(vecDV[i], smallElInfo, largeElInfo, subAssembler, quad, helper[i]);
+ for (size_t i = 0; i < num_rows(vecDV); i++) {
+ if (ot && subAssembler)
+ ot->getVectorAtQPs(vecDV[i], smallElInfo, subAssembler, quad, helper[i]);
+ else if (quad)
+ vecDV[i]->getVecAtQPs(smallElInfo, quad, NULL, helper[i]);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV[i]->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff[i].change_dim(localBasisFct->getNumber());
+ vecDV[i]->getLocalVector(smallElInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ helper[i].change_dim(nBasisFct);
+ for (size_t j = 0; j < nBasisFct; j++)
+ helper[i][j] = localBasisFct->evalUh(*basisFct->getCoords(j), coeff[i]);
+ }
+ }
vec.change_dim(num_rows(helper[0]));
for (size_t iq = 0; iq < num_rows(helper[0]); iq++) {
value_type tmp; resize(tmp, num_rows(vecDV));
@@ -142,58 +217,6 @@ namespace result_of {
};
- struct ValueOf2 : public LazyOperatorTermBase
- {
- typedef double value_type;
-
- DOFVector<value_type>* vecDV;
-
- mutable mtl::dense_vector<value_type> vec;
-
- const int I; // component of problem
-
- ValueOf2(ProblemStat& prob, int I_) : vecDV(prob.getSolution(I_)), I(I_) {}
- ValueOf2(ProblemStat* prob, int I_) : vecDV(prob->getSolution(I_)), I(I_) {}
-
- // for Rosenbrock problems
- ValueOf2(RosenbrockStationary& prob, int I_)
- : vecDV(prob.getStageSolution(I_)), I(I_) {}
- ValueOf2(RosenbrockStationary* prob, int I_)
- : vecDV(prob->getStageSolution(I_)), I(I_) {}
-
- template<typename List>
- inline void insertFeSpaces(List& feSpaces) const
- {
- feSpaces.insert(vecDV->getFeSpace());
- }
-
- inline int getDegree() const
- {
- return vecDV->getFeSpace()->getBasisFcts()->getDegree();
- }
-
- template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
- {
- ot.getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
- }
-
-
- template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
- {
- ot.getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
- }
-
- inline value_type operator()(const int& iq) const { return vec[iq]; }
- inline value_type derivative(const int& iq, int identifier) const
- {
- return static_cast<value_type>(identifier == I);
- }
- };
-
template<typename Vector>
struct ComponentOf : public LazyOperatorTermBase {};
@@ -204,6 +227,7 @@ namespace result_of {
DOFVector<Vector<T> >* vecDV;
mutable mtl::dense_vector<Vector<T> > vec;
+ mutable mtl::dense_vector<Vector<T> > coeff;
int I;
ComponentOf(DOFVector<Vector<T> >& vector, int I_) : vecDV(&vector), I(I_) {}
@@ -221,18 +245,52 @@ namespace result_of {
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* elInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ if (ot && subAssembler)
+ ot->getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ else if (quad)
+ vecDV->getVecAtQPs(elInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+ for (size_t i = 0; i < nBasisFct; i++)
+ vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
+ }
}
template<typename OT>
- inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
- SubAssembler* subAssembler, Quadrature *quad)
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
{
- ot.getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ if (ot && subAssembler)
+ ot->getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ else if (quad)
+ vecDV->getVecAtQPs(smallElInfo, largeElInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(smallElInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+ for (size_t i = 0; i < nBasisFct; i++)
+ vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
+ }
}
inline value_type operator()(const int& iq) const { return vec[iq][I]; }
--
GitLab