#include <GmshReader.hpp>
#include <iostream>
#include <fstream>
#include <set>
#include <rang.hpp>
#include <CellType.hpp>
#include <Connectivity.hpp>
#include <Mesh.hpp>
#include <RefFaceList.hpp>
#include <map>
#include <regex>
#include <iomanip>
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() const;
/**
* 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() const
{
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;
__keywordList["$PhysicalNames"] = PHYSICALNAMES;
__keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;
__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]= true;
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(const 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::
__readPhysicalNames2_2()
{
const int number_of_names = this->_getInteger();
for (int i=0; i<number_of_names; ++i) {
const int physical_dimension = this->_getInteger();
const int physical_number = this->_getInteger();
std::string physical_name;
m_fin >> physical_name;
physical_name = std::regex_replace(physical_name, std::regex("(\")"), "");
PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
auto searched_physical_ref_id = m_physical_ref_map.find(physical_number);
if (searched_physical_ref_id != m_physical_ref_map.end()) {
std::cerr << "Physical reference id '" << physical_ref_id
<< "' already defined as '" << searched_physical_ref_id->second << "'!";
std::exit(0);
}
m_physical_ref_map[physical_number] = physical_ref_id;
}
}
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 << " [" << rang::fg::yellow << "IGNORED" << rang::style::reset << "]";
} 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
__points = Kokkos::View<Point*>("points",
elementNumber[i]);
__points_ref.resize(elementNumber[i]);
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
int& point_number = elementNumber[14];
const int a = __verticesCorrepondance[elementVertices[0]];
__points[point_number] = a;
__points_ref[point_number] = __references[i];
point_number++;
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) {
const size_t nb_cells = (dimension3_mask, elementNumber);
std::vector<CellType> cell_type_vector(nb_cells);
std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
const size_t nb_tetrahedra = __tetrahedra.extent(0);
for (size_t j=0; j<nb_tetrahedra; ++j) {
cell_by_node_vector[j].resize(4);
for (int r=0; r<4; ++r) {
cell_by_node_vector[j][r] = __tetrahedra[j][r];
}
cell_type_vector[j] = CellType::Tetrahedron;
}
const size_t nb_hexahedra = __hexahedra.extent(0);
for (size_t j=0; j<nb_hexahedra; ++j) {
const size_t jh = nb_tetrahedra+j;
cell_by_node_vector[jh].resize(8);
for (int r=0; r<8; ++r) {
cell_by_node_vector[jh][r] = __hexahedra[j][r];
}
cell_type_vector[jh] = CellType::Hexahedron;
}
std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector,
cell_type_vector));
Connectivity3D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int f=0; f<__triangles.extent(0); ++f) {
const unsigned int face_number
= connectivity.getFaceNumber({__triangles[f][0], __triangles[f][1], __triangles[f][2]});
const unsigned int& ref = __triangles_ref[f];
ref_faces_map[ref].push_back(face_number);
}
for (unsigned int f=0; f<__quadrangles.extent(0); ++f) {
const unsigned int face_number
= connectivity.getFaceNumber({__quadrangles[f][0], __quadrangles[f][1],
__quadrangles[f][2], __quadrangles[f][3]});
const unsigned int& ref = __quadrangles_ref[f];
ref_faces_map[ref].push_back(face_number);
}
for (const auto& ref_face_list : ref_faces_map) {
Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
for (size_t j=0; j<ref_face_list.second.size(); ++j) {
face_list[j]=ref_face_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
}
typedef Mesh<Connectivity3D> MeshType;
typedef TinyVector<3, 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];
xr[i][2] = __vertices[i][2];
}
m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
} else if ((dimension2_mask, elementNumber)>0) {
const size_t nb_cells = (dimension2_mask, elementNumber);
std::vector<CellType> cell_type_vector(nb_cells);
std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
const size_t nb_triangles = __triangles.extent(0);
for (size_t j=0; j<nb_triangles; ++j) {
cell_by_node_vector[j].resize(3);
for (int r=0; r<3; ++r) {
cell_by_node_vector[j][r] = __triangles[j][r];
}
cell_type_vector[j] = CellType::Triangle;
}
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_by_node_vector[jq].resize(4);
for (int r=0; r<4; ++r) {
cell_by_node_vector[jq][r] = __quadrangles[j][r];
}
cell_type_vector[jq] = CellType::Quadrangle;
}
std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector,
cell_type_vector));
Connectivity2D& connectivity = *p_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) {
Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
for (size_t j=0; j<ref_face_list.second.size(); ++j) {
face_list[j]=ref_face_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
}
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
for (unsigned int r=0; r<__points.extent(0); ++r) {
const unsigned int point_number = __points[r];
const unsigned int& ref = __points_ref[r];
ref_points_map[ref].push_back(point_number);
}
for (const auto& ref_point_list : ref_points_map) {
Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
for (size_t j=0; j<ref_point_list.second.size(); ++j) {
point_list[j]=ref_point_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_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 = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
} else if ((dimension1_mask, elementNumber)>0) {
const size_t nb_cells = (dimension1_mask, elementNumber);
std::vector<CellType> cell_type_vector(nb_cells);
std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
for (size_t j=0; j<nb_cells; ++j) {
cell_by_node_vector[j].resize(2);
for (int r=0; r<2; ++r) {
cell_by_node_vector[j][r] = __edges[j][r];
}
cell_type_vector[j] = CellType::Line;
}
std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector,
cell_type_vector));
Connectivity1D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
for (unsigned int r=0; r<__points.extent(0); ++r) {
const unsigned int point_number = __points[r];
const unsigned int& ref = __points_ref[r];
ref_points_map[ref].push_back(point_number);
}
for (const auto& ref_point_list : ref_points_map) {
Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
for (size_t j=0; j<ref_point_list.second.size(); ++j) {
point_list[j]=ref_point_list.second[j];
}
const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
}
typedef Mesh<Connectivity1D> MeshType;
typedef TinyVector<1, 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];
}
m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
} else {
std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
std::exit(0);
}
}
GmshReader::Keyword
GmshReader::__nextKeyword()
{
GmshReader::Keyword kw;
std::string aKeyword;
m_fin >> aKeyword;
if (not m_fin) {
kw.second = EndOfFile;
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 = std::make_pair("", Unknown);
while (kw.second != EndOfFile) {
kw = this->__nextKeyword();
switch (kw.second) {
case NODES: {
this->__readVertices();
if (this->__nextKeyword().second != ENDNODES) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': expecting $EndNodes, '"+kw.first+"' was found",
ErrorHandler::normal);
}
break;
}
case ELEMENTS: {
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);
}
break;
}
case PHYSICALNAMES: {
this->__readPhysicalNames2_2();
if (this->__nextKeyword().second != ENDPHYSICALNAMES) {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': expecting $EndNodes, '"+kw.first+"' was found",
ErrorHandler::normal);
}
break;
}
case EndOfFile: {
break;
}
default: {
throw ErrorHandler(__FILE__,__LINE__,
"reading file '"+m_filename
+"': unexpected '"+kw.first+"'",
ErrorHandler::normal);
}
}
}
}