* Forgotten a small file

AMDiS/src/VtkWriter.hh

+#include <list>
+#include <vector>
+#include "DOFVector.h"
+
+namespace AMDiS {
+
+  template<typename T>
+  void VtkWriter::writeFile(T &file)
+  {
+    int nVertices = (*dataCollector)[0]->getNumberVertices();
+    int nElements = (*dataCollector)[0]->getNumberElements();
+    int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);
+
+    if ((dim == 2) && (degree == 2)) {
+      nVertices += (*dataCollector)[0]->getNumberInterpPoints();
+      nElements *= 4;
+    } else if ((dim == 2) && (degree == 3)) {
+      nVertices += (*dataCollector)[0]->getNumberInterpPoints();
+      nElements *= 9;
+    } else if ((dim == 2) && (degree == 4)) {
+      nVertices += (*dataCollector)[0]->getNumberInterpPoints();
+      nElements *= 16;
+    }
+
+    file << "<?xml version=\"1.0\"?>\n";
+    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
+    file << "  <UnstructuredGrid>\n";
+    file << "    <Piece NumberOfPoints=\"" << nVertices 
+	 << "\" NumberOfCells=\"" << nElements << "\">\n";
+    file << "      <Points>\n";
+    file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+    
+    writeVertexCoords(file);
+    
+    file << "        </DataArray>\n";
+    file << "      </Points>\n";
+    file << "      <Cells>\n";
+    file << "        <DataArray type=\"Int32\" Name=\"offsets\">\n";
+    
+    for (int i = 0; i < nElements; i++) {
+      file << " " << (i + 1) * vertices << "\n";
+    }
+
+    file << "        </DataArray>\n";
+    file << "        <DataArray type=\"UInt8\" Name=\"types\">\n";
+
+    for (int i = 0; i < nElements; i++) {
+      switch (vertices) {
+      case 2:
+	file << " 3\n";
+	break;
+      case 3: 
+	file << " 5\n";
+	break;
+      case 4:
+	file << " 10\n";
+	break;		    
+      default:
+	break;
+      }	   
+    }
+
+    file << "        </DataArray>\n";
+    file << "        <DataArray type=\"Int32\" Name=\"connectivity\">\n";
+    
+    writeConnectivity(file);
+    
+    file << "        </DataArray>\n";    
+    file << "      </Cells>\n";
+    file << "      <PointData>\n";
+    
+    for (int i = 0; i < static_cast<int>(dataCollector->size()); i++) {
+      file << "        <DataArray type=\"Float32\" Name=\"value" << i 
+	   << "\" format=\"ascii\">\n";
+      
+      writeVertexValues(file, i);
+      
+      file << "        </DataArray>\n";
+    }
+    
+    file << "      </PointData>\n";
+    file << "    </Piece>\n";
+    file << "  </UnstructuredGrid>\n";
+    file << "</VTKFile>\n";
+  }
+  
+  template<typename T>
+  void writeVertexCoords(T &file)
+  {
+    DOFVector< std::list<VertexInfo> > *vertexInfos = (*dataCollector)[0]->getVertexInfos();
+    DOFVector< std::list<VertexInfo> >::Iterator it(vertexInfos, USED_DOFS);
+    int counter = 0;
+
+    // For all DOFs of vertices, write the coordinates.
+    for (it.reset(); !it.end(); ++it) {
+      // for all vertex infos of this DOF
+      std::list<VertexInfo>::iterator it2;
+      for (it2 = it->begin(); it2 != it->end(); ++it2) {
+	it2->outputIndex = counter++;
+	writeCoord(file, it2->coords);
+      }
+    }
+
+    // For the second dim case, write also the interpolation points.
+    if ((dim == 2) && (degree > 1)) {
+      DOFVector< std::list< WorldVector<double> > > *interpPointCoords = (*dataCollector)[0]->getInterpPointCoords();
+      DOFVector< std::list< WorldVector<double> > >::Iterator pointIt(interpPointCoords, USED_DOFS);
+      
+      for (pointIt.reset(); !pointIt.end(); ++pointIt) {
+	std::list< WorldVector<double> >::iterator it2;
+	for (it2 = pointIt->begin(); it2 != pointIt->end(); ++it2) {
+	  writeCoord(file, *it2);
+	}
+      }
+    }
+  }
+
+  template<typename T>
+  void writeVertexValues(T &file, int componentNo)
+  {
+    //    DOFVector<int> *interpPointInd;
+    //    DOFVector<double> *values;
+    //    DOFVector< std::list<WorldVector<double> > > *dofCoords;
+
+    DOFVector<int> *interpPointInd = (*dataCollector)[componentNo]->getInterpPointInd();
+    DOFVector<double> *values = (*dataCollector)[componentNo]->getValues();
+    DOFVector< std::list<WorldVector<double> > > *dofCoords = (*dataCollector)[componentNo]->getDofCoords();
+
+    /*
+#ifdef _OPENMP
+#pragma omp critical 
+#endif
+    {   
+      interpPointInd = (*dataCollector)[componentNo]->getInterpPointInd();
+      values = (*dataCollector)[componentNo]->getValues();
+      dofCoords = (*dataCollector)[componentNo]->getDofCoords();
+    }
+    */
+    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
+    DOFVector<double>::Iterator valueIt(values, USED_DOFS);
+    DOFVector< std::list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS);
+          
+    // 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 << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\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)[componentNo]->getInterpPointCoords(), USED_DOFS);
+    
+      for (intPointIt.reset(), valueIt.reset(), interpCoordIt.reset();
+	   !intPointIt.end(); 
+	   ++intPointIt, ++valueIt, ++interpCoordIt) {      
+	
+	if (*intPointIt >= 0) {
+	  for (int i = 0; i < static_cast<int>(interpCoordIt->size()); i++) {
+	    file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n";
+	  }
+	}
+      }
+    }    
+  }
+  
+  template<typename T>
+  void writeConnectivity(T &file)
+  {
+    // For the second dim case, and if higher order Lagrange elements are used,
+    // write the connectivity by extra functions.
+    if ((dim == 2) && (degree == 2)) {
+      writeConnectivity_dim2_degree2(file);
+    } else if ((dim == 2) && (degree == 3)) {
+      writeConnectivity_dim2_degree3(file);
+    } else if ((dim == 2) && (degree == 4)) {
+      writeConnectivity_dim2_degree4(file);   
+    } else {
+      std::list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
+      std::list<ElementInfo>::iterator elementIt;
+      int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);
+      
+      for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
+	// for all vertices
+	for (int i = 0; i < vertices; i++) {
+	  file << " " << elementIt->vertexInfo[i]->outputIndex;
+	}
+	file << "\n";
+      } 
+    }
+  }
+  
+  template<typename T>
+  void writeConnectivity_dim2_degree2(T &file)
+  {
+    std::list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
+    std::list<ElementInfo>::iterator elementIt;
+
+    std::vector< std::vector<int> > *interpPoints = (*dataCollector)[0]->getInterpPoints();
+    std::vector< std::vector<int> >::iterator pointIt;
+
+    int nVertices = (*dataCollector)[0]->getNumberVertices();
+
+    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
+	 pointIt != interpPoints->end(); 
+	 ++pointIt, ++elementIt) {
+    
+      file << " " << elementIt->vertexInfo[0]->outputIndex
+	   << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[2] + nVertices << "\n";
+
+      file << " " << elementIt->vertexInfo[2]->outputIndex 
+	   << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[1] + nVertices << "\n";
+
+      file << " " << elementIt->vertexInfo[1]->outputIndex
+	   << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[2] + nVertices << "\n";
+
+      file << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[1] + nVertices 
+	   << " " << (*pointIt)[2] + nVertices << "\n";
+    }
+  }
+
+  template<typename T>
+  void writeConnectivity_dim2_degree3(T &file)
+  {
+    std::list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
+    std::list<ElementInfo>::iterator elementIt;
+
+    std::vector< std::vector<int> > *interpPoints = (*dataCollector)[0]->getInterpPoints();
+    std::vector< std::vector<int> >::iterator pointIt;
+
+    int nVertices = (*dataCollector)[0]->getNumberVertices();
+
+    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
+	 pointIt != interpPoints->end(); 
+	 ++pointIt, ++elementIt) {
+                          
+      file << " " << elementIt->vertexInfo[0]->outputIndex
+	   << " " << (*pointIt)[3] + nVertices
+	   << " " << (*pointIt)[4] + nVertices << "\n";
+
+      file << " " << (*pointIt)[4] + nVertices
+	   << " " << (*pointIt)[5] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << (*pointIt)[3] + nVertices
+	   << " " << (*pointIt)[4] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << (*pointIt)[2] + nVertices
+	   << " " << (*pointIt)[3] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << elementIt->vertexInfo[1]->outputIndex
+	   << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[5] + nVertices << "\n";
+
+      file << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[6] + nVertices
+	   << " " << (*pointIt)[5] + nVertices << "\n";
+
+      file << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[2] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << elementIt->vertexInfo[2]->outputIndex 
+	   << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[2] + nVertices << "\n";
+     
+      
+    }
+  }
+
+  template<typename T>
+  void writeConnectivity_dim2_degree4(T &file)
+  {
+    std::list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
+    std::list<ElementInfo>::iterator elementIt;
+
+    std::vector< std::vector<int> > *interpPoints = (*dataCollector)[0]->getInterpPoints();
+    std::vector< std::vector<int> >::iterator pointIt;
+
+    int nVertices = (*dataCollector)[0]->getNumberVertices();
+
+    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
+	 pointIt != interpPoints->end(); 
+	 ++pointIt, ++elementIt) {
+
+      file << " " << elementIt->vertexInfo[0]->outputIndex 
+	   << " " << (*pointIt)[5] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << (*pointIt)[5] + nVertices
+	   << " " << (*pointIt)[9] + nVertices
+	   << " " << (*pointIt)[6] + nVertices << "\n";
+
+      file << " " << (*pointIt)[6] + nVertices
+	   << " " << (*pointIt)[7] + nVertices
+	   << " " << (*pointIt)[9] + nVertices << "\n";
+
+      file << " " << (*pointIt)[7] + nVertices
+	   << " " << (*pointIt)[9] + nVertices
+	   << " " << (*pointIt)[10] + nVertices << "\n";
+
+      file << " " << (*pointIt)[7] + nVertices
+	   << " " << (*pointIt)[8] + nVertices
+	   << " " << (*pointIt)[10] + nVertices << "\n";
+
+      file << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[8] + nVertices
+	   << " " << (*pointIt)[10] + nVertices << "\n";
+
+      file << " " << elementIt->vertexInfo[1]->outputIndex 
+	   << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[8] + nVertices << "\n";
+
+      file << " " << (*pointIt)[4] + nVertices
+	   << " " << (*pointIt)[5] + nVertices
+	   << " " << (*pointIt)[9] + nVertices << "\n";
+
+      file << " " << (*pointIt)[4] + nVertices
+	   << " " << (*pointIt)[9] + nVertices
+	   << " " << (*pointIt)[11] + nVertices << "\n";
+
+      file << " " << (*pointIt)[9] + nVertices
+	   << " " << (*pointIt)[10] + nVertices
+	   << " " << (*pointIt)[11] + nVertices << "\n";
+
+      file << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[10] + nVertices
+	   << " " << (*pointIt)[11] + nVertices << "\n";
+
+      file << " " << (*pointIt)[0] + nVertices
+	   << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[10] + nVertices << "\n";
+
+      file << " " << (*pointIt)[3] + nVertices
+	   << " " << (*pointIt)[4] + nVertices
+	   << " " << (*pointIt)[11] + nVertices << "\n";
+
+      file << " " << (*pointIt)[2] + nVertices
+	   << " " << (*pointIt)[3] + nVertices
+	   << " " << (*pointIt)[11] + nVertices << "\n";
+
+      file << " " << (*pointIt)[1] + nVertices
+	   << " " << (*pointIt)[2] + nVertices
+	   << " " << (*pointIt)[11] + nVertices << "\n";
+
+      file << " " << elementIt->vertexInfo[2]->outputIndex 
+	   << " " << (*pointIt)[2] + nVertices
+	   << " " << (*pointIt)[3] + nVertices << "\n";
+
+    }
+  }
+
+}