From 57d26e5570b0bbf0c0c504e54e0023903f0a377f Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 8 Feb 2013 12:30:49 +0000
Subject: [PATCH] OperatorTerms updated and BaseClass for Preconditioners added
---
AMDiS/src/FirstOrderTerm.cc | 12 +--
AMDiS/src/FirstOrderTerm.h | 56 ++++++------
AMDiS/src/Functors.h | 86 +++++++++++++------
AMDiS/src/ITL_OEMSolver.h | 1 +
AMDiS/src/ITL_OEMSolver.hh | 36 ++++++++
AMDiS/src/ITL_Preconditioner.h | 8 +-
AMDiS/src/OEMSolver.h | 12 +++
AMDiS/src/io/FileWriter.cc | 8 +-
.../parallel/PetscSolverGlobalBlockMatrix.cc | 9 ++
9 files changed, 165 insertions(+), 63 deletions(-)
diff --git a/AMDiS/src/FirstOrderTerm.cc b/AMDiS/src/FirstOrderTerm.cc
index 2e211892..d88e4a0b 100644
--- a/AMDiS/src/FirstOrderTerm.cc
+++ b/AMDiS/src/FirstOrderTerm.cc
@@ -574,9 +574,9 @@ namespace AMDiS {
}
- // ========== Vec3FctAtQP_FOT ==========
+ // ========== Vec3AtQP_FOT ==========
- Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1,
+ Vec3AtQP_FOT::Vec3AtQP_FOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f_,
@@ -593,7 +593,7 @@ namespace AMDiS {
auxFeSpaces.insert(dv3->getFeSpace());
}
- Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1,
+ Vec3AtQP_FOT::Vec3AtQP_FOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f_,
@@ -610,7 +610,7 @@ namespace AMDiS {
auxFeSpaces.insert(dv3->getFeSpace());
}
- void Vec3FctAtQP_FOT::initElement(const ElInfo* elInfo,
+ void Vec3AtQP_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
@@ -619,7 +619,7 @@ namespace AMDiS {
getVectorAtQPs(vec3, elInfo, subAssembler, quad, vec3AtQPs);
}
- void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo,
+ void Vec3AtQP_FOT::getLb(const ElInfo *elInfo,
vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
@@ -638,7 +638,7 @@ namespace AMDiS {
}
}
- void Vec3FctAtQP_FOT::eval(int nPoints,
+ void Vec3AtQP_FOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
diff --git a/AMDiS/src/FirstOrderTerm.h b/AMDiS/src/FirstOrderTerm.h
index 773bf488..34085b2f 100644
--- a/AMDiS/src/FirstOrderTerm.h
+++ b/AMDiS/src/FirstOrderTerm.h
@@ -132,25 +132,25 @@ namespace AMDiS {
* First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a vector b which
* only consits of entries equal to one.
*/
- class Simple_FOT : public FirstOrderTerm
- {
- public:
- /// Constructor.
- Simple_FOT()
- : FirstOrderTerm(0)
- {}
-
- /// Implements FirstOrderTerm::getLb().
- inline void getLb(const ElInfo *elInfo,
- vector<mtl::dense_vector<double> >& Lb) const
- {
- const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- const int nPoints = static_cast<int>(Lb.size());
-
- for (int iq = 0; iq < nPoints; iq++)
- l1(grdLambda, Lb[iq], 1.0);
- }
- };
+// class Simple_FOT : public FirstOrderTerm
+// {
+// public:
+// /// Constructor.
+// Simple_FOT()
+// : FirstOrderTerm(0)
+// {}
+//
+// /// Implements FirstOrderTerm::getLb().
+// inline void getLb(const ElInfo *elInfo,
+// vector<mtl::dense_vector<double> >& Lb) const
+// {
+// const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+// const int nPoints = static_cast<int>(Lb.size());
+//
+// for (int iq = 0; iq < nPoints; iq++)
+// l1(grdLambda, Lb[iq], 1.0);
+// }
+// };
/**
* \ingroup Assembler
@@ -159,11 +159,11 @@ namespace AMDiS {
* First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a vector b which
* only consits of entries equal to one.
*/
- class FactorSimple_FOT : public FirstOrderTerm
+ class Simple_FOT : public FirstOrderTerm
{
public:
/// Constructor.
- FactorSimple_FOT(double f)
+ Simple_FOT(double f = 1.0)
: FirstOrderTerm(0)
{
factor = new double;
@@ -171,7 +171,7 @@ namespace AMDiS {
}
/// Constructor.
- FactorSimple_FOT(double *fptr)
+ Simple_FOT(double *fptr)
: FirstOrderTerm(0), factor(fptr)
{}
@@ -191,6 +191,9 @@ namespace AMDiS {
double *factor;
};
+ // for compatibility
+ typedef Simple_FOT FactorSimple_FOT;
+
/**
* \ingroup Assembler
*
@@ -703,14 +706,14 @@ namespace AMDiS {
};
- class Vec3FctAtQP_FOT : public FirstOrderTerm
+ class Vec3AtQP_FOT : public FirstOrderTerm
{
public:
- Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
+ Vec3AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f,
WorldVector<double> *b);
- Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
+ Vec3AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f,
int b);
@@ -736,6 +739,9 @@ namespace AMDiS {
WorldVector<double> *b;
};
+
+ // for compatibility
+ typedef Vec3AtQP_FOT Vec3FctAtQP_FOT;
class General_FOT : public FirstOrderTerm
diff --git a/AMDiS/src/Functors.h b/AMDiS/src/Functors.h
index 854532e5..d2f764cf 100644
--- a/AMDiS/src/Functors.h
+++ b/AMDiS/src/Functors.h
@@ -10,52 +10,55 @@ namespace AMDiS {
template<typename T>
struct Id : public AbstractFunction<T,T>
{
+ Id() : AbstractFunction<T,T>(1) {}
T operator()(const T &v) const { return v; }
};
template<typename T1,typename T2>
struct Const : public AbstractFunction<T1,T2>
{
- Const(T1 val_) : val(val_) {}
+ Const(T1 val_) : AbstractFunction<T1,T2>(0), val(val_) {}
T1 operator()(const T2 &v) const { return val; }
private:
T1 val;
};
-template<typename T>
+template<typename T=double>
struct Factor : public AbstractFunction<T,T>
{
- Factor(double fac_) : fac(fac_) {}
+ Factor(double fac_) : AbstractFunction<T,T>(1), fac(fac_) {}
T operator()(const T& x) const { return fac*x; };
private:
double fac;
};
-template<typename T>
+template<typename T=double>
struct Add : public BinaryAbstractFunction<T,T,T>
{
+ Add() : BinaryAbstractFunction<T,T,T>(1) {}
T operator()(const T &v1, const T &v2) const { return v1+v2; }
};
-template<typename T>
+template<typename T=double>
struct AddFactor : public BinaryAbstractFunction<T,T,T>
{
- AddFactor(double factor_ = 1.0) : factor(factor_) {}
+ AddFactor(double factor_ = 1.0) : BinaryAbstractFunction<T,T,T>(1), factor(factor_) {}
T operator()(const T &v1, const T &v2) const { return v1 + factor*v2; }
private:
double factor;
};
-template<typename T>
+template<typename T=double>
struct Subtract : public BinaryAbstractFunction<T,T,T>
{
+ Subtract() : BinaryAbstractFunction<T,T,T>(1) {}
T operator()(const T &v1, const T &v2) const { return v1-v2; }
};
-template<typename T>
+template<typename T=double>
struct AddScal : public AbstractFunction<T,T>
{
- AddScal(T scal_) : scal(scal_) {}
+ AddScal(T scal_) : AbstractFunction<T,T>(1), scal(scal_) {}
T operator()(const T &v) const { return v+scal; }
private:
T scal;
@@ -64,19 +67,21 @@ private:
template<typename T>
struct Mult : public BinaryAbstractFunction<T,T,T>
{
+ Mult() : BinaryAbstractFunction<T,T,T>(2) {}
T operator()(const T &v1, const T &v2) const { return v1*v2; }
};
template<typename T1, typename T2, typename T3>
struct Mult2 : public BinaryAbstractFunction<T1,T2,T3>
{
+ Mult2() : BinaryAbstractFunction<T1,T2,T3>(2) {}
T1 operator()(const T2 &v1, const T3 &v2) const { return v1*v2; }
};
template<typename T>
struct MultScal : public BinaryAbstractFunction<T,T,T>
{
- MultScal(T scal_) : scal(scal_) {}
+ MultScal(T scal_) : BinaryAbstractFunction<T,T,T>(2), scal(scal_) {}
T operator()(const T &v1, const T &v2) const { return v1*v2*scal; }
private:
T scal;
@@ -85,87 +90,120 @@ private:
template<typename T>
struct Max : public BinaryAbstractFunction<T,T,T>
{
+ Max() : BinaryAbstractFunction<T,T,T>(1) {}
T operator()(const T &v1, const T &v2) const { return std::max(v1,v2); }
};
template<typename T>
struct Min : public BinaryAbstractFunction<T,T,T>
{
+ Min() : BinaryAbstractFunction<T,T,T>(1) {}
T operator()(const T &v1, const T &v2) const { return std::min(v1,v2); }
};
template<typename T>
struct Diff : public BinaryAbstractFunction<T,T,T>
{
+ Diff() : BinaryAbstractFunction<T,T,T>(1) {}
T operator()(const T &v1, const T &v2) const { return abs(v1-v2); }
};
-template<typename T>
+template<typename T=double>
struct Abs : public AbstractFunction<T,T>
{
+ Abs() : AbstractFunction<T,T>(1) {}
T operator()(const T &v) const { return abs(v); }
};
-template<typename T>
+template<typename T=double>
struct Signum : public AbstractFunction<T,T>
{
+ Signum() : AbstractFunction<T,T>(0) {}
T operator()(const T &v) const { return (v>0.0?1.0:(v<0.0?-1.0:0.0)); }
};
-template<typename T>
+template<typename T=double>
struct Sqr : public AbstractFunction<T,T>
{
+ Sqr() : AbstractFunction<T, T>(2) {}
T operator()(const T &v) const { return sqr(v); }
};
-template<typename T>
+template<typename T=double>
struct Sqrt : public AbstractFunction<T,T>
{
T operator()(const T &v) const { return sqrt(v); }
};
-template<typename T>
+namespace detail {
+ template<int p>
+ struct Pow
+ {
+ template<typename T>
+ static T eval(const T& v) { return v*Pow<p-1>::eval(v); }
+ };
+
+ template<>
+ struct Pow<1> {
+ template<typename T>
+ static T eval(const T& v) { return v; }
+ };
+
+ template<>
+ struct Pow<0> {
+ template<typename T>
+ static T eval(const T& v) { return 1.0; }
+ };
+}
+
+template<int p, typename T=double>
struct Pow : public AbstractFunction<T,T>
{
- Pow(double p_) : p(p_) {}
- T operator()(const T &v) const { return pow(v,p); }
+ Pow(double factor_=1.0) : AbstractFunction<T,T>(p), factor(factor_) {}
+ T operator()(const T &v) const { return factor * detail::Pow<p>::eval(v); }
private:
- double p;
+ double factor;
};
template<typename T1, typename T2 = ProductType<T1, T1> >
struct Norm2 : public AbstractFunction<T1, T2>
{
+ Norm2() : AbstractFunction<T1, T2>(2) {}
T1 operator()(const T2 &v) const { return sqrt(v*v); }
};
template<typename T1, typename T2>
struct Norm2Sqr : public AbstractFunction<T1, T2>
{
+ Norm2Sqr() : AbstractFunction<T1, T2>(2) {}
T1 operator()(const T2 &v) const { return v*v; }
};
template<typename T>
struct Norm2_comp2 : public BinaryAbstractFunction<T,T,T>
{
+ Norm2_comp2() : BinaryAbstractFunction<T,T,T>(2) {}
T operator()(const T &v1, const T &v2) const { return sqrt(sqr(v1)+sqr(v2)); }
};
template<typename T>
struct Norm2Sqr_comp2 : public BinaryAbstractFunction<T,T,T>
{
+ Norm2Sqr_comp2() : BinaryAbstractFunction<T,T,T>(2) {}
T operator()(const T &v1, const T &v2) const { return sqr(v1)+sqr(v2); }
};
template<typename T>
struct Norm2_comp3 : public TertiaryAbstractFunction<T,T,T,T>
{
+ Norm2_comp3() : TertiaryAbstractFunction<T,T,T,T>(2) {}
T operator()(const T &v1, const T &v2, const T &v3) const { return sqrt(sqr(v1)+sqr(v2)+sqr(v3)); }
};
template<typename T>
struct Norm2Sqr_comp3 : public TertiaryAbstractFunction<T,T,T,T>
{
+ Norm2Sqr_comp3() : TertiaryAbstractFunction<T,T,T,T>(2) {}
T operator()(const T &v1, const T &v2, const T &v3) const { return sqr(v1)+sqr(v2)+sqr(v3); }
};
@@ -208,15 +246,13 @@ struct Vec3WorldVec : public TertiaryAbstractFunction<WorldVector<T>,T,T,T>
return result;
}
};
-
-struct Component : public AbstractFunction<double, WorldVector<double> >
+template<int c, typename T=double>
+struct Component : public AbstractFunction<T, WorldVector<T> >
{
- Component(int comp_) : comp(comp_) {}
- double operator()(const WorldVector<double> &x) const {
- return x[comp];
+ Component() : AbstractFunction<T, WorldVector<T> >(1) {}
+ T operator()(const WorldVector<T> &x) const {
+ return x[c];
}
-private:
- int comp;
};
struct FadeOut : public TertiaryAbstractFunction<double, double, double ,double>
diff --git a/AMDiS/src/ITL_OEMSolver.h b/AMDiS/src/ITL_OEMSolver.h
index 7b126669..8e5d7fe5 100644
--- a/AMDiS/src/ITL_OEMSolver.h
+++ b/AMDiS/src/ITL_OEMSolver.h
@@ -63,6 +63,7 @@ namespace AMDiS {
~ITL_OEMSolver() {}
+
/// Creator class used in the OEMSolverMap.
class Creator : public OEMSolverCreator
diff --git a/AMDiS/src/ITL_OEMSolver.hh b/AMDiS/src/ITL_OEMSolver.hh
index db6e3688..851d8e3e 100644
--- a/AMDiS/src/ITL_OEMSolver.hh
+++ b/AMDiS/src/ITL_OEMSolver.hh
@@ -65,6 +65,11 @@ namespace AMDiS {
void init(const Matrix& A)
{
BaseITL_runner::setPrecon(preconPair, A);
+
+ if (preconPair.l != NULL)
+ preconPair.l->init();
+ if (preconPair.r != NULL)
+ preconPair.r->init();
}
int solve(const Matrix& A, Vector& x, Vector& b)
@@ -88,15 +93,28 @@ namespace AMDiS {
void exit()
{
if (preconPair.l != NULL) {
+ preconPair.l->exit();
delete preconPair.l;
preconPair.l = NULL;
}
if (preconPair.r != NULL) {
+ preconPair.r->exit();
delete preconPair.r;
preconPair.r = NULL;
}
}
+
+ virtual BasePreconditionerType* getLeftPrecon()
+ {
+ return preconPair.l;
+ }
+
+ virtual BasePreconditionerType* getRightPrecon()
+ {
+ return preconPair.r;
+ }
+
private:
PreconPair< Vector > preconPair;
ITLSolver solver;
@@ -112,6 +130,11 @@ namespace AMDiS {
void init(const Matrix& A)
{
BaseITL_runner::setPrecon(preconPair, A);
+
+ if (preconPair.l != NULL)
+ preconPair.l->init();
+ if (preconPair.r != NULL)
+ preconPair.r->init();
}
int solve(const Matrix& A , Vector& x, Vector& b)
@@ -135,15 +158,28 @@ namespace AMDiS {
void exit()
{
if (preconPair.l != NULL) {
+ preconPair.l->exit();
delete preconPair.l;
preconPair.l = NULL;
}
if (preconPair.r != NULL) {
+ preconPair.r->exit();
delete preconPair.r;
preconPair.r = NULL;
}
}
+
+ virtual BasePreconditionerType* getLeftPrecon()
+ {
+ return preconPair.l;
+ }
+
+ virtual BasePreconditionerType* getRightPrecon()
+ {
+ return preconPair.r;
+ }
+
private:
ITLSolver solver;
PreconPair< Vector > preconPair;
diff --git a/AMDiS/src/ITL_Preconditioner.h b/AMDiS/src/ITL_Preconditioner.h
index 12a094ea..2ad84d9f 100644
--- a/AMDiS/src/ITL_Preconditioner.h
+++ b/AMDiS/src/ITL_Preconditioner.h
@@ -34,6 +34,12 @@
namespace AMDiS {
+ struct BasePreconditionerType
+ {
+ virtual void init() {};
+ virtual void exit() {};
+ };
+
/**
* \ingroup Solver
@@ -41,7 +47,7 @@ namespace AMDiS {
* \brief Common base class for wrappers to use ITL preconditioners in AMDiS.
*/
template< class Vec >
- class ITL_BasePreconditioner
+ class ITL_BasePreconditioner : public BasePreconditionerType
{
public:
virtual ~ITL_BasePreconditioner() {}
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index 2eea5d0b..5be6a989 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -185,6 +185,18 @@ namespace AMDiS {
return relative;
}
+
+ virtual BasePreconditionerType* getLeftPrecon()
+ {
+ WARNING("no left preconditioner provided!\n");
+ return NULL;
+ }
+
+ virtual BasePreconditionerType* getRightPrecon()
+ {
+ WARNING("no right preconditioner provided!\n");
+ return NULL;
+ }
/** \} */
diff --git a/AMDiS/src/io/FileWriter.cc b/AMDiS/src/io/FileWriter.cc
index 6d299f08..378762b7 100644
--- a/AMDiS/src/io/FileWriter.cc
+++ b/AMDiS/src/io/FileWriter.cc
@@ -159,16 +159,12 @@ namespace AMDiS {
if (writeParaViewVectorFormat) {
#if HAVE_PARALLEL_DOMAIN_AMDIS
VtkVectorWriter::writeFile(solutionVecs, paraFilename + paraviewFileExt, true, writeAs3dVector);
+ MSG("ParaView file written to %s\n", (paraFilename + paraviewFileExt).c_str());
#else
VtkVectorWriter::writeFile(solutionVecs, fn + paraviewFileExt, true, writeAs3dVector);
+ MSG("ParaView file written to %s\n", (fn + paraviewFileExt).c_str());
#endif
-// VtkVectorWriter::Impl<double> vtkVectorWriter(&dataCollectors, writeAs3dVector);
-// vtkVectorWriter.setCompression(compression);
-// vtkVectorWriter.setWriteAs3dVector(writeAs3dVector);
-// vtkVectorWriter.writeFile(fn + paraviewFileExt);
-
- MSG("ParaView file written to %s\n", (fn + paraviewFileExt).c_str());
}
if (writeParaViewAnimation) {
diff --git a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
index 3bdd3fc3..bdc6cd52 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
@@ -175,6 +175,15 @@ namespace AMDiS {
const FiniteElemSpace *feSpace = componentSpaces[0];
VecDuplicate(getVecRhsInterior(), &petscSolVec);
+
+ for (int i = 0; i < vec.getSize(); i++)
+ {
+ Vec tmp;
+ VecNestGetSubVec(petscSolVec, i, &tmp);
+ setDofVector(tmp, vec.getDOFVector(i));
+ VecAssemblyBegin(tmp);
+ VecAssemblyEnd(tmp);
+ }
// PETSc.
solve(getVecRhsInterior(), petscSolVec);
--
GitLab