general integrate function and getDOFidxAtPoint implemented

AMDiS/src/DOFVector.cc

@@ -892,6 +892,112 @@ namespace AMDiS {
       elInfo = stack.traverseNext(elInfo);
     }
 
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    double localValue = value;
+    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
+#endif
+
+    return value;
+  }
+
+
+
+  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
+		       AbstractFunction<double, std::vector<double> > *fct)
+  {
+    FUNCNAME("integrateGeneral()");
+
+    TEST_EXIT(fct)("No function defined!\n");
+    TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");
+
+    int deg = std::max(fct->getDegree(),
+		       2 * vecs[0].getFeSpace()->getBasisFcts()->getDegree());
+    Quadrature* quad =
+      Quadrature::provideQuadrature(vecs[0].getFeSpace()->getMesh()->getDim(), deg);
+    FastQuadrature *fastQuad =
+      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());
+    for (size_t i = 0; i < vecs.size(); i++)
+      qp[i].change_dim(fastQuad->getNumPoints());
+
+    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);
+    while (elInfo) {
+      for (size_t i = 0; i < vecs.size(); i++)
+	vecs[i].getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
+
+      double tmp = 0.0;
+      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
+	for (size_t i = 0; i < vecs.size(); i++)
+	  qp_local[i] = qp[i][iq];
+ 	tmp += fastQuad->getWeight(iq) * (*fct)(qp_local);
+      }
+      value += tmp * elInfo->getDet();
+
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    double localValue = value;
+    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
+#endif
+
+    return value;
+  }
+
+  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
+				  const std::vector<DOFVector<double>*> &grds,
+				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > *fct)
+  {
+    FUNCNAME("integrateGeneral()");
+
+    TEST_EXIT(fct)("No function defined!\n");
+    TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");
+    TEST_EXIT(grds.size()>0)("No DOFVectors for gradients provided!\n");
+
+    int deg = std::max(fct->getDegree(),
+		       2 * vecs[0].getFeSpace()->getBasisFcts()->getDegree());
+    Quadrature* quad =
+      Quadrature::provideQuadrature(vecs[0].getFeSpace()->getMesh()->getDim(), deg);
+    FastQuadrature *fastQuad =
+      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());
+    std::vector<double> qp_local(vecs.size());
+    std::vector<WorldVector<double> > grd_local(grds.size());
+    for (size_t i = 0; i < vecs.size(); i++)
+      qp[i].change_dim(fastQuad->getNumPoints());
+    for (size_t i = 0; i < grds.size(); i++)
+      qpGrd[i].change_dim(fastQuad->getNumPoints());
+
+    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);
+    while (elInfo) {
+      for (size_t i = 0; i < vecs.size(); 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]);
+
+      double tmp = 0.0;
+      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
+	for (size_t i = 0; i < vecs.size(); i++)
+	  qp_local[i] = qp[i][iq];
+	for (size_t i = 0; i < grds.size(); i++)
+	  grd_local[i] = qpGrd[i][iq];
+ 	tmp += fastQuad->getWeight(iq) * (*fct)(qp_local, grd_local);
+      }
+      value += tmp * elInfo->getDet();
+
+      elInfo = stack.traverseNext(elInfo);
+    }
+
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     double localValue = value;
     MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);

AMDiS/src/DOFVector.h

@@ -578,6 +578,9 @@ namespace AMDiS {
       TEST_EXIT(false)("Please implement your evaluation\n");
     }
 
+    const bool getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo = NULL, bool useOldElInfo = false);
+
+
     /// Writes the data of the DOFVector to an output stream.
     void serialize(std::ostream &out) 
     {
@@ -865,8 +868,13 @@ namespace AMDiS {
 
   double integrate(const FiniteElemSpace* feSpace,
 		   AbstractFunction<double, WorldVector<double> > *fct);
-
-
+  
+  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
+			  AbstractFunction<double, std::vector<double> > *fct);
+  
+  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
+				  const std::vector<DOFVector<double>*> &grds,
+				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > *fct);
 }
 
 #include "DOFVector.hh"

AMDiS/src/DOFVector.hh

@@ -676,6 +676,59 @@ namespace AMDiS {
   }
 
 
