some error corrections

AMDiS/src/DOFVector.cc

@@ -68,10 +68,9 @@ namespace AMDiS {
 
 
   template<>
-  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
-					       ElInfo *oldElInfo) const
+  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
   {
-    FUNCNAME("DOFVector<double>::evalAtPoint()");
+    FUNCNAME("DOFVector<double>::evalAtCoords()");
   
     Mesh *mesh = feSpace->getMesh();
     const BasisFunction *basFcts = feSpace->getBasisFcts();
@@ -79,7 +78,7 @@ namespace AMDiS {
     int dim = mesh->getDim();
     int nBasFcts = basFcts->getNumber();
   
-    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
+    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
     DimVec<double> lambda(dim, NO_INIT);
   
     ElInfo *elInfo = mesh->createNewElInfo();
@@ -87,8 +86,7 @@ namespace AMDiS {
     bool inside = false;
   
     if (oldElInfo && oldElInfo->getMacroElement()) {
-      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, 
-				       oldElInfo->getMacroElement(), NULL, NULL);
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
       delete oldElInfo;
     } else
       inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
@@ -97,8 +95,7 @@ namespace AMDiS {
       oldElInfo = elInfo;
 
     if (inside) {
-      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), 
-			       localIndices);
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
       ElementVector uh(lambda.size());
       for (int i = 0; i < lambda.size(); i++)
         uh[i] = operator[](localIndices[i]);
@@ -106,19 +103,20 @@ namespace AMDiS {
     } else
       throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
 
+    delete [] localIndices;
     if (oldElInfo == NULL)
       delete elInfo;
 
+    if (values != NULL)
+      *values = value;
     return value;
-  }
+  };
 
 
   template<>
-  const WorldVector<double> 
-  DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
-					       ElInfo *oldElInfo) const
+  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
   {
-    FUNCNAME("DOFVector<double>::evalAtPoint()");
+    FUNCNAME("DOFVector<double>::evalAtCoords()");
   
     Mesh *mesh = feSpace->getMesh();
     const BasisFunction *basFcts = feSpace->getBasisFcts();
@@ -126,17 +124,18 @@ namespace AMDiS {
     int dim = mesh->getDim();
     int nBasFcts = basFcts->getNumber();
   
-    vector<DegreeOfFreedom> localIndices(nBasFcts);
+    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
     DimVec<double> lambda(dim, NO_INIT);
   
     ElInfo *elInfo = mesh->createNewElInfo();
 
-    WorldVector<double> value(DEFAULT_VALUE, 0.0);
+    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
+    WorldVector<double> *val = (NULL != values) ? values : &Values;
+
     bool inside = false;
   
     if (oldElInfo && oldElInfo->getMacroElement()) {
-      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, 
-				       oldElInfo->getMacroElement(), NULL, NULL);
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
       delete oldElInfo;
     } else
       inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
@@ -145,20 +144,20 @@ namespace AMDiS {
       oldElInfo = elInfo;
 
     if (inside) {
-      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), 
-			       localIndices);
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
       mtl::dense_vector<WorldVector<double> > uh(lambda.size());
       for (int i = 0; i < lambda.size(); i++)
         uh[i] = operator[](localIndices[i]);
-      value = basFcts->evalUh(lambda, uh);
+      *val = basFcts->evalUh(lambda, uh);
     } else
       throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
 
+    delete [] localIndices;
     if (oldElInfo == NULL)
       delete elInfo;
 
-    return value;
-  }
+    return ((*val));
+  };
 
 
   template<>
