gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
GModelIO_NEU.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributed by Larry Price and Michael Ermakov
7 
8 #include <time.h>
9 #include <algorithm>
10 #include <iterator>
11 #include <limits>
12 #include <unordered_map>
13 #include "GModel.h"
14 #include "OS.h"
15 #include "MTriangle.h"
16 #include "MQuadrangle.h"
17 #include "MTetrahedron.h"
18 #include "MHexahedron.h"
19 #include "MPrism.h"
20 #include "MPyramid.h"
21 #include "GmshVersion.h"
22 
23 namespace {
24  // static const unsigned GAMBIT_TYPE_EDGE = 1;
25  // static const unsigned GAMBIT_TYPE_QUAD = 2;
26  // static const unsigned GAMBIT_TYPE_TRI = 3;
27  static const unsigned GAMBIT_TYPE_HEX = 4;
28  static const unsigned GAMBIT_TYPE_PRI = 5;
29  static const unsigned GAMBIT_TYPE_TET = 6;
30  static const unsigned GAMBIT_TYPE_PYR = 7;
31 
32  // This struct allows us to take advantage of C++11 unordered_map while
33  // maintaining backwards compatibility with C++03
34  template <typename Key, typename Value> struct hashMap {
35  typedef std::unordered_map<Key, Value> _;
36  };
37 
38  const hashMap<std::string, unsigned>::_::value_type rawData[] = {
39  hashMap<std::string, unsigned>::_::value_type("UNSPECIFIED", 0),
40  hashMap<std::string, unsigned>::_::value_type("AXIS", 1),
41  hashMap<std::string, unsigned>::_::value_type("CONJUGATE", 2),
42  hashMap<std::string, unsigned>::_::value_type("CONVECTION", 3),
43  hashMap<std::string, unsigned>::_::value_type("CYCLIC", 4),
44  hashMap<std::string, unsigned>::_::value_type("DEAD", 5),
45  hashMap<std::string, unsigned>::_::value_type("ELEMENT_SID", 6),
46  hashMap<std::string, unsigned>::_::value_type("ESPECIES", 7),
47  hashMap<std::string, unsigned>::_::value_type("EXHAUST_FAN", 8),
48  hashMap<std::string, unsigned>::_::value_type("FAN", 9),
49  hashMap<std::string, unsigned>::_::value_type("FREE_SURFACE", 10),
50  hashMap<std::string, unsigned>::_::value_type("GAP", 11),
51  hashMap<std::string, unsigned>::_::value_type("INFLOW", 12),
52  hashMap<std::string, unsigned>::_::value_type("INLET", 13),
53  hashMap<std::string, unsigned>::_::value_type("INLET_VENT", 14),
54  hashMap<std::string, unsigned>::_::value_type("INTAKE_FAN", 15),
55  hashMap<std::string, unsigned>::_::value_type("INTERFACE", 16),
56  hashMap<std::string, unsigned>::_::value_type("INTERIOR", 17),
57  hashMap<std::string, unsigned>::_::value_type("INTERNAL", 18),
58  hashMap<std::string, unsigned>::_::value_type("LIVE", 19),
59  hashMap<std::string, unsigned>::_::value_type("MASS_FLOW_INLET", 20),
60  hashMap<std::string, unsigned>::_::value_type("MELT", 21),
61  hashMap<std::string, unsigned>::_::value_type("MELT_INTERFACE", 22),
62  hashMap<std::string, unsigned>::_::value_type("MOVING_BOUNDARY", 23),
63  hashMap<std::string, unsigned>::_::value_type("NODE", 24),
64  hashMap<std::string, unsigned>::_::value_type("OUTFLOW", 25),
65  hashMap<std::string, unsigned>::_::value_type("OUTLET", 26),
66  hashMap<std::string, unsigned>::_::value_type("OUTLET_VENT", 27),
67  hashMap<std::string, unsigned>::_::value_type("PERIODIC", 28),
68  hashMap<std::string, unsigned>::_::value_type("PLOT", 29),
69  hashMap<std::string, unsigned>::_::value_type("POROUS", 30),
70  hashMap<std::string, unsigned>::_::value_type("POROUS_JUMP", 31),
71  hashMap<std::string, unsigned>::_::value_type("PRESSURE", 32),
72  hashMap<std::string, unsigned>::_::value_type("PRESSURE_FAR_FIELD", 33),
73  hashMap<std::string, unsigned>::_::value_type("PRESSURE_INFLOW", 34),
74  hashMap<std::string, unsigned>::_::value_type("PRESSURE_INLET", 35),
75  hashMap<std::string, unsigned>::_::value_type("PRESSURE_OUTFLOW", 36),
76  hashMap<std::string, unsigned>::_::value_type("PRESSURE_OUTLET", 37),
77  hashMap<std::string, unsigned>::_::value_type("RADIATION", 38),
78  hashMap<std::string, unsigned>::_::value_type("RADIATOR", 39),
79  hashMap<std::string, unsigned>::_::value_type("RECIRCULATION_INLET", 40),
80  hashMap<std::string, unsigned>::_::value_type("RECIRCULATION_OUTLET", 41),
81  hashMap<std::string, unsigned>::_::value_type("SLIP", 42),
82  hashMap<std::string, unsigned>::_::value_type("SREACTION", 43),
83  hashMap<std::string, unsigned>::_::value_type("SURFACE", 44),
84  hashMap<std::string, unsigned>::_::value_type("SYMMETRY", 45),
85  hashMap<std::string, unsigned>::_::value_type("TRACTION", 46),
86  hashMap<std::string, unsigned>::_::value_type("TRAJECTORY", 47),
87  hashMap<std::string, unsigned>::_::value_type("VELOCITY", 48),
88  hashMap<std::string, unsigned>::_::value_type("VELOCITY_INLET", 49),
89  hashMap<std::string, unsigned>::_::value_type("VENT", 50),
90  hashMap<std::string, unsigned>::_::value_type("WALL", 51),
91  hashMap<std::string, unsigned>::_::value_type("SPRING", 52),
92  };
93  static const hashMap<std::string, unsigned>::_
94  boundaryCodeMap(rawData, rawData + (sizeof rawData / sizeof rawData[0]));
95 
96  // Gambit numbers its faces slightly differently
97  static const unsigned GAMBIT_FACE_TET_MAP[4] = {1, 2, 4, 3};
98  static const unsigned GAMBIT_FACE_PYR_MAP[5] = {5, 2, 4, 3, 1};
99  static const unsigned GAMBIT_FACE_PRI_MAP[5] = {4, 5, 1, 3, 2};
100  static const unsigned GAMBIT_FACE_HEX_MAP[6] = {5, 1, 4, 2, 3, 6};
101 
102  unsigned const gambitBoundaryCode(std::string name)
103  {
104  std::transform(name.begin(), name.end(), name.begin(), ::toupper);
105  auto code = boundaryCodeMap.find(name);
106  return code == boundaryCodeMap.end() ? 0 : code->second;
107  }
108 
109  typedef std::pair<unsigned, unsigned> FacePair;
110  typedef hashMap<unsigned, std::vector<FacePair> >::_ IDFaceMap;
111 
112  bool const sortBCs(FacePair const &lhs, FacePair const &rhs)
113  {
114  return lhs.first < rhs.first;
115  }
116 
117 } // namespace
118 
119 // for a specification of the GAMBIT neutral format:
120 // http://web.stanford.edu/class/me469b/handouts/gambit_write.pdf
121 int GModel::writeNEU(const std::string &name, bool saveAll,
122  double scalingFactor)
123 {
124  FILE *fp = Fopen(name.c_str(), "w");
125  if(!fp) {
126  Msg::Error("Unable to open file '%s'", name.c_str());
127  return 0;
128  }
129 
130  // gather tetrahedra and id normalization information
131  unsigned numTetrahedra = 0;
132  unsigned numHexahedra = 0;
133  unsigned numPrisms = 0;
134  unsigned numPyramids = 0;
135 
136  unsigned lowestId = std::numeric_limits<int>::max();
137  hashMap<unsigned, std::vector<unsigned> >::_ elementGroups;
138 
139  for(auto it = firstRegion(); it != lastRegion(); ++it) {
140  if(saveAll || (*it)->physicals.size()) {
141  numTetrahedra += (*it)->tetrahedra.size();
142  numHexahedra += (*it)->hexahedra.size();
143  numPrisms += (*it)->prisms.size();
144  numPyramids += (*it)->pyramids.size();
145 
146  for(std::size_t phys = 0; phys < (*it)->physicals.size(); ++phys) {
147  for(std::size_t i = 0; i < (*it)->getNumMeshElements(); ++i) {
148  MElement *elem = (*it)->getMeshElement(i);
149  unsigned n = (unsigned)elem->getNum();
150  elementGroups[(*it)->physicals[phys]].push_back(n);
151  if(n < lowestId) lowestId = n - 1;
152  }
153  }
154  }
155  }
156 
157  if(lowestId == std::numeric_limits<int>::max()) {
158  Msg::Warning("No 3d-physicals");
159  return 1;
160  }
161 
162  unsigned numElements = numTetrahedra + numHexahedra + numPrisms + numPyramids;
163 
164  unsigned nVertex = indexMeshVertices(saveAll, 0, false);
165  std::vector<unsigned char> vertex_phys(nVertex);
166  // create association map for vertices and faces
167  hashMap<unsigned, std::vector<unsigned> >::_ vertmap;
168  for(auto it = firstFace(); it != lastFace(); ++it) {
169  for(auto &tri : (*it)->triangles) {
170  for(std::size_t j = 0; j < tri->getNumVertices(); ++j) {
171  unsigned numVertex = tri->getVertex(j)->getNum();
172  vertmap[numVertex].push_back(tri->getNum());
173  vertex_phys[numVertex] = 1;
174  }
175  }
176 
177  for(auto &quad : (*it)->quadrangles) {
178  for(std::size_t j = 0; j < quad->getNumVertices(); ++j) {
179  unsigned numVertex = quad->getVertex(j)->getNum();
180  vertmap[numVertex].push_back(quad->getNum());
181  vertex_phys[numVertex] = 1;
182  }
183  }
184  }
185 
186  // because we use a set_intersection later on, the vectors of vertmap should
187  // be sorted
188  for(auto &it : vertmap) { std::sort(it.second.begin(), it.second.end()); }
189 
190  std::vector<unsigned char> elem_phys(numElements);
191  for(auto it = firstRegion(); it != lastRegion(); ++it) {
192  if(saveAll || (*it)->physicals.size()) {
193  for(std::size_t i = 0; i < (*it)->getNumMeshElements(); ++i) {
194  unsigned num = (*it)->getMeshElement(i)->getNum();
195  unsigned sum = 0;
196  for(unsigned j = 0; j < (*it)->getMeshElement(i)->getNumVertices();
197  ++j) {
198  unsigned vertex = (*it)->getMeshElement(i)->getVertex(j)->getNum();
199 
200  if(vertex_phys[vertex]) { sum++; }
201  }
202 
203  if(sum >= 3) { elem_phys[num - lowestId - 1] = 1; }
204  }
205  }
206  }
207 
208  // determine which faces belong to which element by comparing vertices
209  IDFaceMap facemap;
210  for(auto it = firstRegion(); it != lastRegion(); ++it) {
211  for(std::size_t i = 0; i < (*it)->getNumMeshElements(); ++i) {
212  MElement *element = (*it)->getMeshElement(i);
213  unsigned num = element->getNum();
214  if(elem_phys[num - lowestId - 1]) {
215  unsigned numFaces = element->getNumFaces();
216 
217  for(unsigned faceNum = 0; faceNum < numFaces; ++faceNum) {
218  std::vector<MVertex *> verts;
219  element->getFaceVertices(faceNum, verts);
220 
221  std::vector<unsigned> current = vertmap[verts[0]->getNum()];
222  for(std::size_t j = 1; j < verts.size() && current.size() != 0; ++j) {
223  std::vector<unsigned> common_data;
224  set_intersection(current.begin(), current.end(),
225  vertmap[verts[j]->getNum()].begin(),
226  vertmap[verts[j]->getNum()].end(),
227  std::back_inserter(common_data));
228  current = common_data;
229  }
230 
231  if(current.size() == 1) {
232  unsigned type = element->getType();
233  unsigned face = 0;
234  switch(type) {
235  case TYPE_TET: face = GAMBIT_FACE_TET_MAP[faceNum]; break;
236  case TYPE_HEX: face = GAMBIT_FACE_HEX_MAP[faceNum]; break;
237  case TYPE_PRI: face = GAMBIT_FACE_PRI_MAP[faceNum]; break;
238  case TYPE_PYR: face = GAMBIT_FACE_PYR_MAP[faceNum]; break;
239  default: break;
240  }
241  facemap[current[0]].push_back(FacePair(num, face));
242  }
243  }
244  }
245  }
246  }
247 
248  // return 1;
249 
250  // populate boundary conditions for tetrahedra given triangle physicals
251  IDFaceMap boundaryConditions;
252  for(auto it = firstFace(); it != lastFace(); ++it) {
253  if((*it)->physicals.size()) {
254  for(std::size_t i = 0; i < (*it)->physicals.size(); ++i) {
255  unsigned phys = (*it)->physicals[i];
256 
257  for(std::size_t j = 0; j < (*it)->triangles.size(); ++j) {
258  MTriangle *tri = (*it)->triangles[j];
259  auto tets = facemap.find(tri->getNum());
260  if(tets != facemap.end()) {
261  for(std::size_t tet = 0; tet < tets->second.size(); ++tet) {
262  boundaryConditions[phys].push_back(tets->second[tet]);
263  }
264  }
265  }
266 
267  for(std::size_t j = 0; j < (*it)->quadrangles.size(); ++j) {
268  MQuadrangle *quad = (*it)->quadrangles[j];
269  auto tets = facemap.find(quad->getNum());
270  if(tets != facemap.end()) {
271  for(std::size_t tet = 0; tet < tets->second.size(); ++tet) {
272  boundaryConditions[phys].push_back(tets->second[tet]);
273  }
274  }
275  }
276  }
277  }
278  }
279 
280  // Metadata
281  fprintf(fp, " CONTROL INFO 2.0.0\n");
282  fprintf(fp, "** GAMBIT NEUTRAL FILE\n");
283  fprintf(fp, "Gmsh mesh in GAMBIT neutral file format\n");
284  fprintf(fp, "PROGRAM: Gmsh VERSION: %s\n", GMSH_VERSION);
285 
286  time_t rawtime;
287  time(&rawtime);
288  fprintf(fp, "%s", ctime(&rawtime));
289 
290  fprintf(fp, " NUMNP NELEM NGRPS NBSETS NDFCD NDFVL\n");
291 
292  fprintf(fp, " %9ld %9d %9lu %9lu %9d %9d\n",
293  indexMeshVertices(saveAll, 0, false), numElements,
294  elementGroups.size(), boundaryConditions.size(), getDim(), getDim());
295 
296  fprintf(fp, "ENDOFSECTION\n");
297 
298  // Nodes
299  fprintf(fp, " NODAL COORDINATES 2.0.0\n");
300  std::vector<GEntity *> entities;
301  getEntities(entities);
302  for(std::size_t i = 0; i < entities.size(); ++i) {
303  for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); ++j) {
304  entities[i]->mesh_vertices[j]->writeNEU(fp, getDim(), scalingFactor);
305  }
306  }
307  fprintf(fp, "ENDOFSECTION\n");
308 
309  // Elements
310  fprintf(fp, " ELEMENTS/CELLS 2.0.0\n");
311  for(auto it = firstRegion(); it != lastRegion(); ++it) {
312  std::size_t numPhys = (*it)->physicals.size();
313  if(saveAll || numPhys) {
314  for(auto cell : (*it)->tetrahedra) {
315  cell->writeNEU(fp, GAMBIT_TYPE_TET, lowestId,
316  numPhys ? (*it)->physicals[0] : 0);
317  }
318  for(auto cell : (*it)->hexahedra) {
319  cell->writeNEU(fp, GAMBIT_TYPE_HEX, lowestId,
320  numPhys ? (*it)->physicals[0] : 0);
321  }
322  for(auto cell : (*it)->prisms) {
323  cell->writeNEU(fp, GAMBIT_TYPE_PRI, lowestId,
324  numPhys ? (*it)->physicals[0] : 0);
325  }
326  for(auto cell : (*it)->pyramids) {
327  cell->writeNEU(fp, GAMBIT_TYPE_PYR, lowestId,
328  numPhys ? (*it)->physicals[0] : 0);
329  }
330  }
331  }
332  fprintf(fp, "ENDOFSECTION\n");
333 
334  // Element Groups
335 
336  for(auto it = elementGroups.begin(); it != elementGroups.end(); ++it) {
337  fprintf(fp, " ELEMENT GROUP 2.0.0\n");
338  fprintf(fp,
339  "GROUP: %10d ELEMENTS: %10lu MATERIAL: 0 NFLAGS: %10d\n",
340  it->first, it->second.size(), 1);
341 
342  std::string volumeName = getPhysicalName(3, it->first);
343  if(volumeName.empty()) {
344  char tmp[256];
345  sprintf(tmp, "Material group %d", it->first);
346  volumeName = tmp;
347  }
348  fprintf(fp, "%32s\n", volumeName.c_str());
349 
350  fprintf(fp, " 0");
351 
352  for(unsigned i = 0; i < it->second.size(); ++i) {
353  if(i % 10 == 0) { fprintf(fp, "\n"); }
354  fprintf(fp, "%8d", it->second[i] - lowestId);
355  }
356  fprintf(fp, "\n");
357  fprintf(fp, "ENDOFSECTION\n");
358  }
359 
360  // Boundary Conditions
361 
362  unsigned nphys = getMaxPhysicalNumber(2);
363  for(unsigned iphys = 1; iphys <= nphys; ++iphys) {
364  for(auto it = boundaryConditions.begin(); it != boundaryConditions.end();
365  ++it) {
366  if(it->first == iphys) {
367  fprintf(fp, " BOUNDARY CONDITIONS 2.0.0\n");
368  std::string regionName = getPhysicalName(2, it->first);
369  if(regionName.empty()) {
370  char tmp[256];
371  sprintf(tmp, "PhysicalSurface%d", it->first);
372  regionName = tmp;
373  }
374 
375  fprintf(fp, "%32s%8d%8lu%8d%8d\n", regionName.c_str(), 1,
376  it->second.size(), 0, gambitBoundaryCode(regionName));
377  std::sort(it->second.begin(), it->second.end(), sortBCs);
378  for(auto tfp = it->second.begin(); tfp != it->second.end(); ++tfp) {
379  MElement *element = getMeshElementByTag(tfp->first);
380  unsigned type = element->getType();
381  unsigned gambit_type = 0;
382  switch(type) {
383  case TYPE_TET: gambit_type = GAMBIT_TYPE_TET; break;
384  case TYPE_PYR: gambit_type = GAMBIT_TYPE_PYR; break;
385  case TYPE_PRI: gambit_type = GAMBIT_TYPE_PRI; break;
386  case TYPE_HEX: gambit_type = GAMBIT_TYPE_HEX; break;
387  default: break;
388  }
389 
390  fprintf(fp, "%10d %5d %5d\n", tfp->first - lowestId, gambit_type,
391  tfp->second);
392  }
393  }
394  }
395 
396  fprintf(fp, "ENDOFSECTION\n");
397  }
398 
399  fclose(fp);
400  return 1;
401 }
MElement::getNum
virtual std::size_t getNum() const
Definition: MElement.h:68
MTriangle.h
GModel::writeNEU
int writeNEU(const std::string &name, bool saveAll, double scalingFactor)
Definition: GModelIO_NEU.cpp:121
OS.h
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
transform
static SPoint3 transform(MVertex *vsource, const std::vector< double > &tfo)
Definition: Generator.cpp:1347
TYPE_PRI
#define TYPE_PRI
Definition: GmshDefines.h:70
Fopen
FILE * Fopen(const char *f, const char *mode)
Definition: OS.cpp:273
GModel::lastFace
fiter lastFace()
Definition: GModel.h:359
GModel::getPhysicalName
std::string getPhysicalName(int dim, int num) const
Definition: GModel.cpp:961
MHexahedron.h
GModel::getMaxPhysicalNumber
int getMaxPhysicalNumber(int dim)
Definition: GModel.cpp:917
GModel::getDim
int getDim() const
Definition: GModel.cpp:989
MPyramid.h
GModel::firstFace
fiter firstFace()
Definition: GModel.h:355
MElement
Definition: MElement.h:30
GModel::getMeshElementByTag
MElement * getMeshElementByTag(int n)
Definition: GModel.h:538
GModel::getEntities
void getEntities(std::vector< GEntity * > &entities, int dim=-1) const
Definition: GModel.cpp:651
element
Definition: shapeFunctions.h:12
MTriangle
Definition: MTriangle.h:26
TYPE_PYR
#define TYPE_PYR
Definition: GmshDefines.h:69
GModel::firstRegion
riter firstRegion()
Definition: GModel.h:354
GModel::lastRegion
riter lastRegion()
Definition: GModel.h:358
MTetrahedron.h
MQuadrangle.h
TYPE_HEX
#define TYPE_HEX
Definition: GmshDefines.h:71
MPrism.h
TYPE_TET
#define TYPE_TET
Definition: GmshDefines.h:68
GModel.h
MQuadrangle
Definition: MQuadrangle.h:26
GModel::indexMeshVertices
std::size_t indexMeshVertices(bool all, int singlePartition=0, bool renumber=true)
Definition: GModel.cpp:2135
GModel::current
static GModel * current(int index=-1)
Definition: GModel.cpp:136
elem
Definition: OctreeInternals.h:17