Fixed a problem with virtual function in ProblemInstat, rewritten BasisFct::interpol

AMDiS/src/BasisFunction.h

@@ -181,7 +181,7 @@ namespace AMDiS {
     virtual DimVec<double> *getCoords(int i) const = 0;
 
     /** \brief
-     * Returns a pointer to a const vector with interpolation coefficients of the
+     * Fills a vector with interpolation coefficients of the
      * function f; if indices is a pointer to NULL, the coefficient for all 
      * basis functions are calculated and the i-th entry in the vector is the 
      * coefficient of the i-th basis function; if indices is non NULL, only the 
@@ -197,16 +197,18 @@ namespace AMDiS {
      * during mesh traversal.
      * Must be implemented by sub classes.
      */
-    virtual const double* interpol(const ElInfo *el_info, 
-				   int n, const int *indices, 
-				   AbstractFunction<double, WorldVector<double> > *f,
-				   double *coeff) = 0;
+    virtual void interpol(const ElInfo *el_info, 
+			  int n, 
+			  const int *indices, 
+			  AbstractFunction<double, WorldVector<double> > *f,
+			  mtl::dense_vector<double> &coeff) const = 0;
 
     /// WorldVector<double> valued interpol function.
-    virtual const WorldVector<double>* interpol(const ElInfo *el_info, int no, 
-						const int *b_no,
-						AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
-						WorldVector<double> *vec) = 0;
+    virtual void interpol(const ElInfo *el_info, 
+			  int no, 
+			  const int *b_no,
+			  AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
+			  mtl::dense_vector<WorldVector<double> >& coeff) const = 0;
 
     /// Returns the i-th local basis function
     inline BasFctType *getPhi(int i) const 

AMDiS/src/DOFVector.h

@@ -572,9 +572,6 @@ namespace AMDiS {
     /// Returns the average value of the DOFVector.
     T average() const;
 
-    /// Used by interpol while mesh traversal.
-    void interpolFct(ElInfo* elinfo);
-
     /// Prints \ref vec to stdout.
     void print() const; 
 
@@ -641,12 +638,6 @@ namespace AMDiS {
  
     /// Specifies what operation should be performed after coarsening
     CoarsenOperation coarsenOperation;
-
-    /// Used by \ref interpol
-    AbstractFunction<T, WorldVector<double> > *interFct;
-
-    /// Used for mesh traversal
-    static DOFVector<T> *traverseVector;
   }; 
 
 
@@ -795,7 +786,7 @@ namespace AMDiS {
 
   // y = a*x + y
   template<typename T>
-  void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);
+  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
 
   // matrix vector product
   template<typename T>
@@ -856,6 +847,12 @@ namespace AMDiS {
     return vec->getUsedSize();
   }
 
+  template<typename T>
+  inline int size(const DOFVector<T>& vec) 
+  {
+    return vec.getUsedSize();
+  }
+
   template<typename T>
   inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
   {

AMDiS/src/DOFVector.hh

@@ -132,9 +132,6 @@ namespace AMDiS {
     vec.clear();
   }
 
-  template<typename T>
-  DOFVector<T> * DOFVector<T>::traverseVector = NULL;
-
   template<typename T>
   void DOFVectorBase<T>::addElementVector(T factor, 
 					  const ElementVector &elVec, 
@@ -450,8 +447,6 @@ namespace AMDiS {
 
     TEST_EXIT_DBG(fct)("No function to interpolate!\n");
 
-    interFct = fct;
-
     if (!this->getFeSpace()) {
       MSG("no dof admin in vec %s, skipping interpolation\n",
 	  this->getName().c_str());
@@ -476,36 +471,27 @@ namespace AMDiS {
       return;
     }
 
-    traverseVector = this;
+    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
+    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
+    int nBasFcts = basFct->getNumber();
+    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
+    mtl::dense_vector<double> fctInterpolValues(nBasFcts);
 
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
 					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
     while (elInfo) {
-      interpolFct(elInfo);
+      basFct->interpol(elInfo, 0, NULL, fct, fctInterpolValues);
+      basFct->getLocalIndices(const_cast<Element*>(elInfo->getElement()), 
+			      admin, myLocalIndices);
+      for (int i = 0; i < nBasFcts; i++)
+	vec[myLocalIndices[i]] = fctInterpolValues[i];
+
       elInfo = stack.traverseNext(elInfo);
     }
   }
 
 
-  template<typename T>
-  void DOFVector<T>::interpolFct(ElInfo* elinfo)
-  {
-    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
-    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
-    const T *inter_val = 
-      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
-						   traverseVector->interFct, NULL);
-
-    int nBasFcts = basFct->getNumber();
-    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
-    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), 
-			    admin, myLocalIndices);
-    for (int i = 0; i < nBasFcts; i++)
-      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
-  }
-
-
   template<typename T>
   double DOFVector<T>::Int(Quadrature* q) const
   {
@@ -1047,7 +1033,6 @@ namespace AMDiS {
     this->nBasFcts = rhs.nBasFcts;
     vec = rhs.vec;
     this->elementVector.change_dim(this->nBasFcts);
-    interFct = rhs.interFct;
     coarsenOperation = rhs.coarsenOperation;
     this->operators = rhs.operators;
     this->operatorFactor = rhs.operatorFactor;

AMDiS/src/Lagrange.cc

@@ -811,43 +811,21 @@ namespace AMDiS {
   }
 
 
-  const double* Lagrange::interpol(const ElInfo *el_info, 
-				   int no, const int *b_no, 
-				   AbstractFunction<double, WorldVector<double> > *f, 
-				   double *vec)
+  void Lagrange::interpol(const ElInfo *el_info, 
+			  int no, 
+			  const int *b_no, 
+			  AbstractFunction<double, WorldVector<double> > *f, 
+			  mtl::dense_vector<double> &rvec) const
   {
     FUNCNAME("Lagrange::interpol()");
 
-    static double* localVec = NULL;
-    static int localVecSize = 0;
-
-    double *rvec = NULL;
-
-    if (vec) {
-      rvec = vec;
-    } else {
-      if (localVec && nBasFcts > localVecSize)  {
-	delete [] localVec;
-	localVec = new double[nBasFcts]; 
-      }
-      if (!localVec)
-	localVec = new double[nBasFcts]; 
-
-      localVecSize = nBasFcts;
-      rvec = localVec;
-    } 
-
     WorldVector<double> x;
 
     el_info->testFlag(Mesh::FILL_COORDS);
  
     if (b_no) {
-      if (no <= 0  ||  no > getNumber()) {
-	ERROR("something is wrong, doing nothing\n");
-	rvec[0] = 0.0;
-	return(const_cast<const double *>(rvec));
-      }
-	
+      TEST_EXIT_DBG(no >= 0 && no < getNumber())("Something is wrong!\n");
+
       for (int i = 0; i < no; i++) {
 	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
 	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
@@ -865,23 +843,17 @@ namespace AMDiS {
 	rvec[i] = (*f)(x);
       }
     }  
-    
-    return(const_cast<const double *>( rvec));
   }
 
 
-  const WorldVector<double>* 
-  Lagrange::interpol(const ElInfo *el_info, 
-		     int no, 
-		     const int *b_no,
-		     AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
-		     WorldVector<double> *vec)
+  void Lagrange::interpol(const ElInfo *el_info, 
+			  int no, 
+			  const int *b_no,
+			  AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
+			  mtl::dense_vector<WorldVector<double> > &rvec) const
   {
-    FUNCNAME("*Lagrange::interpol_d()");
+    FUNCNAME("*Lagrange::interpol()");
 
-    static WorldVector<double> *inter;
-    WorldVector<double> *rvec = 
-      vec ? vec : (inter ? inter : inter = new WorldVector<double>[getNumber()]);
     WorldVector<double> x;
 
     el_info->testFlag(Mesh::FILL_COORDS);
@@ -889,12 +861,8 @@ namespace AMDiS {
     int vertices = Global::getGeo(VERTEX, dim);
 
     if (b_no) {
-      if (no <= 0  ||  no > getNumber()) {
-	ERROR("something is wrong, doing nothing\n");
-	for (int i = 0; i < dow; i++)
-	  rvec[0][i] = 0.0;
-	return(const_cast<const WorldVector<double> *>( rvec));
-      }
+      TEST_EXIT_DBG(no >= 0 && no < getNumber())("Something is wrong!\n");
+
       for (int i = 0; i < no; i++) {
 	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
 	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
@@ -911,8 +879,6 @@ namespace AMDiS {
 	rvec[i] = (*f)(x);
       }
     }  
-    
-    return(const_cast<const WorldVector<double> *>( rvec));
   }
 
 

AMDiS/src/Lagrange.h

@@ -64,16 +64,15 @@ namespace AMDiS {
     static Lagrange* getLagrange(int dim, int degree);
 
     /// Implements BasisFunction::interpol
-    const double *interpol(const ElInfo *, int, const int *, 
-			   AbstractFunction<double, WorldVector<double> >*, 
-			   double *);
+    void interpol(const ElInfo *, int, const int *, 
+		  AbstractFunction<double, WorldVector<double> >*, 
+		  mtl::dense_vector<double>&) const;
 
     /// Implements BasisFunction::interpol
-    const WorldVector<double> *interpol(const ElInfo *, int, 
-					const int *b_no,
-					AbstractFunction<WorldVector<double>, 
-					WorldVector<double> >*, 
-					WorldVector<double> *);
+    void interpol(const ElInfo *, int, 
+		  const int *b_no,
+		  AbstractFunction<WorldVector<double>, WorldVector<double> >*, 
+		  mtl::dense_vector<WorldVector<double> >&) const;
 
     /// Returns the barycentric coordinates of the i-th basis function.
     DimVec<double> *getCoords(int i) const;

AMDiS/src/ProblemInstat.h

@@ -164,13 +164,13 @@ namespace AMDiS {
 		    Flag adoptFlag = INIT_NOTHING);
 
     /// Used in \ref initialize().
-    void createUhOld();
+    virtual void createUhOld();
 
     /// Implementation of ProblemTimeInterface::initTimestep().
-    void initTimestep(AdaptInfo *adaptInfo);
+    virtual void initTimestep(AdaptInfo *adaptInfo);
 
     /// Implementation of ProblemTimeInterface::closeTimestep().
-    void closeTimestep(AdaptInfo *adaptInfo);
+    virtual void closeTimestep(AdaptInfo *adaptInfo);
   
     /// Returns \ref problemStat.
     inline ProblemStatSeq* getStatProblem() 
@@ -190,11 +190,11 @@ namespace AMDiS {
     }
 
     /// Used by \ref problemInitial
-    void transferInitialSolution(AdaptInfo *adaptInfo);  
+    virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
 
-    void serialize(std::ostream &out) {}
+    virtual void serialize(std::ostream &out) {}
 
-    void deserialize(std::istream &in) {}
+    virtual void deserialize(std::istream &in) {}
 
   protected:
     /// Space problem solved in each timestep.

AMDiS/src/ProblemStat.cc

@@ -1065,7 +1065,7 @@ namespace AMDiS {
   
 
   void ProblemStatSeq::dualAssemble(AdaptInfo *adaptInfo, Flag flag, 
-				bool asmMatrix, bool asmVector)
+				    bool asmMatrix, bool asmVector)
   {
     FUNCNAME("ProblemStat::dualAssemble()");
 
@@ -1631,7 +1631,8 @@ namespace AMDiS {
 
   void ProblemStatSeq::assembleOnOneMesh(FiniteElemSpace *feSpace, 
 					 Flag assembleFlag,
-					 DOFMatrix *matrix, DOFVector<double> *vector)
+					 DOFMatrix *matrix, 
+					 DOFVector<double> *vector)
   {
     FUNCNAME("ProblemStat::assembleOnOnMesh()");