diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 9eaaea15fad8d0ea1f19b3a7f9bb037acbdb2eb3..ca443f084613785b3823475924ee4d87b800354b 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -2,12 +2,817 @@
#include <iostream>
#include <fstream>
+#include <set>
+#include <rang.hpp>
+
+template<typename T>
+inline std::string stringify(const T & t)
+{
+ std::ostringstream oss;
+ oss << t;
+ return oss.str();
+}
+
+template<>
+inline std::string stringify<std::string>(const std::string& t)
+{
+ return t;
+}
+
+class ErrorHandler
+{
+public:
+ enum Type {
+ asked, /**< execution request by the user*/
+ compilation, /**< syntax error in a language */
+ normal, /**< normal error due to a bad use of ff3d */
+ unexpected /**< Unexpected execution error */
+ };
+
+private:
+ const std::string __filename; /**< The source file name where the error occured */
+ const size_t __lineNumber; /**< The line number where exception was raised */
+ const std::string __errorMessage; /**< The reporting message */
+
+ const Type __type; /**< the type of the error */
+public:
+ /**
+ * Prints the error message
+ *
+ */
+ virtual void writeErrorMessage();
+
+ /**
+ * The copy constructor
+ *
+ * @param e an handled error
+ */
+ ErrorHandler(const ErrorHandler& e)
+ : __filename(e.__filename),
+ __lineNumber(e.__lineNumber),
+ __errorMessage(e.__errorMessage),
+ __type(e.__type)
+ {
+ ;
+ }
+
+ /**
+ * Constructor
+ *
+ * @param filename the source file name
+ * @param lineNumber the source line
+ * @param errorMessage the reported message
+ * @param type the type of the error
+ */
+ ErrorHandler(const std::string& filename,
+ const size_t& lineNumber,
+ const std::string& errorMessage,
+ const Type& type);
+
+ /**
+ * The destructor
+ *
+ */
+ virtual ~ErrorHandler()
+ {
+ ;
+ }
+};
+void ErrorHandler::writeErrorMessage()
+{
+ switch(__type) {
+ case asked: {
+ std::cerr << "\nremark: exit command explicitly called\n";
+ }
+ case normal: {
+ std::cerr << '\n' << __filename << ':' << __lineNumber
+ << ":remark: emitted the following message\n";
+ std::cerr << "error: " << __errorMessage << '\n';
+ break;
+ }
+ case compilation: {
+ std::cerr << "\nline " << __lineNumber << ':' << __errorMessage << '\n';
+ break;
+ }
+ case unexpected: {
+ std::cerr << '\n' << __filename << ':' << __lineNumber << ":\n" << __errorMessage << '\n';
+ std::cerr << "\nUNEXPECTED ERROR: this should not occure, please report it\n";
+ std::cerr << "\nBUG REPORT: Please send bug reports to:\n"
+ << " ff3d-dev@nongnu.org or freefem@ann.jussieu.fr\n"
+ << "or better, use the Bug Tracking System:\n"
+ << " http://savannah.nongnu.org/bugs/?group=ff3d\n";
+ break;
+ }
+ default: {
+ std::cerr << __filename << ':' << __lineNumber << ": " << __errorMessage << '\n';
+ std::cerr << __FILE__ << ':' << __LINE__ << ":remark: error type not implemented!\n";
+ }
+ }
+}
+
+ErrorHandler::
+ErrorHandler(const std::string& filename,
+ const size_t& lineNumber,
+ const std::string& errorMessage,
+ const Type& type)
+ : __filename(filename),
+ __lineNumber(lineNumber),
+ __errorMessage(errorMessage),
+ __type(type)
+{
+ ;
+}
+
GmshReader::GmshReader(const std::string& filename)
+ : m_filename(filename)
+{
+ m_fin.open(m_filename);
+ if (not m_fin) {
+ std::cerr << "cannot read file '" << m_filename << "'\n";
+ std::exit(0);
+ }
+
+ // gmsh 2.2 format keywords
+ __keywordList["$MeshFormat"] = MESHFORMAT;
+ __keywordList["$EndMeshFormat"] = ENDMESHFORMAT;
+ __keywordList["$Nodes"] = NODES;
+ __keywordList["$EndNodes"] = ENDNODES;
+ __keywordList["$Elements"] = ELEMENTS;
+ __keywordList["$EndElements"] = ENDELEMENTS;
+
+ __numberOfPrimitiveNodes.resize(16);
+ __numberOfPrimitiveNodes[ 0] = 2; // edge
+ __numberOfPrimitiveNodes[ 1] = 3; // triangle
+ __numberOfPrimitiveNodes[ 2] = 4; // quadrangle
+ __numberOfPrimitiveNodes[ 3] = 4; // Tetrahedron
+ __numberOfPrimitiveNodes[ 4] = 8; // Hexaredron
+ __numberOfPrimitiveNodes[ 5] = 6; // Prism
+ __numberOfPrimitiveNodes[ 6] = 5; // Pyramid
+ __numberOfPrimitiveNodes[ 7] = 3; // second order edge
+ __numberOfPrimitiveNodes[ 8] = 6; // second order triangle
+ __numberOfPrimitiveNodes[ 9] = 9; // second order quadrangle
+ __numberOfPrimitiveNodes[10] = 10; // second order tetrahedron
+ __numberOfPrimitiveNodes[11] = 27; // second order hexahedron
+ __numberOfPrimitiveNodes[12] = 18; // second order prism
+ __numberOfPrimitiveNodes[13] = 14; // second order pyramid
+ __numberOfPrimitiveNodes[14] = 1; // point
+
+ __primitivesNames[0] = "edges";
+ __supportedPrimitives[0] = false;
+ __primitivesNames[1] = "triangles";
+ __supportedPrimitives[1] = true;
+ __primitivesNames[2] = "quadrangles";
+ __supportedPrimitives[2] = true;
+ __primitivesNames[3] = "tetrahedra";
+ __supportedPrimitives[3] = true;
+ __primitivesNames[4] = "hexahedra";
+ __supportedPrimitives[4] = true;
+ __primitivesNames[5] = "prisms";
+ __supportedPrimitives[5] = false;
+ __primitivesNames[6] = "pyramids";
+ __supportedPrimitives[6] = false;
+ __primitivesNames[7] = "second order edges";
+ __supportedPrimitives[7] = false;
+ __primitivesNames[8] = "second order triangles";
+ __supportedPrimitives[8] = false;
+ __primitivesNames[9] = "second order quadrangles";
+ __supportedPrimitives[9] = false;
+ __primitivesNames[10] = "second order tetrahedra";
+ __supportedPrimitives[10]= false;
+ __primitivesNames[11] = "second order hexahedra";
+ __supportedPrimitives[11]= false;
+ __primitivesNames[12] = "second order prisms";
+ __supportedPrimitives[12]= false;
+ __primitivesNames[13] = "second order pyramids";
+ __supportedPrimitives[13]= false;
+ __primitivesNames[14] = "point";
+ __supportedPrimitives[14]= false;
+
+ std::cout << "Reading file '" << m_filename << "'\n";
+
+ // Getting vertices list
+ GmshReader::Keyword kw = this->__nextKeyword();
+ switch(kw.second) {
+ // case NOD: {
+ // this->__readGmsh1();
+ // break;
+ // }
+ case MESHFORMAT: {
+ double fileVersion = this->_getReal();
+ if (fileVersion != 2.2) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "Cannot read Gmsh format '"+stringify(fileVersion)+"'",
+ ErrorHandler::normal);
+ }
+ int fileType = this->_getInteger();
+ __binary = (fileType == 1);
+ if ((fileType < 0) or (fileType > 1)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "Cannot read Gmsh file type '"+stringify(fileType)+"'",
+ ErrorHandler::normal);
+ }
+
+ int dataSize = this->_getInteger();
+ if (dataSize != sizeof(double)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "Data size not supported '"+stringify(dataSize)+"'",
+ ErrorHandler::normal);
+ }
+
+ if (__binary) {
+// int one=0;
+
+// fseek(__ifh,1L,SEEK_CUR);
+// fread(reinterpret_cast<char*>(&one),sizeof(int),1,__ifh);
+
+// if (one == 1) {
+// __convertEndian = false;
+// } else {
+// invertEndianess(one);
+
+// if (one == 1) {
+// __convertEndian = true;
+// } else {
+// throw ErrorHandler(__FILE__,__LINE__,
+// "Cannot determine endianess",
+// ErrorHandler::normal);
+// }
+// }
+
+// std::cout << "- Binary format: ";
+// #ifdef WORDS_BIGENDIAN
+// if (not __convertEndian) {
+// std::cout << "big endian\n";
+// } else {
+// std::cout << "little endian\n";
+// }
+// #else // WORDS_BIGENDIAN
+// if (not __convertEndian) {
+// std::cout << "little endian\n";
+// } else {
+// std::cout << "big endian\n";
+// }
+// #endif // WORDS_BIGENDIAN
+ }
+
+ kw = this->__nextKeyword();
+ if (kw.second != ENDMESHFORMAT) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': expecting $EndMeshFormat, '"+kw.first+"' was found",
+ ErrorHandler::normal);
+ }
+
+ this->__readGmshFormat2_2();
+
+ break;
+ }
+ default: {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "cannot determine format version of '"+m_filename+"'",
+ ErrorHandler::normal);
+ }
+ }
+
+ this->__proceedData();
+ // this->__createMesh();
+}
+
+void GmshReader::__readVertices()
+{
+ const int numberOfVerices = this->_getInteger();
+ std::cout << "- Number Of Vertices: " << numberOfVerices << '\n';
+ if (numberOfVerices<0) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+this->m_filename
+ +"': number of vertices is negative",
+ ErrorHandler::normal);
+ }
+
+ __verticesNumbers.resize(numberOfVerices);
+ __vertices = Kokkos::View<R3*>("vertices",numberOfVerices);
+
+ if (not __binary) {
+ for (int i=0; i<numberOfVerices; ++i) {
+ __verticesNumbers[i] = this->_getInteger();
+ const double x = this->_getReal();
+ const double y = this->_getReal();
+ const double z = this->_getReal();
+ __vertices[i] = TinyVector<3, double>(x,y,z);
+ }
+ } else {
+ // fseek(m_fin.file()__ifh,1L,SEEK_CUR);
+ // for (int i=0; i<numberOfVerices; ++i) {
+ // int number = 0;
+ // fread(reinterpret_cast<char*>(&number),sizeof(int),1,__ifh);
+ // __verticesNumbers[i] = number;
+ // TinyVector<3,double> x;
+ // fread(reinterpret_cast<char*>(&(x[0])),sizeof(double),3,__ifh);
+ // for (size_t j=0; j<3; ++j) {
+ // __vertices[i][j] = x[j];
+ // }
+ // }
+
+ // if (__convertEndian) {
+ // for (int i=0; i<numberOfVerices; ++i) {
+ // invertEndianess(__verticesNumbers[i]);
+ // for (size_t j=0; j<3; ++j) {
+ // invertEndianess(vertices[i][j]);
+ // }
+ // }
+ // }
+ }
+}
+
+// void
+// GmshReader::__readElements1()
+// {
+// const int numberOfElements = this->_getInteger();
+// std::cout << "- Number Of Elements: " << numberOfElements << '\n';
+// if (numberOfElements<0) {
+// throw ErrorHandler(__FILE__,__LINE__,
+// "reading file '"+m_filename
+// +"': number of elements is negative",
+// ErrorHandler::normal);
+// }
+
+// __elementType.resize(numberOfElements);
+// __references.resize(numberOfElements);
+// __elementVertices.resize(numberOfElements);
+
+// for (int i=0; i<numberOfElements; ++i) {
+// this->_getInteger(); // drops element number
+// const int elementType = this->_getInteger()-1;
+
+// if ((elementType < 0) or (elementType > 14)) {
+// throw ErrorHandler(__FILE__,__LINE__,
+// "reading file '"+m_filename
+// +"': unknown element type '"+stringify(elementType)+"'",
+// ErrorHandler::normal);
+// }
+// __elementType[i] = elementType;
+// __references[i] = this->_getInteger(); // physical reference
+// this->_getInteger(); // drops entity reference
+
+// const int numberOfVertices = this->_getInteger();
+// if (numberOfVertices != __numberOfPrimitiveNodes[elementType]) {
+// std::stringstream errorMsg;
+// errorMsg << "reading file '" <<m_filename
+// << "':number of vertices (" << numberOfVertices
+// << ") does not match expected ("
+// << __numberOfPrimitiveNodes[elementType] << ")" << std::ends;
+
+// throw ErrorHandler(__FILE__,__LINE__,
+// errorMsg.str(),
+// ErrorHandler::normal);
+// }
+// __elementVertices[i].resize(numberOfVertices);
+// for (int j=0; j<numberOfVertices; ++j) {
+// __elementVertices[i][j] = this->_getInteger();
+// }
+// }
+// }
+
+void
+GmshReader::__readElements2_2()
+{
+ const int numberOfElements = this->_getInteger();
+ std::cout << "- Number Of Elements: " << numberOfElements << '\n';
+ if (numberOfElements<0) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': number of elements is negative",
+ ErrorHandler::normal);
+ }
+
+ __elementType.resize(numberOfElements);
+ __references.resize(numberOfElements);
+ __elementVertices.resize(numberOfElements);
+
+ if (not __binary) {
+ for (int i=0; i<numberOfElements; ++i) {
+ this->_getInteger(); // drops element number
+ const int elementType = this->_getInteger()-1;
+
+ if ((elementType < 0) or (elementType > 14)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': unknown element type '"+stringify(elementType)+"'",
+ ErrorHandler::normal);
+ }
+ __elementType[i] = elementType;
+ const int numberOfTags = this->_getInteger();
+ __references[i] = this->_getInteger(); // physical reference
+ for (int tag=1; tag<numberOfTags; ++tag) {
+ this->_getInteger(); // drops remaining tags
+ }
+
+ const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
+ __elementVertices[i].resize(numberOfVertices);
+ for (int j=0; j<numberOfVertices; ++j) {
+ __elementVertices[i][j] = this->_getInteger();
+ }
+ }
+ } else {
+ // fseek(__ifh,1L,SEEK_CUR);
+ // int i=0;
+ // for (;i<numberOfElements;) {
+ // int elementType = 0;
+ // fread(reinterpret_cast<char*>(&elementType),sizeof(int),1,__ifh);
+ // if (__convertEndian) {
+ // invertEndianess(elementType);
+ // }
+ // --elementType;
+
+ // if ((elementType < 0) or (elementType > 14)) {
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "reading file '"+m_filename
+ // +"': unknown element type '"+stringify(elementType)+"'",
+ // ErrorHandler::normal);
+ // }
+
+ // int elementTypeNumber = 0;
+ // fread(reinterpret_cast<char*>(&elementTypeNumber),sizeof(int),1,__ifh);
+ // if (__convertEndian) {
+ // invertEndianess(elementTypeNumber);
+ // }
+
+ // int numberOfTags = 0;
+ // fread(reinterpret_cast<char*>(&numberOfTags),sizeof(int),1,__ifh);
+ // if (__convertEndian) {
+ // invertEndianess(numberOfTags);
+ // }
+
+ // for (int k=0; k<elementTypeNumber; ++k) {
+ // __elementType[i] = elementType;
+
+ // int elementNumber = 0;
+ // fread(reinterpret_cast<char*>(&elementNumber),sizeof(int),1,__ifh);
+
+ // int reference = 0;
+ // fread(reinterpret_cast<char*>(&reference),sizeof(int),1,__ifh);
+
+ // __references[i] = reference; // physical reference
+ // for (int tag=1; tag<numberOfTags; ++tag) {
+ // int t;
+ // fread(reinterpret_cast<char*>(&t),sizeof(int),1,__ifh);
+ // }
+
+ // const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
+ // __elementVertices[i].resize(numberOfVertices);
+ // fread(reinterpret_cast<char*>(&(__elementVertices[i][0])),sizeof(int),numberOfVertices,__ifh);
+ // ++i;
+ // }
+ // }
+
+ // if (__convertEndian) {
+ // for (size_t i=0; i<__references.size(); ++i) {
+ // invertEndianess(__references[i]);
+ // }
+ // for (size_t i=0; i<__elementVertices.size(); ++i) {
+ // for (size_t j=0; j<__elementVertices[i].size(); ++i) {
+ // invertEndianess(__elementVertices[i][j]);
+ // }
+ // }
+ // }
+ }
+}
+
+void
+GmshReader::__proceedData()
{
- std::ifstream fin(filename);
- if (not fin) {
- std::cerr << "cannot read file '" << filename << "'\n";
+ TinyVector<15,int> elementNumber = zero;
+ std::vector<std::set<size_t> > elementReferences(15);
+
+ for (size_t i=0; i<__elementType.size(); ++i) {
+ const int elementType = __elementType[i];
+ elementNumber[elementType]++;
+ elementReferences[elementType].insert(__references[i]);
+ }
+
+ for (size_t i=0; i<elementNumber.dimension(); ++i) {
+ if (elementNumber[i]>0) {
+ std::cout << " - Number of "
+ << __primitivesNames[i]
+ << ": " << elementNumber[i];
+ if (not(__supportedPrimitives[i])) {
+ std::cout << " [IGNORED]";
+ } else {
+ std::cout << " references: ";
+ for(std::set<size_t>::const_iterator iref
+ = elementReferences[i].begin();
+ iref != elementReferences[i].end(); ++iref) {
+ if (iref != elementReferences[i].begin()) std::cout << ',';
+ std::cout << *iref;
+ }
+
+ switch (i) {
+ // Supported entities
+ case 1: { // triangles
+ __triangles
+ = Kokkos::View<Triangle*>("triangles",
+ elementNumber[i]);
+ __triangles_ref.resize(elementNumber[i]);
+ break;
+ }
+ case 2: { // quadrangles
+ __quadrangles
+ = Kokkos::View<Quadrangle*>("quadrangles",
+ elementNumber[i]);
+ __quadrangles_ref.resize(elementNumber[i]);
+ break;
+ }
+ case 3: { // tetrahedra
+ __tetrahedra
+ = Kokkos::View<Tetrahedron*>("tetrahedra",
+ elementNumber[i]);
+ __tetrahedra_ref.resize(elementNumber[i]);
+ break;
+ }
+ case 4: {// hexahedra
+ __hexahedra
+ = Kokkos::View<Hexahedron*>("hexahedra",
+ elementNumber[i]);
+ __hexahedra_ref.resize(elementNumber[i]);
+ }
+ // Ignored entities
+ case 0: // edge
+ case 14: {// point
+ break;
+ }
+ // Unsupported entities
+ case 5: // prism
+ case 6: // pyramid
+ case 7: // second order edge
+ case 8: // second order triangle
+ case 9: // second order quadrangle
+ case 10: // second order tetrahedron
+ case 11: // second order hexahedron
+ case 12: // second order prism
+ case 13: // second order pyramid
+ default: {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': ff3d cannot read this kind of element",
+ ErrorHandler::normal);
+ }
+ }
+ }
+ std::cout << '\n';
+ }
+ }
+
+ std::cout << "- Building correspondance table\n";
+ int minNumber = std::numeric_limits<int>::max();
+ int maxNumber = std::numeric_limits<int>::min();
+ for (size_t i=0; i<__verticesNumbers.size(); ++i) {
+ const int& vertexNumber = __verticesNumbers[i];
+ minNumber = std::min(minNumber,vertexNumber);
+ maxNumber = std::max(maxNumber,vertexNumber);
+ }
+
+ if (minNumber<0) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': ff3d does not implement negative vertices number",
+ ErrorHandler::normal);
+ }
+
+ // A value of -1 means that the vertex is unknown
+ __verticesCorrepondance.resize(maxNumber+1,-1);
+
+ for (size_t i=0; i<__verticesNumbers.size(); ++i) {
+
+ __verticesCorrepondance[__verticesNumbers[i]] = i;
+ }
+
+ // reset element number to count them while filling structures
+ elementNumber = zero;
+
+ // Element structures filling
+ for (size_t i=0; i<__elementType.size(); ++i) {
+ std::vector<int>& elementVertices = __elementVertices[i];
+ switch (__elementType[i]) {
+ // Supported entities
+ case 1: { // triangles
+ int& triangleNumber = elementNumber[1];
+
+ const int a = __verticesCorrepondance[elementVertices[0]];
+ const int b = __verticesCorrepondance[elementVertices[1]];
+ const int c = __verticesCorrepondance[elementVertices[2]];
+ if ((a<0) or (b<0) or (c<0)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': error reading element "+stringify(i)
+ +" [bad vertices definition]",
+ ErrorHandler::normal);
+ }
+ __triangles[triangleNumber]
+ = Triangle(a, b, c);
+ __triangles_ref[triangleNumber] = __references[i];
+ triangleNumber++; // one more triangle
+ break;
+ }
+ case 2: { // quadrangle
+ int& quadrilateralNumber = elementNumber[2];
+
+ const int a = __verticesCorrepondance[elementVertices[0]];
+ const int b = __verticesCorrepondance[elementVertices[1]];
+ const int c = __verticesCorrepondance[elementVertices[2]];
+ const int d = __verticesCorrepondance[elementVertices[3]];
+ if ((a<0) or (b<0) or (c<0) or (d<0)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': error reading element "+stringify(i)
+ +" [bad vertices definition]",
+ ErrorHandler::normal);
+ }
+ __quadrangles[quadrilateralNumber] = Quadrangle(a,b,c,d);
+ __quadrangles_ref[quadrilateralNumber] = __references[i];
+ quadrilateralNumber++; // one more quadrangle
+ break;
+ }
+ case 3: { // tetrahedra
+ int& tetrahedronNumber = elementNumber[3];
+
+ const int a = __verticesCorrepondance[elementVertices[0]];
+ const int b = __verticesCorrepondance[elementVertices[1]];
+ const int c = __verticesCorrepondance[elementVertices[2]];
+ const int d = __verticesCorrepondance[elementVertices[3]];
+ if ((a<0) or (b<0) or (c<0) or (d<0)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': error reading element "+stringify(i)
+ +" [bad vertices definition]",
+ ErrorHandler::normal);
+ }
+ __tetrahedra[tetrahedronNumber] = Tetrahedron(a, b, c, d);
+ __tetrahedra_ref[tetrahedronNumber] = __references[i];
+ tetrahedronNumber++; // one more tetrahedron
+ break;
+ }
+ case 4: { // hexaredron
+ int& hexahedronNumber = elementNumber[4];
+
+ const int a = __verticesCorrepondance[elementVertices[0]];
+ const int b = __verticesCorrepondance[elementVertices[1]];
+ const int c = __verticesCorrepondance[elementVertices[2]];
+ const int d = __verticesCorrepondance[elementVertices[3]];
+ const int e = __verticesCorrepondance[elementVertices[4]];
+ const int f = __verticesCorrepondance[elementVertices[5]];
+ const int g = __verticesCorrepondance[elementVertices[6]];
+ const int h = __verticesCorrepondance[elementVertices[7]];
+ if ((a<0) or (b<0) or (c<0) or (d<0) or (e<0) or (f<0) or (g<0) or (h<0)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': error reading element "+stringify(i)
+ +" [bad vertices definition]",
+ ErrorHandler::normal);
+ }
+ __hexahedra[hexahedronNumber]
+ = Hexahedron(a, b, c, d,
+ e, f, g, h);
+ __hexahedra_ref[hexahedronNumber] = __references[i];
+ hexahedronNumber++; // one more hexahedron
+ break; }
+ // Unsupported entities
+ case 0: // edge
+ case 14:{// point
+ break;
+ }
+ case 5: // prism
+ case 6: // pyramid
+ case 7: // second order edge
+ case 8: // second order triangle
+ case 9: // second order quadrangle
+ case 10: // second order tetrahedron
+ case 11: // second order hexahedron
+ case 12: // second order prism
+ case 13:{// second order pyramid
+ }
+ default: {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': ff3d cannot read this kind of element",
+ ErrorHandler::normal);
+ }
+ }
+ }
+
+ TinyVector<15,int> dimension0_mask = zero;
+ dimension0_mask[14]=1;
+ TinyVector<15,int> dimension1_mask = zero;
+ dimension1_mask[0]=1;
+ dimension1_mask[7]=1;
+ TinyVector<15,int> dimension2_mask = zero;
+ dimension2_mask[1]=1;
+ dimension2_mask[2]=1;
+ dimension2_mask[8]=1;
+ dimension2_mask[9]=1;
+ TinyVector<15,int> dimension3_mask = zero;
+ dimension3_mask[3]=1;
+ dimension3_mask[4]=1;
+ dimension3_mask[5]=1;
+ dimension3_mask[6]=1;
+ dimension3_mask[10]=1;
+ dimension3_mask[11]=1;
+ dimension3_mask[12]=1;
+ dimension3_mask[13]=1;
+
+ std::cout << "- dimension 0 entities : " << (dimension0_mask, elementNumber) << '\n';
+ std::cout << "- dimension 1 entities : " << (dimension1_mask, elementNumber) << '\n';
+ std::cout << "- dimension 2 entities : " << (dimension2_mask, elementNumber) << '\n';
+ std::cout << "- dimension 3 entities : " << (dimension3_mask, elementNumber) << '\n';
+ if ((dimension3_mask, elementNumber)>0) {
+ std::cerr << "*** using a 3d mesh (NIY)\n";
+ std::exit(0);
+ } else if ((dimension2_mask, elementNumber)>0) {
+ std::cerr << "*** using a 2d mesh (NIY)\n";
+ std::exit(0);
+ } else if ((dimension1_mask, elementNumber)>0) {
+ std::cerr << "*** using a 1d mesh (NIY)\n";
std::exit(0);
+ } else {
+ std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
+ std::exit(0);
+ }
+}
+
+GmshReader::Keyword
+GmshReader::__nextKeyword()
+{
+ GmshReader::Keyword kw;
+
+ std::cerr << "warning: " << rang::fg::red << __PRETTY_FUNCTION__ << rang::fg::reset << " keyword validity not checked!\n";
+ int retval = 1;
+ std::string aKeyword;
+ m_fin >> aKeyword;
+ if (retval < 0) {
+ kw.second = EndOfFile;
+ return kw;
+ } else if (retval == 0) {
+ kw.second = Unknown;
+ return kw;
+ }
+
+ KeywordList::iterator i = __keywordList.find(aKeyword.c_str());
+
+ if (i != __keywordList.end()) {
+ kw.first = (*i).first;
+ kw.second = (*i).second;
+ return kw;
+ }
+
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "reading file '"+m_filename
+ // +"': unknown keyword '"+aKeyword+"'",
+ // ErrorHandler::normal);
+
+ kw.first = aKeyword;
+ kw.second = Unknown;
+ return kw;
+}
+
+void GmshReader::
+__readGmshFormat2_2()
+{
+ std::cout << "- Reading Gmsh format 2.2\n";
+ GmshReader::Keyword kw;
+ kw = this->__nextKeyword();
+ if (kw.second != NODES) {
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "reading file '"+m_filename
+ // +"': expecting $Nodes, '"+kw.first+"' was found",
+ // ErrorHandler::normal);
+ }
+
+ this->__readVertices();
+
+ kw = this->__nextKeyword();
+ if (kw.second != ENDNODES) {
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "reading file '"+m_filename
+ // +"': expecting $EndNodes, '"+kw.first+"' was found",
+ // ErrorHandler::normal);
+ }
+
+ // Getting elements list
+ kw = this->__nextKeyword();
+ if (kw.second != ELEMENTS) {
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "reading file '"+m_filename
+ // +"': expecting $Elements, '"+kw.first+"' was found",
+ // ErrorHandler::normal);
+ }
+
+ this->__readElements2_2();
+
+ kw = this->__nextKeyword();
+ if (kw.second != ENDELEMENTS) {
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "reading file '"+m_filename
+ // +"': expecting $EndElements, '"+kw.first+"' was found",
+ // ErrorHandler::normal);
}
}
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index eafacfeaceaf055f5344542165beec55f9473623..61aba7d4f24d6e55afad6fdf9936b9737f34a750 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -2,9 +2,182 @@
#define GMSH_READER_HPP
#include <string>
+#include <array>
+#include <vector>
+#include <map>
+#include <fstream>
+#include <Kokkos_Core.hpp>
+#include <TinyVector.hpp>
class GmshReader
{
+public:
+ typedef TinyVector<3, double> R3;
+private:
+ std::ifstream m_fin;
+ const std::string m_filename;
+
+ int _getInteger()
+ {
+ int i;
+ m_fin >> i;
+
+ return std::move(i);
+ }
+
+ double _getReal()
+ {
+ double d;
+ m_fin >> d;
+
+ return std::move(d);
+ }
+
+ /**
+ * Gmsh format provides a numbered, none ordrered and none dense
+ * vertices list, this stores the number of read vertices.
+ */
+ std::vector<int> __verticesNumbers;
+
+ Kokkos::View<R3*> __vertices;
+ typedef TinyVector<3,unsigned int> Triangle;
+ Kokkos::View<Triangle*> __triangles;
+ std::vector<int> __triangles_ref;
+ typedef TinyVector<4,unsigned int> Quadrangle;
+ Kokkos::View<Quadrangle*> __quadrangles;
+ std::vector<int> __quadrangles_ref;
+ typedef TinyVector<4,unsigned int> Tetrahedron;
+ Kokkos::View<Tetrahedron*> __tetrahedra;
+ std::vector<int> __tetrahedra_ref;
+ typedef TinyVector<8,unsigned int> Hexahedron;
+ Kokkos::View<Hexahedron*> __hexahedra;
+ std::vector<int> __hexahedra_ref;
+
+ /**
+ * Gmsh format provides a numbered, none ordrered and none dense
+ * vertices list, this provides vertices renumbering correspondance
+ */
+ std::vector<int> __verticesCorrepondance;
+
+ /**
+ * elements types
+ */
+ std::vector<short> __elementType;
+
+ /**
+ * References
+ */
+ std::vector<int> __references;
+
+ /**
+ * References
+ */
+ std::vector<std::vector<int> > __elementVertices;
+
+ /**
+ * Stores the number of nodes associated to each primitive
+ *
+ */
+ std::vector<int> __numberOfPrimitiveNodes;
+
+ /**
+ * Array of boolean describing gmsh supported primitives by ff3d
+ *
+ */
+ std::array<bool, 15> __supportedPrimitives;
+
+ /**
+ * Primitives names according to there number
+ *
+ */
+ std::map<int, std::string> __primitivesNames;
+
+ bool __binary; /**< true if binary format */
+ bool __convertEndian; /**< true if needs to adapt endianess */
+
+ /**
+ * Adapts gmsh data to ff3d's structures
+ *
+ */
+ void __proceedData();
+
+ /**
+ * Reads data in format 2.2
+ *
+ */
+ void __readGmshFormat2_2();
+
+ /**
+ * List of allowed keyword in mesh file
+ *
+ */
+ enum KeywordType {
+ // gmsh 2.2
+ MESHFORMAT,
+ ENDMESHFORMAT,
+ NODES,
+ ENDNODES,
+ ELEMENTS,
+ ENDELEMENTS,
+
+ Unknown,
+ EndOfFile
+ };
+
+ /**
+ * List of known keywords
+ *
+ */
+ typedef std::map<std::string, int> KeywordList;
+
+ KeywordList __keywordList; /**< The keyword list */
+
+ /**
+ * Type for keyword
+ *
+ */
+ typedef std::pair<std::string, int> Keyword;
+
+ /**
+ * Skips next comments if exists
+ *
+ */
+ void __skipComments();
+
+ /**
+ * Reads the next keyword and returns its token
+ *
+ * @return KeywordToken
+ */
+ Keyword __nextKeyword();
+
+ /**
+ * get list of vertices
+ *
+ */
+ void __getVertices();
+
+ /**
+ * Read list of vertices
+ *
+ */
+ void __readVertices();
+
+ /**
+ * Read all elements in format 2.0
+ *
+ */
+ void __readElements2_2();
+
+ // /**
+ // * Common interface for writing references
+ // *
+ // * @param references the set of computed references
+ // * @param objectName the type of refernces
+ // */
+ // void __writeReferences(const std::set<size_t>& references,
+ // std::string objectName);
+
public:
GmshReader(const std::string& filename);
~GmshReader() = default;