+
+  template<typename T>
+  bool DOFVector<T>::getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo, bool useOldElInfo)
+  { FUNCNAME("DOFVector::getDOFidxAtPoint()");
+
+    Mesh *mesh = feSpace->getMesh();
+    const BasisFunction *basFcts = feSpace->getBasisFcts();
+
+    int dim = mesh->getDim();
+    int numBasFcts = basFcts->getNumber();
+
+    std::vector<DegreeOfFreedom> localIndices(numBasFcts);
+    DimVec<double> lambda(dim, NO_INIT);
+
+    ElInfo *elInfo = mesh->createNewElInfo();
+    idx = 0;
+    inside = false;
+
+    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
+      delete oldElInfo;
+    } else {
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
+    }
+    
+    if (oldElInfo)
+      oldElInfo = elInfo;
+    
+    if (!inside)
+      return false;
+    
+    basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+    FixVec<WorldVector<double>, VERTEX> coords = elInfo->getCoords();
+    int minIdx = -1;
+    double minDist = 1.e15;
+    
+    for (int i = 0; i < coords.size(); i++) {
+      WorldVector<double> dist = coords[i] - p;
+      double newDist = norm(dist);
+      if (newDist < minDist) {
+	minIdx = i;
+	minDist = newDist;
+      }
+    }
+    
+    if (minIdx >= 0)
+      idx = localIndices[minIdx];
+    
+    if(!oldElInfo) delete elInfo;
+    return inside;
+  };
+
+
   template<typename T>
   void DOFVector<T>::compressDOFIndexed(int first, int last, 
 					std::vector<DegreeOfFreedom> &newDOF)

AMDiS/src/Initfile.h

@@ -168,18 +168,18 @@ namespace AMDiS {
       // accepted brackets and delimiters for vector input
       std::string begBrackets= "{[(";
       std::string endBrackets= "}])";
-      std::string delims= ",;";
+      std::string delims= ";,";
 
       c.clear();
-      std::string val = trim(val_); 
+      std::string val = trim(val_);
+      bool hasBrackets = true;
       size_t pos = begBrackets.find(val[0]);
       if (pos == std::string::npos)
-        throw WrongVectorFormat("cannot convert "
-            "'" + val + "' into a list. No leading bracket found!");
-      if (val[val.length() - 1] != endBrackets[pos])
+        hasBrackets = false;
+      if (hasBrackets && val[val.length() - 1] != endBrackets[pos])
         throw WrongVectorFormat("begin and end bracket are different in"
             " value '" + val + "'");
-      size_t oldPos = 1;
+      size_t oldPos = (hasBrackets ? 1 : 0);
       size_t curDelim = 0;
       typedef typename Container::value_type ValueType;
       ValueType swap;
@@ -194,11 +194,11 @@ namespace AMDiS {
           pos = val.find(delims[curDelim], oldPos);
         }
         //last entry
-        std::string curWord = val.substr(oldPos, val.length() - 1 - oldPos);
+        std::string curWord = val.substr(oldPos, val.length() - (hasBrackets ? 1 : 0) - oldPos);
         convert(curWord, swap);
         c.push_back(swap);
       } catch (NoDelim nd) {
-        std::string curWord = val.substr(1, val.length() - 2);
+        std::string curWord = val.substr(oldPos, val.length() - (hasBrackets ? 2 : 0));
         curWord = trim(curWord);
         if (curWord.length() > 0) {
           // container with one entry

AMDiS/src/SystemVector.h

@@ -49,13 +49,13 @@ namespace AMDiS {
 
     /// Copy Constructor.
     SystemVector(const SystemVector& rhs)
-      : name(rhs.name),
-	feSpace(rhs.feSpace),
-	vectors(rhs.vectors.getSize())
+      : name(rhs.getName()),
+	feSpace(rhs.getFeSpaces()),
+	vectors(rhs.getNumVectors())
       
     {
       for (int i = 0; i < vectors.getSize(); i++)
-	vectors[i] = new DOFVector<double>(*rhs.vectors[i]);
+	vectors[i] = new DOFVector<double>(*rhs.getDOFVector(i));
     }
 
     ~SystemVector() 
@@ -91,6 +91,11 @@ namespace AMDiS {
       return vectors[index];
     }
 
+    inline std::string getName() const
+    {
+      return name;
+    }
+
     /// Returns sum of used vector sizes.
     inline int getUsedSize() const 
     {