corrections in VtkVectorWriter

AMDiS/src/CouplingProblemStat.h

@@ -84,11 +84,6 @@ namespace AMDiS {
 
       // create Meshes and FeSpaces
 
-      // one global macro-mesh for all problems
-      string meshName("");
-      Parameters::get(name + "->mesh", meshName);
-      TEST_EXIT(meshName != "")("No mesh name specified for \"%s->mesh\"!\n", 
-                                name.c_str());
 
       // all problems must have the same dimension (?)
       int dim = 0;
@@ -96,59 +91,82 @@ namespace AMDiS {
       TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
                     name.c_str());
 
-      map<int, Mesh*> meshForRefinementSet;
-      vector< std::set<Mesh*> > meshesForProblems(problems.size());
-      map<pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
+      std::map<std::string, Mesh*> meshByName;
+      std::vector< std::set<Mesh*> > meshesForProblems(problems.size());
+      std::map<std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
       for (size_t i = 0; i < problems.size(); ++i) {
         TEST_EXIT(problems[i])("problem[%d] does not exist!\n",i);
         nComponents += problems[i]->getNumComponents();
         for (size_t j = 0; j < problems[i]->getNumComponents(); j++) {
           // mesh
-          int refSet = -1;
-          Parameters::get(problems[i]->getName() + "->" +
-                          "refinement set[" + boost::lexical_cast<string>(j) + "]", 
-                          refSet);
-          if (refSet < 0)
-            refSet = 0;
-
-          if (meshForRefinementSet.find(refSet) == meshForRefinementSet.end()) {
+
+	  string meshName("");
+	  Parameters::get(problems[i]->getName() + "->mesh", meshName);
+	  TEST_EXIT(meshName != "")("No mesh name specified for \"%s->mesh\"!\n",
+				    problems[i]->getName().c_str());
+
+          if (meshByName.find(meshName) == meshByName.end()) {
             Mesh *newMesh = new Mesh(meshName, dim);
-            meshForRefinementSet[refSet] = newMesh;
+            meshByName[meshName] = newMesh;
             meshes.push_back(newMesh);
             meshesForProblems[i].insert(newMesh);
             nMeshes++;
           } else
-            meshesForProblems[i].insert(meshForRefinementSet[refSet]);
+            meshesForProblems[i].insert(meshByName[meshName]);
 
-          problems[i]->setComponentMesh(j, meshForRefinementSet[refSet]);
+          problems[i]->setComponentMesh(j, meshByName[meshName]);
             
           // feSpace
           int degree = 1;
           Parameters::get(problems[i]->getName() + "->polynomial degree[" + 
                           boost::lexical_cast<string>(j) + "]", degree);
 
-          if (feSpaceMap[pair<Mesh*, int>(meshForRefinementSet[refSet], degree)] == NULL) {
+          if (feSpaceMap[pair<Mesh*, int>(meshByName[meshName], degree)] == NULL) {
             stringstream s;
             s << problems[i]->getName() << "->feSpace[" << j << "]";
     
             FiniteElemSpace *newFeSpace = 
               FiniteElemSpace::provideFeSpace(NULL, Lagrange::getLagrange(dim, degree),
-                                              meshForRefinementSet[refSet], s.str());
-            feSpaceMap[pair<Mesh*, int>(meshForRefinementSet[refSet], degree)] = newFeSpace;
+                                              meshByName[meshName], s.str());
+            feSpaceMap[pair<Mesh*, int>(meshByName[meshName], degree)] = newFeSpace;
             feSpaces.push_back(newFeSpace);
           }
     
-//           problems[i].setComponentSpace(i, feSpaceMap[pair<Mesh*, int>(meshForRefinementSet[refSet], degree)]);
+//           problems[i]->setComponentSpace(j, feSpaceMap[pair<Mesh*, int>(meshByName[meshName], degree)]);
+
         }
       }
 
-      for (size_t i = 0; i < problems.size(); ++i) {
+      for (size_t i = 0; i < problems.size(); i++) {
         vector<Mesh*> problemMeshes(meshesForProblems[i].begin(), meshesForProblems[i].end());
         problems[i]->setMeshes(problemMeshes);
         problems[i]->setRefinementManager(refinementManager);
         problems[i]->setCoarseningManager(coarseningManager);
-        problems[i]->initialize(INIT_ALL);
+        problems[i]->initialize(INIT_ALL - INIT_MESH);
+      }
+
+      for (size_t i = 0; i < meshes.size(); i++) {
+	int globalRefinements = 0;
+
+	// If AMDiS is compiled for parallel computations, the global refinements are
+	// ignored here. Later, each rank will add the global refinements to its
+	// private mesh.
+#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
+	Parameters::get(meshes[i]->getName() + "->global refinements",
+			globalRefinements);
+#endif
+
+	bool initMesh = initFlag.isSet(INIT_MESH);
+
+	// Initialize the meshes if there is no serialization file.
+	if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
+	  meshes[i]->initialize();
+
+	// do global refinements
+	if (initMesh && meshes[i])
+	  refinementManager->globalRefine(meshes[i], globalRefinements);
       }
+      
     }
 
     void createRefCoarseManager() 

AMDiS/src/DOFVector.h

@@ -904,6 +904,10 @@ 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);
+
+  template<typename T>
+  T integrate_VecAndCoords(const DOFVector<double> &vec,
+			   BinaryAbstractFunction<T, double, WorldVector<double> > *fct);
 }
 
 #include "DOFVector.hh"

