Vector VTK-Writer added

AMDiS/src/DOFIterator.h

@@ -312,7 +312,199 @@ namespace AMDiS {
     /// Iterator for \ref iteratedObject
     typename std::vector<T>::iterator it;
   };
+  
+  
+
+  /** \ingroup DOFAdministration
+    * \brief
+    * Implements a DOFIterator for a const DOFIndexed<T> object
+    */
+  template<typename T>
+  class DOFConstIterator : public DOFIteratorBase
+  {
+  public:
+    /// Constructs a DOFIterator for cont of type t
+    DOFConstIterator(const DOFIndexed<T> *obj, DOFIteratorType t)
+      : DOFIteratorBase(dynamic_cast<DOFAdmin*>(obj->getFeSpace()->getAdmin()), t),
+	iteratedObject(obj)
+    {}
+
+    /// Constructs a DOFIterator for cont of type t
+    DOFConstIterator(DOFAdmin *admin,
+		const DOFIndexed<T> *obj,
+		DOFIteratorType t)
+      : DOFIteratorBase(admin, t),
+	iteratedObject(obj)
+    {}
+
+    /// Dereference operator
+    inline const T& operator*()
+    {
+      return *it;
+    }
+
+    /// Dereference operator
+    inline const T* operator->()
+    {
+      return &(*it);
+    }
+
+    inline bool operator!=(const DOFIterator<T>& rhs)
+    {
+      if (this->iteratedObject != rhs.iteratedObject)
+	return true;
+
+      if (this->it != rhs.it)
+	return true;
+
+      return false;
+    }
+
+    inline bool operator==(const DOFIterator<T>& rhs)
+    {
+      return !(this->operator==(rhs));
+    }
 
+  protected:
+    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
+    inline void goToBeginOfIteratedObject()
+    {
+      it = iteratedObject->begin();
+    }
+
+    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
+    inline void goToEndOfIteratedObject()
+    {
+      it = iteratedObject->end();
+    }
+
+    /// Implementation of DOFIteratorBase::incObjectIterator()
+    inline void incObjectIterator()
+    {
+      ++it;
+    }
+
+    /// Implementation of DOFIteratorBase::incObjectIterator()
+    inline void decObjectIterator()
+    {
+      --it;
+    }
+
+  protected:
+    /// Object that is iterated
+    const DOFIndexed<T> *iteratedObject;
+
+    /// Iterator for \ref iteratedObject
+    typename std::vector<T>::const_iterator it;
+  };
+
+
+
+  /** \ingroup DOFAdministration
+    * \brief
+    * Implements a DOFIterator for a vector of DOFVector<T> objects
+    */
+  template<typename T>
+  class DOFVectorIterator : public DOFIteratorBase
+  {
+  public:
+    typedef typename std::vector<T>::iterator VectorIterator;
+    
+    /// Constructs a DOFIterator for cont of type t
+    DOFVectorIterator(std::vector<DOFVector<T>*> &obj, DOFIteratorType t) 
+      : DOFIteratorBase(dynamic_cast<DOFAdmin*>(obj[0]->getFeSpace()->getAdmin()), t),
+	iteratedObject(obj)
+    {
+      for (size_t i = 0; i < obj.size(); i++)
+	it.push_back(new VectorIterator);
+    }
+
+    /// Constructs a DOFIterator for cont of type t
+    DOFVectorIterator(DOFAdmin *admin,
+		std::vector<DOFVector<T>*> &obj, 
+		DOFIteratorType t) 
+      : DOFIteratorBase(admin, t),
+	iteratedObject(obj)
+    {
+      for (size_t i = 0; i < obj.size(); i++)
+	it.push_back(new VectorIterator);
+    }
+    
+    ~DOFVectorIterator()
+    {
+      for (size_t i = 0; i < it.size(); i++)
+	delete it[i];
+    }
+
+    /// Dereference operator
+    inline std::vector<T> operator*() 
+    { 
+      std::vector<T> result(it.size());
+      for (size_t i = 0; i < it.size(); i++)
+	result[i] = *(*it[i]);
+      return result;
+    }
+
+    /// Dereference operator
+    inline T* operator->() 
+    { 
+      throw std::runtime_error("operator-> not available for DOFVectorIterator!");
+      return &(*(*it[0])); 
+    }
+
+    inline bool operator!=(const DOFVectorIterator<T>& rhs) 
+    {
+      if (this->iteratedObject != rhs.iteratedObject) 
+	return true;
+
+      if (this->it != rhs.it) 
+	return true;
+
+      return false;
+    }
+
+    inline bool operator==(const DOFVectorIterator<T>& rhs) 
+    {
+      return !(this->operator==(rhs));
+    }
+
+  protected:
+    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
+    inline void goToBeginOfIteratedObject() 
+    { 
+      for (size_t i = 0; i < it.size(); i++)
+	*(it[i]) = iteratedObject[i]->begin(); 
+    }
+
+    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
+    inline void goToEndOfIteratedObject() 
+    { 
+      for (size_t i = 0; i < it.size(); i++)
+	*(it[i]) = iteratedObject[i]->end(); 
+    }
+
+    /// Implementation of DOFIteratorBase::incObjectIterator()
+    inline void incObjectIterator() 
+    { 
+      for (size_t i = 0; i < it.size(); i++)
+	++(*(it[i])); 
+    }
+
+    /// Implementation of DOFIteratorBase::incObjectIterator()
+    inline void decObjectIterator() 
+    { 
+      for (size_t i = 0; i < it.size(); i++)
+	--(*(it[i]));
+    }
+
+  protected:
+    /// Object that is iterated
+    std::vector<DOFVector<T>*> iteratedObject;
+
+    /// Iterator for \ref iteratedObject
+    std::vector<typename std::vector<T>::iterator*> it;
+  };
+  
 }
 
 #endif // AMDIS_DOFITERATOR_H

