// // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology // All rights reserved. // Authors: Simon Vey, Thomas Witkowski et al. // // This file is part of AMDiS // // See also license.opensource.txt in the distribution. /** \file VtkVectorWriter.hh */ #ifndef AMDIS_VTKVECTORWRITER_HH #define AMDIS_VTKVECTORWRITER_HH #include <list> #include <vector> #include "DOFVector.h" #include <stdio.h> #include <string> #include <fstream> #include <sstream> #include <cmath> #include <boost/filesystem/operations.hpp> #include <boost/filesystem/convenience.hpp> #include <boost/lexical_cast.hpp> #ifdef HAVE_PARALLEL_DOMAIN_AMDIS #include <mpi.h> #endif #include "DataCollector.h" #include "SurfaceRegion_ED.h" #include "ElementRegion_ED.h" namespace AMDiS { template<typename S> int VtkVectorWriter::Impl<S>::writeFile(std::string name) { FUNCNAME("VtkVectorWriter::Impl<S>::writeFile()"); boost::iostreams::filtering_ostream file; switch (compress) { case GZIP: file.push(boost::iostreams::gzip_compressor()); name.append(".gz"); break; case BZIP2: file.push(boost::iostreams::bzip2_compressor()); name.append(".bz2"); break; default: break; } { std::ofstream swapfile(name.c_str(), std::ios::out | std::ios::trunc); TEST_EXIT(swapfile.is_open()) ("Cannot open file %s for writing!\n", name.c_str()); swapfile.close(); } file.push(boost::iostreams::file_descriptor_sink(name, std::ios::trunc)); writeFileToStream(file); return 0; } template<typename S> void VtkVectorWriter::Impl<S>::writeParallelFile(std::string name, int nRanks, std::string fnPrefix, std::string fnPostfix) { FUNCNAME("VtkVectorWriter::Impl<S>::writeParallelFile()"); boost::iostreams::filtering_ostream file; { std::ofstream swapfile(name.c_str(), std::ios::out | std::ios::trunc); TEST_EXIT(swapfile.is_open()) ("Cannot open file %s for writing!\n", name.c_str()); swapfile.close(); } file.push(boost::iostreams::file_descriptor_sink(name, std::ios::trunc)); file << "<?xml version=\"1.0\"?>\n"; file << "<VTKFile type=\"PUnstructuredGrid\">\n"; file << " <PUnstructuredGrid GhostLevel=\"0\">\n"; file << " <PPoints>\n" << " <PDataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\"/>\n" << " </PPoints>\n"; file << " <PCells>\n" << " <PDataArray type=\"Int32\" Name=\"offsets\"/>\n" << " <PDataArray type=\"UInt8\" Name=\"types\"/>\n" << " <PDataArray type=\"Int32\" Name=\"connectivity\"/>\n" << " </PCells>\n"; file << " <PPointData>\n"; for (int i = 0; i < static_cast<int>(dataCollector->size()); i++) file << " <PDataArray type=\"Float32\" Name=\"value" << i << "\" format=\"ascii\"/>\n"; file << " </PPointData>\n"; for (int i = 0; i < nRanks; i++) { std::stringstream oss; oss << fnPrefix << "-p" << i << "-" << fnPostfix; boost::filesystem::path filepath(oss.str()); file << " <Piece Source=\"" << boost::filesystem::basename(filepath) << boost::filesystem::extension(filepath) << "\"/>\n"; } file << " </PUnstructuredGrid>\n"; file << "</VTKFile>\n"; } template<typename S> template<typename T> void VtkVectorWriter::Impl<S>::writeFileToStream(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++) { S temp = (*((*dataCollector)[i]->getValues()))[0]; int numComponent = num_rows(temp); file << " <DataArray type=\"Float32\" Name=\"value" << i << "\" format=\"ascii\"" << (numComponent > 1 ? " NumberOfComponents=\"" + boost::lexical_cast<std::string>(numComponent) + "\"" : "") << ">\n"; writeVertexValues(file, i); file << " </DataArray>\n"; } file << " </PointData>\n"; file << " </Piece>\n"; file << " </UnstructuredGrid>\n"; file << "</VTKFile>\n"; } // template<typename S> template<typename T> void VtkVectorWriter::Impl<S>::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 for (std::list<VertexInfo>::iterator 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); counter = 0; for (pointIt.reset(); !pointIt.end(); ++pointIt) { for (std::list< WorldVector<double> >::iterator it2 = pointIt->begin(); it2 != pointIt->end(); ++it2) { counter++; writeCoord(file, *it2); } } } } template<typename S> template<typename T> void VtkVectorWriter::Impl<S>::writeVertexValues(T &file, int componentNo) { DOFVector<int> *interpPointInd = (*dataCollector)[componentNo]->getInterpPointInd(); DOFVector<S> *values = (*dataCollector)[componentNo]->getValues(); DOFVector< std::list<WorldVector<double> > > *dofCoords = (*dataCollector)[componentNo]->getDofCoords(); DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS); typename DOFVector<S>::Iterator valueIt(values, USED_DOFS); DOFVector< std::list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS); file << std::scientific; file.precision(5); // 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)[componentNo]->getInterpPointCoords(), USED_DOFS); for (intPointIt.reset(), valueIt.reset(), interpCoordIt.reset(); !intPointIt.end(); ++intPointIt, ++valueIt, ++interpCoordIt) { if (*intPointIt >= 0) { for (unsigned int i = 0; i < interpCoordIt->size(); i++) { file << " "; print(*valueIt,file); file << "\n"; } } } } } template<typename S> template<typename T> void VtkVectorWriter::Impl<S>::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 S> template<typename T> void VtkVectorWriter::Impl<S>::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 S> template<typename T> void VtkVectorWriter::Impl<S>::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 S> template<typename T> void VtkVectorWriter::Impl<S>::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"; } } } #endif // AMDIS_VTKVECTORWRITER_HH