From 30debe8f18725a439c472a7858a4ab6392e71d0f Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 24 Feb 2012 23:13:34 +0000
Subject: [PATCH] several changes in integrate routines of DOFVectors
---
AMDiS/CMakeLists.txt | 4 +-
AMDiS/src/DOFVector.cc | 187 ------------------------
AMDiS/src/DOFVector.h | 93 ++++++++----
AMDiS/src/DOFVector.hh | 258 ++++++++++++++++++++++++++++-----
AMDiS/src/FixVec.h | 26 ++++
AMDiS/src/parallel/MpiHelper.h | 15 ++
6 files changed, 335 insertions(+), 248 deletions(-)
diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 9ebae66d..1f795cf0 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -8,8 +8,8 @@ SET(SOURCE_DIR ${AMDIS_SOURCE_DIR}/src)
#TODO: use the cmake build type
SET(MTL_INCLUDE_DIR ${LIB_DIR}/mtl4/ CACHE PATH "mtl4 directory")
-set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -g -Wall -DDEBUG=0")
-set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -DDEBUG=1 -O0")
+set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -g -Wall -Wno-unused-but-set-variable -DDEBUG=0")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wno-unused-but-set-variable -DDEBUG=1 -O0")
if(NOT CurrentRevision)
find_package(Subversion)
if(Subversion_FOUND)
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 9de37947..3eaf62de 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -714,193 +714,6 @@ namespace AMDiS {
}
- double integrate(const DOFVector<double> &vec1,
- const DOFVector<double> &vec2,
- BinaryAbstractFunction<double, double, double> *fct)
- {
- if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
- return intSingleMesh(vec1, vec2, fct);
- else
- return intDualMesh(vec1, vec2, fct);
- }
-
-
- double intSingleMesh(const DOFVector<double> &vec1,
- const DOFVector<double> &vec2,
- BinaryAbstractFunction<double, double, double> *fct)
- {
- FUNCNAME("intDualmesh()");
-
- TEST_EXIT(fct)("No function defined!\n");
-
- int deg = std::max(fct->getDegree(),
- 2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
- Quadrature* quad =
- Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
- FastQuadrature *fastQuad =
- FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
-
- mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
- mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
-
- double value = 0.0;
- Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
- while (elInfo) {
- vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
- vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);
-
- double tmp = 0.0;
- for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
- tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
- value += tmp * elInfo->getDet();
-
- elInfo = stack.traverseNext(elInfo);
- }
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- double localValue = value;
- MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
-#endif
-
- return value;
- }
-
-
- double intDualMesh(const DOFVector<double> &vec1,
- const DOFVector<double> &vec2,
- BinaryAbstractFunction<double, double, double> *fct)
- {
- FUNCNAME("intDualmesh()");
-
- TEST_EXIT(fct)("No function defined!\n");
-
- int deg = std::max(fct->getDegree(),
- 2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
- Quadrature* quad =
- Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
- FastQuadrature *fastQuad =
- FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
-
- mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
- mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
-
- double value = 0.0;
- Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
- DualTraverse dualTraverse;
- DualElInfo dualElInfo;
- bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
- vec2.getFeSpace()->getMesh(),
- -1, -1, traverseFlag, traverseFlag,
- dualElInfo);
- while (cont) {
- vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
- vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);
-
- double tmp = 0.0;
- for (int iq = 0; iq < quad->getNumPoints(); iq++)
- tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
- value += tmp * dualElInfo.smallElInfo->getDet();
-
- cont = dualTraverse.traverseNext(dualElInfo);
- }
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- double localValue = value;
- MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
-#endif
-
- return value;
- }
-
-
- double integrate(const DOFVector<double> &vec,
- AbstractFunction<double, WorldVector<double> > *fct)
- {
- FUNCNAME("integrate()");
-
- TEST_EXIT(fct)("No function defined!\n");
-
- const FiniteElemSpace *feSpace = vec.getFeSpace();
- Mesh *mesh = feSpace->getMesh();
-
- int deg = std::max(fct->getDegree(), 2 * feSpace->getBasisFcts()->getDegree());
- Quadrature* quad = Quadrature::provideQuadrature(mesh->getDim(), deg);
- FastQuadrature *fastQuad =
- FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);
-
- mtl::dense_vector<double> qp(fastQuad->getNumPoints());
- WorldVector<double> coords;
-
- double value = 0.0;
- Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
- while (elInfo) {
- vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
-
- double tmp = 0.0;
- for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
- elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
- tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
- }
-
- value += tmp * elInfo->getDet();
-
- elInfo = stack.traverseNext(elInfo);
- }
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- double localValue = value;
- MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
-#endif
-
- return value;
- }
-
-
- double integrate(const FiniteElemSpace* feSpace,
- AbstractFunction<double, WorldVector<double> > *fct)
- {
- FUNCNAME("integrate()");
-
- TEST_EXIT(fct)("No function defined!\n");
-
- int deg = std::max(fct->getDegree(), feSpace->getBasisFcts()->getDegree());
- Quadrature* quad =
- Quadrature::provideQuadrature(feSpace->getMesh()->getDim(), deg);
- FastQuadrature *fastQuad =
- FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);
-
- WorldVector<double> coords;
-
- double value = 0.0;
- Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, traverseFlag);
- while (elInfo) {
- double tmp = 0.0;
- for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
- elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
- tmp += fastQuad->getWeight(iq) * (*fct)(coords);
- }
-
- value += tmp * elInfo->getDet();
-
- elInfo = stack.traverseNext(elInfo);
- }
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- double localValue = value;
- MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
-#endif
-
- return value;
- }
-
-
-
double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
AbstractFunction<double, std::vector<double> > *fct)
{
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 4b29a429..56e219db 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -872,46 +872,89 @@ namespace AMDiS {
template<typename T>
std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
std::vector<DOFVector<double>*> *res);
+
+
+ /// Computes the integral: \f$ \int f(\vec{x})\,\text{d}\vec{x}\f$
+ template<typename TOut>
+ TOut integrate_Coords(const FiniteElemSpace* feSpace,
+ AbstractFunction<TOut, WorldVector<double> > *fct);
+
+ template<typename TOut>
+ TOut integrate(const FiniteElemSpace* feSpace,
+ AbstractFunction<TOut, WorldVector<double> > *fct)
+ {
+ return integrate_Coords(feSpace, fct);
+ }
+
+ /// Computes the integral: \f$ \int f(v(\vec{x}))\,\text{d}\vec{x}\f$
+ template<typename TOut, typename T>
+ TOut integrate_Vec(const DOFVector<T> &vec,
+ AbstractFunction<TOut, T> *fct);
+
+ template<typename TOut, typename T>
+ TOut integrate(const DOFVector<T> &vec,
+ AbstractFunction<TOut, T> *fct)
+ {
+ return integrate_Vec(vec, fct);
+ }
/** \brief
- * Computes the integral of a function that includes two different DOFVectors. This
- * function works also for the case that the DOFVectors are defined on two different
+ * Computes the integral of a function that includes two different DOFVectors:
+ * \f$ \int f(v(\vec{x}), w(\vec{x}))\,\text{d}\vec{x}\f$
+ * This function works also for the case that the DOFVectors are defined on two different
* meshes.
*/
- double integrate(const DOFVector<double> &vec1,
- const DOFVector<double> &vec2,
- BinaryAbstractFunction<double, double, double> *fct);
-
+ template<typename TOut, typename T1, typename T2>
+ TOut integrate_Vec2(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct);
+
+ template<typename TOut, typename T1, typename T2>
+ TOut integrate(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct)
+ {
+ return integrate_Vec2(vec1, vec2, fct);
+ }
+
/// Computes the integral of a function with two DOFVectors defined on the same mesh.
- double intSingleMesh(const DOFVector<double> &vec1,
- const DOFVector<double> &vec2,
- BinaryAbstractFunction<double, double, double> *fct);
+ template<typename TOut, typename T1, typename T2>
+ TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct);
/// Computes the integral of a function with two DOFVectors defined on different meshes.
- double intDualMesh(const DOFVector<double> &vec1,
- const DOFVector<double> &vec2,
- BinaryAbstractFunction<double, double, double> *fct);
-
- double integrate(const DOFVector<double> &vec,
- AbstractFunction<double, WorldVector<double> > *fct);
-
- double integrate(const FiniteElemSpace* feSpace,
- AbstractFunction<double, WorldVector<double> > *fct);
+ template<typename TOut, typename T1, typename T2>
+ TOut int_Vec2_DualMesh(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct);
+
+ /// Computes the integral: \f$ \int v(\vec{x}) f(\vec{x})\,\text{d}\vec{x}\f$
+ template<typename T1, typename T2>
+ typename ProductType<T1,T2>::type integrate_VecTimesCoords(const DOFVector<T1> &vec,
+ AbstractFunction<T2, WorldVector<double> > *fct);
+
+ /// Computes the integral: \f$ \int f(v(\vec{x}), \vec{x})\,\text{d}\vec{x}\f$
+ template<typename TOut, typename T>
+ TOut integrate_VecAndCoords(const DOFVector<T> &vec,
+ BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct);
+ template<typename TOut, typename T>
+ TOut integrate_VecAndCoords(const DOFVector<T> &vec,
+ BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
+ {
+ return integrate(vec, fct);
+ }
+
+ /// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
AbstractFunction<double, std::vector<double> > *fct);
+ /// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i,\{\nabla w_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
const std::vector<DOFVector<double>*> &grds,
BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct);
- template<typename T>
- T integrate_VecAndCoords(const DOFVector<double> &vec,
- BinaryAbstractFunction<T, double, WorldVector<double> > *fct);
-
- template<typename T1, typename T2>
- T2 integrate_Vec(const DOFVector<T1> &vec,
- AbstractFunction<T2, T1> *fct);
}
#include "DOFVector.hh"
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 9091a4dc..c81f285f 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -34,6 +34,11 @@
#include "Operator.h"
#include "Initfile.h"
#include "Traverse.h"
+#include "DualTraverse.h"
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#include "MpiHelper.h"
+#endif
// Defining the interface for MTL4
namespace mtl {
@@ -516,7 +521,7 @@ namespace AMDiS {
FastQuadrature *quadFast =
FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
- double result = 0.0;
+ T result; nullify(result);
int nPoints = quadFast->getNumPoints();
mtl::dense_vector<T> uh_vec(nPoints);
TraverseStack stack;
@@ -524,7 +529,7 @@ namespace AMDiS {
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
while (elInfo) {
double det = elInfo->getDet();
- double normT = 0.0;
+ T normT; nullify(normT);
this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
for (int iq = 0; iq < nPoints; iq++)
normT += quadFast->getWeight(iq) * (uh_vec[iq]);
@@ -534,18 +539,58 @@ namespace AMDiS {
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- double localResult = result;
- MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
+ mpi::globalAdd(result);
#endif
return result;
}
- template<typename T>
- T integrate_VecAndCoords(const DOFVector<double> &vec,
- BinaryAbstractFunction<T, double, WorldVector<double> > *fct)
+
+ template<typename TOut>
+ TOut integrate_Coords(const FiniteElemSpace* feSpace,
+ AbstractFunction<TOut, WorldVector<double> > *fct)
{
- FUNCNAME("integrate_VecAndCoords()");
+ FUNCNAME("integrate_Coords()");
+
+ TEST_EXIT(fct)("No function defined!\n");
+
+ int deg = std::max(fct->getDegree(), feSpace->getBasisFcts()->getDegree());
+ Quadrature* quad =
+ Quadrature::provideQuadrature(feSpace->getMesh()->getDim(), deg);
+ FastQuadrature *fastQuad =
+ FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);
+
+ WorldVector<double> coords;
+
+ TOut value; nullify(value);
+ Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, traverseFlag);
+ while (elInfo) {
+ TOut tmp; nullify(tmp);
+ for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
+ elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
+ tmp += fastQuad->getWeight(iq) * (*fct)(coords);
+ }
+
+ value += tmp * elInfo->getDet();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ mpi::globalAdd(value);
+#endif
+
+ return value;
+ }
+
+
+ template<typename TOut, typename T>
+ TOut integrate_Vec(const DOFVector<T> &vec,
+ AbstractFunction<TOut, T> *fct)
+ {
+ FUNCNAME("integrate_Vec()");
TEST_EXIT(fct)("No function defined!\n");
@@ -556,39 +601,184 @@ namespace AMDiS {
FastQuadrature *fastQuad =
FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
- mtl::dense_vector<double> qp(fastQuad->getNumPoints());
- WorldVector<double> coordsAtQP;
+ mtl::dense_vector<T> qp(fastQuad->getNumPoints());
+
+ TOut value; nullify(value);
- T value;
- nullify(value);
-
Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
while (elInfo) {
vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
- T tmp;
- nullify(tmp);
+ TOut tmp; nullify(tmp);
for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
- elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
- tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
+ tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
}
value += tmp * elInfo->getDet();
elInfo = stack.traverseNext(elInfo);
}
-// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-// double localValue = value;
-// MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
-// #endif
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ mpi::globalAdd(value);
+#endif
return value;
}
+
+
+ template<typename TOut, typename T1, typename T2>
+ TOut integrate_Vec2(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct)
+ {
+ if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
+ return int_Vec2_SingleMesh(vec1, vec2, fct);
+ else
+ return int_Vec2_DualMesh(vec1, vec2, fct);
+ }
+
+ template<typename TOut, typename T1, typename T2>
+ TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct)
+ {
+ FUNCNAME("intDualmesh()");
+
+ TEST_EXIT(fct)("No function defined!\n");
+
+ int deg = std::max(fct->getDegree(),
+ 2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
+ Quadrature* quad =
+ Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
+ FastQuadrature *fastQuad =
+ FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
+
+ mtl::dense_vector<T1> qp1(fastQuad->getNumPoints());
+ mtl::dense_vector<T2> qp2(fastQuad->getNumPoints());
+
+ TOut value; nullify(value);
+ Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
+ while (elInfo) {
+ vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
+ vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);
+
+ TOut tmp; nullify(tmp);
+ for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
+ tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
+ value += tmp * elInfo->getDet();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ mpi::globalAdd(value);
+#endif
+
+ return value;
+ }
+
+
+ template<typename TOut, typename T1, typename T2>
+ TOut int_Vec2_DualMesh(const DOFVector<T1> &vec1,
+ const DOFVector<T2> &vec2,
+ BinaryAbstractFunction<TOut, T1, T2> *fct)
+ {
+ FUNCNAME("intDualmesh()");
+
+ TEST_EXIT(fct)("No function defined!\n");
+
+ int deg = std::max(fct->getDegree(),
+ 2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
+ Quadrature* quad =
+ Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
+ FastQuadrature *fastQuad =
+ FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
+
+ mtl::dense_vector<T1> qp1(fastQuad->getNumPoints());
+ mtl::dense_vector<T2> qp2(fastQuad->getNumPoints());
+
+ TOut value; nullify(value);
+ Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
+ DualTraverse dualTraverse;
+ DualElInfo dualElInfo;
+ bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
+ vec2.getFeSpace()->getMesh(),
+ -1, -1, traverseFlag, traverseFlag,
+ dualElInfo);
+ while (cont) {
+ vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
+ vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);
+
+ TOut tmp; nullify(tmp);
+ for (int iq = 0; iq < quad->getNumPoints(); iq++)
+ tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
+ value += tmp * dualElInfo.smallElInfo->getDet();
+
+ cont = dualTraverse.traverseNext(dualElInfo);
+ }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ mpi::globalAdd(value);
+#endif
+
+ return value;
+ }
+
+
template<typename T1, typename T2>
- T2 integrate_Vec(const DOFVector<T1> &vec,
- AbstractFunction<T2, T1> *fct)
+ typename ProductType<T1,T2>::type integrate_VecTimesCoords(const DOFVector<T1> &vec,
+ AbstractFunction<T2, WorldVector<double> > *fct)
+ {
+ FUNCNAME("integrate_VecTimesCoords()");
+
+ TEST_EXIT(fct)("No function defined!\n");
+
+ typedef typename ProductType<T1,T2>::type TOut;
+
+ const FiniteElemSpace *feSpace = vec.getFeSpace();
+ Mesh *mesh = feSpace->getMesh();
+
+ int deg = std::max(fct->getDegree(), 2 * feSpace->getBasisFcts()->getDegree());
+ Quadrature* quad = Quadrature::provideQuadrature(mesh->getDim(), deg);
+ FastQuadrature *fastQuad =
+ FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);
+
+ mtl::dense_vector<T1> qp(fastQuad->getNumPoints());
+ WorldVector<double> coords;
+
+ TOut value; nullify(value);
+ Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
+ while (elInfo) {
+ vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
+
+ TOut tmp; nullify(tmp);
+ for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
+ elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
+ tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
+ }
+
+ value += tmp * elInfo->getDet();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ mpi::globalAdd(value);
+#endif
+
+ return value;
+ }
+
+
+ template<typename TOut, typename T>
+ TOut integrate_VecAndCoords(const DOFVector<T> &vec,
+ BinaryAbstractFunction<T, T, WorldVector<double> > *fct)
{
FUNCNAME("integrate_VecAndCoords()");
@@ -601,33 +791,33 @@ namespace AMDiS {
FastQuadrature *fastQuad =
FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
- mtl::dense_vector<T1> qp(fastQuad->getNumPoints());
-
- T2 value;
- nullify(value);
+ mtl::dense_vector<T> qp(fastQuad->getNumPoints());
+ WorldVector<double> coordsAtQP;
+ TOut value; nullify(value);
+
Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
while (elInfo) {
vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
- T2 tmp;
- nullify(tmp);
+ TOut tmp; nullify(tmp);
for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
- tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
+ elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
+ tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
}
value += tmp * elInfo->getDet();
elInfo = stack.traverseNext(elInfo);
}
-// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-// double localValue = value;
-// MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
-// #endif
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ mpi::globalAdd(value);
+#endif
return value;
}
+
template<typename T>
double DOFVector<T>::L1Norm(Quadrature* q) const
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index d8395de5..33c6e873 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -645,7 +645,33 @@ namespace AMDiS {
}
};
+ /// defines of wich type the product T1*T2 is
+ template<typename T1, typename T2>
+ struct ProductType
+ {
+ typedef T1 type;
+ };
+
+ // scalar times vector
+ template<typename T1>
+ struct ProductType<T1, WorldVector<T1> >
+ {
+ typedef WorldVector<typename ProductType<T1,T1>::type > type;
+ };
+ // vector times scalar
+ template<typename T1>
+ struct ProductType<WorldVector<T1>, T1>
+ {
+ typedef WorldVector<typename ProductType<T1,T1>::type > type;
+ };
+
+ // vector times vector
+ template<typename T1, typename T2>
+ struct ProductType<WorldVector<T1>, WorldVector<T2> >
+ {
+ typedef typename ProductType<T1,T2>::type type;
+ };
}
diff --git a/AMDiS/src/parallel/MpiHelper.h b/AMDiS/src/parallel/MpiHelper.h
index 9c5e88a1..81fb5ddd 100644
--- a/AMDiS/src/parallel/MpiHelper.h
+++ b/AMDiS/src/parallel/MpiHelper.h
@@ -35,15 +35,30 @@ namespace AMDiS {
void globalAdd(double &value);
void globalAdd(int &value);
+
+ template<typename T>
+ void globalAdd(T &value) {
+ WARNING("Unknown type for globalAdd. Can not sum up the values of all processors!\n");
+ }
void globalMin(double &value);
void globalMin(int &value);
+ template<typename T>
+ void globalMin(T &value) {
+ WARNING("Unknown type for globalMin. Can not determine minimal value of all processors!\n");
+ }
+
void globalMax(double &value);
void globalMax(int &value);
+ template<typename T>
+ void globalMax(T &value) {
+ WARNING("Unknown type for globalMax. Can not determine maximal value of all processors!\n");
+ }
+
inline void startRand()
{
srand(time(NULL) * (MPI::COMM_WORLD.Get_rank() + 1));
--
GitLab