AMDiS/src/Debug.cc

@@ -69,6 +69,7 @@ namespace AMDiS {
       if (rank == -1 || myRank == rank) {
 	DOFVector<double> tmp(feSpace, "tmp");
 	io::VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu", 
+				 false,
 			     false);
       }
     }

AMDiS/src/ProblemStat.cc

@@ -639,6 +639,39 @@ namespace AMDiS {
 					     componentMeshes[i], 
 					     solution->getDOFVector(i)));      
     }
+    
+    // create a filewrite for groups of components to write vector-values output
+    int nVectors = 0;
+    Parameters::get(name + "->output->num vectors", nVectors);
+    if (nVectors > 0) {
+      for (int j = 0; j < nVectors; j++) {
+	numberedName = name + "->output->vector[" + boost::lexical_cast<string>(j) + "]";
+	
+	filename = "";
+	Parameters::get(numberedName + "->filename", filename);
+	std::string componentName = "";
+	Parameters::get(numberedName + "->name", componentName);
+	std::vector<std::string> names; 
+	if (name != "")
+	  names.push_back(componentName);
+	std::vector<int> comp;
+	Parameters::get(numberedName + "->components", comp);
+	
+	if (filename != "" && comp.size() > 0) {
+	  // Create own filewriters for each component of the problem
+	  std::vector<DOFVector<double>*> vectors;
+	  for (size_t i = 0; i < comp.size(); i++)
+	    vectors.push_back(solution->getDOFVector(comp[i]));
+	  
+	  fileWriters.push_back(new FileWriter(numberedName, 
+					      componentMeshes[comp[0]], 
+					      vectors,
+					      names
+					      ));      
+	  
+	}
+      }
+    }
 
     int writeSerialization = 0;
     Parameters::get(name + "->output->write serialization", writeSerialization);

AMDiS/src/io/FileWriter.cc

