#include <GmshReader.hpp>
#include <iostream>
#include <fstream>
#include <set>
#include <rang.hpp>
#include <Connectivity2D.hpp>
#include <Mesh.hpp>
#include <map>
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)
{
try {
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] = true;
__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();
}
catch(ErrorHandler e) {
e.writeErrorMessage();
std::exit(0);
}
}
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()
{
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 0: { // edges
__edges
= Kokkos::View<Edge*>("edges",
elementNumber[i]);
__edges_ref.resize(elementNumber[i]);
break;
}
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 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 0: { // edge
int& edgeNumber = elementNumber[0];
const int a = __verticesCorrepondance[elementVertices[0]];
const int b = __verticesCorrepondance[elementVertices[1]];
if ((a<0) or (b<0)) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': error reading element "+stringify(i)
+" [bad vertices definition]",
ErrorHandler::normal);
}
__edges[edgeNumber]
= Edge(a, b);
__edges_ref[edgeNumber] = __references[i];
edgeNumber++; // one more edge
break;
}
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 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) {
const size_t nb_cells = (dimension2_mask, elementNumber);
size_t max_nb_node_per_cell=3;
if (elementNumber[2] > 0) {
max_nb_node_per_cell = 4;
}
const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
const size_t nb_triangles = __triangles.extent(0);
for (size_t j=0; j<nb_triangles; ++j) {
cell_nodes(j,0) = __triangles[j][0];
cell_nodes(j,1) = __triangles[j][1];
cell_nodes(j,2) = __triangles[j][2];
}
const size_t nb_quadrangles = __quadrangles.extent(0);
for (size_t j=0; j<nb_quadrangles; ++j) {
const size_t jq = j+nb_triangles;
cell_nodes(jq,0) = __quadrangles[j][0];
cell_nodes(jq,1) = __quadrangles[j][1];
cell_nodes(jq,2) = __quadrangles[j][2];
cell_nodes(jq,3) = __quadrangles[j][3];
}
const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
for (size_t j=0; j<nb_triangles; ++j) {
cell_nb_nodes[j] = 3;
}
for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
cell_nb_nodes[j] = 4;
}
m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
Connectivity2D& connectivity = *m_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int e=0; e<__edges.extent(0); ++e) {
const unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
const unsigned int ref = __edges_ref[e];
ref_faces_map[ref].push_back(edge_number);
}
for (const auto& ref_face_list : ref_faces_map) {
connectivity.addFaceBoundary(ref_face_list.first,
ref_face_list.second);
}
const Kokkos::View<const unsigned int**> face_nodes = connectivity.faceNodes();
for (size_t i=0; i<connectivity.numberOfFaceBoundaries(); ++i) {
const Connectivity2D::FacesBoundary& faces_boundary
= connectivity.facesBoundary(i);
const unsigned int ref = faces_boundary.first;
std::set<unsigned int> node_id_set;
const Kokkos::View<const unsigned int*> face_ids
= faces_boundary.second;
for (unsigned int l=0; l<face_ids.extent(0); ++l) {
for (unsigned short r=0; r<2; ++r) {
node_id_set.insert(face_nodes(face_ids[l],r));
}
}
Kokkos::View<unsigned int*> node_id_list("node_ids", node_id_set.size());
{
int r=0;
for (auto node_id : node_id_set) {
node_id_list[r] = node_id;
++r;
}
}
connectivity.addNodeBoundary(ref, node_id_list);
}
typedef Mesh<Connectivity2D> MeshType;
typedef TinyVector<2, double> Rd;
Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
for (size_t i=0; i<__vertices.extent(0); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
}
m_mesh = new MeshType(connectivity, xr);
MeshType& mesh = *m_mesh;
std::ofstream gnuplot("mesh.gnu");
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
for (int r=0; r<mesh.connectivity().cellNbNodes()[j]; ++r) {
const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,r)];
gnuplot << x[0] << ' ' << x[1] << '\n';
}
const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,0)];
gnuplot << x[0] << ' ' << x[1] << "\n\n";
}
gnuplot.close();
} 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);
}
}