From 2d8a6968f878f0fb90abbcfd83caf549b98073b5 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 26 May 2010 08:29:11 +0000
Subject: [PATCH] Added new functions for integrating abstract function with
two DOFVectors.
---
AMDiS/src/AMDiS.h | 3 +-
AMDiS/src/AMDiS_fwd.h | 13 ++++++
AMDiS/src/DOFVector.cc | 92 ++++++++++++++++++++++++++++++++++++++++++
AMDiS/src/DOFVector.h | 20 +++++++++
4 files changed, 127 insertions(+), 1 deletion(-)
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index a39c5a0a..2d9d4949 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -1,6 +1,6 @@
#ifndef AMDIS_H
#define AMDIS_H
-#include "MathFunctions.h"
+
#include "AbstractFunction.h"
#include "AdaptInfo.h"
#include "AdaptInstationary.h"
@@ -49,6 +49,7 @@
#include "MacroElement.h"
#include "MacroWriter.h"
#include "Marker.h"
+#include "MathFunctions.h"
#include "MatrixVector.h"
#include "Mesh.h"
#include "MeshStructure.h"
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index 272c8904..686225f7 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -101,6 +101,19 @@ namespace AMDiS {
struct AtomicBoundary;
template<typename ReturnType, typename ArgumentType> class AbstractFunction;
+ template<typename ReturnType,
+ typename ArgumentType1,
+ typename ArgumentType2> class BinaryAbstractFunction;
+ template<typename ReturnType,
+ typename ArgumentType1,
+ typename ArgumentType2,
+ typename ArgumentType3> class TertiaryAbstractFunction;
+ template<typename ReturnType,
+ typename ArgumentType1,
+ typename ArgumentType2,
+ typename ArgumentType3,
+ typename ArgumentType4> class QuartAbstractFunction;
+
template<typename T> class DOFIndexed;
template<typename T> class DOFVectorBase;
template<typename T> class DOFVector;
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index cb978648..ab89835f 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -865,5 +865,97 @@ 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 = 2 * vec1.getFESpace()->getBasisFcts()->getDegree();
+ Quadrature* quad =
+ Quadrature::provideQuadrature(vec1.getFESpace()->getMesh()->getDim(), deg);
+ FastQuadrature *fastQuad =
+ FastQuadrature::provideFastQuadrature(vec1.getFESpace()->getBasisFcts(), *quad, INIT_PHI);
+
+ std::vector<double> qp1(fastQuad->getNumPoints());
+ std::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) {
+ double det = elInfo->getDet();
+ double tmp = 0.0;
+ vec1.getVecAtQPs(elInfo, quad, fastQuad, &(qp1[0]));
+ vec2.getVecAtQPs(elInfo, quad, fastQuad, &(qp2[0]));
+ for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
+ tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
+
+ value += tmp * det;
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ 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 = 2 * vec1.getFESpace()->getBasisFcts()->getDegree();
+ Quadrature* quad =
+ Quadrature::provideQuadrature(vec1.getFESpace()->getMesh()->getDim(), deg);
+ FastQuadrature *fastQuad =
+ FastQuadrature::provideFastQuadrature(vec1.getFESpace()->getBasisFcts(), *quad, INIT_PHI);
+
+ std::vector<double> qp1(fastQuad->getNumPoints());
+ std::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) {
+ double det = dualElInfo.smallElInfo->getDet();
+ double tmp = 0.0;
+ vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo,
+ quad, fastQuad, &(qp1[0]));
+ vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo,
+ quad, fastQuad, &(qp2[0]));
+ for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
+ tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
+
+ value += tmp * det;
+
+ cont = dualTraverse.traverseNext(dualElInfo);
+ }
+
+ return value;
+ }
+
}
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 4ba37da1..d8716809 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -789,6 +789,26 @@ namespace AMDiS {
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *result);
+
+ /** \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
+ * meshes.
+ */
+ double integrate(const DOFVector<double> &vec1,
+ const DOFVector<double> &vec2,
+ BinaryAbstractFunction<double, double, double> *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);
+
+ /// 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);
+
}
#include "DOFVector.hh"
--
GitLab