@@ -67,6 +67,9 @@ namespace AMDiS
       solutionVecs.resize(vecs->getSize());
       for (int i = 0; i < static_cast<int>(vecs->getSize()); i++)
 	solutionVecs[i] = vecs->getDOFVector(i);
+      
+      for (size_t i = 0; i < solutionVecs.size(); i++)
+	solutionNames.push_back(solutionVecs[i]->getName());
     }
 
 
@@ -171,22 +174,25 @@ namespace AMDiS
       }
 
       if (writeParaViewFormat) {
-	VtkWriter::Aux vtkWriter(&dataCollectors, VtkWriter::Vtuformat(paraViewMode), (paraViewPrecision == 1));
+	VtkWriter::Aux vtkWriter(&dataCollectors, 
+				 solutionNames,
+				 VtkWriter::Vtuformat(paraViewMode), (paraViewPrecision == 1), writeParaViewVectorFormat);
 	vtkWriter.writeFile(fn + paraviewFileExt);
 
   #if HAVE_PARALLEL_DOMAIN_AMDIS
 	if (MPI::COMM_WORLD.Get_rank() == 0) {
-	  vector<string> componentNames;
-	  for (unsigned int i = 0; i < dataCollectors.size(); i++)
-	    componentNames.push_back(dataCollectors[i]->getValues()->getName());
+// 	  vector<string> componentNames;
+// 	  for (unsigned int i = 0; i < dataCollectors.size(); i++)
+// 	    componentNames.push_back(dataCollectors[i]->getValues()->getName());
 
 	  VtkWriter::detail::writeParallelFile(paraFilename + paraviewParallelFileExt,
 				      MPI::COMM_WORLD.Get_size(),
 				      filename, 
 				      postfix,
-				      componentNames,
+				      solutionNames,
 				      VtkWriter::Vtuformat(paraViewMode),
-				      (paraViewPrecision == 1));
+				      (paraViewPrecision == 1),
+				      writeParaViewVectorFormat);
 	}
   #endif
 
@@ -194,7 +200,7 @@ namespace AMDiS
       }
 
       // write vtu-vector files
