OperatorTerms updated and BaseClass for Preconditioners added

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,

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

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>

AMDiS/src/ITL_OEMSolver.h

@@ -63,6 +63,7 @@ namespace AMDiS {
 
   
     ~ITL_OEMSolver() {}
+    
 
     /// Creator class used in the OEMSolverMap.
     class Creator : public OEMSolverCreator

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;

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() {}

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;
+    }
 
     /** \} */ 
 

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) {

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);