From 2982e0f319a7df9b6030a150ee0b2eeacd6938bd Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 30 Mar 2010 16:18:19 +0000
Subject: [PATCH] Some small changes in Operators. Thanks to Simon ...
---
AMDiS/src/Operator.cc | 108 +++++++++++++-----------------------------
AMDiS/src/Operator.h | 67 +++++---------------------
2 files changed, 47 insertions(+), 128 deletions(-)
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 354b97ff..b125ef2e 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -767,8 +767,9 @@ namespace AMDiS {
((*b) * grdUhAtQP[iq]);
}
+
FctVecAtQP_FOT::FctVecAtQP_FOT(DOFVectorBase<double> *dv,
- AbstractFunction<double, WorldVector<double> > *f_,
+ BinaryAbstractFunction<double, WorldVector<double>, double> *f_,
WorldVector<double> *b_)
: FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_)
{
@@ -791,9 +792,9 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
if (b)
- lb(Lambda, *b, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
+ lb(Lambda, *b, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
else
- l1(Lambda, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
+ l1(Lambda, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
}
}
@@ -807,21 +808,29 @@ namespace AMDiS {
if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++)
result[iq] +=
- fac * (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
+ fac * (*f)(coordsAtQPs[iq], vecAtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
}
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af,
WorldVector<double> *b_)
- : FirstOrderTerm(0), vec1(dv1), vec2(dv2), b(b_)
+ : FirstOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af), b(b_)
{
+ TEST_EXIT(dv1)("No first vector!\n");
+ TEST_EXIT(dv2)("No second vector!\n");
+
auxFeSpaces.insert(dv1->getFESpace());
auxFeSpaces.insert(dv2->getFESpace());
}
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af,
int bIdx)
- : FirstOrderTerm(0), vec1(dv1), vec2(dv2)
+ : FirstOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
{
+ TEST_EXIT(dv1)("No first vector!\n");
+ TEST_EXIT(dv2)("No second vector!\n");
+
bOne = bIdx;
auxFeSpaces.insert(dv1->getFESpace());
auxFeSpaces.insert(dv2->getFESpace());
@@ -842,13 +851,13 @@ namespace AMDiS {
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(Lambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
+ lb_one(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
+ lb(Lambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
- l1(Lambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
+ l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
}
}
@@ -861,11 +870,14 @@ namespace AMDiS {
{
if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++)
- result[iq] += fac * vec1AtQPs[iq] * vec2AtQPs[iq] * ((*b) * grdUhAtQP[iq]);
+ result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
}
- Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
- BinaryAbstractFunction<double, double, double> *f_,
+
+ Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ DOFVectorBase<double> *dv3,
+ TertiaryAbstractFunction<double, double, double, double> *f_,
WorldVector<double> *bvec)
: FirstOrderTerm(f_->getDegree()),
vec1(dv1),
@@ -879,8 +891,11 @@ namespace AMDiS {
auxFeSpaces.insert(dv3->getFESpace());
}
- Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
- BinaryAbstractFunction<double, double, double> *f_,
+
+ Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ DOFVectorBase<double> *dv3,
+ TertiaryAbstractFunction<double, double, double, double> *f_,
int b)
: FirstOrderTerm(f_->getDegree()),
vec1(dv1),
@@ -910,17 +925,18 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(Lambda, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
+ lb_one(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++) {
if (b)
- lb(Lambda, *b, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
+ lb(Lambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
else
- l1(Lambda, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
+ l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
}
}
}
+
void Vec3FctAtQP_FOT::eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -930,8 +946,7 @@ namespace AMDiS {
{
if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++)
- result[iq] += fac *
- (vec1AtQPs[iq]) * (*f)(vec2AtQPs[iq] ,vec3AtQPs[iq]) *
+ result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq] ,vec3AtQPs[iq]) *
((*b) * grdUhAtQP[iq]);
}
@@ -1053,15 +1068,6 @@ namespace AMDiS {
auxFeSpaces.insert(dv->getFESpace());
}
- VecAndGradientAtQP_SOT::VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
- BinaryAbstractFunction<double, double, WorldVector<double> > *af)
- : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
- {
- TEST_EXIT(dv)("No vector!\n");
-
- auxFeSpaces.insert(dv->getFESpace());
- }
-
VecMatrixGradientAtQP_SOT::VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
@@ -1325,14 +1331,6 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void VecAndGradientAtQP_SOT::initElement(const ElInfo* elInfo,
- SubAssembler* subAssembler,
- Quadrature *quad)
- {
- vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
- gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
- }
-
void VecMatrixGradientAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -1581,14 +1579,6 @@ namespace AMDiS {
l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
}
- void VecAndGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
- {
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
- }
-
void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
@@ -2156,36 +2146,6 @@ namespace AMDiS {
}
}
- void VecAndGradientAtQP_SOT::eval(int nPoints,
- const double *uhAtQP,
- const WorldVector<double> *grdUhAtQP,
- const WorldMatrix<double> *D2UhAtQP,
- double *result,
- double fac)
- {
- int dow = Global::getGeo(WORLD);
-
- if (D2UhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
- double resultQP = 0.0;
- for (int i = 0; i < dow; i++)
- resultQP += D2UhAtQP[iq][i][i];
- result[iq] += fac * factor * resultQP;
- }
- }
- }
-
- void VecAndGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result)
- {
- int nPoints = grdUhAtQP.size();
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
- }
-
void VecMatrixGradientAtQP_SOT::eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 54fb1667..624c5a7d 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -845,54 +845,6 @@ namespace AMDiS {
};
- /**
- * \ingroup Assembler
- *
- * \brief
- * Laplace multiplied with a function which maps the gradient of a DOFVector
- * at each quadrature point to a double:
- * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
- */
- class VecAndGradientAtQP_SOT : public SecondOrderTerm
- {
- public:
- /// Constructor.
- VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
- BinaryAbstractFunction<double, double, WorldVector<double> > *af);
-
- /// Implementation of \ref OperatorTerm::initElement().
- void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = NULL);
-
- /// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
-
- /// Implenetation of SecondOrderTerm::eval().
- void eval(int nPoints,
- const double *uhAtQP,
- const WorldVector<double> *grdUhAtQP,
- const WorldMatrix<double> *D2UhAtQP,
- double *result,
- double factor);
-
- /// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result);
-
- protected:
- DOFVectorBase<double>* vec;
-
- /// Function wich maps \ref gradAtQPs to a double.
- BinaryAbstractFunction<double, double, WorldVector<double> > *f;
-
- /** \brief
- * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
- * at each quadrature point.
- */
- double* vecAtQPs;
- WorldVector<double>* gradAtQPs;
- };
-
/**
* \ingroup Assembler
*
@@ -1819,7 +1771,7 @@ namespace AMDiS {
{
public:
FctVecAtQP_FOT(DOFVectorBase<double> *dv,
- AbstractFunction<double, WorldVector<double> > *f_,
+ BinaryAbstractFunction<double, WorldVector<double>, double> *f_,
WorldVector<double> *b_);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
@@ -1839,7 +1791,7 @@ namespace AMDiS {
DOFVectorBase<double>* vec;
double *vecAtQPs;
WorldVector<double> *coordsAtQPs;
- AbstractFunction<double, WorldVector<double> > *f;
+ BinaryAbstractFunction<double, WorldVector<double>, double> *f;
WorldVector<double> *b;
};
@@ -1848,9 +1800,11 @@ namespace AMDiS {
{
public:
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af,
WorldVector<double> *b_);
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af,
int bIdx);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
@@ -1871,6 +1825,8 @@ namespace AMDiS {
DOFVectorBase<double>* vec2;
double *vec1AtQPs;
double *vec2AtQPs;
+ /// Function f.
+ BinaryAbstractFunction<double, double, double> *f;
WorldVector<double> *b;
};
@@ -1879,11 +1835,11 @@ namespace AMDiS {
{
public:
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
- BinaryAbstractFunction<double, double, double> *f,
+ TertiaryAbstractFunction<double, double, double, double> *f,
WorldVector<double> *b);
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
- BinaryAbstractFunction<double, double, double> *f,
+ TertiaryAbstractFunction<double, double, double, double> *f,
int b);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
@@ -1903,7 +1859,7 @@ namespace AMDiS {
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
DOFVectorBase<double>* vec3;
- BinaryAbstractFunction<double, double, double>* f;
+ TertiaryAbstractFunction<double, double, double, double>* f;
double *vec1AtQPs;
double *vec2AtQPs;
double *vec3AtQPs;
@@ -2507,6 +2463,9 @@ namespace AMDiS {
* \ingroup Assembler
*
* \brief
+ * Laplace multiplied with a function which maps the gradient of a DOFVector
+ * at each quadrature point to a double:
+ * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
*/
class VecAndGradAtQP_SOT : public SecondOrderTerm
{
@@ -2521,7 +2480,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
/// Implements SecondOrderTerm::eval().
void eval(int nPoints,
--
GitLab