-      if (writeParaViewVectorFormat) {
+      if (writeParaViewVectorFormat && !writeParaViewFormat) {
 	VtkVectorWriter::writeFile(solutionVecs, fn_ + paraviewFileExt, true, writeAs3dVector);
 	MSG("ParaView file written to %s\n", (fn_ + paraviewFileExt).c_str());
       }

AMDiS/src/io/FileWriter.h

@@ -109,8 +109,10 @@ namespace AMDiS {
 
       /// Constructor for a filewriter with more than one data component.
       FileWriter(std::string name,
-			  Mesh *mesh,
-			  std::vector< DOFVector<T>* > vecs);
+		  Mesh *mesh,
+		  std::vector< DOFVector<T>* > vecs,
+		  std::vector< std::string > componentNames = std::vector< std::string >()
+		);
 
       /// Constructor for a filewriter with more than one data component.
       FileWriter(std::string name,
@@ -239,6 +241,9 @@ namespace AMDiS {
 
       /// Pointers to the std::vectors which store the solution.
       std::vector<DOFVector<T>*> solutionVecs;
+      
+      /// Names of the DOFVectors
+      std::vector<std::string> solutionNames;
 
       /** \brief
       * Stores the number of temporal solutions std::vectors, which have been created

AMDiS/src/io/FileWriter.hh

@@ -63,15 +63,19 @@ namespace AMDiS
 
       solutionVecs.resize(1);
       solutionVecs[0] = vec;
+      
+      solutionNames.push_back(vec->getName());
     }
 
 
     template<typename T>
     FileWriter<T>::FileWriter(std::string name_,
 			  Mesh *mesh_, 
-			  std::vector<DOFVector<T>*> vecs)
+			  std::vector<DOFVector<T>*> vecs,
+			  std::vector< std::string > componentNames)
       : name(name_),
-	mesh(mesh_)
+	mesh(mesh_),
+	solutionNames(componentNames)
     {
       initialize();
 
@@ -82,6 +86,11 @@ namespace AMDiS
   */
       feSpace = vecs[0]->getFeSpace();
       solutionVecs = vecs;
+      
+      if (solutionNames.size() < vecs.size()) {
+	for (size_t i = solutionNames.size(); i < vecs.size(); i++)
+	  solutionNames.push_back(vecs[i]->getName());
+      }
     }
 
 

AMDiS/src/io/VtkWriter.cc

@@ -139,7 +139,10 @@ namespace AMDiS
 		     bool writeParallel
 		    )
       {
-	::AMDiS::io::VtkWriter::Aux writer(&dcList, format, highPrecision);
+	vector<string> componentNames;
+	for (size_t i = 0; i < dcList.size(); i++)
+	  componentNames.push_back(dcList[i]->getValues()->getName());
+	::AMDiS::io::VtkWriter::Aux writer(&dcList, componentNames, format, highPrecision);
 	
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
 	if (writeParallel) {
@@ -150,9 +153,6 @@ namespace AMDiS
 	  string name = filename.substr(0, sPos);
 	  
 	  if (MPI::COMM_WORLD.Get_rank() == 0) {
-	    vector<string> componentNames;
-	    for (size_t i = 0; i < dcList.size(); i++)
-	      componentNames.push_back(dcList[i]->getValues()->getName());
 
 	    detail::writeParallelFile(name + ".pvtu", MPI::COMM_WORLD.Get_size(), 
 				      name, ".vtu", componentNames, format, highPrecision);

AMDiS/src/io/detail/VtkWriter.cc

@@ -175,7 +175,9 @@ namespace AMDiS  { namespace io {
 			      string fnPrefix, string fnPostfix,
 			      vector<string> &componentNames,
 			      ::AMDiS::io::VtkWriter::Vtuformat format,
-			      bool highPrecision)
+			      bool highPrecision,
+			      bool writeAsVector
+			      )
 	{
 	  FUNCNAME("writeParallelFile()");
 
@@ -222,8 +224,8 @@ namespace AMDiS  { namespace io {
 	      << "      <PDataArray type=\"Int32\" Name=\"connectivity\"/>\n"
 	      << "    </PCells>\n";
 	  file << "    <PPointData>\n";
-
-	  for (size_t i = 0; i < componentNames.size(); i++) {
+	  
+	  for (size_t i = 0; i < writeAsVector ? 1 : componentNames.size(); i++) {
 	    if(highPrecision && format != ::AMDiS::io::VtkWriter::ASCII)
 	      file << "      <PDataArray type=\"Float64\" Name=\"";
 	    else 
@@ -232,11 +234,15 @@ namespace AMDiS  { namespace io {
 	    file << componentNames[i];
 	    
 	    if(format == ::AMDiS::io::VtkWriter::ASCII)
-	      file << "\" format=\"ascii\"/>\n";
+	      file << "\" format=\"ascii\"";
 	    else
-	      file << "\" format=\"appended\"/>\n"; 
+	      file << "\" format=\"appended\"";
+	    if (writeAsVector && componentNames.size() > 1)
+	      file << " NumberOfComponents=\"" << std::max(3,static_cast<int>(componentNames.size())) << "\"";
+	    file << "/>\n"; 
 	  }
 	  
+	  
 	  file << "    </PPointData>\n";
 
 	  for (int i = 0; i < nRanks; i++) {

AMDiS/src/io/detail/VtkWriter.h

@@ -109,9 +109,15 @@ namespace AMDiS  { namespace io {
 	    bool highPrecision;
 	};
 
-	Aux(std::vector<DataCollector<>*> *dc, Vtuformat f = ASCII, bool hp = false)
+	Aux(std::vector<DataCollector<>*> *dc, 
+	    std::vector<std::string>& names_, 
+	    Vtuformat f = ASCII, 
+	    bool hp = false, 
+	    bool writeAsVector_ = false)
 	  : dataCollector(dc),
-	    bstream(hp)
+	    componentNames(names_),
+	    bstream(hp),
+	    writeAsVector(writeAsVector_)
 	{
 #ifndef HAVE_COMPRESSION
 	  FUNCNAME("VtkWriter::Aux::Aux()");
@@ -151,16 +157,50 @@ namespace AMDiS  { namespace io {
 	template<typename T>
 	void writeVertexValues(T &file, int componentNo);
 
+	/// Writes all values of vertices and interpolation point as vector to an output file.
+	template<typename T>
+	void writeVertexValues(T &file);
+
 	/// Writes the connectivity of all simplices to an output file.
 	template<typename OutputStream>
 	void writeConnectivity(OutputStream &file);
-
-      private:
 	
+	template<typename Stream>
+	Stream& print(const WorldVector<double> &l, Stream &o) 
+	{
+	  for (size_t i = 0; i < static_cast<size_t>(l.getSize()); i++) {
+	    o << (fabs(l[i]) < 1e-40 ? 0.0 : l[i]) << " ";
+	  }
+	  for (int i = l.getSize(); i < 3; i++)
+	    o << 0.0 << " ";
+	  return o;
+	}
+	
+	template<typename Stream>
+	Stream& print(const std::vector<double> &l, Stream &o) 
+	{
+	  for (size_t i = 0; i < l.size(); i++) {
+	    o << (fabs(l[i]) < 1e-40 ? 0.0 : l[i]) << " ";
+	  }
+	  for (int i = l.size(); i < 3; i++)
+	    o << 0.0 << " ";
+	  return o;
+	}
+
+	template<typename T, typename Stream>
+	Stream& print(const T &value, Stream &o) 
+	{
+	  o << (fabs(value) < 1e-40 ? 0.0 : value);
+	  return o;
+	}
+
+      private:	
 	std::string getStreamData();
       
 	/// List of DataCollectors, for each component of the problem one.
 	std::vector<DataCollector<>*> *dataCollector;
+	
+	std::vector<std::string> componentNames;
 
 	/// Degree of the basis function of the problem.
 	int degree;
@@ -173,6 +213,9 @@ namespace AMDiS  { namespace io {
 	
 	/// The binary stream for filtering data. Only used in append and appended_compressed mode.
 	BinaryStream bstream; 
+	
+	/// write data-vectors as vector-valued component
+	bool writeAsVector;
       };
       
       template <class T> 
@@ -249,7 +292,9 @@ namespace AMDiS  { namespace io {
 			      std::string fnPrefix, std::string fnPostfix,
 			      std::vector<std::string> &componentNames,
 			      ::AMDiS::io::VtkWriter::Vtuformat format = ::AMDiS::io::VtkWriter::ASCII,
-			      bool highPrecision = false);
+			      bool highPrecision = false,
+			      bool writeAsVector = false
+			      );
 
 	/// Writes the connectivity for the case dim = 2 and degree = 2 to an output file.
 	template<typename S, typename OutputStream>

AMDiS/src/io/detail/VtkWriter.hh

@@ -96,10 +96,15 @@ namespace AMDiS  { namespace io {
 	file << "      </Cells>\n";
 	file << "      <PointData>\n";
 	
-	for (int i = 0; i < static_cast<int>(dataCollector->size()); i++) {
+	int nValues = static_cast<int>(dataCollector->size());
+	nValues = writeAsVector ? std::min(nValues,1) : nValues;	
+	for (int i = 0; i < nValues; i++) {
 	  file << "        <DataArray type=\"Float32\" Name=\"" 
-	      << (*dataCollector)[i]->getValues()->getName() 
-	      << "\" format=\"ascii\">\n";
+	      << componentNames[i]
+	      << "\" format=\"ascii\"";
+	  if (writeAsVector && dataCollector->size() > 1)
+	    file << " NumberOfComponents=\"" << std::max(3,static_cast<int>(dataCollector->size())) << "\"";
+	  file << ">\n"; 
 
 	  writeVertexValues(file, i);
 	  
@@ -119,7 +124,8 @@ namespace AMDiS  { namespace io {
 	
 	bstream.str("");
 	
-	int valueSize = static_cast<int>(dataCollector->size());
+	int valueSize = writeAsVector ? 1 : static_cast<int>(dataCollector->size());
+	valueSize = std::min(valueSize, static_cast<int>(dataCollector->size()));
 	int nDataArray = 4 + valueSize;
 	vector<string> dataBase64Vec(nDataArray, "");
 	vector<size_t> dataOffsetVec(nDataArray, 0);
@@ -202,8 +208,11 @@ namespace AMDiS  { namespace io {
 	  else
 	    file << "        <DataArray type=\"Float32\" Name=\"";
 	  
-	  file << (*dataCollector)[i]->getValues()->getName() 
-	      << "\" format=\"appended\" offset=\"" << dataOffsetVec[i] << "\"/>\n";
+	  file << componentNames[i] 
+	      << "\" format=\"appended\" offset=\"" << dataOffsetVec[i] << "\"";
+	  if (writeAsVector && dataCollector->size() > 1)
+	    file << " NumberOfComponents=\"" << std::max(3,static_cast<int>(dataCollector->size())) << "\"";
+	  file << "/>\n"; 
 	}
 	
 	file << "      </PointData>\n";
@@ -268,6 +277,11 @@ namespace AMDiS  { namespace io {
       template<typename T>
       void Aux::writeVertexValues(T &file, int componentNo)
       {
+	if (writeAsVector) {
+	  writeVertexValues(file);
+	  return;
+	}
+	
 	DOFVector<int> *interpPointInd = (*dataCollector)[componentNo]->getInterpPointInd();
 	DOFVector<double> *values = (*dataCollector)[componentNo]->getValues();
 	DOFVector<std::list<WorldVector<double> > > *dofCoords = (*dataCollector)[componentNo]->getDofCoords();
@@ -305,6 +319,59 @@ namespace AMDiS  { namespace io {
 	  }
 	}    
       }
+      
+      
+
+
+      template<typename T>
+      void Aux::writeVertexValues(T &file)
+      {	
+	DOFVector<int> *interpPointInd = (*dataCollector)[0]->getInterpPointInd();
+	std::vector<DOFVector<double>*> values(dataCollector->size());	
+	DOFVector<std::list<WorldVector<double> > > *dofCoords = (*dataCollector)[0]->getDofCoords();
+	
+	for (size_t i = 0; i < dataCollector->size(); i++)
+	  values[i] = (*dataCollector)[i]->getValues();
+	
+	DOFIterator<int> intPointIt(interpPointInd, USED_DOFS);
+	DOFVectorIterator<double> valueIt(values, USED_DOFS);
+	DOFIterator<std::list<WorldVector<double> > > coordIt(dofCoords, USED_DOFS);
+	
+	file << std::scientific;
+	file.precision(15);
+	    
+	// Write the values for all vertex DOFs.
+	for (intPointIt.reset(), valueIt.reset(), coordIt.reset(); 
+	     !intPointIt.end(); 
+	     ++intPointIt, ++valueIt, ++coordIt) {
+
+	  if (*intPointIt == -2) 
+	    for (int i = 0; i < static_cast<int>(coordIt->size()); i++) {
+	      file << " " ;
+	      print(*valueIt, file);
+	      file << "\n";
+	    }
+	}        
+
+	// For the second dim case, write also the values of the interpolation points.
+	if ((dim == 2) && (degree > 1)) {
+	  DOFVector< std::list<WorldVector<double> > >::Iterator 
+	    interpCoordIt((*dataCollector)[0]->getInterpPointCoords(), USED_DOFS);
+