Commit 2982e0f3 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Some small changes in Operators. Thanks to Simon ...

parent e32cd0e6
...@@ -767,8 +767,9 @@ namespace AMDiS { ...@@ -767,8 +767,9 @@ namespace AMDiS {
((*b) * grdUhAtQP[iq]); ((*b) * grdUhAtQP[iq]);
} }
FctVecAtQP_FOT::FctVecAtQP_FOT(DOFVectorBase<double> *dv, FctVecAtQP_FOT::FctVecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, WorldVector<double> > *f_, BinaryAbstractFunction<double, WorldVector<double>, double> *f_,
WorldVector<double> *b_) WorldVector<double> *b_)
: FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_) : FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_)
{ {
...@@ -791,9 +792,9 @@ namespace AMDiS { ...@@ -791,9 +792,9 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
if (b) if (b)
lb(Lambda, *b, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq])); lb(Lambda, *b, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
else 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 { ...@@ -807,21 +808,29 @@ namespace AMDiS {
if (grdUhAtQP) if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++) for (int iq = 0; iq < nPoints; iq++)
result[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, Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af,
WorldVector<double> *b_) 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(dv1->getFESpace());
auxFeSpaces.insert(dv2->getFESpace()); auxFeSpaces.insert(dv2->getFESpace());
} }
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af,
int bIdx) 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; bOne = bIdx;
auxFeSpaces.insert(dv1->getFESpace()); auxFeSpaces.insert(dv1->getFESpace());
auxFeSpaces.insert(dv2->getFESpace()); auxFeSpaces.insert(dv2->getFESpace());
...@@ -842,13 +851,13 @@ namespace AMDiS { ...@@ -842,13 +851,13 @@ namespace AMDiS {
if (bOne > -1) { if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++) 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) { } else if (b) {
for (int iq = 0; iq < nPoints; iq++) 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 { } else {
for (int iq = 0; iq < nPoints; iq++) 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 { ...@@ -861,11 +870,14 @@ namespace AMDiS {
{ {
if (grdUhAtQP) if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++) 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) WorldVector<double> *bvec)
: FirstOrderTerm(f_->getDegree()), : FirstOrderTerm(f_->getDegree()),
vec1(dv1), vec1(dv1),
...@@ -879,8 +891,11 @@ namespace AMDiS { ...@@ -879,8 +891,11 @@ namespace AMDiS {
auxFeSpaces.insert(dv3->getFESpace()); 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) int b)
: FirstOrderTerm(f_->getDegree()), : FirstOrderTerm(f_->getDegree()),
vec1(dv1), vec1(dv1),
...@@ -910,17 +925,18 @@ namespace AMDiS { ...@@ -910,17 +925,18 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
if (bOne > -1) { if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++) 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 { } else {
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
if (b) 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 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, void Vec3FctAtQP_FOT::eval(int nPoints,
const double *uhAtQP, const double *uhAtQP,
const WorldVector<double> *grdUhAtQP, const WorldVector<double> *grdUhAtQP,
...@@ -930,8 +946,7 @@ namespace AMDiS { ...@@ -930,8 +946,7 @@ namespace AMDiS {
{ {
if (grdUhAtQP) if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++) for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq] ,vec3AtQPs[iq]) *
(vec1AtQPs[iq]) * (*f)(vec2AtQPs[iq] ,vec3AtQPs[iq]) *
((*b) * grdUhAtQP[iq]); ((*b) * grdUhAtQP[iq]);
} }
...@@ -1053,15 +1068,6 @@ namespace AMDiS { ...@@ -1053,15 +1068,6 @@ namespace AMDiS {
auxFeSpaces.insert(dv->getFESpace()); 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, VecMatrixGradientAtQP_SOT::VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af, BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf, AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
...@@ -1325,14 +1331,6 @@ namespace AMDiS { ...@@ -1325,14 +1331,6 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad); 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, void VecMatrixGradientAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler, SubAssembler* subAssembler,
Quadrature *quad) Quadrature *quad)
...@@ -1581,14 +1579,6 @@ namespace AMDiS { ...@@ -1581,14 +1579,6 @@ namespace AMDiS {
l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq])); 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, void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const DimMat<double> **LALt) const
{ {
...@@ -2156,36 +2146,6 @@ namespace AMDiS { ...@@ -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, void VecMatrixGradientAtQP_SOT::eval(int nPoints,
const double *uhAtQP, const double *uhAtQP,
const WorldVector<double> *grdUhAtQP, const WorldVector<double> *grdUhAtQP,
......
...@@ -845,54 +845,6 @@ namespace AMDiS { ...@@ -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 * \ingroup Assembler
* *
...@@ -1819,7 +1771,7 @@ namespace AMDiS { ...@@ -1819,7 +1771,7 @@ namespace AMDiS {
{ {
public: public:
FctVecAtQP_FOT(DOFVectorBase<double> *dv, FctVecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, WorldVector<double> > *f_, BinaryAbstractFunction<double, WorldVector<double>, double> *f_,
WorldVector<double> *b_); WorldVector<double> *b_);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
...@@ -1839,7 +1791,7 @@ namespace AMDiS { ...@@ -1839,7 +1791,7 @@ namespace AMDiS {
DOFVectorBase<double>* vec; DOFVectorBase<double>* vec;
double *vecAtQPs; double *vecAtQPs;
WorldVector<double> *coordsAtQPs; WorldVector<double> *coordsAtQPs;
AbstractFunction<double, WorldVector<double> > *f; BinaryAbstractFunction<double, WorldVector<double>, double> *f;
WorldVector<double> *b; WorldVector<double> *b;
}; };
...@@ -1848,9 +1800,11 @@ namespace AMDiS { ...@@ -1848,9 +1800,11 @@ namespace AMDiS {
{ {
public: public:
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af,
WorldVector<double> *b_); WorldVector<double> *b_);
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af,
int bIdx); int bIdx);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
...@@ -1871,6 +1825,8 @@ namespace AMDiS { ...@@ -1871,6 +1825,8 @@ namespace AMDiS {
DOFVectorBase<double>* vec2; DOFVectorBase<double>* vec2;
double *vec1AtQPs; double *vec1AtQPs;
double *vec2AtQPs; double *vec2AtQPs;
/// Function f.
BinaryAbstractFunction<double, double, double> *f;
WorldVector<double> *b; WorldVector<double> *b;
}; };
...@@ -1879,11 +1835,11 @@ namespace AMDiS { ...@@ -1879,11 +1835,11 @@ namespace AMDiS {
{ {
public: public:
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3, 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); WorldVector<double> *b);
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3, Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f, TertiaryAbstractFunction<double, double, double, double> *f,
int b); int b);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
...@@ -1903,7 +1859,7 @@ namespace AMDiS { ...@@ -1903,7 +1859,7 @@ namespace AMDiS {
DOFVectorBase<double>* vec1; DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2; DOFVectorBase<double>* vec2;
DOFVectorBase<double>* vec3; DOFVectorBase<double>* vec3;
BinaryAbstractFunction<double, double, double>* f; TertiaryAbstractFunction<double, double, double, double>* f;
double *vec1AtQPs; double *vec1AtQPs;
double *vec2AtQPs; double *vec2AtQPs;
double *vec3AtQPs; double *vec3AtQPs;
...@@ -2507,6 +2463,9 @@ namespace AMDiS { ...@@ -2507,6 +2463,9 @@ namespace AMDiS {
* \ingroup Assembler * \ingroup Assembler
* *
* \brief * \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 class VecAndGradAtQP_SOT : public SecondOrderTerm
{ {
...@@ -2521,7 +2480,7 @@ namespace AMDiS { ...@@ -2521,7 +2480,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt(). /// 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(). /// Implements SecondOrderTerm::eval().
void eval(int nPoints, void eval(int nPoints,
......
Supports Markdown
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