AMDiS/src/DOFVector.hh

@@ -541,7 +541,51 @@ namespace AMDiS {
     return result;
   }
 
+  template<typename T>
+  T integrate_VecAndCoords(const DOFVector<double> &vec,
+			   BinaryAbstractFunction<T, double, WorldVector<double> > *fct)
+  {
+    FUNCNAME("integrate_VecAndCoords()");
+
+    TEST_EXIT(fct)("No function defined!\n");
+
+    int deg = std::max(fct->getDegree(),
+		       2 * vec.getFeSpace()->getBasisFcts()->getDegree());
+    Quadrature* quad =
+      Quadrature::provideQuadrature(vec.getFeSpace()->getMesh()->getDim(), deg);
+    FastQuadrature *fastQuad =
+      FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
+
+    mtl::dense_vector<double> qp(fastQuad->getNumPoints());
+    WorldVector<double> coordsAtQP;
+
+    T value;
+    nullify(value);
+    
+    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
+    while (elInfo) {
+      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
+      T tmp;
+      nullify(tmp);
+      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
+	elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
+ 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
+      }
+      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;
+  }
+  
   template<typename T>
   double DOFVector<T>::L1Norm(Quadrature* q) const
   {

AMDiS/src/ProblemStat.h

@@ -433,6 +433,16 @@ namespace AMDiS {
       feSpaces[comp] = feSpace;
     }
 
+    void setComponentSpace(int comp, FiniteElemSpace *feSpace)
+    {
+      if (static_cast<int>(componentSpaces.size()) < nComponents)
+        componentSpaces.resize(nComponents);
+      TEST_EXIT(comp >= 0 && comp < nComponents)
+        ("Component number not in feasable range!");
+
+      componentSpaces[comp] = feSpace;
+    }
+    
     /// Sets \ref estimator.
     inline void setEstimator(Estimator* est, int comp = 0) 
     { 

AMDiS/src/io/FileWriter.cc

@@ -140,6 +140,7 @@ namespace AMDiS {
     writeAMDiSFormat = 0;
     writeParaViewFormat = 0;
     writeParaViewVectorFormat = 0;
+    writeAs3dVector = false;
     writeParaViewAnimation = 0;
     writePeriodicFormat = 0;
     writePngFormat = 0;
@@ -169,6 +170,7 @@ namespace AMDiS {
     Parameters::get(name + "->AMDiS data ext", amdisDataExt);
     Parameters::get(name + "->ParaView format", writeParaViewFormat);
     Parameters::get(name + "->ParaView vector format", writeParaViewVectorFormat);
+    Parameters::get(name + "->write vector as 3d vector", writeAs3dVector);
     Parameters::get(name + "->ParaView animation", writeParaViewAnimation);
     Parameters::get(name + "->ParaView ext", paraviewFileExt);    
     Parameters::get(name + "->Periodic format", writePeriodicFormat);
@@ -288,12 +290,7 @@ namespace AMDiS {
     }
 
     if (writeParaViewVectorFormat) {
-      VtkVectorWriter::writeFile(solutionVecs, fn + paraviewFileExt);
-
-#if HAVE_PARALLEL_DOMAIN_AMDIS
-      if (MPI::COMM_WORLD.Get_rank() == 0)
-	VtkVectorWriter::writeFile(solutionVecs, fn + paraviewFileExt, true);
-#endif
+      VtkVectorWriter::writeFile(solutionVecs, fn + paraviewFileExt, true, writeAs3dVector);
 
       MSG("ParaView file written to %s\n", (fn + paraviewFileExt).c_str());
     }

AMDiS/src/io/FileWriter.h

@@ -170,6 +170,9 @@ namespace AMDiS {
     /// 0: Don't write ParaView vector files; 1: Write ParaView vector files.
     int writeParaViewVectorFormat;
 
+    /// 1: extend number of component to 3, so that paraview can display the vector as worldvector
+    bool writeAs3dVector;
+
     /// 0: Don't write ParaView animation file; 1: Write ParaView animation file.
     int writeParaViewAnimation;
 

AMDiS/src/io/VtkVectorWriter.cc

@@ -75,7 +75,8 @@ namespace AMDiS {
 
   void VtkVectorWriter::writeFile(std::vector<DOFVector<double>* > &values,
 			    std::string filename,
-			    bool writeParallel)
+			    bool writeParallel,
+			    bool writeAs3dVector)
   {
     DOFVector<std::vector<double> > *newValues = new DOFVector<std::vector<double> >(values[0]->getFeSpace(), "values");
     std::vector<DOFIterator<double>* > iterators;
@@ -88,13 +89,16 @@ namespace AMDiS {
     for(resultIter.reset(); !resultIter.end(); resultIter++)
     {
       std::vector<double> val;
-      for (size_t i = 0; i < iterators.size(); i++)
+      for (size_t i = 0; i < static_cast<size_t>(iterators.size()); i++)
 	val.push_back(*(*(iterators[i])));
 
       *resultIter = val;
+      
+      for (size_t i = 0; i < static_cast<size_t>(iterators.size()); i++)
+	(*(iterators[i]))++;
     }
 
-    writeFile(newValues, filename, writeParallel);
+    writeFile(newValues, filename, writeParallel, writeAs3dVector);
     for (size_t i = 0; i < iterators.size(); i++)
       delete iterators[i];
     delete newValues;
@@ -103,7 +107,8 @@ namespace AMDiS {
 
   void VtkVectorWriter::writeFile(WorldVector<DOFVector<double>* > &values,
 			    std::string filename,
-			    bool writeParallel)
+			    bool writeParallel,
+			    bool writeAs3dVector)
   {
     DOFVector<WorldVector<double> > *newValues = new DOFVector<WorldVector<double> >(values[0]->getFeSpace(), "values");
     WorldVector<DOFIterator<double>* > iterators;
@@ -117,9 +122,12 @@ namespace AMDiS {
     {
       for (size_t i = 0; i < static_cast<size_t>(iterators.getSize()); i++)
 	(*resultIter)[i] = *(*(iterators[i]));
+
+      for (size_t i = 0; i < static_cast<size_t>(iterators.size()); i++)
+	(*(iterators[i]))++;
     }
 
-    writeFile(newValues, filename, writeParallel);
+    writeFile(newValues, filename, writeParallel, writeAs3dVector);
     for (size_t i = 0; i < static_cast<size_t>(iterators.getSize()); i++)
       delete iterators[i];
     delete newValues;
@@ -128,7 +136,8 @@ namespace AMDiS {
 
   void VtkVectorWriter::writeFile(SystemVector *values,
 			    std::string filename,
-			    bool writeParallel)
+			    bool writeParallel,
+			    bool writeAs3dVector)
   {
     DOFVector<std::vector<double> > *newValues = new DOFVector<std::vector<double> >(values->getDOFVector(0)->getFeSpace(), "values");
     std::vector<DOFIterator<double>* > iterators;
@@ -141,13 +150,16 @@ namespace AMDiS {
     for(resultIter.reset(); !resultIter.end(); resultIter++)
     {
       std::vector<double> val;
-      for (size_t i = 0; i < iterators.size(); i++)
+      for (size_t i = 0; i < static_cast<size_t>(iterators.size()); i++)
 	val.push_back(*(*(iterators[i])));
 
       *resultIter = val;
+
+      for (size_t i = 0; i < static_cast<size_t>(iterators.size()); i++)
+	(*(iterators[i]))++;
     }
 
-    writeFile(newValues, filename, writeParallel);
+    writeFile(newValues, filename, writeParallel, writeAs3dVector);
     for (size_t i = 0; i < iterators.size(); i++)
       delete iterators[i];
     delete newValues;

AMDiS/src/io/VtkVectorWriter.h

@@ -41,10 +41,10 @@ namespace AMDiS {
   {
   template<typename S>  
   struct Impl {
-     Impl(std::vector<DataCollector<S>*> *dc)
+     Impl(std::vector<DataCollector<S>*> *dc, bool writeAs3dVector_=false)
       : dataCollector(dc),
 	compress(NONE),
-	writeAs3dVector(false),
+	writeAs3dVector(writeAs3dVector_),
 	three(3)
     {
       degree = (*dataCollector)[0]->getFeSpace()->getBasisFcts()->getDegree();
@@ -131,6 +131,10 @@ namespace AMDiS {
 	o << (fabs(*it) < 1e-40 ? 0.0 : *it) << (i < l.size() - 1 ? " " : "");
 	++it;
       }
+      if (writeAs3dVector) {
+	for (int i = l.size(); i < 3; i++)
+	  o << 0.0 << " ";
+      }
       return o;
     }
 
@@ -198,41 +202,48 @@ namespace AMDiS {
     template<typename T>
     static void writeFile(DOFVector<T> *values,
 			  std::string filename, 
-			  bool writeParallel = true)
+			  bool writeParallel = true,
+			  bool writeAs3dVector = false)
     {
       std::vector<DataCollector<T>*> dcList(0);
       dcList.push_back(new DataCollector<T>(values->getFeSpace(), values));
-      writeFile(dcList, filename, writeParallel);
+      writeFile(dcList, filename, writeParallel, writeAs3dVector);
       delete dcList[0];
     }
 
     /// May be used to simply write ParaView files.
-    static void writeFile(DOFVector<double> &values, 
+    template<typename T>
+    static void writeFile(DOFVector<T> &values,
 			  std::string filename,
-			  bool writeParallel = true)
+			  bool writeParallel = true,
+			  bool writeAs3dVector = false)
     {
-      writeFile(&values, filename, writeParallel);
+      writeFile(&values, filename, writeParallel, writeAs3dVector);
     }
 
     /// May be used to simply write ParaView files with a list of values.
     static void writeFile(std::vector<DOFVector<double>*> &values,
 			  std::string filename, 
-			  bool writeParallel = true);
+			  bool writeParallel = true,
+			  bool writeAs3dVector = false);
 
     static void writeFile(WorldVector<DOFVector<double>*> &values,
 			  std::string filename, 
-			  bool writeParallel = true);
+			  bool writeParallel = true,
+			  bool writeAs3dVector = false);
 
     static void writeFile(SystemVector *values,
 			  std::string filename, 
-			  bool writeParallel = true);
+			  bool writeParallel = true,
+			  bool writeAs3dVector = false);
 
     template<typename T>
     static void writeFile(std::vector<DataCollector<T>*> &dcList,
 			  std::string filename,
-			  bool writeParallel = true)
+			  bool writeParallel = true,
+			  bool writeAs3dVector = false)
     {
-      Impl<T> writer(&dcList);
+      Impl<T> writer(&dcList, writeAs3dVector);
 
   #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
       if (writeParallel) {

AMDiS/src/io/VtkVectorWriter.hh

@@ -185,7 +185,8 @@ namespace AMDiS {
     for (int i = 0; i < static_cast<int>(dataCollector->size()); i++) {
       S temp = (*((*dataCollector)[i]->getValues()))[0];
       int numComponent = num_rows(temp);
-      file << "        <DataArray type=\"Float32\" Name=\"value" << i 
+      file << "        <DataArray type=\"Float32\" Name=\""
+	   << (*dataCollector)[i]->getValues()->getName()
 	   << "\" format=\"ascii\"" << (numComponent > 1 ? " NumberOfComponents=\"" + boost::lexical_cast<std::string>(numComponent) + "\"" : "") << ">\n";
 
       writeVertexValues(file, i);