@@ -911,11 +910,11 @@ namespace AMDiS {
     TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");
 
     int deg = std::max(fct->getDegree(),
-		       2 * vecs[0].getFeSpace()->getBasisFcts()->getDegree());
+		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
     Quadrature* quad =
-      Quadrature::provideQuadrature(vecs[0].getFeSpace()->getMesh()->getDim(), deg);
+      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
     FastQuadrature *fastQuad =
-      FastQuadrature::provideFastQuadrature(vecs[0].getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
+      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
 
     std::vector<mtl::dense_vector<double> > qp(vecs.size());
     std::vector<double> qp_local(vecs.size());
@@ -925,10 +924,10 @@ namespace AMDiS {
     double value = 0.0;
     Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
+    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
     while (elInfo) {
       for (size_t i = 0; i < vecs.size(); i++)
-	vecs[i].getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
+	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
 
       double tmp = 0.0;
       for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
@@ -951,7 +950,7 @@ namespace AMDiS {
 
   double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
 				  const std::vector<DOFVector<double>*> &grds,
-				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > *fct)
+				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct)
   {
     FUNCNAME("integrateGeneral()");
 
@@ -960,11 +959,11 @@ namespace AMDiS {
     TEST_EXIT(grds.size()>0)("No DOFVectors for gradients provided!\n");
 
     int deg = std::max(fct->getDegree(),
-		       2 * vecs[0].getFeSpace()->getBasisFcts()->getDegree());
+		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
     Quadrature* quad =
-      Quadrature::provideQuadrature(vecs[0].getFeSpace()->getMesh()->getDim(), deg);
+      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
     FastQuadrature *fastQuad =
-      FastQuadrature::provideFastQuadrature(vecs[0].getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
+      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
 
     std::vector<mtl::dense_vector<double> > qp(vecs.size());
     std::vector<mtl::dense_vector<WorldVector<double> > > qpGrd(vecs.size());
@@ -978,12 +977,12 @@ namespace AMDiS {
     double value = 0.0;
     Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
+    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
     while (elInfo) {
       for (size_t i = 0; i < vecs.size(); i++)
-	vecs[i].getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
+	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
       for (size_t i = 0; i < grds.size(); i++)
-	grds[i].getGradientAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
+	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
 
       double tmp = 0.0;
       for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {

AMDiS/src/DOFVector.h

@@ -65,13 +65,17 @@ namespace AMDiS {
 
     virtual ~DOFVectorBase();
 
-    /// For the given element, this function returns an array of all DOFs of this
-    /// DOFVector that are defined on this element.
+    /** \brief 
+     * For the given element, this function returns an array of all DOFs of this
+     * DOFVector that are defined on this element.
+     */
     virtual void getLocalVector(const Element *el, 
 				mtl::dense_vector<T>& localVec) const;
 
-    /// Evaluates the DOF vector at a set of quadrature points defined on the 
-    /// given element.
+    /** \brief
+     * Evaluates the DOF vector at a set of quadrature points defined on the 
+     * given element.
+     */
     void getVecAtQPs(const ElInfo *elInfo, 
 		     const Quadrature *quad,
 		     const FastQuadrature *quadFast,
@@ -83,8 +87,10 @@ namespace AMDiS {
 		     const FastQuadrature *quadFast,
 		     mtl::dense_vector<T>& vecAtQPs) const;
 
-    /// Evaluates the gradient of a DOF vector at a set of quadrature points 
-    /// defined on the given element.
+    /** \brief
+     * Evaluates the gradient of a DOF vector at a set of quadrature points defined on the
+     * given element.
+     */
     void getGrdAtQPs( const ElInfo *elInfo,
 		      const Quadrature *quad,
 		      const FastQuadrature *quadFast,
@@ -96,8 +102,10 @@ namespace AMDiS {
 		      const FastQuadrature *quadFast,
 		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;
 
-    /// Evaluates the comp'th component of the derivative of a DOF vector at a 
-    /// set of quadrature points defined on the given element.
+    /** \brief
+     * Evaluates the comp'th component of the derivative of a DOF vector at a set of quadrature points defined on the
+     * given element.
+     */
     void getDerivativeAtQPs( const ElInfo *elInfo,
 		      const Quadrature *quad,
 		      const FastQuadrature *quadFast,
@@ -111,8 +119,10 @@ namespace AMDiS {
 		      int comp,
 		      mtl::dense_vector<T> &derivativeAtQPs) const;
 
-    /// Evaluates the jacobian of a DOF vector at a set of quadrature points 
-    /// defined on the given element.
+    /** \brief
+     * Evaluates the jacobian of a DOF vector at a set of quadrature points defined on the
+     * given element.
+     */
     void getD2AtQPs(const ElInfo *elInfo,
 		    const Quadrature *quad,
 		    const FastQuadrature *quadFast,
@@ -124,8 +134,10 @@ namespace AMDiS {
       return feSpace;
     }
 
-    /// Assembles the element vector for the given ellement and adds the element
-    /// matrix to the current DOF vector.
+    /** \brief
+     * Assembles the element vector for the given ellement and adds the
+     * element matrix to the current DOF vector.
+     */ 
     void assemble(T factor, ElInfo *elInfo,			    
 		  const BoundaryType *bound, 
 		  Operator *op = NULL);
@@ -142,9 +154,11 @@ namespace AMDiS {
 			  ElInfo *elInfo,
 			  bool add = true); 
 
-    /// That function must be called after the matrix assembling has been 
-    /// finished. This makes it possible to start some cleanup or matrix data 
-    /// compressing procedures.
+    /* \brief
+     * That function must be called after the matrix assembling has been finished. 
+     * This makes it possible to start some cleanup or matrix data compressing 
+     * procedures.
+     */
     void finishAssembling();
  
     inline void addOperator(Operator* op, 
@@ -377,8 +391,10 @@ namespace AMDiS {
       return vec.end(); 
     }
   
-    /// Used by DOFAdmin to compress this DOFVector. Implementation of
-    /// DOFIndexedBase::compress()
+    /** \brief
+     * Used by DOFAdmin to compress this DOFVector. Implementation of
+     * DOFIndexedBase::compress()
+     */
     virtual void compressDOFIndexed(int first, int last,
 				    std::vector<DegreeOfFreedom> &newDof);
 
@@ -457,18 +473,22 @@ namespace AMDiS {
     /// Calculates Integral of this DOFVector
     double Int(Quadrature* q = NULL) const;
 
-    /// Calculates Integral of this DOFVector over parts of the domain
-    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
+    /** \brief
+     * Calculates Integral of this DOFVector over parts of the domain
+     * boundary, indicated by boundaryType. Implemented for DOFVector<double>
+     **/
     T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
     {
       FUNCNAME("DOFVector::IntOnBoundary())");
       TEST_EXIT(false)("Please implement your integration\n");
       return 0.0;
     }
-    
-    /// Calculates Integral of this DOFVector times normal vector over parts 
-    /// of the domain boundary, indicated by boundaryType. Implemented for 
-    /// DOFVector<WorldVector<double> >
+
+    /** \brief
+     * Calculates Integral of this DOFVector times normal vector over parts 
+     * of the domain boundary, indicated by boundaryType. Implemented for 
+     * DOFVector<WorldVector<double> >
+     **/
     double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
     {
       FUNCNAME("DOFVector::IntOnBoundaryNormal())");
@@ -561,25 +581,31 @@ namespace AMDiS {
     ///
     int calcMemoryUsage() const;
 
-    /// Computes the coefficients of the interpolant of the function fct and
-    /// stores these in the DOFVector
+    /** \brief
+     * Computes the coefficients of the interpolant of the function fct and
+     * stores these in the DOFVector
+     */
     void interpol(AbstractFunction<T, WorldVector<double> > *fct);
 
     void interpol(DOFVector<T> *v, double factor = 1.0);
 
 
-    /// Eval DOFVector at given point p. If oldElInfo != NULL the search for 
-    /// the element, where p is inside, starts from oldElInfo. implemented for: 
-    /// double, WorldVector<double>
-    inline const T evalAtPoint(WorldVector<double> &p, 
-			       ElInfo *oldElInfo = NULL) const 
+    /** eval DOFVector at given point p. If oldElInfo != NULL the search for the element, where p is inside,
+     *  starts from oldElInfo. 
+     *  implemented for: double, WorldVector< double >
+     */
+    inline const T& evalAtPoint(
+      WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const 
     {
       FUNCNAME("DOFVector::evalAtPoint())");
       TEST_EXIT(false)("Please implement your evaluation\n");
+      return *value;
     }
 
-    const bool getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo = NULL, bool useOldElInfo = false);
-
+    /** determine the DegreeOfFreedom that has coords with minimal euclidean distance to WorldVector p
+     *  return true if DOF is found, and false otherwise
+     */
+    const bool getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo = NULL, bool useOldElInfo = false) const;
 
     /// Writes the data of the DOFVector to an output stream.
     void serialize(std::ostream &out) 
@@ -598,11 +624,13 @@ namespace AMDiS {
       in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
     }
 
+//     DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
     DOFVector<typename GradientType<T>::type>*
       getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
 
     WorldVector<DOFVector<T>*> *getGradient(WorldVector<DOFVector<T>*> *grad) const;
 
+//     DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
     DOFVector<typename GradientType<T>::type>*
       getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
 
@@ -631,13 +659,12 @@ namespace AMDiS {
     BoundaryType boundaryType, Quadrature* q) const;
 
   template<>
-  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
-					      ElInfo *oldElInfo) const;
+  const double& DOFVector<double>::evalAtPoint(
+    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
 
   template<>
-  const WorldVector<double> 
-    DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p,
- 						 ElInfo *oldElInfo) const;
+  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
+    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
 
   template<>
   void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
@@ -664,8 +691,10 @@ namespace AMDiS {
 		       public DOFContainer
   {
   public:  
-    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
-    /// as DOFContainer at DOFAdmin
+    /** \brief
+     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
+     * as DOFContainer at DOFAdmin
+     */
     DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
       : DOFVector<DegreeOfFreedom>(feSpace_, name_)
     {
@@ -678,8 +707,10 @@ namespace AMDiS {
       feSpace->getAdmin()->removeDOFContainer(this);
     }
 
-    /// Implements DOFContainer::operator[]() by calling 
-    /// DOFVector<DegreeOfFreedom>::operator[]()
+    /** \brief
+     * Implements DOFContainer::operator[]() by calling 
+     * DOFVector<DegreeOfFreedom>::operator[]()
+     */
     DegreeOfFreedom& operator[](DegreeOfFreedom i) 
     {
       return DOFVector<DegreeOfFreedom>::operator[](i);
@@ -843,22 +874,20 @@ namespace AMDiS {
 					     std::vector<DOFVector<double>*> *res);
   
   /** \brief
-   * Computes the integral of a function that includes two different DOFVectors.
-   * This function works also for the case that the DOFVectors are defined on
-   * two different meshes.
+   * Computes the integral of a function that includes two different DOFVectors. This
+   * function works also for the case that the DOFVectors are defined on two different
+   * meshes.
    */
   double integrate(const DOFVector<double> &vec1,
 		   const DOFVector<double> &vec2,
 		   BinaryAbstractFunction<double, double, double> *fct);
 
-  /// Computes the integral of a function with two DOFVectors defined on the 
-  /// same mesh.
+  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
   double intSingleMesh(const DOFVector<double> &vec1,
 		       const DOFVector<double> &vec2,
 		       BinaryAbstractFunction<double, double, double> *fct);
 
-  /// Computes the integral of a function with two DOFVectors defined on 
-  /// different meshes.
+  /// Computes the integral of a function with two DOFVectors defined on different meshes.
   double intDualMesh(const DOFVector<double> &vec1,
 		     const DOFVector<double> &vec2,
 		     BinaryAbstractFunction<double, double, double> *fct);
@@ -874,7 +903,7 @@ namespace AMDiS {
   
   double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
 				  const std::vector<DOFVector<double>*> &grds,
-				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double