Commit 2d8a6968 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Added new functions for integrating abstract function with two DOFVectors.

parent dc85616d
#ifndef AMDIS_H #ifndef AMDIS_H
#define AMDIS_H #define AMDIS_H
#include "MathFunctions.h"
#include "AbstractFunction.h" #include "AbstractFunction.h"
#include "AdaptInfo.h" #include "AdaptInfo.h"
#include "AdaptInstationary.h" #include "AdaptInstationary.h"
...@@ -49,6 +49,7 @@ ...@@ -49,6 +49,7 @@
#include "MacroElement.h" #include "MacroElement.h"
#include "MacroWriter.h" #include "MacroWriter.h"
#include "Marker.h" #include "Marker.h"
#include "MathFunctions.h"
#include "MatrixVector.h" #include "MatrixVector.h"
#include "Mesh.h" #include "Mesh.h"
#include "MeshStructure.h" #include "MeshStructure.h"
......
...@@ -101,6 +101,19 @@ namespace AMDiS { ...@@ -101,6 +101,19 @@ namespace AMDiS {
struct AtomicBoundary; struct AtomicBoundary;
template<typename ReturnType, typename ArgumentType> class AbstractFunction; 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 DOFIndexed;
template<typename T> class DOFVectorBase; template<typename T> class DOFVectorBase;
template<typename T> class DOFVector; template<typename T> class DOFVector;
......
...@@ -865,5 +865,97 @@ namespace AMDiS { ...@@ -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;
}
} }
...@@ -789,6 +789,26 @@ namespace AMDiS { ...@@ -789,6 +789,26 @@ namespace AMDiS {
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec, WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *result); 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" #include "DOFVector.hh"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment