Reading and writing of FPMA format is done.

Unstructured_mesh/UnstructuredMesh/MeshIO/MeshReader/FPMAMeshReader.h

+#ifndef FPMAMESHREADER_H
+#define FPMAMESHREADER_H
+#include "MeshReader.h"
+#include "../../MeshElements/MeshElement.h"
+#include <iostream>
+#include <unordered_map>
+
+template <unsigned int MeshDimension>
+class FPMAMeshReader : public MeshReader<MeshDimension>{
+};
+
+template <>
+class FPMAMeshReader<3> : public MeshReader<3> {
+public:
+
+
+
+
+    template<typename IndexType, typename Real, unsigned int ...Reserve>
+    void loadVertices(std::istream& ist,MeshElements<3, IndexType, Real, Reserve...>& mesh){
+        IndexType numVert;
+        ist >> numVert;
+
+        mesh.getVertices().resize(numVert);
+
+        for (IndexType i = 0; i < numVert; i++) {
+
+            typename MeshElements<3, IndexType, Real, Reserve...>::Vertex& vert = mesh.getVertices().at(i);
+            vert.setIndex(i);
+            ist >> vert[0];
+            ist >> vert[1];
+            ist >> vert[2];
+        }
+    }
+
+
+    template<typename IndexType, typename Real, unsigned int ...Reserve>
+    void loadFaces(std::istream& ist,MeshElements<3, IndexType, Real, Reserve...>& mesh){
+        // map of constructed edges
+        std::unordered_map<std::string, IndexType> edges;
+
+        IndexType numFace;
+        ist >> numFace;
+
+        mesh.getFaces().resize(numFace);
+
+        // read face verts
+        for (IndexType faceIndex = 0; faceIndex < numFace; faceIndex++) {
+            mesh.getFaces().at(faceIndex).setIndex(faceIndex);
+            IndexType numVert;
+            ist >> numVert;
+
+            std::vector<IndexType> vertices(numVert);
+            for(IndexType j = 0; j < numVert; j++){
+                ist >> vertices.at(j);
+            }
+
+
+            for(IndexType j = 0; j < numVert; j++){
+
+                IndexType iA = vertices.at(j), iB = vertices.at((j+1)%numVert);
+                std::string edgeKey = iA < iB ? std::to_string(iA) +";"+ std::to_string(iB) : std::to_string(iB) +";"+ std::to_string(iA);
+                typename std::unordered_map<std::string, IndexType>::iterator edgeIt = edges.find(edgeKey);
+
+                IndexType edgeIndex = IndexType();
+
+                if (edgeIt == edges.end()){
+
+                    edgeIndex = mesh.getEdges().size();
+                    mesh.getEdges().push_back({});
+                    mesh.getEdges().at(edgeIndex).setVertexAIndex(iA);
+                    mesh.getEdges().at(edgeIndex).setVertexBIndex(iB);
+                    mesh.getEdges().at(edgeIndex).setIndex(edgeIndex);
+
+                    edges[edgeKey] = edgeIndex;
+                } else {
+                    edgeIndex = edgeIt->second;
+                }
+                try {
+                    mesh.getFaces().at(faceIndex).getSubelements().addSubelement(edgeIndex, true);
+                } catch (std::runtime_error& err) {
+                    throw std::runtime_error(std::string("The number of edges has overflew the prealocated memory ") +
+                                             std::to_string(mesh.getFaces().at(faceIndex).getSubelements().size()) +
+                                             " in face number " + std::to_string(faceIndex) + " while adding edge number " +
+                                             std::to_string(edgeIndex) + " ( vertA : " +std::to_string(iA) +
+                                             "vertB : " + std::to_string(iB) + ")\n" + err.what());
+                }
+            }
+
+        }
+        mesh.getEdges().shrink_to_fit();
+    }
+
+
+    template<typename IndexType, typename Real, unsigned int ...Reserve>
+    void loadCells(std::istream& ist,MeshElements<3, IndexType, Real, Reserve...>& mesh){
+        IndexType numCells;
+        ist >> numCells;
+
+        mesh.getCells().resize(numCells);
+        // read cells
+        for (IndexType cellIndex = 0; cellIndex < numCells; cellIndex++) {
+            mesh.getCells().at(cellIndex).setIndex(cellIndex);
+
+            // read cell faces
+            IndexType numFaces;
+            ist >> numFaces;
+
+
+            IndexType prevFace = INVALID_INDEX(IndexType);
+            for(IndexType j = 0; j < numFaces; j++){
+                IndexType faceIndex;
+                ist >> faceIndex;
+
+                if (j == 0){
+                    mesh.getCells().at(cellIndex).setBoundaryElementIndex(faceIndex);
+                }
+
+                if (prevFace != INVALID_INDEX(IndexType)){
+                    mesh.getFaces().at(prevFace).setNextBElem(faceIndex, cellIndex);
+
+                }
+
+                if (j == numFaces -1) {
+                    mesh.getFaces().at(faceIndex).setNextBElem(mesh.getCells().at(cellIndex).getBoundaryElementIndex(), cellIndex);
+                }
+                prevFace = faceIndex;
+            }
+
+
+        }
+
+    }
+
+
+
+    template<typename IndexType, typename Real, unsigned int ...Reserve>
+    void loadFromStream(std::istream& ist,MeshElements<3, IndexType, Real, Reserve...>& mesh){
+        ist.seekg(ist.beg);
+
+        loadVertices(ist, mesh);
+
+        loadFaces(ist, mesh);
+
+        loadCells(ist, mesh);
+
+    }
+
+
+};
+
+
+#endif // FPMAMESHREADER_H

Unstructured_mesh/UnstructuredMesh/MeshIO/MeshReader/VTKMeshReader.h

@@ -344,7 +344,9 @@ public:
 
             IndexType prevFaceIndex = INVALID_INDEX(IndexType);
             for (IndexType fi = 0; fi < faceOrder.size(); fi++) {
+
                 std::vector<int>& f = faceOrder.at(fi);
+
                 std::vector<IndexType> faceEdges;
                 for (int& index : f) {
                     faceEdges.push_back(edgeIndexes.at(index).first);

Unstructured_mesh/UnstructuredMesh/MeshIO/MeshWriter/FPMAMeshWriter.h

+#ifndef FPMAMESHWRITER_H
+#define FPMAMESHWRITER_H
+
+#include "MeshWriter.h"
+#include "../../MeshElements/MeshElement.h"
+#include "../../MeshDataContainer/MeshDataContainer.h"
+#include "../../MeshFunctions/MeshFunctions.h"
+#include <map>
+#include <ostream>
+
+template<unsigned int MeshDimension, typename IndexType = size_t, typename Real = double>
+class FPMAMeshWriter : public MeshWriter<MeshDimension>{
+public:
+    FPMAMeshWriter() = default;
+    template<unsigned int ...Reserve>
+    FPMAMeshWriter(const MeshElements<MeshDimension, IndexType, Real, Reserve...>&){}
+};
+
+
+template<typename IndexType, typename Real>
+class FPMAMeshWriter<3, IndexType, Real> : public MeshWriter<3>{
+
+     using writer = MeshWriter<3>;
+    size_t lastHash;
+    /**
+     * @brief faceVert<HR>
+     * correctly erdered vertices for all faces
+     */
+    MeshDataContainer<std::vector<IndexType>, 2> faceVert;
+
+    /**
+     * @brief cellFace<HR>
+     * original order of faces
+     */
+    MeshDataContainer<std::vector<IndexType>, 3> cellFace;
+
+    /**
+     * @brief indexFace<HR>
+     * This funcion return indexes of vertices of a face in correct order
+     * with respect to edge orientation
+     * in output vector verticesIndexed
+     * @param mesh the structure of the mesh
+     * @param face the face of the mesh to be indexed
+     * @param cell the cell which are the vertices being indexed to
+     * @param faceEdgeOri the orientation of the edges to the faces
+     * @param verticesIndexed output vector of indexed vertices
+     */
+    template<unsigned int ...Reserve>
+    void indexFace(MeshElements<3, IndexType, Real, Reserve...>& mesh,
+                   typename MeshElements<3, IndexType, Real, Reserve...>::Face& face,
+                   std::vector<IndexType>& verticesIndexed){
+
+        // export the face in "left" direction see VTK export
+        IndexType startVertex = mesh.getEdges().at(face.getSubelements()[0].index).getVertexBIndex();
+        IndexType nextVertex = mesh.getEdges().at(face.getSubelements()[0].index).getVertexAIndex();
+
+
+        verticesIndexed.push_back(startVertex);
+
+        IndexType lastWrittenEdge = face.getSubelements()[0].index;
+        while (startVertex != nextVertex){
+            for (auto& sube : face.getSubelements()) {
+                auto &edge = mesh.getEdges().at(sube.index);
+
+                if (edge.getIndex() != lastWrittenEdge) {
+                    if (edge.getVertexAIndex() == nextVertex) {
+                        lastWrittenEdge = edge.getIndex();
+                        verticesIndexed.push_back(edge.getVertexAIndex());
+                        nextVertex = edge.getVertexBIndex();
+                    } else if (edge.getVertexBIndex() == nextVertex) {
+                        lastWrittenEdge = edge.getIndex();
+                        verticesIndexed.push_back(edge.getVertexBIndex());
+                        nextVertex = edge.getVertexAIndex();
+                    }
+                }
+            }
+        }
+        verticesIndexed.shrink_to_fit();
+    }
+
+
+    /**
+     * @brief indexMesh<HR>
+     * This function creates vector of indexes of vertices for each cell in order suitable
+     * for VTK output. Moreover in case of elements type different from
+     * any VTK cell type, the cell is split into tetrahedrons. The tetrahedrons are
+     * made for each edge of every faces. The tetrahedrons is made of both edges vertices
+     * and cell and face center.
+     * @param mesh Mesh to be indexed
+     * @param cellTypes Vector of known cell types. If a cell type is not known then a method
+     * of splitting into tetrahedrons is used.
+     */
+    template<unsigned int ...Reserve>
+    void indexMesh(MeshElements<3, IndexType, Real, Reserve...>& mesh){
+
+
+        faceVert.template getDataByPos<0>().clear();
+        faceVert.alocateData(mesh);
+DBGMSG("indexing mesh");
+        // write cells of the mesh
+        // prepare connections
+        auto cellVert = MeshConnections<3,0>::connections(mesh);
+
+        for (typename MeshElements<3, IndexType, Real, Reserve...>::Face& face : mesh.getFaces()){
+            indexFace(mesh, face, faceVert.at(face));
+        }
+
+        cellFace = MeshConnections<3,2, Order::ORDER_ORIGINAL>::connections(mesh);
+    }
+
+public:
+    /**
+     * @brief writeToStream
+     * Exports the mesh to the output stream in VTK format.
+     * If the mesh is written for the first time, it is indexed.
+     * @see indexMesh
+     * @param ost
+     * @param mesh
+     * @param cellTypes the types of cells in NativeType format
+     */
+    template<unsigned int ...Reserve>
+    void writeToStream(std::ostream& ost,
+                       MeshElements<3, IndexType, Real, Reserve...>& mesh){
+        // create index of mesh if the mesh has changed
+        size_t curHash = writer::computeHash(mesh);
+
+        // if the mesh is not the same as it was,
+        // then update the index
+        if (lastHash != curHash){
+            indexMesh(mesh);
+        }
+        lastHash = curHash;
+
+        // first write verices
+        ost << mesh.getVertices().size() << std::endl;
+
+        for(auto& vert : mesh.getVertices()) {
+            ost << vert[0] << ' ' << vert[1] << ' ' << vert[2] <<"\n";
+        }
+
+        ost << std::endl;
+
+        // write faces
+        ost << mesh.getFaces().size() << std::endl;
+        for(auto& face : mesh.getFaces()) {
+            ost << faceVert.at(face).size() << ' ';
+            for (IndexType& vertIndex : faceVert.at(face)) {
+                ost << vertIndex << ' ';
+            }
+            ost << std::endl;
+        }
+
+        // write cells of the mesh
+        ost << mesh.getCells().size() << std::endl;
+        for(auto& cell : mesh.getCells()) {
+            ost << cellFace.at(cell).size() << ' ';
+            for (IndexType& faceIndex : cellFace.at(cell)) {
+                ost << faceIndex << ' ';
+            }
+            ost << std::endl;
+        }
+        ost << 0;
+    }
+};
+
+
+
+
+
+
+#endif // FPMAMESHWRITER_H

Unstructured_mesh/Unstructured_mesh.pro

@@ -21,8 +21,10 @@ HEADERS += \
     UnstructuredMesh/MeshElements/MeshElement.h \
     UnstructuredMesh/MeshFunctions/MeshFunctions.h \
     UnstructuredMesh/MeshIO/MeshNativeType.h \
+    UnstructuredMesh/MeshIO/MeshReader/FPMAMeshReader.h \
     UnstructuredMesh/MeshIO/MeshReader/MeshReader.h \
     UnstructuredMesh/MeshIO/MeshReader/VTKMeshReader.h \
+    UnstructuredMesh/MeshIO/MeshWriter/FPMAMeshWriter.h \
     UnstructuredMesh/MeshIO/MeshWriter/MeshWriter.h \
     UnstructuredMesh/MeshIO/MeshWriter/VTKMeshWriter.h \
     UnstructuredMesh/UnstructedMeshDefine.h \

Unstructured_mesh/Unstructured_mesh.pro.user

 <?xml version="1.0" encoding="UTF-8"?>
 <!DOCTYPE QtCreatorProject>
-<!-- Written by QtCreator 4.10.0, 2019-10-22T14:43:15. -->
+<!-- Written by QtCreator 4.10.0, 2019-10-25T12:44:25. -->
 <qtcreator>
  <data>
   <variable>EnvironmentId</variable>

Unstructured_mesh/main.cpp

@@ -5,6 +5,10 @@
 #include "UnstructuredMesh/MeshFunctions/MeshFunctions.h"
 #include "UnstructuredMesh/MeshIO/MeshReader/VTKMeshReader.h"
 #include "UnstructuredMesh/MeshIO/MeshWriter/VTKMeshWriter.h"
+
+#include "UnstructuredMesh/MeshIO/MeshReader/FPMAMeshReader.h"
+#include "UnstructuredMesh/MeshIO/MeshWriter/FPMAMeshWriter.h"
+
 #include "UnstructuredMesh/MeshDataContainer/MemberApproach.h"
 #include <fstream>
 #include <list>
@@ -817,11 +821,80 @@ DBGVAR(mesh.getVertices().size(),mesh.getEdges().size(), mesh.getFaces().size(),
     ofstream out3D("3D_test_mesh_output.vtk");
     writer.writeHeader(out3D, "test data");
     writer.writeToStream(out3D, mesh, reader.getCellTypes());
+
+}
+
+
+void testFPMARW(){
+    UnstructuredMesh<3, size_t, double, 8> mesh;
+    FPMAMeshReader<3> reader;
+    ifstream file("Poly_simple.fpma");
+    reader.loadFromStream(file, mesh);
+
+    DBGVAR(mesh.getCells().size(), mesh.getVertices().size());
+
+    auto faceVert = MeshConnections<2,0,ORDER_ORIGINAL>::connections(mesh);
+    DBGVAR(faceVert.getDataByPos<0>());
+
+    mesh.initializeCenters();
+    VTKMeshWriter<3, size_t, double> writer;
+    ofstream ofile("Poly_simple.vtk");
+    writer.writeHeader(ofile, "fpma_output_test");
+    writer.writeToStream(ofile, mesh, MeshDataContainer<MeshNativeType<3>::ElementType, 3>(mesh, MeshNativeType<3>::POLYHEDRON));
+
+    ofile << "CELL_DATA " << writer.cellVert.getDataByPos<0>().size() << endl;
+    ofile << "SCALARS cell_wrt_vertex_colour double 1\nLOOKUP_TABLE default" << endl;
+    size_t realIndex = 0;
+    for (size_t i = 0; i < writer.cellVert.getDataByPos<0>().size(); i++) {
+        auto iterator = writer.backwardCellIndexMapping.find(i);
+        if (iterator == writer.backwardCellIndexMapping.end()){
+            ofile << realIndex << ' ';
+            realIndex++;
+        } else {
+            ofile << iterator->second << ' ';
+            realIndex = iterator->second;
+        }
+    }
+    ofile.close();
+
+    FPMAMeshWriter<3,size_t, double> writerFPMA;
+    ofstream ofile1("Poly_simple_test.fpma");
+    writerFPMA.writeToStream(ofile1, mesh);
+
+
 }
 
 
+void testFPMA_poly(){
+    UnstructuredMesh<3, size_t, double, 12> mesh;
+    FPMAMeshReader<3> reader;
+    ifstream file("Spark_mesh.fpma");
+    reader.loadFromStream(file, mesh);
+
+    DBGVAR(mesh.getCells().size(), mesh.getVertices().size());
 
+    mesh.initializeCenters();
+    VTKMeshWriter<3, size_t, double> writer;
+    ofstream ofile("Spark_mesh.vtk");
+    writer.writeHeader(ofile, "fpma_output_test");
+    writer.writeToStream(ofile, mesh, MeshDataContainer<MeshNativeType<3>::ElementType, 3>(mesh, MeshNativeType<3>::POLYHEDRON));
 
+    ofile << "CELL_DATA " << writer.cellVert.getDataByPos<0>().size() << endl;
+    ofile << "SCALARS cell_wrt_vertex_colour double 1\nLOOKUP_TABLE default" << endl;
+    size_t realIndex = 0;
+    for (size_t i = 0; i < writer.cellVert.getDataByPos<0>().size(); i++) {
+        auto iterator = writer.backwardCellIndexMapping.find(i);
+        if (iterator == writer.backwardCellIndexMapping.end()){
+            ofile << realIndex << ' ';
+            realIndex++;
+        } else {
+            ofile << iterator->second << ' ';
+            realIndex = iterator->second;
+        }
+    }
+    ofile.close();
+
+}
 
 
 int main()
@@ -836,4 +909,5 @@ int main()
     //m.ComputeElementMeasures();
     //test3DMeshLoad();
 
+    testFPMA_poly();
 }