1 // Gmsh - Copyright (C) 1997-2021 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 #include <sstream>
7 #include <iomanip>
8 #include <ctime>
9 #include "GModel.h"
10 #include "OS.h"
11 #include "Context.h"
12 #include "GmshMessage.h"
13 #include "MElement.h"
14 #include "MPoint.h"
15 #include "MLine.h"
16 #include "MTriangle.h"
17 #include "MQuadrangle.h"
18 #include "MTetrahedron.h"
19 #include "MHexahedron.h"
20 #include "MPrism.h"
21 #include "MPyramid.h"
22 #include "MTrihedron.h"
23 #include "StringUtils.h"
24 #include "discreteVertex.h"
25 #include "discreteEdge.h"
26 #include "discreteFace.h"
27 #include "discreteRegion.h"
28 
readMSHPhysicals(FILE * fp,GEntity * ge)29 static int readMSHPhysicals(FILE *fp, GEntity *ge)
30 {
31   int nump;
32   if(fscanf(fp, "%d", &nump) != 1) return 0;
33   for(int i = 0; i < nump; i++) {
34     int tag;
35     if(fscanf(fp, "%d", &tag) != 1) return 0;
36     ge->physicals.push_back(tag);
37   }
38   return 1;
39 }
40 
readMSHEntities(FILE * fp,GModel * gm)41 static void readMSHEntities(FILE *fp, GModel *gm)
42 {
43   int nv, ne, nf, nr;
44   int tag;
45   if(fscanf(fp, "%d %d %d %d", &nv, &ne, &nf, &nr) != 4) return;
46   for(int i = 0; i < nv; i++) {
47     if(fscanf(fp, "%d", &tag) != 1) return;
48     GVertex *gv = gm->getVertexByTag(tag);
49     if(!gv) {
50       gv = new discreteVertex(gm, tag);
51       gm->add(gv);
52     }
53     if(!readMSHPhysicals(fp, gv)) return;
54   }
55   for(int i = 0; i < ne; i++) {
56     int n;
57     if(fscanf(fp, "%d %d", &tag, &n) != 2) return;
58     GEdge *ge = gm->getEdgeByTag(tag);
59     if(!ge) {
60       GVertex *v1 = nullptr, *v2 = nullptr;
61       for(int j = 0; j < n; j++) {
62         int tagv;
63         if(fscanf(fp, "%d", &tagv) != 1) {
64           delete ge;
65           return;
66         }
67         GVertex *v = gm->getVertexByTag(tagv);
68         if(!v) Msg::Error("Unknown GVertex %d", tagv);
69         if(j == 0) v1 = v;
70         if(j == 1) v2 = v;
71       }
72       ge = new discreteEdge(gm, tag, v1, v2);
73       gm->add(ge);
74     }
75     if(!readMSHPhysicals(fp, ge)) return;
76   }
77   for(int i = 0; i < nf; i++) {
78     int n;
79     if(fscanf(fp, "%d %d", &tag, &n) != 2) return;
80     GFace *gf = gm->getFaceByTag(tag);
81     if(!gf) {
82       discreteFace *df = new discreteFace(gm, tag);
83       std::vector<int> edges, signs;
84       for(int j = 0; j < n; j++) {
85         int tage;
86         if(fscanf(fp, "%d", &tage) != 1) {
87           delete df;
88           return;
89         }
90         edges.push_back(std::abs(tage));
91         int sign = tage > 0 ? 1 : -1;
92         signs.push_back(sign);
93       }
94       df->setBoundEdges(edges, signs);
95       gm->add(df);
96       gf = df;
97     }
98     if(!readMSHPhysicals(fp, gf)) return;
99   }
100   for(int i = 0; i < nr; i++) {
101     int n;
102     if(fscanf(fp, "%d %d", &tag, &n) != 2) return;
103     GRegion *gr = gm->getRegionByTag(tag);
104     if(!gr) {
105       discreteRegion *dr = new discreteRegion(gm, tag);
106       std::vector<int> faces, signs;
107       for(int j = 0; j < n; j++) {
108         int tagf;
109         if(fscanf(fp, "%d", &tagf) != 1) {
110           delete dr;
111           return;
112         }
113         faces.push_back(std::abs(tagf));
114         int sign = tagf > 0 ? 1 : -1;
115         signs.push_back(sign);
116       }
117       dr->setBoundFaces(faces, signs);
118       gm->add(dr);
119       gr = dr;
120     }
121     if(!readMSHPhysicals(fp, gr)) return;
122   }
123 }
124 
readMSHPeriodicNodes(FILE * fp,GModel * gm)125 void readMSHPeriodicNodes(FILE *fp, GModel *gm) // also used in MSH2
126 {
127   int count;
128   if(fscanf(fp, "%d", &count) != 1) return;
129   for(int i = 0; i < count; i++) {
130     int dim, slave, master;
131 
132     if(fscanf(fp, "%d %d %d", &dim, &slave, &master) != 3) continue;
133 
134     GEntity *s = nullptr, *m = nullptr;
135     switch(dim) {
136     case 0:
137       s = gm->getVertexByTag(slave);
138       m = gm->getVertexByTag(master);
139       break;
140     case 1:
141       s = gm->getEdgeByTag(slave);
142       m = gm->getEdgeByTag(master);
143       break;
144     case 2:
145       s = gm->getFaceByTag(slave);
146       m = gm->getFaceByTag(master);
147       break;
148     }
149 
150     // we need to continue parsing, otherwise we end up reading on the wrong
151     // position
152 
153     bool completePer = s && m;
154 
155     char token[7];
156     fpos_t pos;
157     fgetpos(fp, &pos);
158     if(fscanf(fp, "%s", token) != 1) return;
159     if(strcmp(token, "Affine") == 0) {
160       std::vector<double> tfo(16);
161       for(int i = 0; i < 16; i++) {
162         if(fscanf(fp, "%lf", &tfo[i]) != 1) return;
163       }
164       if(completePer) s->setMeshMaster(m, tfo);
165     }
166     else {
167       fsetpos(fp, &pos);
168       if(completePer) s->setMeshMaster(m);
169     }
170     int numv;
171     if(fscanf(fp, "%d", &numv) != 1) numv = 0;
172     for(int j = 0; j < numv; j++) {
173       int v1, v2;
174       if(fscanf(fp, "%d %d", &v1, &v2) != 2) continue;
175       MVertex *mv1 = gm->getMeshVertexByTag(v1);
176       MVertex *mv2 = gm->getMeshVertexByTag(v2);
177       if(completePer) s->correspondingVertices[mv1] = mv2;
178     }
179     if(!completePer) {
180       if(!s)
181         Msg::Info("Could not find periodic slave entity %d of dimension %d",
182                   slave, dim);
183       if(!m)
184         Msg::Info("Could not find periodic master entity %d of dimension %d",
185                   master, dim);
186     }
187   }
188 }
189 
_readMSH3(const std::string & name)190 int GModel::_readMSH3(const std::string &name)
191 {
192   FILE *fp = Fopen(name.c_str(), "rb");
193   if(!fp) {
194     Msg::Error("Unable to open file '%s'", name.c_str());
195     return 0;
196   }
197 
198   char str[256] = "";
199   double version = 0.;
200   bool binary = false, swap = false, postpro = false;
201   int minVertex = 0;
202   std::map<int, std::vector<MElement *> > elements[11];
203   std::size_t oldNumPartitions = getNumPartitions();
204 
205   while(1) {
206     while(str[0] != '$') {
207       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
208     }
209 
210     if(feof(fp)) break;
211 
212     // $MeshFormat section
213     if(!strncmp(&str[1], "MeshFormat", 10)) {
214       if(!fgets(str, sizeof(str), fp)) {
215         fclose(fp);
216         return 0;
217       }
218       int format, size;
219       if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) {
220         fclose(fp);
221         return 0;
222       }
223       if(version < 3. || version >= 4.) {
224         Msg::Error("Wrong MSH file version %g for MSH3 reader", version);
225         fclose(fp);
226         return 0;
227       }
228       if(format) {
229         binary = true;
230         Msg::Debug("Mesh is in binary format");
231         int one;
232         if(fread(&one, sizeof(int), 1, fp) != 1) {
233           fclose(fp);
234           return 0;
235         }
236         if(one != 1) {
237           swap = true;
238           Msg::Debug("Swapping bytes from binary file");
239         }
240       }
241     }
242 
243     // $PhysicalNames section
244     else if(!strncmp(&str[1], "PhysicalNames", 13)) {
245       if(!fgets(str, sizeof(str), fp)) {
246         fclose(fp);
247         return 0;
248       }
249       int numNames;
250       if(sscanf(str, "%d", &numNames) != 1) {
251         fclose(fp);
252         return 0;
253       }
254       for(int i = 0; i < numNames; i++) {
255         int dim, num;
256         if(fscanf(fp, "%d", &dim) != 1) {
257           fclose(fp);
258           return 0;
259         }
260         if(fscanf(fp, "%d", &num) != 1) {
261           fclose(fp);
262           return 0;
263         }
264         if(!fgets(str, sizeof(str), fp)) {
265           fclose(fp);
266           return 0;
267         }
268         std::string name = ExtractDoubleQuotedString(str, 256);
269         if(name.size()) setPhysicalName(name, dim, num);
270       }
271     }
272 
273     // $Entities section
274     else if(!strncmp(&str[1], "Entities", 8)) {
275       readMSHEntities(fp, this);
276     }
277 
278     // $Nodes section
279     else if(!strncmp(&str[1], "Nodes", 5)) {
280       if(!fgets(str, sizeof(str), fp)) {
281         fclose(fp);
282         return 0;
283       }
284       int numVertices;
285       if(sscanf(str, "%d", &numVertices) != 1) {
286         fclose(fp);
287         return 0;
288       }
289       Msg::Info("%d nodes", numVertices);
290       Msg::StartProgressMeter(numVertices);
291       _vertexMapCache.clear();
292       _vertexVectorCache.clear();
293       int maxVertex = -1;
294       minVertex = numVertices + 1;
295       for(int i = 0; i < numVertices; i++) {
296         int num, entity, dim;
297         double xyz[3];
298         MVertex *vertex = nullptr;
299         if(!binary) {
300           if(fscanf(fp, "%d %lf %lf %lf %d", &num, &xyz[0], &xyz[1], &xyz[2],
301                     &entity) != 5) {
302             fclose(fp);
303             return 0;
304           }
305         }
306         else {
307           if(fread(&num, sizeof(int), 1, fp) != 1) {
308             fclose(fp);
309             return 0;
310           }
311           if(swap) SwapBytes((char *)&num, sizeof(int), 1);
312           if(fread(xyz, sizeof(double), 3, fp) != 3) {
313             fclose(fp);
314             return 0;
315           }
316           if(swap) SwapBytes((char *)xyz, sizeof(double), 3);
317           if(fread(&entity, sizeof(int), 1, fp) != 1) {
318             fclose(fp);
319             return 0;
320           }
321           if(swap) SwapBytes((char *)&entity, sizeof(int), 1);
322         }
323         if(!entity) {
324           vertex = new MVertex(xyz[0], xyz[1], xyz[2], nullptr, num);
325         }
326         else {
327           if(!binary) {
328             if(fscanf(fp, "%d", &dim) != 1) {
329               fclose(fp);
330               return 0;
331             }
332           }
333           else {
334             if(fread(&dim, sizeof(int), 1, fp) != 1) {
335               fclose(fp);
336               return 0;
337             }
338             if(swap) SwapBytes((char *)&dim, sizeof(int), 1);
339           }
340           switch(dim) {
341           case 0: {
342             GVertex *gv = getVertexByTag(entity);
343             // FIXME -- cannot call this: it destroys _vertexMapCache
344             // if(gv) gv->deleteMesh();
345             vertex = new MVertex(xyz[0], xyz[1], xyz[2], gv, num);
346           } break;
347           case 1: {
348             GEdge *ge = getEdgeByTag(entity);
349             double u;
350             if(!binary) {
351               if(fscanf(fp, "%lf", &u) != 1) {
352                 fclose(fp);
353                 return 0;
354               }
355             }
356             else {
357               if(fread(&u, sizeof(double), 1, fp) != 1) {
358                 fclose(fp);
359                 return 0;
360               }
361               if(swap) SwapBytes((char *)&u, sizeof(double), 1);
362             }
363             vertex = new MEdgeVertex(xyz[0], xyz[1], xyz[2], ge, u, num);
364           } break;
365           case 2: {
366             GFace *gf = getFaceByTag(entity);
367             double uv[2];
368             if(!binary) {
369               if(fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) {
370                 fclose(fp);
371                 return 0;
372               }
373             }
374             else {
375               if(fread(uv, sizeof(double), 2, fp) != 2) {
376                 fclose(fp);
377                 return 0;
378               }
379               if(swap) SwapBytes((char *)uv, sizeof(double), 2);
380             }
381             vertex =
382               new MFaceVertex(xyz[0], xyz[1], xyz[2], gf, uv[0], uv[1], num);
383           } break;
384           case 3: {
385             GRegion *gr = getRegionByTag(entity);
386             double uvw[3];
387             if(!binary) {
388               if(fscanf(fp, "%lf %lf %lf", &uvw[0], &uvw[1], &uvw[2]) != 3) {
389                 fclose(fp);
390                 return 0;
391               }
392             }
393             else {
394               if(fread(uvw, sizeof(double), 3, fp) != 3) {
395                 fclose(fp);
396                 return 0;
397               }
398               if(swap) SwapBytes((char *)uvw, sizeof(double), 3);
399             }
400             vertex = new MVertex(xyz[0], xyz[1], xyz[2], gr, num);
401           } break;
402           default:
403             Msg::Error("Wrong entity dimension for node %d", num);
404             fclose(fp);
405             return 0;
406           }
407         }
408         minVertex = std::min(minVertex, num);
409         maxVertex = std::max(maxVertex, num);
410         if(_vertexMapCache.count(num))
411           Msg::Warning("Skipping duplicate node %d", num);
412         _vertexMapCache[num] = vertex;
413         if(numVertices > 100000)
414           Msg::ProgressMeter(i + 1, true, "Reading nodes");
415       }
416       Msg::StopProgressMeter();
417       // if the vertex numbering is dense, transfer the map into a vector to
418       // speed up element creation
419       if((int)_vertexMapCache.size() == numVertices &&
420          ((minVertex == 1 && maxVertex == numVertices) ||
421           (minVertex == 0 && maxVertex == numVertices - 1))) {
422         Msg::Debug("Vertex numbering is dense");
423         _vertexVectorCache.resize(_vertexMapCache.size() + 1);
424         if(minVertex == 1)
425           _vertexVectorCache[0] = nullptr;
426         else
427           _vertexVectorCache[numVertices] = nullptr;
428         for(auto it = _vertexMapCache.begin(); it != _vertexMapCache.end(); ++it)
429           _vertexVectorCache[it->first] = it->second;
430         _vertexMapCache.clear();
431       }
432     }
433 
434     // $Elements section
435     else if(!strncmp(&str[1], "Elements", 8)) {
436       if(!fgets(str, sizeof(str), fp)) {
437         fclose(fp);
438         return 0;
439       }
440       int numElements;
441       if(sscanf(str, "%d", &numElements) != 1) {
442         fclose(fp);
443         return 0;
444       }
445       Msg::Info("%d elements", numElements);
446       Msg::StartProgressMeter(numElements);
447       for(int i = 0; i < numElements; i++) {
448         int num, type, entity, numData;
449         if(!binary) {
450           if(fscanf(fp, "%d %d %d %d", &num, &type, &entity, &numData) != 4) {
451             fclose(fp);
452             return 0;
453           }
454         }
455         else {
456           if(fread(&num, sizeof(int), 1, fp) != 1) {
457             fclose(fp);
458             return 0;
459           }
460           if(swap) SwapBytes((char *)&num, sizeof(int), 1);
461           if(fread(&type, sizeof(int), 1, fp) != 1) {
462             fclose(fp);
463             return 0;
464           }
465           if(swap) SwapBytes((char *)&type, sizeof(int), 1);
466           if(fread(&entity, sizeof(int), 1, fp) != 1) {
467             fclose(fp);
468             return 0;
469           }
470           if(swap) SwapBytes((char *)&entity, sizeof(int), 1);
471           if(fread(&numData, sizeof(int), 1, fp) != 1) {
472             fclose(fp);
473             return 0;
474           }
475           if(swap) SwapBytes((char *)&numData, sizeof(int), 1);
476         }
477         std::vector<int> data;
478         if(numData > 0) {
479           data.resize(numData);
480           if(!binary) {
481             for(int j = 0; j < numData; j++) {
482               if(fscanf(fp, "%d", &data[j]) != 1) {
483                 fclose(fp);
484                 return 0;
485               }
486             }
487           }
488           else {
489             if((int)fread(&data[0], sizeof(int), numData, fp) != numData) {
490               fclose(fp);
491               return 0;
492             }
493             if(swap) SwapBytes((char *)&data[0], sizeof(int), numData);
494           }
495         }
496         MElementFactory f;
497         MElement *element = f.create(num, type, data, this);
498         if(!element) {
499           fclose(fp);
500           return 0;
501         }
502         switch(element->getType()) {
503         case TYPE_PNT: elements[0][entity].push_back(element); break;
504         case TYPE_LIN: elements[1][entity].push_back(element); break;
505         case TYPE_TRI: elements[2][entity].push_back(element); break;
506         case TYPE_QUA: elements[3][entity].push_back(element); break;
507         case TYPE_TET: elements[4][entity].push_back(element); break;
508         case TYPE_HEX: elements[5][entity].push_back(element); break;
509         case TYPE_PRI: elements[6][entity].push_back(element); break;
510         case TYPE_PYR: elements[7][entity].push_back(element); break;
511         case TYPE_TRIH: elements[8][entity].push_back(element); break;
512         case TYPE_POLYG: elements[9][entity].push_back(element); break;
513         case TYPE_POLYH: elements[10][entity].push_back(element); break;
514         }
515         if(numElements > 100000)
516           Msg::ProgressMeter(i + 1, true, "Reading elements");
517       }
518       Msg::StopProgressMeter();
519     }
520 
521     // $Periodical section
522     else if(!strncmp(&str[1], "Periodic", 8)) {
523       readMSHPeriodicNodes(fp, this);
524     }
525 
526     // Post-processing sections
527     else if(!strncmp(&str[1], "NodeData", 8) ||
528             !strncmp(&str[1], "ElementData", 11) ||
529             !strncmp(&str[1], "ElementNodeData", 15)) {
530       postpro = true;
531       break;
532     }
533 
534     do {
535       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
536     } while(str[0] != '$');
537   }
538 
539   // store the elements in their associated elementary entity. If the
540   // entity does not exist, create a new (discrete) one.
541   for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
542     _storeElementsInEntities(elements[i]);
543 
544   // associate the correct geometrical entity with each mesh vertex
545   _associateEntityWithMeshVertices();
546 
547   // store the vertices in their associated geometrical entity
548   if(_vertexVectorCache.size())
549     _storeVerticesInEntities(_vertexVectorCache);
550   else
551     _storeVerticesInEntities(_vertexMapCache);
552 
553   for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
554     _storeParentsInSubElements(elements[i]);
555 
556   fclose(fp);
557 
558   // create new partition entities if the mesh is partitioned
559   if(CTX::instance()->mesh.partitionConvertMsh2 &&
560      getNumPartitions() > oldNumPartitions)
561     convertOldPartitioningToNewOne();
562 
563   return postpro ? 2 : 1;
564 }
565 
writeMSHPhysicals(FILE * fp,GEntity * ge)566 static void writeMSHPhysicals(FILE *fp, GEntity *ge)
567 {
568   std::vector<int> phys = ge->physicals;
569   fprintf(fp, "%d ", (int)phys.size());
570   for(auto itp = phys.begin(); itp != phys.end(); itp++)
571     fprintf(fp, "%d ", *itp);
572 }
573 
writeMSHEntities(FILE * fp,GModel * gm)574 void writeMSHEntities(FILE *fp, GModel *gm) // also used in MSH2
575 {
576   fprintf(fp, "$Entities\n");
577   fprintf(fp, "%lu %lu %lu %lu\n", gm->getNumVertices(), gm->getNumEdges(),
578           gm->getNumFaces(), gm->getNumRegions());
579   for(auto it = gm->firstVertex(); it != gm->lastVertex(); ++it) {
580     fprintf(fp, "%d ", (*it)->tag());
581     writeMSHPhysicals(fp, *it);
582     fprintf(fp, "\n");
583   }
584   for(auto it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
585     std::list<GVertex *> vertices;
586     if((*it)->getBeginVertex()) vertices.push_back((*it)->getBeginVertex());
587     if((*it)->getEndVertex()) vertices.push_back((*it)->getEndVertex());
588     fprintf(fp, "%d %d ", (*it)->tag(), (int)vertices.size());
589     for(auto itv = vertices.begin(); itv != vertices.end(); itv++) {
590       fprintf(fp, "%d ", (*itv)->tag());
591     }
592     writeMSHPhysicals(fp, *it);
593     fprintf(fp, "\n");
594   }
595   for(auto it = gm->firstFace(); it != gm->lastFace(); ++it) {
596     std::vector<GEdge *> const &edges = (*it)->edges();
597     std::vector<int> const &ori = (*it)->edgeOrientations();
598     fprintf(fp, "%d %d ", (*it)->tag(), (int)edges.size());
599     std::vector<int> tags;
600     for(auto ite = edges.begin(); ite != edges.end(); ite++)
601       tags.push_back((*ite)->tag());
602 
603     std::vector<int> signs(ori.begin(), ori.end());
604 
605     if(tags.size() == signs.size()) {
606       for(std::size_t i = 0; i < tags.size(); i++)
607         tags[i] *= (signs[i] > 0 ? 1 : -1);
608     }
609     for(std::size_t i = 0; i < tags.size(); i++) fprintf(fp, "%d ", tags[i]);
610     writeMSHPhysicals(fp, *it);
611     fprintf(fp, "\n");
612   }
613   for(auto it = gm->firstRegion(); it != gm->lastRegion(); ++it) {
614     std::vector<GFace *> faces = (*it)->faces();
615     std::vector<int> ori = (*it)->faceOrientations();
616     fprintf(fp, "%d %d ", (*it)->tag(), (int)faces.size());
617     std::vector<int> tags, signs;
618     for(auto itf = faces.begin(); itf != faces.end(); itf++)
619       tags.push_back((*itf)->tag());
620     for(auto itf = ori.begin(); itf != ori.end(); itf++)
621       signs.push_back(*itf);
622     if(tags.size() == signs.size()) {
623       for(std::size_t i = 0; i < tags.size(); i++)
624         tags[i] *= (signs[i] > 0 ? 1 : -1);
625     }
626     for(std::size_t i = 0; i < tags.size(); i++) fprintf(fp, "%d ", tags[i]);
627     writeMSHPhysicals(fp, *it);
628     fprintf(fp, "\n");
629   }
630   fprintf(fp, "$EndEntities\n");
631 }
632 
getNumElementsMSH(GEntity * ge,bool saveAll,int saveSinglePartition)633 static int getNumElementsMSH(GEntity *ge, bool saveAll, int saveSinglePartition)
634 {
635   int n = 0;
636   if(saveAll || ge->physicals.size()) {
637     if(saveSinglePartition <= 0)
638       n = ge->getNumMeshElements();
639     else
640       for(std::size_t i = 0; i < ge->getNumMeshElements(); i++)
641         if(ge->getMeshElement(i)->getPartition() == saveSinglePartition) n++;
642   }
643   return n;
644 }
645 
writeElementMSH(FILE * fp,GModel * model,MElement * ele,bool binary,int elementary)646 static void writeElementMSH(FILE *fp, GModel *model, MElement *ele, bool binary,
647                             int elementary)
648 {
649   if(model->getGhostCells().size()) {
650     std::vector<short> ghosts;
651     auto itp = model->getGhostCells().equal_range(ele);
652     for(auto it = itp.first; it != itp.second; it++)
653       ghosts.push_back(it->second);
654     ele->writeMSH3(fp, binary, elementary, &ghosts);
655   }
656   else
657     ele->writeMSH3(fp, binary, elementary);
658 }
659 
660 template <class T>
writeElementsMSH(FILE * fp,GModel * model,GEntity * ge,std::vector<T * > & ele,bool saveAll,int saveSinglePartition,bool binary)661 static void writeElementsMSH(FILE *fp, GModel *model, GEntity *ge,
662                              std::vector<T *> &ele, bool saveAll,
663                              int saveSinglePartition, bool binary)
664 {
665   if(saveAll || ge->physicals.size()) {
666     for(std::size_t i = 0; i < ele.size(); i++) {
667       if(saveSinglePartition && ele[i]->getPartition() != saveSinglePartition)
668         continue;
669       writeElementMSH(fp, model, ele[i], binary, ge->tag());
670     }
671   }
672 }
673 
writeMSHPeriodicNodes(FILE * fp,std::vector<GEntity * > & entities,bool renumber,bool saveAll)674 void writeMSHPeriodicNodes(FILE *fp, std::vector<GEntity *> &entities,
675                            bool renumber, bool saveAll) // also used in MSH2
676 {
677   int count = 0;
678   for(std::size_t i = 0; i < entities.size(); i++) {
679     if(entities[i]->getMeshMaster() != entities[i] &&
680        (saveAll || (entities[i]->physicals.size() &&
681                     entities[i]->getMeshMaster()->physicals.size()))) {
682       count++;
683     }
684   }
685   if(!count) return;
686   fprintf(fp, "$Periodic\n");
687   fprintf(fp, "%d\n", count);
688   for(std::size_t i = 0; i < entities.size(); i++) {
689     GEntity *g_slave = entities[i];
690     GEntity *g_master = g_slave->getMeshMaster();
691     if(g_slave != g_master &&
692        (saveAll || (entities[i]->physicals.size() &&
693                     entities[i]->getMeshMaster()->physicals.size()))) {
694       fprintf(fp, "%d %d %d\n", g_slave->dim(), g_slave->tag(),
695               g_master->tag());
696 
697       if(g_slave->affineTransform.size() == 16) {
698         fprintf(fp, "Affine");
699         for(int i = 0; i < 16; i++)
700           fprintf(fp, " %.16g", g_slave->affineTransform[i]);
701         fprintf(fp, "\n");
702       }
703 
704       std::map<MVertex *, MVertex *, MVertexPtrLessThan> corrVert = g_slave->correspondingVertices;
705       if(CTX::instance()->mesh.hoSavePeriodic)
706         corrVert.insert(g_slave->correspondingHighOrderVertices.begin(),
707                         g_slave->correspondingHighOrderVertices.end());
708 
709       fprintf(fp, "%d\n", (int)corrVert.size());
710       for(auto it = corrVert.begin(); it != corrVert.end(); it++) {
711         MVertex *v1 = it->first;
712         MVertex *v2 = it->second;
713         if(renumber)
714           fprintf(fp, "%ld %ld\n", v1->getIndex(), v2->getIndex());
715         else
716           fprintf(fp, "%lu %lu\n", v1->getNum(), v2->getNum());
717       }
718     }
719   }
720   fprintf(fp, "$EndPeriodic\n");
721 }
722 
_writeMSH3(const std::string & name,double version,bool binary,bool saveAll,bool saveParametric,double scalingFactor,int elementStartNum,int saveSinglePartition,bool append)723 int GModel::_writeMSH3(const std::string &name, double version, bool binary,
724                        bool saveAll, bool saveParametric, double scalingFactor,
725                        int elementStartNum, int saveSinglePartition,
726                        bool append)
727 {
728   if(version < 3. || version >= 4.) {
729     Msg::Error("Wrong MSH file version %g for MSH3 writer", version);
730     return 0;
731   }
732 
733   FILE *fp;
734   if(append)
735     fp = Fopen(name.c_str(), binary ? "ab" : "a");
736   else
737     fp = Fopen(name.c_str(), binary ? "wb" : "w");
738   if(!fp) {
739     Msg::Error("Unable to open file '%s'", name.c_str());
740     return 0;
741   }
742 
743   // should make this available to users, and should allow to renumber elements,
744   // too. Renumbering should be disabled by default.
745   bool renumber = false;
746 
747   // if there are no physicals we save all the elements
748   if(noPhysicalGroups()) saveAll = true;
749 
750   // get the number of vertices and index the vertices in a continuous
751   // sequence
752   int numVertices = indexMeshVertices(saveAll, saveSinglePartition, renumber);
753 
754   // get the number of elements we need to save
755   std::vector<GEntity *> entities;
756   getEntities(entities);
757   int numElements = 0;
758   for(std::size_t i = 0; i < entities.size(); i++)
759     numElements += getNumElementsMSH(entities[i], saveAll, saveSinglePartition);
760 
761   fprintf(fp, "$MeshFormat\n");
762   fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
763   if(binary) {
764     int one = 1;
765     fwrite(&one, sizeof(int), 1, fp);
766     fprintf(fp, "\n");
767   }
768   fprintf(fp, "$EndMeshFormat\n");
769 
770   if(numPhysicalNames()) {
771     fprintf(fp, "$PhysicalNames\n");
772     fprintf(fp, "%d\n", numPhysicalNames());
773     for(auto it = firstPhysicalName(); it != lastPhysicalName(); it++) {
774       std::string name = it->second;
775       if(name.size() > 254) name.resize(254);
776       fprintf(fp, "%d %d \"%s\"\n", it->first.first, it->first.second,
777               name.c_str());
778     }
779     fprintf(fp, "$EndPhysicalNames\n");
780   }
781 
782   writeMSHEntities(fp, this);
783 
784   fprintf(fp, "$Nodes\n");
785   fprintf(fp, "%d\n", numVertices);
786   for(std::size_t i = 0; i < entities.size(); i++)
787     for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
788       entities[i]->mesh_vertices[j]->writeMSH(fp, binary, saveParametric,
789                                               scalingFactor);
790 
791   if(binary) fprintf(fp, "\n");
792   fprintf(fp, "$EndNodes\n");
793 
794   fprintf(fp, "$Elements\n");
795   fprintf(fp, "%d\n", numElements);
796 
797   _elementIndexCache.clear();
798 
799   for(auto it = firstRegion(); it != lastRegion(); ++it)
800     writeElementsMSH(fp, this, *it, (*it)->tetrahedra, saveAll,
801                      saveSinglePartition, binary);
802   for(auto it = firstRegion(); it != lastRegion(); ++it)
803     writeElementsMSH(fp, this, *it, (*it)->hexahedra, saveAll,
804                      saveSinglePartition, binary);
805   for(auto it = firstRegion(); it != lastRegion(); ++it)
806     writeElementsMSH(fp, this, *it, (*it)->prisms, saveAll, saveSinglePartition,
807                      binary);
808   for(auto it = firstRegion(); it != lastRegion(); ++it)
809     writeElementsMSH(fp, this, *it, (*it)->pyramids, saveAll,
810                      saveSinglePartition, binary);
811   for(auto it = firstRegion(); it != lastRegion(); ++it)
812     writeElementsMSH(fp, this, *it, (*it)->trihedra, saveAll,
813                      saveSinglePartition, binary);
814   for(auto it = firstFace(); it != lastFace(); ++it)
815     writeElementsMSH(fp, this, *it, (*it)->triangles, saveAll,
816                      saveSinglePartition, binary);
817   for(auto it = firstFace(); it != lastFace(); ++it)
818     writeElementsMSH(fp, this, *it, (*it)->quadrangles, saveAll,
819                      saveSinglePartition, binary);
820   for(auto it = firstEdge(); it != lastEdge(); ++it)
821     writeElementsMSH(fp, this, *it, (*it)->lines, saveAll, saveSinglePartition,
822                      binary);
823   for(auto it = firstVertex(); it != lastVertex(); ++it)
824     writeElementsMSH(fp, this, *it, (*it)->points, saveAll, saveSinglePartition,
825                      binary);
826 
827   if(binary) fprintf(fp, "\n");
828 
829   fprintf(fp, "$EndElements\n");
830 
831   // save periodic nodes
832   writeMSHPeriodicNodes(fp, entities, renumber, saveAll);
833 
834   fclose(fp);
835 
836   return 1;
837 }
838 
_writePartitionedMSH3(const std::string & baseName,double version,bool binary,bool saveAll,bool saveParametric,double scalingFactor)839 int GModel::_writePartitionedMSH3(const std::string &baseName, double version,
840                                   bool binary, bool saveAll,
841                                   bool saveParametric, double scalingFactor)
842 {
843   if(version < 3. || version >= 4.) {
844     Msg::Error("Wrong MSH file version %g for MSH3 writer", version);
845     return 0;
846   }
847 
848   for(std::size_t partition = 0; partition < getNumPartitions(); partition++) {
849     std::ostringstream sstream;
850     sstream << baseName << "_" << std::setw(6) << std::setfill('0')
851             << partition;
852     Msg::Info("Writing partition %d in file '%s'", partition,
853               sstream.str().c_str());
854     _writeMSH3(sstream.str(), version, binary, saveAll, saveParametric,
855                scalingFactor, 0, partition, false);
856   }
857   return 1;
858 }
859