diff --git a/AMDiS/src/VtkWriter.hh b/AMDiS/src/VtkWriter.hh new file mode 100644 index 0000000000000000000000000000000000000000..f5557991778a95562650d736c1ab62ea84726e45 --- /dev/null +++ b/AMDiS/src/VtkWriter.hh @@ -0,0 +1,369 @@ +#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"; + + } + } + +}