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 // Contributor(s):
7 //   Anthony Royer
8 
9 #include <cstdio>
10 #include <fstream>
11 #include <vector>
12 #include <map>
13 #include <algorithm>
14 #include <sstream>
15 #include <string>
16 #include <cstdlib>
17 #include <limits>
18 
19 #include "GmshDefines.h"
20 #include "OS.h"
21 #include "Context.h"
22 #include "GModel.h"
23 #include "GEntity.h"
24 #include "partitionRegion.h"
25 #include "partitionFace.h"
26 #include "partitionEdge.h"
27 #include "partitionVertex.h"
28 #include "discreteEdge.h"
29 #include "ghostFace.h"
30 #include "ghostEdge.h"
31 #include "ghostRegion.h"
32 #include "MPoint.h"
33 #include "MLine.h"
34 #include "MTriangle.h"
35 #include "MQuadrangle.h"
36 #include "MTetrahedron.h"
37 #include "MHexahedron.h"
38 #include "MPrism.h"
39 #include "MPyramid.h"
40 #include "MTrihedron.h"
41 #include "StringUtils.h"
42 
readMSH4Physicals(GModel * const model,FILE * fp,GEntity * const entity,bool binary,char * str,bool swap)43 static bool readMSH4Physicals(GModel *const model, FILE *fp,
44                               GEntity *const entity, bool binary, char *str,
45                               bool swap)
46 {
47   std::size_t numPhy = 0;
48   if(binary) {
49     if(fread(&numPhy, sizeof(std::size_t), 1, fp) != 1) { return false; }
50     if(swap) SwapBytes((char *)&numPhy, sizeof(std::size_t), 1);
51 
52     std::vector<int> phyTags(numPhy);
53     if(fread(&phyTags[0], sizeof(int), numPhy, fp) != numPhy) { return false; }
54     if(swap) SwapBytes((char *)&phyTags[0], sizeof(int), numPhy);
55 
56     for(std::size_t i = 0; i < numPhy; i++) {
57       entity->addPhysicalEntity(phyTags[i]);
58     }
59   }
60   else {
61     sscanf(str, "%lu %[0-9- ]", &numPhy, str);
62     for(std::size_t i = 0; i < numPhy; i++) {
63       int phyTag = 0;
64 
65       if(i == numPhy - 1 && entity->dim() == 0) {
66         if(sscanf(str, "%d", &phyTag) != 1) { return false; }
67       }
68       else {
69         if(sscanf(str, "%d %[0-9- ]", &phyTag, str) != 2) { return false; }
70       }
71 
72       entity->addPhysicalEntity(phyTag);
73     }
74   }
75   return true;
76 }
77 
readMSH4BoundingEntities(GModel * const model,FILE * fp,GEntity * const entity,bool binary,char * str,bool swap)78 static bool readMSH4BoundingEntities(GModel *const model, FILE *fp,
79                                      GEntity *const entity, bool binary,
80                                      char *str, bool swap)
81 {
82   std::size_t numBrep = 0;
83   std::vector<GEntity *> boundingEntities;
84   std::vector<int> boundingSign;
85 
86   if(binary) {
87     if(fread(&numBrep, sizeof(std::size_t), 1, fp) != 1) { return false; }
88     if(swap) SwapBytes((char *)&numBrep, sizeof(std::size_t), 1);
89 
90     std::vector<int> brepTags(numBrep);
91     if(fread(&brepTags[0], sizeof(int), numBrep, fp) != numBrep) {
92       return false;
93     }
94     if(swap) SwapBytes((char *)&brepTags[0], sizeof(int), numBrep);
95 
96     for(std::size_t i = 0; i < numBrep; i++) {
97       GEntity *brep =
98         model->getEntityByTag(entity->dim() - 1, std::abs(brepTags[i]));
99       if(!brep) {
100         Msg::Warning("Entity %d not found in the Brep of entity %d",
101                      brepTags[i], entity->tag());
102       }
103       else {
104         boundingEntities.push_back(brep);
105         boundingSign.push_back((std::abs(brepTags[i]) == brepTags[i] ? 1 : -1));
106       }
107     }
108   }
109   else {
110     sscanf(str, "%lu %[0-9- ]", &numBrep, str);
111     for(std::size_t i = 0; i < numBrep; i++) {
112       int entityTag = 0;
113 
114       if(i != numBrep - 1) {
115         if(sscanf(str, "%d %[0-9- ]", &entityTag, str) != 2) { return false; }
116       }
117       else {
118         if(sscanf(str, "%d", &entityTag) != 1) { return false; }
119       }
120 
121       GEntity *brep =
122         model->getEntityByTag(entity->dim() - 1, std::abs(entityTag));
123       if(!brep) {
124         Msg::Warning("Entity %d not found in the Brep of entity %d", entityTag,
125                      entity->tag());
126       }
127       else {
128         boundingEntities.push_back(brep);
129         boundingSign.push_back((std::abs(entityTag) == entityTag ? 1 : -1));
130       }
131     }
132   }
133 
134   switch(entity->dim()) {
135   case 1:
136     if(boundingEntities.size() == 2) {
137       reinterpret_cast<GEdge *>(entity)->setBeginVertex(
138         reinterpret_cast<GVertex *>(boundingEntities[0]));
139       reinterpret_cast<GEdge *>(entity)->setEndVertex(
140         reinterpret_cast<GVertex *>(boundingEntities[1]));
141     }
142     else if(boundingEntities.size() == 1) {
143       if(boundingSign[0] == 1) {
144         reinterpret_cast<GEdge *>(entity)->setBeginVertex(
145           reinterpret_cast<GVertex *>(boundingEntities[0]));
146       }
147       else {
148         reinterpret_cast<GEdge *>(entity)->setEndVertex(
149           reinterpret_cast<GVertex *>(boundingEntities[0]));
150       }
151     }
152     break;
153   case 2: {
154     std::vector<int> tags(boundingEntities.size());
155     for(std::size_t i = 0; i < boundingEntities.size(); i++)
156       tags[i] = std::abs(boundingEntities[i]->tag());
157     reinterpret_cast<GFace *>(entity)->setBoundEdges(tags, boundingSign);
158   } break;
159   case 3: {
160     std::vector<int> tags(boundingEntities.size());
161     for(std::size_t i = 0; i < boundingEntities.size(); i++)
162       tags[i] = std::abs(boundingEntities[i]->tag());
163     reinterpret_cast<GRegion *>(entity)->setBoundFaces(tags, boundingSign);
164   } break;
165   default: break;
166   }
167   return true;
168 }
169 
readMSH4EntityInfo(FILE * fp,bool binary,char * str,int sizeofstr,bool swap,double version,bool partition,int dim,int & tag,int & parentDim,int & parentTag,std::vector<int> & partitions,double & minX,double & minY,double & minZ,double & maxX,double & maxY,double & maxZ)170 static bool readMSH4EntityInfo(FILE *fp, bool binary, char *str, int sizeofstr,
171                                bool swap, double version, bool partition,
172                                int dim, int &tag, int &parentDim,
173                                int &parentTag, std::vector<int> &partitions,
174                                double &minX, double &minY, double &minZ,
175                                double &maxX, double &maxY, double &maxZ)
176 {
177   if(partition) {
178     if(binary) {
179       int dataInt[3];
180       if(fread(dataInt, sizeof(int), 3, fp) != 3) { return false; }
181       if(swap) SwapBytes((char *)dataInt, sizeof(int), 3);
182       tag = dataInt[0];
183       parentDim = dataInt[1];
184       parentTag = dataInt[2];
185       std::size_t numPart;
186       if(fread(&numPart, sizeof(std::size_t), 1, fp) != 1) { return false; }
187       partitions.resize(numPart, 0);
188       if(fread(&partitions[0], sizeof(int), numPart, fp) != numPart) {
189         return false;
190       }
191       if(swap) SwapBytes((char *)&partitions[0], sizeof(int), numPart);
192       double dataDouble[6];
193       const std::size_t nbb = (dim > 0) ? 6 : 3;
194       if(fread(dataDouble, sizeof(double), nbb, fp) != nbb) { return false; }
195       if(swap) SwapBytes((char *)dataDouble, sizeof(double), nbb);
196       minX = dataDouble[0];
197       minY = dataDouble[1];
198       minZ = dataDouble[2];
199       maxX = dataDouble[(nbb == 6) ? 3 : 0];
200       maxY = dataDouble[(nbb == 6) ? 4 : 1];
201       maxZ = dataDouble[(nbb == 6) ? 5 : 2];
202     }
203     else {
204       std::size_t numPart = 0;
205       if(fscanf(fp, "%d %d %d %lu", &tag, &parentDim, &parentTag, &numPart) !=
206          4) {
207         return false;
208       }
209       partitions.resize(numPart, 0);
210       for(std::size_t i = 0; i < numPart; i++) {
211         if(fscanf(fp, "%d", &partitions[i]) != 1) { return false; }
212       }
213       if(version < 4.1 || dim > 0) {
214         if(fscanf(fp, "%lf %lf %lf %lf %lf %lf", &minX, &minY, &minZ, &maxX,
215                   &maxY, &maxZ) != 6) {
216           return false;
217         }
218       }
219       else {
220         if(fscanf(fp, "%lf %lf %lf", &minX, &minY, &minZ) != 3) {
221           return false;
222         }
223         maxX = minX;
224         maxY = minY;
225         maxZ = minZ;
226       }
227       if(!fgets(str, sizeofstr, fp)) { return false; }
228     }
229   }
230   else {
231     if(binary) {
232       int dataInt;
233       if(fread(&dataInt, sizeof(int), 1, fp) != 1) { return false; }
234       if(swap) SwapBytes((char *)&dataInt, sizeof(int), 1);
235       double dataDouble[6];
236       const std::size_t nbb = (dim > 0) ? 6 : 3;
237       if(fread(dataDouble, sizeof(double), nbb, fp) != nbb) { return false; }
238       if(swap) SwapBytes((char *)dataDouble, sizeof(double), nbb);
239       tag = dataInt;
240       minX = dataDouble[0];
241       minY = dataDouble[1];
242       minZ = dataDouble[2];
243       maxX = dataDouble[(nbb == 6) ? 3 : 0];
244       maxY = dataDouble[(nbb == 6) ? 4 : 1];
245       maxZ = dataDouble[(nbb == 6) ? 5 : 2];
246     }
247     else {
248       if(version < 4.1 || dim > 0) {
249         if(fscanf(fp, "%d %lf %lf %lf %lf %lf %lf", &tag, &minX, &minY, &minZ,
250                   &maxX, &maxY, &maxZ) != 7) {
251           return false;
252         }
253       }
254       else {
255         if(fscanf(fp, "%d %lf %lf %lf", &tag, &minX, &minY, &minZ) != 4) {
256           return false;
257         }
258         maxX = minX;
259         maxY = minY;
260         maxZ = minZ;
261       }
262       if(!fgets(str, sizeofstr, fp)) { return false; }
263     }
264   }
265   return true;
266 }
267 
readMSH4Entities(GModel * const model,FILE * fp,bool partition,bool binary,bool swap,double version)268 static bool readMSH4Entities(GModel *const model, FILE *fp, bool partition,
269                              bool binary, bool swap, double version)
270 {
271   if(partition) {
272     std::size_t numPartitions = 0;
273     std::size_t ghostSize = 0;
274     std::vector<int> ghostTags;
275 
276     if(binary) {
277       if(fread(&numPartitions, sizeof(std::size_t), 1, fp) != 1) {
278         return false;
279       }
280       if(swap) SwapBytes((char *)&numPartitions, sizeof(std::size_t), 1);
281 
282       if(fread(&ghostSize, sizeof(std::size_t), 1, fp) != 1) { return false; }
283       if(swap) SwapBytes((char *)&ghostSize, sizeof(std::size_t), 1);
284       if(ghostSize) {
285         ghostTags.resize(2 * ghostSize);
286         if(fread(&ghostTags[0], sizeof(int), 2 * ghostSize, fp) !=
287            2 * ghostSize) {
288           return false;
289         }
290         if(swap) SwapBytes((char *)&ghostTags[0], sizeof(int), 2 * ghostSize);
291       }
292     }
293     else {
294       if(fscanf(fp, "%lu", &numPartitions) != 1) { return false; }
295       if(fscanf(fp, "%lu", &ghostSize) != 1) { return false; }
296       if(ghostSize) {
297         ghostTags.resize(2 * ghostSize);
298         for(std::size_t i = 0; i < 2 * ghostSize; i += 2) {
299           if(fscanf(fp, "%d %d", &ghostTags[i], &ghostTags[i + 1]) != 2) {
300             return false;
301           }
302         }
303       }
304     }
305 
306     model->setNumPartitions(numPartitions);
307     Msg::Info("%lu partitions", model->getNumPartitions());
308     for(std::size_t i = 0; i < 2 * ghostSize; i += 2) {
309       switch(model->getDim()) {
310       case 1: {
311         ghostEdge *ghostEntities =
312           new ghostEdge(model, ghostTags[i], ghostTags[i + 1]);
313         model->add(ghostEntities);
314       } break;
315       case 2: {
316         ghostFace *ghostEntities =
317           new ghostFace(model, ghostTags[i], ghostTags[i + 1]);
318         model->add(ghostEntities);
319       } break;
320       case 3: {
321         ghostRegion *ghostEntities =
322           new ghostRegion(model, ghostTags[i], ghostTags[i + 1]);
323         model->add(ghostEntities);
324       } break;
325       default: break;
326       }
327     }
328   }
329 
330   std::size_t numEntities[4] = {0, 0, 0, 0};
331   if(binary) {
332     if(fread(numEntities, sizeof(std::size_t), 4, fp) != 4) { return false; }
333     if(swap) SwapBytes((char *)numEntities, sizeof(std::size_t), 4);
334   }
335   else {
336     if(fscanf(fp, "%lu %lu %lu %lu", &numEntities[0], &numEntities[1],
337               &numEntities[2], &numEntities[3]) != 4) {
338       return false;
339     }
340   }
341 
342   // max length of line for ascii input file (should be large enough to handle
343   // entities with many entities on their boundary)
344   int nume = numEntities[0] + numEntities[1] + numEntities[2] + numEntities[3];
345   std::size_t strl = std::max(4096, 128 * nume);
346   char *str = new char[strl];
347 
348   if(partition)
349     Msg::Info("%d partition entit%s", nume, nume > 1 ? "ies" : "y");
350   else
351     Msg::Info("%d entit%s", nume, nume > 1 ? "ies" : "y");
352 
353   for(int dim = 0; dim < 4; dim++) {
354     for(std::size_t i = 0; i < numEntities[dim]; i++) {
355       int tag = 0, parentDim = 0, parentTag = 0;
356       std::vector<int> partitions;
357       double minX = 0., minY = 0., minZ = 0., maxX = 0., maxY = 0., maxZ = 0.;
358       if(!readMSH4EntityInfo(fp, binary, str, strl, swap, version, partition,
359                              dim, tag, parentDim, parentTag, partitions, minX,
360                              minY, minZ, maxX, maxY, maxZ)) {
361         delete[] str;
362         return false;
363       }
364 
365       switch(dim) {
366       case 0: {
367         GVertex *gv = model->getVertexByTag(tag);
368         if(!gv) {
369           if(partition) {
370             gv = new partitionVertex(model, tag, partitions);
371             if(parentTag)
372               static_cast<partitionVertex *>(gv)->setParentEntity(
373                 model->getEntityByTag(parentDim, parentTag));
374           }
375           else {
376             gv = new discreteVertex(model, tag, minX, minY, minZ);
377           }
378           model->add(gv);
379         }
380         if(!readMSH4Physicals(model, fp, gv, binary, str, swap)) {
381           delete[] str;
382           return false;
383         }
384       } break;
385       case 1: {
386         GEdge *ge = model->getEdgeByTag(tag);
387         if(!ge) {
388           if(partition) {
389             ge = new partitionEdge(model, tag, nullptr, nullptr, partitions);
390             if(parentTag)
391               static_cast<partitionEdge *>(ge)->setParentEntity(
392                 model->getEntityByTag(parentDim, parentTag));
393           }
394           else {
395             ge = new discreteEdge(model, tag, nullptr, nullptr);
396           }
397           model->add(ge);
398         }
399         if(!readMSH4Physicals(model, fp, ge, binary, str, swap)) {
400           delete[] str;
401           return false;
402         }
403         if(!readMSH4BoundingEntities(model, fp, ge, binary, str, swap)) {
404           delete[] str;
405           return false;
406         }
407       } break;
408       case 2: {
409         GFace *gf = model->getFaceByTag(tag);
410         if(!gf) {
411           if(partition) {
412             gf = new partitionFace(model, tag, partitions);
413             if(parentTag)
414               static_cast<partitionFace *>(gf)->setParentEntity(
415                 model->getEntityByTag(parentDim, parentTag));
416           }
417           else {
418             gf = new discreteFace(model, tag);
419           }
420           model->add(gf);
421         }
422         if(!readMSH4Physicals(model, fp, gf, binary, str, swap)) {
423           delete[] str;
424           return false;
425         }
426         if(!readMSH4BoundingEntities(model, fp, gf, binary, str, swap)) {
427           delete[] str;
428           return false;
429         }
430       } break;
431       case 3: {
432         GRegion *gr = model->getRegionByTag(tag);
433         if(!gr) {
434           if(partition) {
435             gr = new partitionRegion(model, tag, partitions);
436             if(parentTag)
437               static_cast<partitionRegion *>(gr)->setParentEntity(
438                 model->getEntityByTag(parentDim, parentTag));
439           }
440           else {
441             gr = new discreteRegion(model, tag);
442           }
443           model->add(gr);
444         }
445         if(!readMSH4Physicals(model, fp, gr, binary, str, swap)) {
446           delete[] str;
447           return false;
448         }
449         if(!readMSH4BoundingEntities(model, fp, gr, binary, str, swap)) {
450           delete[] str;
451           return false;
452         }
453       } break;
454       }
455     }
456   }
457   delete[] str;
458   return true;
459 }
460 
461 static std::pair<std::size_t, MVertex *> *
readMSH4Nodes(GModel * const model,FILE * fp,bool binary,bool & dense,std::size_t & totalNumNodes,std::size_t & maxNodeNum,bool swap,double version)462 readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
463               std::size_t &totalNumNodes, std::size_t &maxNodeNum, bool swap,
464               double version)
465 {
466   std::size_t numBlock = 0, minTag = 0, maxTag = 0;
467   totalNumNodes = 0;
468   maxNodeNum = 0;
469 
470   if(binary) {
471     std::size_t data[4];
472     if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return nullptr; }
473     if(swap) SwapBytes((char *)data, sizeof(std::size_t), 4);
474     numBlock = data[0];
475     totalNumNodes = data[1];
476     minTag = data[2];
477     maxTag = data[3];
478   }
479   else {
480     if(version >= 4.1) {
481       if(fscanf(fp, "%lu %lu %lu %lu", &numBlock, &totalNumNodes, &minTag,
482                 &maxTag) != 4) {
483         return nullptr;
484       }
485     }
486     else {
487       if(fscanf(fp, "%lu %lu", &numBlock, &totalNumNodes) != 2) {
488         return nullptr;
489       }
490     }
491   }
492 
493   std::size_t nodeRead = 0;
494   std::size_t minNodeNum = std::numeric_limits<std::size_t>::max();
495 
496   std::pair<std::size_t, MVertex *> *vertexCache =
497     new std::pair<std::size_t, MVertex *>[totalNumNodes];
498 
499   Msg::Info("%lu node%s", totalNumNodes, totalNumNodes > 1 ? "s" : "");
500   Msg::StartProgressMeter(totalNumNodes);
501 
502   for(std::size_t i = 0; i < numBlock; i++) {
503     int parametric = 0;
504     int entityTag = 0, entityDim = 0;
505     std::size_t numNodes = 0;
506 
507     if(binary) {
508       int data[3];
509       if(fread(data, sizeof(int), 3, fp) != 3) {
510         delete[] vertexCache;
511         return nullptr;
512       }
513       if(swap) SwapBytes((char *)data, sizeof(int), 3);
514       entityDim = data[0];
515       entityTag = data[1];
516       parametric = data[2];
517 
518       if(fread(&numNodes, sizeof(std::size_t), 1, fp) != 1) {
519         delete[] vertexCache;
520         return nullptr;
521       }
522       if(swap) SwapBytes((char *)&numNodes, sizeof(std::size_t), 1);
523     }
524     else {
525       if(version >= 4.1) {
526         if(fscanf(fp, "%d %d %d %lu", &entityDim, &entityTag, &parametric,
527                   &numNodes) != 4) {
528           delete[] vertexCache;
529           return nullptr;
530         }
531       }
532       else {
533         if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, &parametric,
534                   &numNodes) != 4) {
535           delete[] vertexCache;
536           return nullptr;
537         }
538       }
539     }
540 
541     GEntity *entity = model->getEntityByTag(entityDim, entityTag);
542     if(!entity) {
543       switch(entityDim) {
544       case 0: {
545         Msg::Info("Creating discrete point %d", entityTag);
546         GVertex *gv = new discreteVertex(model, entityTag);
547         GModel::current()->add(gv);
548         entity = gv;
549         break;
550       }
551       case 1: {
552         Msg::Info("Creating discrete curve %d", entityTag);
553         GEdge *ge = new discreteEdge(model, entityTag, nullptr, nullptr);
554         GModel::current()->add(ge);
555         entity = ge;
556         break;
557       }
558       case 2: {
559         Msg::Info("Creating discrete surface %d", entityTag);
560         GFace *gf = new discreteFace(model, entityTag);
561         GModel::current()->add(gf);
562         entity = gf;
563         break;
564       }
565       case 3: {
566         Msg::Info("Creating discrete volume %d", entityTag);
567         GRegion *gr = new discreteRegion(model, entityTag);
568         GModel::current()->add(gr);
569         entity = gr;
570         break;
571       }
572       default: {
573         Msg::Error("Invalid dimension %d to create discrete entity", entityDim);
574         delete[] vertexCache;
575         return nullptr;
576       }
577       }
578     }
579 
580     std::size_t n = 3;
581     if(parametric) n += entityDim;
582 
583     std::vector<std::size_t> tags(numNodes);
584     if(binary) {
585       if(fread(&tags[0], sizeof(std::size_t), numNodes, fp) != numNodes) {
586         delete[] vertexCache;
587         return nullptr;
588       }
589       if(swap) SwapBytes((char *)&tags[0], sizeof(std::size_t), numNodes);
590       std::vector<double> coord(n * numNodes);
591       if(fread(&coord[0], sizeof(double), n * numNodes, fp) != n * numNodes) {
592         delete[] vertexCache;
593         return nullptr;
594       }
595       if(swap) SwapBytes((char *)&coord[0], sizeof(double), n * numNodes);
596       std::size_t k = 0;
597       for(std::size_t j = 0; j < numNodes; j++) {
598         MVertex *mv = nullptr;
599         std::size_t tagNode = tags[j];
600         if(n == 5) {
601           mv = new MFaceVertex(coord[k], coord[k + 1], coord[k + 2], entity,
602                                coord[k + 3], coord[k + 4], tagNode);
603         }
604         else if(n == 4) {
605           mv = new MEdgeVertex(coord[k], coord[k + 1], coord[k + 2], entity,
606                                coord[k + 3], tagNode);
607         }
608         else {
609           mv =
610             new MVertex(coord[k], coord[k + 1], coord[k + 2], entity, tagNode);
611         }
612         k += n;
613         entity->addMeshVertex(mv);
614         mv->setEntity(entity);
615         minNodeNum = std::min(minNodeNum, tagNode);
616         maxNodeNum = std::max(maxNodeNum, tagNode);
617         vertexCache[nodeRead] = std::make_pair(tagNode, mv);
618         nodeRead++;
619         if(totalNumNodes > 100000)
620           Msg::ProgressMeter(nodeRead, true, "Reading nodes");
621       }
622     }
623     else {
624       if(version >= 4.1) {
625         for(std::size_t j = 0; j < numNodes; j++) {
626           if(fscanf(fp, "%lu", &tags[j]) != 1) {
627             delete[] vertexCache;
628             return nullptr;
629           }
630         }
631       }
632       for(std::size_t j = 0; j < numNodes; j++) {
633         std::size_t tagNode = 0;
634         if(version >= 4.1) { tagNode = tags[j]; }
635         else {
636           if(fscanf(fp, "%lu", &tagNode) != 1) {
637             delete[] vertexCache;
638             return nullptr;
639           }
640         }
641         MVertex *mv = nullptr;
642         if(n == 5) {
643           double x, y, z, u, v;
644           if(fscanf(fp, "%lf %lf %lf %lf %lf", &x, &y, &z, &u, &v) != 5) {
645             delete[] vertexCache;
646             return nullptr;
647           }
648           mv = new MFaceVertex(x, y, z, entity, u, v, tagNode);
649         }
650         else if(n == 4) {
651           double x, y, z, u;
652           if(fscanf(fp, "%lf %lf %lf %lf", &x, &y, &z, &u) != 4) {
653             delete[] vertexCache;
654             return nullptr;
655           }
656           mv = new MEdgeVertex(x, y, z, entity, u, tagNode);
657         }
658         else {
659           double x, y, z;
660           if(fscanf(fp, "%lf %lf %lf", &x, &y, &z) != 3) {
661             delete[] vertexCache;
662             return nullptr;
663           }
664           // discard extra parametric coordinates, as Gmsh does not use them
665           for(std::size_t k = 3; k < n; k++) {
666             double dummy;
667             if(fscanf(fp, "%lf", &dummy) != 1) {
668               delete[] vertexCache;
669               return nullptr;
670             }
671           }
672           mv = new MVertex(x, y, z, entity, tagNode);
673         }
674         entity->addMeshVertex(mv);
675         mv->setEntity(entity);
676         minNodeNum = std::min(minNodeNum, tagNode);
677         maxNodeNum = std::max(maxNodeNum, tagNode);
678         vertexCache[nodeRead] = std::make_pair(tagNode, mv);
679         nodeRead++;
680         if(totalNumNodes > 100000)
681           Msg::ProgressMeter(nodeRead, true, "Reading nodes");
682       }
683     }
684   }
685 
686   if(version >= 4.1) { // consistency check
687     if(minTag != minNodeNum || maxTag != maxNodeNum)
688       Msg::Warning("Min/Max node tags reported in section header are wrong: "
689                    "(%d/%d) != (%d/%d)",
690                    minTag, maxTag, minNodeNum, maxNodeNum);
691   }
692 
693   // if the vertex numbering is (fairly) dense, we fill the vector cache,
694   // otherwise we fill the map cache
695   if(minNodeNum == 1 && maxNodeNum == totalNumNodes) {
696     Msg::Debug("Vertex numbering is dense");
697     dense = true;
698   }
699   else if(maxNodeNum < 10 * totalNumNodes) {
700     Msg::Debug(
701       "Vertex numbering is fairly dense - still caching with a vector");
702     dense = true;
703   }
704   else {
705     Msg::Debug("Vertex numbering is not dense");
706     dense = false;
707   }
708 
709   return vertexCache;
710 }
711 
712 static std::pair<std::size_t, std::pair<MElement *, int> > *
readMSH4Elements(GModel * const model,FILE * fp,bool binary,bool & dense,std::size_t & totalNumElements,std::size_t & maxElementNum,bool swap,double version)713 readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
714                  std::size_t &totalNumElements, std::size_t &maxElementNum,
715                  bool swap, double version)
716 {
717   char str[10000]; // 1000 nodes for order 9 hex, 10 digits each
718   std::size_t numBlock = 0, minTag = 0, maxTag = 0;
719   totalNumElements = 0;
720   maxElementNum = 0;
721 
722   if(binary) {
723     std::size_t data[4];
724     if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return nullptr; }
725     if(swap) SwapBytes((char *)data, sizeof(std::size_t), 4);
726     numBlock = data[0];
727     totalNumElements = data[1];
728     minTag = data[2];
729     maxTag = data[3];
730   }
731   else {
732     if(version >= 4.1) {
733       if(fscanf(fp, "%lu %lu %lu %lu", &numBlock, &totalNumElements, &minTag,
734                 &maxTag) != 4) {
735         return nullptr;
736       }
737     }
738     else {
739       if(fscanf(fp, "%lu %lu", &numBlock, &totalNumElements) != 2) {
740         return nullptr;
741       }
742     }
743   }
744 
745   std::size_t elementRead = 0;
746   std::size_t minElementNum = std::numeric_limits<std::size_t>::max();
747 
748   std::pair<std::size_t, std::pair<MElement *, int> > *elementCache =
749     new std::pair<std::size_t, std::pair<MElement *, int> >[totalNumElements];
750   Msg::Info("%lu element%s", totalNumElements, totalNumElements > 1 ? "s" : "");
751   Msg::StartProgressMeter(totalNumElements);
752 
753   for(std::size_t i = 0; i < numBlock; i++) {
754     int entityTag = 0, entityDim = 0, elmType = 0;
755     std::size_t numElements = 0;
756 
757     if(binary) {
758       int data[3];
759       if(fread(data, sizeof(int), 3, fp) != 3) {
760         delete[] elementCache;
761         return nullptr;
762       }
763       if(swap) SwapBytes((char *)data, sizeof(int), 3);
764       entityDim = data[0];
765       entityTag = data[1];
766       elmType = data[2];
767 
768       if(fread(&numElements, sizeof(std::size_t), 1, fp) != 1) {
769         delete[] elementCache;
770         return nullptr;
771       }
772       if(swap) SwapBytes((char *)&numElements, sizeof(std::size_t), 1);
773     }
774     else {
775       if(version >= 4.1) {
776         if(fscanf(fp, "%d %d %d %lu", &entityDim, &entityTag, &elmType,
777                   &numElements) != 4) {
778           delete[] elementCache;
779           return nullptr;
780         }
781       }
782       else {
783         if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, &elmType,
784                   &numElements) != 4) {
785           delete[] elementCache;
786           return nullptr;
787         }
788       }
789     }
790 
791     GEntity *entity = model->getEntityByTag(entityDim, entityTag);
792     if(!entity) {
793       Msg::Error("Unknown entity %d of dimension %d", entityTag, entityDim);
794       delete[] elementCache;
795       return nullptr;
796     }
797     if(entity->geomType() == GEntity::GhostCurve) {
798       static_cast<ghostEdge *>(entity)->haveMesh(true);
799     }
800     else if(entity->geomType() == GEntity::GhostSurface) {
801       static_cast<ghostFace *>(entity)->haveMesh(true);
802     }
803     else if(entity->geomType() == GEntity::GhostVolume) {
804       static_cast<ghostRegion *>(entity)->haveMesh(true);
805     }
806 
807     const int numVertPerElm = MElement::getInfoMSH(elmType);
808     if(binary) {
809       std::size_t n = 1 + numVertPerElm;
810       std::vector<std::size_t> data(numElements * n);
811       if(fread(&data[0], sizeof(std::size_t), numElements * n, fp) !=
812          numElements * n) {
813         delete[] elementCache;
814         return nullptr;
815       }
816       if(swap)
817         SwapBytes((char *)&data[0], sizeof(std::size_t), numElements * n);
818 
819       std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)nullptr);
820       for(std::size_t j = 0; j < numElements * n; j += n) {
821         for(int k = 0; k < numVertPerElm; k++) {
822           vertices[k] = model->getMeshVertexByTag(data[j + k + 1]);
823           if(!vertices[k]) {
824             Msg::Error("Unknown node %lu in element %lu", data[j + k + 1],
825                        data[j]);
826             delete[] elementCache;
827             return nullptr;
828           }
829         }
830 
831         MElementFactory elementFactory;
832         MElement *element = elementFactory.create(
833           elmType, vertices, data[j], 0, false, 0, nullptr, nullptr, nullptr);
834         if(!element) {
835           Msg::Error("Could not create element %lu of type %d", data[j],
836                      elmType);
837           delete[] elementCache;
838           return nullptr;
839         }
840         if(entity->geomType() != GEntity::GhostCurve &&
841            entity->geomType() != GEntity::GhostSurface &&
842            entity->geomType() != GEntity::GhostVolume) {
843           entity->addElement(element->getType(), element);
844         }
845 
846         minElementNum = std::min(minElementNum, data[j]);
847         maxElementNum = std::max(maxElementNum, data[j]);
848 
849         elementCache[elementRead] =
850           std::make_pair(data[j], std::make_pair(element, entityTag));
851         elementRead++;
852 
853         if(totalNumElements > 100000)
854           Msg::ProgressMeter(elementRead, true, "Reading elements");
855       }
856     }
857     else {
858       for(std::size_t j = 0; j < numElements; j++) {
859         std::size_t elmTag = 0;
860         if(fscanf(fp, "%lu", &elmTag) != 1) {
861           delete[] elementCache;
862           return nullptr;
863         }
864         if(!fgets(str, sizeof(str), fp)) {
865           delete[] elementCache;
866           return nullptr;
867         }
868 
869         std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)nullptr);
870 
871         for(int k = 0; k < numVertPerElm; k++) {
872           std::size_t vertexTag = 0;
873           if(k != numVertPerElm - 1) {
874             if(sscanf(str, "%lu %[0-9- ]", &vertexTag, str) != 2) {
875               delete[] elementCache;
876               return nullptr;
877             }
878           }
879           else {
880             if(sscanf(str, "%lu", &vertexTag) != 1) {
881               delete[] elementCache;
882               return nullptr;
883             }
884           }
885 
886           vertices[k] = model->getMeshVertexByTag(vertexTag);
887           if(!vertices[k]) {
888             Msg::Error("Unknown node %lu in element %lu", vertexTag, elmTag);
889             delete[] elementCache;
890             return nullptr;
891           }
892         }
893 
894         MElementFactory elementFactory;
895         MElement *element = elementFactory.create(
896           elmType, vertices, elmTag, 0, false, 0, nullptr, nullptr, nullptr);
897         if(!element) {
898           Msg::Error("Could not create element %lu of type %d", elmTag,
899                      elmType);
900           delete[] elementCache;
901           return nullptr;
902         }
903         if(entity->geomType() != GEntity::GhostCurve &&
904            entity->geomType() != GEntity::GhostSurface &&
905            entity->geomType() != GEntity::GhostVolume) {
906           entity->addElement(element->getType(), element);
907         }
908 
909         minElementNum = std::min(minElementNum, elmTag);
910         maxElementNum = std::max(maxElementNum, elmTag);
911 
912         elementCache[elementRead] =
913           std::make_pair(elmTag, std::make_pair(element, entityTag));
914         elementRead++;
915 
916         if(totalNumElements > 100000)
917           Msg::ProgressMeter(elementRead, true, "Reading elements");
918       }
919     }
920   }
921   // if the vertex numbering is dense, we fill the vector cache, otherwise we
922   // fill the map cache
923   if(minElementNum == 1 && maxElementNum == totalNumElements) {
924     Msg::Debug("Element numbering is dense");
925     dense = true;
926   }
927   else if(maxElementNum < 10 * totalNumElements) {
928     Msg::Debug(
929       "Element numbering is fairly dense - still caching with a vector");
930     dense = true;
931   }
932   else {
933     Msg::Debug("Element numbering is not dense");
934     dense = false;
935   }
936 
937   return elementCache;
938 }
939 
readMSH4PeriodicNodes(GModel * const model,FILE * fp,bool binary,bool swap,double version)940 static bool readMSH4PeriodicNodes(GModel *const model, FILE *fp, bool binary,
941                                   bool swap, double version)
942 {
943   std::size_t numPeriodicLinks = 0;
944   if(binary) {
945     if(fread(&numPeriodicLinks, sizeof(std::size_t), 1, fp) != 1) {
946       return false;
947     }
948     if(swap) SwapBytes((char *)&numPeriodicLinks, sizeof(std::size_t), 1);
949   }
950   else {
951     if(fscanf(fp, "%lu", &numPeriodicLinks) != 1) { return false; }
952   }
953 
954   for(std::size_t i = 0; i < numPeriodicLinks; i++) {
955     int slaveDim = 0, slaveTag = 0, masterTag = 0;
956 
957     if(binary) {
958       int data[3];
959       if(fread(data, sizeof(int), 3, fp) != 3) { return false; }
960       if(swap) SwapBytes((char *)data, sizeof(int), 3);
961       slaveDim = data[0];
962       slaveTag = data[1];
963       masterTag = data[2];
964     }
965     else {
966       if(fscanf(fp, "%d %d %d", &slaveDim, &slaveTag, &masterTag) != 3) {
967         return false;
968       }
969     }
970 
971     GEntity *slave = nullptr, *master = nullptr;
972     switch(slaveDim) {
973     case 0:
974       slave = model->getVertexByTag(slaveTag);
975       master = model->getVertexByTag(masterTag);
976       break;
977     case 1:
978       slave = model->getEdgeByTag(slaveTag);
979       master = model->getEdgeByTag(masterTag);
980       break;
981     case 2:
982       slave = model->getFaceByTag(slaveTag);
983       master = model->getFaceByTag(masterTag);
984       break;
985     }
986 
987     if(!slave) {
988       Msg::Info("Could not find periodic entity %d of dimension %d", slaveTag,
989                 slaveDim);
990       continue;
991     }
992     if(!master) {
993       Msg::Info("Could not find periodic source entity %d of dimension %d",
994                 masterTag, slaveDim);
995       continue;
996     }
997 
998     std::size_t correspondingVertexSize = 0;
999     if(binary) {
1000       std::size_t numAffine;
1001       if(fread(&numAffine, sizeof(std::size_t), 1, fp) != 1) { return false; }
1002       if(swap) SwapBytes((char *)&numAffine, sizeof(std::size_t), 1);
1003       if(numAffine) {
1004         std::vector<double> tfo(numAffine);
1005         if(fread(&tfo[0], sizeof(double), numAffine, fp) != numAffine) {
1006           return false;
1007         }
1008         if(swap) SwapBytes((char *)&tfo[0], sizeof(double), numAffine);
1009         slave->setMeshMaster(master, tfo);
1010       }
1011       else {
1012         slave->setMeshMaster(master);
1013       }
1014       if(fread(&correspondingVertexSize, sizeof(std::size_t), 1, fp) != 1) {
1015         return false;
1016       }
1017       if(swap)
1018         SwapBytes((char *)&correspondingVertexSize, sizeof(std::size_t), 1);
1019     }
1020     else {
1021       if(version >= 4.1) {
1022         std::size_t numAffine;
1023         if(!fscanf(fp, "%lu", &numAffine)) { return false; }
1024         if(numAffine) {
1025           std::vector<double> tfo(numAffine);
1026           for(std::size_t i = 0; i < numAffine; i++) {
1027             if(fscanf(fp, "%lf", &tfo[i]) != 1) { return false; }
1028           }
1029           slave->setMeshMaster(master, tfo);
1030         }
1031         else {
1032           slave->setMeshMaster(master);
1033         }
1034         if(fscanf(fp, "%lu", &correspondingVertexSize) != 1) { return false; }
1035       }
1036       else {
1037         char affine[256];
1038         if(!fscanf(fp, "%s", affine)) { return false; }
1039         if(!strncmp(affine, "Affine", 6)) {
1040           if(!fgets(affine, sizeof(affine), fp)) { return false; }
1041           std::vector<double> tfo(16);
1042           if(sscanf(affine,
1043                     "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf "
1044                     "%lf %lf %lf %lf",
1045                     &tfo[0], &tfo[1], &tfo[2], &tfo[3], &tfo[4], &tfo[5],
1046                     &tfo[6], &tfo[7], &tfo[8], &tfo[9], &tfo[10], &tfo[11],
1047                     &tfo[12], &tfo[13], &tfo[14], &tfo[15]) != 16) {
1048             return false;
1049           }
1050           slave->setMeshMaster(master, tfo);
1051           if(fscanf(fp, "%lu", &correspondingVertexSize) != 1) { return false; }
1052         }
1053         else {
1054           slave->setMeshMaster(master);
1055           if(sscanf(affine, "%lu", &correspondingVertexSize) != 1) {
1056             return false;
1057           }
1058         }
1059       }
1060     }
1061 
1062     for(std::size_t j = 0; j < correspondingVertexSize; j++) {
1063       std::size_t v1 = 0, v2 = 0;
1064       if(binary) {
1065         std::size_t data[2];
1066         if(fread(data, sizeof(std::size_t), 2, fp) != 2) { return false; }
1067         if(swap) SwapBytes((char *)data, sizeof(std::size_t), 2);
1068         v1 = data[0];
1069         v2 = data[1];
1070       }
1071       else {
1072         if(fscanf(fp, "%lu %lu", &v1, &v2) != 2) { return false; }
1073       }
1074       MVertex *mv1 = model->getMeshVertexByTag(v1);
1075       MVertex *mv2 = model->getMeshVertexByTag(v2);
1076 
1077       if(mv1 && mv2) { slave->correspondingVertices[mv1] = mv2; }
1078       else {
1079         if(!mv1) { Msg::Info("Could not find periodic node %d", v1); }
1080         else if(!mv2) {
1081           Msg::Info("Could not find periodic node %d", v2);
1082         }
1083       }
1084     }
1085   }
1086   return true;
1087 }
1088 
readMSH4GhostElements(GModel * const model,FILE * fp,bool binary,bool swap)1089 static bool readMSH4GhostElements(GModel *const model, FILE *fp, bool binary,
1090                                   bool swap)
1091 {
1092   std::size_t numGhostCells = 0;
1093   if(binary) {
1094     if(fread(&numGhostCells, sizeof(std::size_t), 1, fp) != 1) { return false; }
1095     if(swap) SwapBytes((char *)&numGhostCells, sizeof(std::size_t), 1);
1096   }
1097   else {
1098     if(fscanf(fp, "%lu", &numGhostCells) != 1) { return false; }
1099   }
1100 
1101   std::multimap<std::pair<MElement *, int>, int> ghostCells;
1102   for(std::size_t i = 0; i < numGhostCells; i++) {
1103     std::size_t elmTag = 0;
1104     int partNum = 0;
1105     std::size_t numGhostPartitions = 0;
1106     char str[1024];
1107 
1108     if(binary) {
1109       if(fread(&elmTag, sizeof(std::size_t), 1, fp) != 1) { return false; }
1110       if(swap) SwapBytes((char *)&elmTag, sizeof(std::size_t), 1);
1111       if(fread(&partNum, sizeof(int), 1, fp) != 1) { return false; }
1112       if(swap) SwapBytes((char *)&partNum, sizeof(int), 1);
1113       if(fread(&numGhostPartitions, sizeof(std::size_t), 1, fp) != 1) {
1114         return false;
1115       }
1116       if(swap) SwapBytes((char *)&numGhostPartitions, sizeof(std::size_t), 1);
1117     }
1118     else {
1119       if(fscanf(fp, "%lu %d %lu", &elmTag, &partNum, &numGhostPartitions) !=
1120          3) {
1121         return false;
1122       }
1123       if(!fgets(str, sizeof(str), fp)) { return false; }
1124     }
1125 
1126     MElement *elm = model->getMeshElementByTag(elmTag);
1127     if(!elm) {
1128       Msg::Error("No element with tag %lu", elmTag);
1129       continue;
1130     }
1131 
1132     for(std::size_t j = 0; j < numGhostPartitions; j++) {
1133       int ghostPartition = 0;
1134 
1135       if(binary) {
1136         if(fread(&ghostPartition, sizeof(int), 1, fp) != 1) { return false; }
1137         if(swap) SwapBytes((char *)&ghostPartition, sizeof(int), 1);
1138       }
1139       else {
1140         if(j == numGhostPartitions - 1) {
1141           if(sscanf(str, "%d", &ghostPartition) != 1) { return false; }
1142         }
1143         else {
1144           if(sscanf(str, "%d %[0-9- ]", &ghostPartition, str) != 2) {
1145             return false;
1146           }
1147         }
1148       }
1149 
1150       ghostCells.insert(std::make_pair(std::make_pair(elm, partNum),
1151                                        ghostPartition));
1152     }
1153   }
1154 
1155   std::vector<GEntity *> ghostEntities(model->getNumPartitions() + 1, nullptr);
1156   std::vector<GEntity *> entities;
1157   model->getEntities(entities);
1158   for(std::size_t i = 0; i < entities.size(); i++) {
1159     GEntity *ge = entities[i];
1160     int partNum = -1;
1161     if(ge->geomType() == GEntity::GhostCurve)
1162       partNum = static_cast<ghostEdge *>(ge)->getPartition();
1163     else if(ge->geomType() == GEntity::GhostSurface)
1164       partNum = static_cast<ghostFace *>(ge)->getPartition();
1165     else if(ge->geomType() == GEntity::GhostVolume)
1166       partNum = static_cast<ghostRegion *>(ge)->getPartition();
1167     if(partNum >= 0 && partNum < (int)ghostEntities.size())
1168       ghostEntities[partNum] = ge;
1169   }
1170 
1171   for(auto it = ghostCells.begin(); it != ghostCells.end(); ++it) {
1172     if(it->second >= (int)ghostEntities.size()) {
1173       Msg::Error("Invalid partition %d in ghost elements", it->second);
1174       return false;
1175     }
1176     GEntity *ge = ghostEntities[it->second];
1177     if(!ge) {
1178       Msg::Warning("Missing ghost entity on partition %d", it->second);
1179     }
1180     else if(ge->geomType() == GEntity::GhostCurve) {
1181       static_cast<ghostEdge *>(ge)->addElement(
1182         it->first.first->getType(), it->first.first, it->first.second);
1183     }
1184     else if(ge->geomType() == GEntity::GhostSurface) {
1185       static_cast<ghostFace *>(ge)->addElement(
1186         it->first.first->getType(), it->first.first, it->first.second);
1187     }
1188     else if(ge->geomType() == GEntity::GhostVolume) {
1189       static_cast<ghostRegion *>(ge)->addElement(
1190         it->first.first->getType(), it->first.first, it->first.second);
1191     }
1192   }
1193   return true;
1194 }
1195 
readMSH4Parametrizations(GModel * const model,FILE * fp,bool binary)1196 static bool readMSH4Parametrizations(GModel *const model, FILE *fp, bool binary)
1197 {
1198   if(CTX::instance()->mesh.ignoreParametrizationMsh4) return true;
1199 
1200   std::size_t nParamE, nParamF;
1201 
1202   if(binary) {
1203     if(fread(&nParamE, sizeof(std::size_t), 1, fp) != 1) { return false; }
1204     if(fread(&nParamF, sizeof(std::size_t), 1, fp) != 1) { return false; }
1205   }
1206   else {
1207     if(fscanf(fp, "%lu %lu", &nParamE, &nParamF) != 2) { return false; }
1208   }
1209 
1210   // only report surface parametrizations
1211   Msg::Info("%lu parametrizations", nParamF);
1212   Msg::StartProgressMeter(nParamF);
1213 
1214   for(std::size_t edge = 0; edge < nParamE; edge++) {
1215     int tag;
1216     if(binary) {
1217       if(fread(&tag, sizeof(int), 1, fp) != 1) { return false; }
1218     }
1219     else {
1220       if(fscanf(fp, "%d", &tag) != 1) { return false; }
1221     }
1222     GEdge *ge = model->getEdgeByTag(tag);
1223     if(ge) {
1224       discreteEdge *de = dynamic_cast<discreteEdge *>(ge);
1225       if(de) {
1226         if(!de->readParametrization(fp, binary)) return false;
1227       }
1228     }
1229   }
1230 
1231   for(std::size_t face = 0; face < nParamF; face++) {
1232     int tag;
1233     if(binary) {
1234       if(fread(&tag, sizeof(int), 1, fp) != 1) { return false; }
1235     }
1236     else {
1237       if(fscanf(fp, "%d", &tag) != 1) { return false; }
1238     }
1239     GFace *gf = model->getFaceByTag(tag);
1240     if(gf) {
1241       discreteFace *df = dynamic_cast<discreteFace *>(gf);
1242       if(df) {
1243         if(!df->readParametrization(fp, binary)) return false;
1244       }
1245     }
1246     Msg::ProgressMeter(face, true, "Processing parametrizations");
1247   }
1248 
1249   Msg::StopProgressMeter();
1250 
1251   return true;
1252 }
1253 
_readMSH4(const std::string & name)1254 int GModel::_readMSH4(const std::string &name)
1255 {
1256   bool partitioned = false;
1257   FILE *fp = Fopen(name.c_str(), "rb");
1258   if(!fp) {
1259     Msg::Error("Unable to open file '%s'", name.c_str());
1260     return 0;
1261   }
1262 
1263   char str[1024] = "x";
1264   double version = 1.0;
1265   bool binary = false, swap = false, postpro = false;
1266 
1267   while(1) {
1268     while(str[0] != '$') {
1269       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
1270     }
1271 
1272     std::string sectionName(&str[1]);
1273     std::string endSectionName = "End" + sectionName;
1274 
1275     if(feof(fp)) break;
1276 
1277     if(!strncmp(&str[1], "MeshFormat", 10)) {
1278       if(!fgets(str, sizeof(str), fp) || feof(fp)) {
1279         fclose(fp);
1280         return 0;
1281       }
1282 
1283       int format;
1284       std::size_t size;
1285       if(sscanf(str, "%lf %d %lu", &version, &format, &size) != 3) {
1286         fclose(fp);
1287         return 0;
1288       }
1289       if(format) {
1290         binary = true;
1291         Msg::Debug("Mesh is in binary format");
1292         int one;
1293         if(fread(&one, sizeof(int), 1, fp) != 1) {
1294           fclose(fp);
1295           return 0;
1296         }
1297         if(one != 1) {
1298           swap = true;
1299           Msg::Debug("Swapping bytes from binary file");
1300         }
1301       }
1302 
1303       if(binary && size != sizeof(std::size_t)) {
1304         Msg::Error("Binary file has sizeof(size_t) = %d, not matching "
1305                    "machine sizeof(size_t) = %d",
1306                    size, sizeof(std::size_t));
1307         return false;
1308       }
1309       if(binary && version < 4.1) {
1310         Msg::Error("Can only read MSH 4.0 format in ASCII mode");
1311         return false;
1312       }
1313     }
1314     else if(!strncmp(&str[1], "PhysicalNames", 13)) {
1315       if(!fgets(str, sizeof(str), fp) || feof(fp)) {
1316         fclose(fp);
1317         return 0;
1318       }
1319       int numPhysicalNames = 0;
1320       if(sscanf(str, "%d", &numPhysicalNames) != 1) {
1321         fclose(fp);
1322         return 0;
1323       }
1324       std::vector<GModel::piter> iterators;
1325       getInnerPhysicalNamesIterators(iterators);
1326       for(int i = 0; i < numPhysicalNames; i++) {
1327         int dim = 0, tag = 0;
1328         if(fscanf(fp, "%d %d", &dim, &tag) != 2) {
1329           fclose(fp);
1330           return 0;
1331         }
1332         char name[256];
1333         if(!fgets(name, sizeof(name), fp)) {
1334           fclose(fp);
1335           return 0;
1336         }
1337         std::string physicalName = ExtractDoubleQuotedString(name, 256);
1338         if(physicalName.size())
1339           iterators[dim] =
1340             setPhysicalName(iterators[dim], physicalName, dim, tag);
1341       }
1342     }
1343     else if(!strncmp(&str[1], "Entities", 8)) {
1344       if(!readMSH4Entities(this, fp, false, binary, swap, version)) {
1345         Msg::Error("Could not read entities");
1346         fclose(fp);
1347         return 0;
1348       }
1349     }
1350     else if(!strncmp(&str[1], "PartitionedEntities", 19)) {
1351       if(!readMSH4Entities(this, fp, true, binary, swap, version)) {
1352         Msg::Error("Could not read partitioned entities");
1353         fclose(fp);
1354         return 0;
1355       }
1356       partitioned = true;
1357     }
1358     else if(!strncmp(&str[1], "Nodes", 5)) {
1359       _vertexVectorCache.clear();
1360       _vertexMapCache.clear();
1361       bool dense = false;
1362       std::size_t totalNumNodes = 0, maxNodeNum;
1363       std::pair<std::size_t, MVertex *> *vertexCache = readMSH4Nodes(
1364         this, fp, binary, dense, totalNumNodes, maxNodeNum, swap, version);
1365       Msg::StopProgressMeter();
1366       if(!vertexCache) {
1367         Msg::Error("Could not read nodes");
1368         fclose(fp);
1369         return false;
1370       }
1371       if(dense) {
1372         _vertexVectorCache.resize(maxNodeNum + 1, nullptr);
1373         for(std::size_t i = 0; i < totalNumNodes; i++) {
1374           if(!_vertexVectorCache[vertexCache[i].first]) {
1375             _vertexVectorCache[vertexCache[i].first] = vertexCache[i].second;
1376           }
1377           else {
1378             Msg::Info("Skipping duplicate node %d", vertexCache[i].first);
1379           }
1380         }
1381       }
1382       else {
1383         for(std::size_t i = 0; i < totalNumNodes; i++) {
1384           if(_vertexMapCache.count(vertexCache[i].first) == 0) {
1385             _vertexMapCache[vertexCache[i].first] = vertexCache[i].second;
1386           }
1387           else {
1388             Msg::Info("Skipping duplicate node %d", vertexCache[i].first);
1389           }
1390         }
1391       }
1392       delete[] vertexCache;
1393     }
1394     else if(!strncmp(&str[1], "Elements", 8)) {
1395       bool dense = false;
1396       std::size_t totalNumElements = 0, maxElementNum = 0;
1397       std::pair<std::size_t, std::pair<MElement *, int> > *elementCache =
1398         readMSH4Elements(this, fp, binary, dense, totalNumElements,
1399                          maxElementNum, swap, version);
1400       Msg::StopProgressMeter();
1401       if(!elementCache) {
1402         Msg::Error("Could not read elements");
1403         fclose(fp);
1404         return 0;
1405       }
1406       if(dense) {
1407         _elementVectorCache.resize(maxElementNum + 1, std::make_pair(nullptr, 0));
1408         for(std::size_t i = 0; i < totalNumElements; i++) {
1409           if(!_elementVectorCache[elementCache[i].first].first) {
1410             _elementVectorCache[elementCache[i].first] = elementCache[i].second;
1411           }
1412           else {
1413             Msg::Info("Skipping duplicate element %d", elementCache[i].first);
1414           }
1415         }
1416       }
1417       else {
1418         for(std::size_t i = 0; i < totalNumElements; i++) {
1419           if(_elementMapCache.count(elementCache[i].first) == 0) {
1420             _elementMapCache[elementCache[i].first] = elementCache[i].second;
1421           }
1422           else {
1423             Msg::Info("Skipping duplicate element %d", elementCache[i].first);
1424           }
1425         }
1426       }
1427       delete[] elementCache;
1428     }
1429     else if(!strncmp(&str[1], "Periodic", 8)) {
1430       if(!readMSH4PeriodicNodes(this, fp, binary, swap, version)) {
1431         Msg::Error("Could not read periodic section");
1432         fclose(fp);
1433         return 0;
1434       }
1435     }
1436     else if(!strncmp(&str[1], "GhostElements", 13)) {
1437       if(!readMSH4GhostElements(this, fp, binary, swap)) {
1438         Msg::Error("Could not read ghost elements");
1439         fclose(fp);
1440         return 0;
1441       }
1442     }
1443     else if(!strncmp(&str[1], "Parametrizations", 16)) {
1444       if(!readMSH4Parametrizations(this, fp, binary)) {
1445         Msg::Error("Could not read parametrizations");
1446         fclose(fp);
1447         return 0;
1448       }
1449     }
1450     else if(!strncmp(&str[1], "NodeData", 8) ||
1451             !strncmp(&str[1], "ElementData", 11) ||
1452             !strncmp(&str[1], "ElementNodeData", 15)) {
1453       postpro = true;
1454       break;
1455     }
1456 
1457     while(strncmp(&str[1], endSectionName.c_str(), endSectionName.size())) {
1458       if(!fgets(str, sizeof(str), fp) || feof(fp)) { break; }
1459     }
1460     str[0] = 'a';
1461   }
1462 
1463   fclose(fp);
1464 
1465   if(partitioned) {
1466     // This part is added to ensure the compatibility between the new
1467     // partitioning and the old one.
1468     std::vector<GEntity *> entities;
1469     getEntities(entities);
1470     for(std::size_t i = 0; i < entities.size(); i++) {
1471       if(entities[i]->geomType() == GEntity::PartitionPoint) {
1472         partitionVertex *pv = static_cast<partitionVertex *>(entities[i]);
1473         if(pv->numPartitions() == 1) {
1474           const int part = pv->getPartition(0);
1475           for(std::size_t j = 0; j < pv->getNumMeshElements(); j++) {
1476             pv->getMeshElement(j)->setPartition(part);
1477           }
1478         }
1479       }
1480       else if(entities[i]->geomType() == GEntity::PartitionCurve) {
1481         partitionEdge *pe = static_cast<partitionEdge *>(entities[i]);
1482         if(pe->numPartitions() == 1) {
1483           const int part = pe->getPartition(0);
1484           for(std::size_t j = 0; j < pe->getNumMeshElements(); j++) {
1485             pe->getMeshElement(j)->setPartition(part);
1486           }
1487         }
1488       }
1489       else if(entities[i]->geomType() == GEntity::PartitionSurface) {
1490         partitionFace *pf = static_cast<partitionFace *>(entities[i]);
1491         if(pf->numPartitions() == 1) {
1492           const int part = pf->getPartition(0);
1493           for(std::size_t j = 0; j < pf->getNumMeshElements(); j++) {
1494             pf->getMeshElement(j)->setPartition(part);
1495           }
1496         }
1497       }
1498       else if(entities[i]->geomType() == GEntity::PartitionVolume) {
1499         partitionRegion *pr = static_cast<partitionRegion *>(entities[i]);
1500         if(pr->numPartitions() == 1) {
1501           const int part = pr->getPartition(0);
1502           for(std::size_t j = 0; j < pr->getNumMeshElements(); j++) {
1503             pr->getMeshElement(j)->setPartition(part);
1504           }
1505         }
1506       }
1507     }
1508   }
1509 
1510   return postpro ? 2 : 1;
1511 }
1512 
writeMSH4Physicals(FILE * fp,GEntity * const entity,bool binary)1513 static void writeMSH4Physicals(FILE *fp, GEntity *const entity, bool binary)
1514 {
1515   if(binary) {
1516     std::vector<int> phys = entity->getPhysicalEntities();
1517     std::size_t phySize = phys.size();
1518     fwrite(&phySize, sizeof(std::size_t), 1, fp);
1519     for(std::size_t i = 0; i < phys.size(); i++) {
1520       int phy = phys[i];
1521       fwrite(&phy, sizeof(int), 1, fp);
1522     }
1523   }
1524   else {
1525     std::vector<int> phys = entity->getPhysicalEntities();
1526     fprintf(fp, "%lu", phys.size());
1527     for(std::size_t i = 0; i < phys.size(); i++) {
1528       fprintf(fp, " %d", phys[i]);
1529     }
1530     fprintf(fp, " ");
1531   }
1532 }
1533 
writeMSH4BoundingBox(SBoundingBox3d boundBox,FILE * fp,double scalingFactor,bool binary,int dim,double version)1534 static void writeMSH4BoundingBox(SBoundingBox3d boundBox, FILE *fp,
1535                                  double scalingFactor, bool binary, int dim,
1536                                  double version)
1537 {
1538   double bb[6] = {0., 0., 0., 0., 0., 0.};
1539   if(!boundBox.empty()) {
1540     boundBox *= scalingFactor;
1541     bb[0] = boundBox.min().x();
1542     bb[1] = boundBox.min().y();
1543     bb[2] = boundBox.min().z();
1544     bb[3] = boundBox.max().x();
1545     bb[4] = boundBox.max().y();
1546     bb[5] = boundBox.max().z();
1547   }
1548   if(binary) { fwrite(bb, sizeof(double), (dim > 0) ? 6 : 3, fp); }
1549   else {
1550     std::size_t n = (version < 4.1 || dim > 0) ? 6 : 3;
1551     for(std::size_t i = 0; i < n; i++) fprintf(fp, "%.16g ", bb[i]);
1552   }
1553 }
1554 
writeMSH4Entities(GModel * const model,FILE * fp,bool partition,bool binary,double scalingFactor,double version)1555 static void writeMSH4Entities(GModel *const model, FILE *fp, bool partition,
1556                               bool binary, double scalingFactor, double version)
1557 {
1558   std::set<GEntity *, GEntityPtrFullLessThan> ghost;
1559   std::set<GRegion *, GEntityPtrLessThan> regions;
1560   std::set<GFace *, GEntityPtrLessThan> faces;
1561   std::set<GEdge *, GEntityPtrLessThan> edges;
1562   std::set<GVertex *, GEntityPtrLessThan> vertices;
1563 
1564   if(partition) {
1565     for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
1566       if((*it)->geomType() == GEntity::PartitionPoint) vertices.insert(*it);
1567     }
1568     for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
1569       if((*it)->geomType() == GEntity::PartitionCurve) edges.insert(*it);
1570       if((*it)->geomType() == GEntity::GhostCurve) ghost.insert(*it);
1571     }
1572     for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
1573       if((*it)->geomType() == GEntity::PartitionSurface) faces.insert(*it);
1574       if((*it)->geomType() == GEntity::GhostSurface) ghost.insert(*it);
1575     }
1576     for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
1577       if((*it)->geomType() == GEntity::PartitionVolume) regions.insert(*it);
1578       if((*it)->geomType() == GEntity::GhostVolume) ghost.insert(*it);
1579     }
1580   }
1581   else {
1582     for(auto it = model->firstVertex(); it != model->lastVertex(); ++it)
1583       if((*it)->geomType() != GEntity::PartitionPoint) vertices.insert(*it);
1584     for(auto it = model->firstEdge(); it != model->lastEdge(); ++it)
1585       if((*it)->geomType() != GEntity::PartitionCurve &&
1586          (*it)->geomType() != GEntity::GhostCurve)
1587         edges.insert(*it);
1588     for(auto it = model->firstFace(); it != model->lastFace(); ++it)
1589       if((*it)->geomType() != GEntity::PartitionSurface &&
1590          (*it)->geomType() != GEntity::GhostSurface)
1591         faces.insert(*it);
1592     for(auto it = model->firstRegion(); it != model->lastRegion(); ++it)
1593       if((*it)->geomType() != GEntity::PartitionVolume &&
1594          (*it)->geomType() != GEntity::GhostVolume)
1595         regions.insert(*it);
1596   }
1597 
1598   if(partition)
1599     fprintf(fp, "$PartitionedEntities\n");
1600   else
1601     fprintf(fp, "$Entities\n");
1602 
1603   if(binary) {
1604     if(partition) {
1605       std::size_t numPartitions = model->getNumPartitions();
1606       fwrite(&numPartitions, sizeof(std::size_t), 1, fp);
1607 
1608       // write the ghostentities' tag
1609       std::size_t ghostSize = ghost.size();
1610       std::vector<int> tags;
1611       if(ghostSize) {
1612         tags.resize(2 * ghostSize);
1613         int index = 0;
1614         for(auto it = ghost.begin(); it != ghost.end(); ++it) {
1615           if((*it)->geomType() == GEntity::GhostCurve) {
1616             tags[index] = (*it)->tag();
1617             tags[++index] = static_cast<ghostEdge *>(*it)->getPartition();
1618           }
1619           else if((*it)->geomType() == GEntity::GhostSurface) {
1620             tags[index] = (*it)->tag();
1621             tags[++index] = static_cast<ghostFace *>(*it)->getPartition();
1622           }
1623           else if((*it)->geomType() == GEntity::GhostVolume) {
1624             tags[index] = (*it)->tag();
1625             tags[++index] = static_cast<ghostRegion *>(*it)->getPartition();
1626           }
1627           index++;
1628         }
1629       }
1630       fwrite(&ghostSize, sizeof(std::size_t), 1, fp);
1631       if(ghostSize) { fwrite(&tags[0], sizeof(int), 2 * ghostSize, fp); }
1632     }
1633     std::size_t verticesSize = vertices.size();
1634     std::size_t edgesSize = edges.size();
1635     std::size_t facesSize = faces.size();
1636     std::size_t regionsSize = regions.size();
1637     fwrite(&verticesSize, sizeof(std::size_t), 1, fp);
1638     fwrite(&edgesSize, sizeof(std::size_t), 1, fp);
1639     fwrite(&facesSize, sizeof(std::size_t), 1, fp);
1640     fwrite(&regionsSize, sizeof(std::size_t), 1, fp);
1641 
1642     for(auto it = vertices.begin(); it != vertices.end(); ++it) {
1643       int entityTag = (*it)->tag();
1644       fwrite(&entityTag, sizeof(int), 1, fp);
1645       if(partition) {
1646         partitionVertex *pv = static_cast<partitionVertex *>(*it);
1647         int parentEntityDim = 0, parentEntityTag = 0;
1648         if(pv->getParentEntity()) {
1649           parentEntityDim = pv->getParentEntity()->dim();
1650           parentEntityTag = pv->getParentEntity()->tag();
1651         }
1652         fwrite(&parentEntityDim, sizeof(int), 1, fp);
1653         fwrite(&parentEntityTag, sizeof(int), 1, fp);
1654         std::vector<int> partitions(pv->getPartitions().begin(),
1655                                     pv->getPartitions().end()); // FIXME
1656         std::size_t numPart = partitions.size();
1657         fwrite(&numPart, sizeof(std::size_t), 1, fp);
1658         fwrite(&partitions[0], sizeof(int), partitions.size(), fp);
1659       }
1660       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 0,
1661                            version);
1662       writeMSH4Physicals(fp, *it, binary);
1663     }
1664 
1665     for(auto it = edges.begin(); it != edges.end(); ++it) {
1666       std::vector<GVertex *> vertices;
1667       std::vector<int> ori;
1668       if((*it)->getBeginVertex()) {
1669         vertices.push_back((*it)->getBeginVertex());
1670         ori.push_back(1);
1671       }
1672       if((*it)->getEndVertex()) { // convention that the end vertex is negative
1673         vertices.push_back((*it)->getEndVertex());
1674         ori.push_back(-1);
1675       }
1676       std::size_t verticesSize = vertices.size();
1677       int entityTag = (*it)->tag();
1678       fwrite(&entityTag, sizeof(int), 1, fp);
1679       if(partition) {
1680         partitionEdge *pe = static_cast<partitionEdge *>(*it);
1681         int parentEntityDim = 0, parentEntityTag = 0;
1682         if(pe->getParentEntity()) {
1683           parentEntityDim = pe->getParentEntity()->dim();
1684           parentEntityTag = pe->getParentEntity()->tag();
1685         }
1686         fwrite(&parentEntityDim, sizeof(int), 1, fp);
1687         fwrite(&parentEntityTag, sizeof(int), 1, fp);
1688         std::vector<int> partitions(pe->getPartitions().begin(),
1689                                     pe->getPartitions().end()); // FIXME
1690         std::size_t numPart = partitions.size();
1691         fwrite(&numPart, sizeof(std::size_t), 1, fp);
1692         fwrite(&partitions[0], sizeof(int), partitions.size(), fp);
1693       }
1694       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 1,
1695                            version);
1696       writeMSH4Physicals(fp, *it, binary);
1697       fwrite(&verticesSize, sizeof(std::size_t), 1, fp);
1698       int oriI = 0;
1699       for(auto itv = vertices.begin(); itv != vertices.end(); itv++) {
1700         int brepTag = ori[oriI] * (*itv)->tag();
1701         fwrite(&brepTag, sizeof(int), 1, fp);
1702         oriI++;
1703       }
1704     }
1705 
1706     for(auto it = faces.begin(); it != faces.end(); ++it) {
1707       std::vector<GEdge *> const &edges = (*it)->edges();
1708       std::vector<int> const &ori = (*it)->edgeOrientations();
1709       std::size_t edgesSize = edges.size();
1710       int entityTag = (*it)->tag();
1711       fwrite(&entityTag, sizeof(int), 1, fp);
1712       if(partition) {
1713         partitionFace *pf = static_cast<partitionFace *>(*it);
1714         int parentEntityDim = 0, parentEntityTag = 0;
1715         if(pf->getParentEntity()) {
1716           parentEntityDim = pf->getParentEntity()->dim();
1717           parentEntityTag = pf->getParentEntity()->tag();
1718         }
1719         fwrite(&parentEntityDim, sizeof(int), 1, fp);
1720         fwrite(&parentEntityTag, sizeof(int), 1, fp);
1721         std::vector<int> partitions(pf->getPartitions().begin(),
1722                                     pf->getPartitions().end()); // FIXME
1723         std::size_t numPart = partitions.size();
1724         fwrite(&numPart, sizeof(std::size_t), 1, fp);
1725         fwrite(&partitions[0], sizeof(int), partitions.size(), fp);
1726       }
1727       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 2,
1728                            version);
1729       writeMSH4Physicals(fp, *it, binary);
1730       fwrite(&edgesSize, sizeof(std::size_t), 1, fp);
1731       std::vector<int> tags, signs;
1732       for(auto ite = edges.begin(); ite != edges.end(); ite++)
1733         tags.push_back((*ite)->tag());
1734 
1735       signs.insert(signs.end(), ori.begin(), ori.end());
1736 
1737       if(tags.size() == signs.size()) {
1738         for(std::size_t i = 0; i < tags.size(); i++)
1739           tags[i] *= (signs[i] > 0 ? 1 : -1);
1740       }
1741       for(std::size_t i = 0; i < tags.size(); i++) {
1742         int brepTag = tags[i];
1743         fwrite(&brepTag, sizeof(int), 1, fp);
1744       }
1745     }
1746 
1747     for(auto it = regions.begin(); it != regions.end(); ++it) {
1748       std::vector<GFace *> faces = (*it)->faces();
1749       std::vector<int> const &ori = (*it)->faceOrientations();
1750       std::size_t facesSize = faces.size();
1751       int entityTag = (*it)->tag();
1752       fwrite(&entityTag, sizeof(int), 1, fp);
1753       if(partition) {
1754         partitionRegion *pr = static_cast<partitionRegion *>(*it);
1755         int parentEntityDim = 0, parentEntityTag = 0;
1756         if(pr->getParentEntity()) {
1757           parentEntityDim = pr->getParentEntity()->dim();
1758           parentEntityTag = pr->getParentEntity()->tag();
1759         }
1760         fwrite(&parentEntityDim, sizeof(int), 1, fp);
1761         fwrite(&parentEntityTag, sizeof(int), 1, fp);
1762         std::vector<int> partitions(pr->getPartitions().begin(),
1763                                     pr->getPartitions().end()); // FIXME
1764         std::size_t numPart = partitions.size();
1765         fwrite(&numPart, sizeof(std::size_t), 1, fp);
1766         fwrite(&partitions[0], sizeof(int), partitions.size(), fp);
1767       }
1768       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 3,
1769                            version);
1770       writeMSH4Physicals(fp, *it, binary);
1771       fwrite(&facesSize, sizeof(std::size_t), 1, fp);
1772       std::vector<int> tags, signs;
1773       for(auto itf = faces.begin(); itf != faces.end(); itf++)
1774         tags.push_back((*itf)->tag());
1775       for(auto itf = ori.begin(); itf != ori.end(); itf++)
1776         signs.push_back(*itf);
1777       if(tags.size() == signs.size()) {
1778         for(std::size_t i = 0; i < tags.size(); i++)
1779           tags[i] *= (signs[i] > 0 ? 1 : -1);
1780       }
1781       for(std::size_t i = 0; i < tags.size(); i++) {
1782         int brepTag = tags[i];
1783         fwrite(&brepTag, sizeof(int), 1, fp);
1784       }
1785     }
1786     fprintf(fp, "\n");
1787   }
1788   else {
1789     if(partition) {
1790       fprintf(fp, "%lu\n", model->getNumPartitions());
1791 
1792       // write the ghostentities' tag
1793       std::size_t ghostSize = ghost.size();
1794       std::vector<int> tags;
1795       if(ghostSize) {
1796         tags.resize(2 * ghostSize);
1797         int index = 0;
1798         for(auto it = ghost.begin(); it != ghost.end(); ++it) {
1799           if((*it)->geomType() == GEntity::GhostCurve) {
1800             tags[index] = (*it)->tag();
1801             tags[++index] = static_cast<ghostEdge *>(*it)->getPartition();
1802           }
1803           else if((*it)->geomType() == GEntity::GhostSurface) {
1804             tags[index] = (*it)->tag();
1805             tags[++index] = static_cast<ghostFace *>(*it)->getPartition();
1806           }
1807           else if((*it)->geomType() == GEntity::GhostVolume) {
1808             tags[index] = (*it)->tag();
1809             tags[++index] = static_cast<ghostRegion *>(*it)->getPartition();
1810           }
1811           index++;
1812         }
1813       }
1814       fprintf(fp, "%lu\n", ghostSize);
1815       if(ghostSize) {
1816         for(std::size_t i = 0; i < 2 * ghostSize; i += 2) {
1817           fprintf(fp, "%d %d\n", tags[i], tags[i + 1]);
1818         }
1819       }
1820     }
1821     fprintf(fp, "%lu %lu %lu %lu\n", vertices.size(), edges.size(),
1822             faces.size(), regions.size());
1823 
1824     for(auto it = vertices.begin(); it != vertices.end(); ++it) {
1825       fprintf(fp, "%d ", (*it)->tag());
1826       if(partition) {
1827         partitionVertex *pv = static_cast<partitionVertex *>(*it);
1828         int parentEntityDim = 0, parentEntityTag = 0;
1829         if(pv->getParentEntity()) {
1830           parentEntityDim = pv->getParentEntity()->dim();
1831           parentEntityTag = pv->getParentEntity()->tag();
1832         }
1833         fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag);
1834         std::vector<int> partitions(pv->getPartitions().begin(),
1835                                     pv->getPartitions().end()); // FIXME
1836         fprintf(fp, "%lu ", partitions.size());
1837         for(std::size_t i = 0; i < partitions.size(); i++)
1838           fprintf(fp, "%d ", partitions[i]);
1839       }
1840       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 0,
1841                            version);
1842       writeMSH4Physicals(fp, *it, binary);
1843       fprintf(fp, "\n");
1844     }
1845 
1846     for(auto it = edges.begin(); it != edges.end(); ++it) {
1847       std::vector<GVertex *> vertices;
1848       std::vector<int> ori;
1849       if((*it)->getBeginVertex()) {
1850         vertices.push_back((*it)->getBeginVertex());
1851         ori.push_back(1);
1852       }
1853       if((*it)->getEndVertex()) { // I use the convention that the end vertex is
1854                                   // negative
1855         vertices.push_back((*it)->getEndVertex());
1856         ori.push_back(-1);
1857       }
1858       fprintf(fp, "%d ", (*it)->tag());
1859       if(partition) {
1860         partitionEdge *pe = static_cast<partitionEdge *>(*it);
1861         int parentEntityDim = 0, parentEntityTag = 0;
1862         if(pe->getParentEntity()) {
1863           parentEntityDim = pe->getParentEntity()->dim();
1864           parentEntityTag = pe->getParentEntity()->tag();
1865         }
1866         fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag);
1867         std::vector<int> partitions(pe->getPartitions().begin(),
1868                                     pe->getPartitions().end()); // FIXME
1869         fprintf(fp, "%lu ", partitions.size());
1870         for(std::size_t i = 0; i < partitions.size(); i++)
1871           fprintf(fp, "%d ", partitions[i]);
1872       }
1873       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 1,
1874                            version);
1875       writeMSH4Physicals(fp, *it, binary);
1876       fprintf(fp, "%lu ", vertices.size());
1877       int oriI = 0;
1878       for(auto itv = vertices.begin(); itv != vertices.end(); itv++) {
1879         fprintf(fp, "%d ", ori[oriI] * (*itv)->tag());
1880         oriI++;
1881       }
1882       fprintf(fp, "\n");
1883     }
1884 
1885     for(auto it = faces.begin(); it != faces.end(); ++it) {
1886       std::vector<GEdge *> const &edges = (*it)->edges();
1887       std::vector<int> const &ori = (*it)->edgeOrientations();
1888       fprintf(fp, "%d ", (*it)->tag());
1889       if(partition) {
1890         partitionFace *pf = static_cast<partitionFace *>(*it);
1891         int parentEntityDim = 0, parentEntityTag = 0;
1892         if(pf->getParentEntity()) {
1893           parentEntityDim = pf->getParentEntity()->dim();
1894           parentEntityTag = pf->getParentEntity()->tag();
1895         }
1896         fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag);
1897         std::vector<int> partitions(pf->getPartitions().begin(),
1898                                     pf->getPartitions().end()); // FIXME
1899         fprintf(fp, "%lu ", partitions.size());
1900         for(std::size_t i = 0; i < partitions.size(); i++)
1901           fprintf(fp, "%d ", partitions[i]);
1902       }
1903       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 2,
1904                            version);
1905       writeMSH4Physicals(fp, *it, binary);
1906       fprintf(fp, "%lu ", edges.size());
1907       std::vector<int> tags, signs;
1908       for(auto ite = edges.begin(); ite != edges.end(); ite++)
1909         tags.push_back((*ite)->tag());
1910       for(auto ite = ori.begin(); ite != ori.end(); ite++)
1911         signs.push_back(*ite);
1912       if(tags.size() == signs.size()) {
1913         for(std::size_t i = 0; i < tags.size(); i++)
1914           tags[i] *= (signs[i] > 0 ? 1 : -1);
1915       }
1916       for(std::size_t i = 0; i < tags.size(); i++) fprintf(fp, "%d ", tags[i]);
1917       fprintf(fp, "\n");
1918     }
1919 
1920     for(auto it = regions.begin(); it != regions.end(); ++it) {
1921       std::vector<GFace *> const &faces = (*it)->faces();
1922       std::vector<int> const &ori = (*it)->faceOrientations();
1923       fprintf(fp, "%d ", (*it)->tag());
1924       if(partition) {
1925         partitionRegion *pr = static_cast<partitionRegion *>(*it);
1926         int parentEntityDim = 0, parentEntityTag = 0;
1927         if(pr->getParentEntity()) {
1928           parentEntityDim = pr->getParentEntity()->dim();
1929           parentEntityTag = pr->getParentEntity()->tag();
1930         }
1931         fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag);
1932 
1933         fprintf(fp, "%lu ", pr->getPartitions().size());
1934 
1935         for(auto const partition : pr->getPartitions()) {
1936           fprintf(fp, "%d ", partition);
1937         }
1938       }
1939       writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 3,
1940                            version);
1941       writeMSH4Physicals(fp, *it, binary);
1942       fprintf(fp, "%lu ", faces.size());
1943 
1944       std::vector<int> tags(faces.size());
1945       std::transform(begin(faces), end(faces), begin(tags),
1946                      [](const GFace *const face) { return face->tag(); });
1947 
1948       std::vector<int> signs(ori.begin(), ori.end());
1949 
1950       if(tags.size() == signs.size()) {
1951         for(std::size_t i = 0; i < tags.size(); i++)
1952           tags[i] *= (signs[i] > 0 ? 1 : -1);
1953       }
1954 
1955       for(auto const tag : tags) { fprintf(fp, "%d ", tag); }
1956       fprintf(fp, "\n");
1957     }
1958   }
1959 
1960   if(partition)
1961     fprintf(fp, "$EndPartitionedEntities\n");
1962   else
1963     fprintf(fp, "$EndEntities\n");
1964 }
1965 
writeMSH4EntityNodes(GEntity * ge,FILE * fp,bool binary,int saveParametric,double scalingFactor,double version)1966 static void writeMSH4EntityNodes(GEntity *ge, FILE *fp, bool binary,
1967                                  int saveParametric, double scalingFactor,
1968                                  double version)
1969 {
1970   int parametric = saveParametric;
1971   if(ge->dim() != 1 && ge->dim() != 2)
1972     parametric = 0; // Gmsh only stores parametric coordinates for dim 1 and 2
1973 
1974   if(binary) {
1975     int entityDim = ge->dim();
1976     int entityTag = ge->tag();
1977     std::size_t numVerts = ge->getNumMeshVertices();
1978     fwrite(&entityDim, sizeof(int), 1, fp);
1979     fwrite(&entityTag, sizeof(int), 1, fp);
1980     fwrite(&parametric, sizeof(int), 1, fp);
1981     fwrite(&numVerts, sizeof(std::size_t), 1, fp);
1982   }
1983   else {
1984     fprintf(fp, "%d %d %d %lu\n", (version >= 4.1) ? ge->dim() : ge->tag(),
1985             (version >= 4.1) ? ge->tag() : ge->dim(), parametric,
1986             ge->getNumMeshVertices());
1987   }
1988 
1989   std::size_t N = ge->getNumMeshVertices();
1990   std::size_t n = 3;
1991   if(parametric) n += ge->dim();
1992 
1993   if(binary) {
1994     std::vector<size_t> tags(N);
1995     for(std::size_t i = 0; i < N; i++) tags[i] = ge->getMeshVertex(i)->getNum();
1996     fwrite(&tags[0], sizeof(std::size_t), N, fp);
1997     std::vector<double> coord(n * N);
1998     std::size_t j = 0;
1999     for(std::size_t i = 0; i < N; i++) {
2000       MVertex *mv = ge->getMeshVertex(i);
2001       coord[j++] = mv->x() * scalingFactor;
2002       coord[j++] = mv->y() * scalingFactor;
2003       coord[j++] = mv->z() * scalingFactor;
2004       if(n >= 4) mv->getParameter(0, coord[j++]);
2005       if(n == 5) mv->getParameter(1, coord[j++]);
2006     }
2007     fwrite(&coord[0], sizeof(double), n * N, fp);
2008   }
2009   else {
2010     if(version >= 4.1) {
2011       for(std::size_t i = 0; i < N; i++)
2012         fprintf(fp, "%lu\n", ge->getMeshVertex(i)->getNum());
2013     }
2014     for(std::size_t i = 0; i < N; i++) {
2015       MVertex *mv = ge->getMeshVertex(i);
2016       double x = mv->x() * scalingFactor;
2017       double y = mv->y() * scalingFactor;
2018       double z = mv->z() * scalingFactor;
2019       if(version < 4.1) fprintf(fp, "%lu ", mv->getNum());
2020       if(n == 5) {
2021         double u, v;
2022         mv->getParameter(0, u);
2023         mv->getParameter(1, v);
2024         fprintf(fp, "%.16g %.16g %.16g %.16g %.16g\n", x, y, z, u, v);
2025       }
2026       else if(n == 4) {
2027         double u;
2028         mv->getParameter(0, u);
2029         fprintf(fp, "%.16g %.16g %.16g %.16g\n", x, y, z, u);
2030       }
2031       else {
2032         fprintf(fp, "%.16g %.16g %.16g\n", x, y, z);
2033       }
2034     }
2035   }
2036 }
2037 
2038 static std::size_t
getAdditionalEntities(std::set<GRegion *,GEntityPtrLessThan> & regions,std::set<GFace *,GEntityPtrLessThan> & faces,std::set<GEdge *,GEntityPtrLessThan> & edges,std::set<GVertex *,GEntityPtrLessThan> & vertices)2039 getAdditionalEntities(std::set<GRegion *, GEntityPtrLessThan> &regions,
2040                       std::set<GFace *, GEntityPtrLessThan> &faces,
2041                       std::set<GEdge *, GEntityPtrLessThan> &edges,
2042                       std::set<GVertex *, GEntityPtrLessThan> &vertices)
2043 {
2044   std::size_t numVertices = 0;
2045 
2046   for(auto it = vertices.begin(); it != vertices.end(); ++it) {
2047     numVertices += (*it)->getNumMeshVertices();
2048   }
2049 
2050   for(auto it = edges.begin(); it != edges.end(); ++it) {
2051     numVertices += (*it)->getNumMeshVertices();
2052     for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) {
2053       for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices();
2054           j++) {
2055         if((*it)->getMeshElement(i)->getVertex(j)->onWhat() != (*it)) {
2056           GEntity *entity = (*it)->getMeshElement(i)->getVertex(j)->onWhat();
2057 
2058           switch(entity->dim()) {
2059           case 0:
2060             if(vertices.find(static_cast<GVertex *>(entity)) ==
2061                vertices.end()) {
2062               vertices.insert(static_cast<GVertex *>(entity));
2063               numVertices += entity->getNumMeshVertices();
2064             }
2065             break;
2066           case 1:
2067             if(edges.find(static_cast<GEdge *>(entity)) == edges.end()) {
2068               edges.insert(static_cast<GEdge *>(entity));
2069               numVertices += entity->getNumMeshVertices();
2070             }
2071             break;
2072           case 2:
2073             if(faces.find(static_cast<GFace *>(entity)) == faces.end()) {
2074               faces.insert(static_cast<GFace *>(entity));
2075               numVertices += entity->getNumMeshVertices();
2076             }
2077             break;
2078           case 3:
2079             if(regions.find(static_cast<GRegion *>(entity)) == regions.end()) {
2080               regions.insert(static_cast<GRegion *>(entity));
2081               numVertices += entity->getNumMeshVertices();
2082             }
2083             break;
2084           default: break;
2085           }
2086         }
2087       }
2088     }
2089   }
2090 
2091   for(auto it = faces.begin(); it != faces.end(); ++it) {
2092     numVertices += (*it)->getNumMeshVertices();
2093     for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) {
2094       for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices();
2095           j++) {
2096         if((*it)->getMeshElement(i)->getVertex(j)->onWhat() != (*it)) {
2097           GEntity *entity = (*it)->getMeshElement(i)->getVertex(j)->onWhat();
2098 
2099           switch(entity->dim()) {
2100           case 0:
2101             if(vertices.find(static_cast<GVertex *>(entity)) ==
2102                vertices.end()) {
2103               vertices.insert(static_cast<GVertex *>(entity));
2104               numVertices += entity->getNumMeshVertices();
2105             }
2106             break;
2107           case 1:
2108             if(edges.find(static_cast<GEdge *>(entity)) == edges.end()) {
2109               edges.insert(static_cast<GEdge *>(entity));
2110               numVertices += entity->getNumMeshVertices();
2111             }
2112             break;
2113           case 2:
2114             if(faces.find(static_cast<GFace *>(entity)) == faces.end()) {
2115               faces.insert(static_cast<GFace *>(entity));
2116               numVertices += entity->getNumMeshVertices();
2117             }
2118             break;
2119           case 3:
2120             if(regions.find(static_cast<GRegion *>(entity)) == regions.end()) {
2121               regions.insert(static_cast<GRegion *>(entity));
2122               numVertices += entity->getNumMeshVertices();
2123             }
2124             break;
2125           default: break;
2126           }
2127         }
2128       }
2129     }
2130   }
2131 
2132   for(auto it = regions.begin(); it != regions.end(); ++it) {
2133     numVertices += (*it)->getNumMeshVertices();
2134     for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) {
2135       for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices();
2136           j++) {
2137         if((*it)->getMeshElement(i)->getVertex(j)->onWhat() != (*it)) {
2138           GEntity *entity = (*it)->getMeshElement(i)->getVertex(j)->onWhat();
2139 
2140           switch(entity->dim()) {
2141           case 0:
2142             if(vertices.find(static_cast<GVertex *>(entity)) ==
2143                vertices.end()) {
2144               vertices.insert(static_cast<GVertex *>(entity));
2145               numVertices += entity->getNumMeshVertices();
2146             }
2147             break;
2148           case 1:
2149             if(edges.find(static_cast<GEdge *>(entity)) == edges.end()) {
2150               edges.insert(static_cast<GEdge *>(entity));
2151               numVertices += entity->getNumMeshVertices();
2152             }
2153             break;
2154           case 2:
2155             if(faces.find(static_cast<GFace *>(entity)) == faces.end()) {
2156               faces.insert(static_cast<GFace *>(entity));
2157               numVertices += entity->getNumMeshVertices();
2158             }
2159             break;
2160           case 3:
2161             if(regions.find(static_cast<GRegion *>(entity)) == regions.end()) {
2162               regions.insert(static_cast<GRegion *>(entity));
2163               numVertices += entity->getNumMeshVertices();
2164             }
2165             break;
2166           default: break;
2167           }
2168         }
2169       }
2170     }
2171   }
2172 
2173   return numVertices;
2174 }
2175 
2176 static void
getEntitiesToSave(GModel * const model,bool partitioned,int partitionToSave,bool saveAll,std::set<GRegion *,GEntityPtrLessThan> & regions,std::set<GFace *,GEntityPtrLessThan> & faces,std::set<GEdge *,GEntityPtrLessThan> & edges,std::set<GVertex *,GEntityPtrLessThan> & vertices)2177 getEntitiesToSave(GModel *const model, bool partitioned,
2178                   int partitionToSave, bool saveAll,
2179                   std::set<GRegion *, GEntityPtrLessThan> &regions,
2180                   std::set<GFace *, GEntityPtrLessThan> &faces,
2181                   std::set<GEdge *, GEntityPtrLessThan> &edges,
2182                   std::set<GVertex *, GEntityPtrLessThan> &vertices)
2183 {
2184   if(partitioned) {
2185     for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
2186       if((*it)->geomType() == GEntity::PartitionPoint) {
2187         partitionVertex *pv = static_cast<partitionVertex *>(*it);
2188         if(!partitionToSave ||
2189            std::find(pv->getPartitions().begin(), pv->getPartitions().end(),
2190                      partitionToSave) != pv->getPartitions().end())
2191           vertices.insert(pv);
2192       }
2193     }
2194     for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
2195       if((*it)->geomType() == GEntity::PartitionCurve) {
2196         partitionEdge *pe = static_cast<partitionEdge *>(*it);
2197         if(!partitionToSave ||
2198            std::find(pe->getPartitions().begin(), pe->getPartitions().end(),
2199                      partitionToSave) != pe->getPartitions().end())
2200           edges.insert(pe);
2201       }
2202       else if((*it)->geomType() == GEntity::GhostCurve) {
2203         ghostEdge *ge = static_cast<ghostEdge *>(*it);
2204         if(ge->getPartition() == partitionToSave)
2205           edges.insert(ge);
2206       }
2207     }
2208     for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
2209       if((*it)->geomType() == GEntity::PartitionSurface) {
2210         partitionFace *pf = static_cast<partitionFace *>(*it);
2211         if(!partitionToSave ||
2212            std::find(pf->getPartitions().begin(), pf->getPartitions().end(),
2213                      partitionToSave) != pf->getPartitions().end())
2214           faces.insert(pf);
2215       }
2216       else if((*it)->geomType() == GEntity::GhostSurface) {
2217         ghostFace *gf = static_cast<ghostFace *>(*it);
2218         if(gf->getPartition() == partitionToSave)
2219           faces.insert(gf);
2220       }
2221     }
2222     for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
2223       if((*it)->geomType() == GEntity::PartitionVolume) {
2224         partitionRegion *pr = static_cast<partitionRegion *>(*it);
2225         if(!partitionToSave ||
2226            std::find(pr->getPartitions().begin(), pr->getPartitions().end(),
2227                      partitionToSave) != pr->getPartitions().end())
2228           regions.insert(pr);
2229       }
2230       else if((*it)->geomType() == GEntity::GhostVolume) {
2231         ghostRegion *gr = static_cast<ghostRegion *>(*it);
2232         if(gr->getPartition() == partitionToSave)
2233           regions.insert(gr);
2234       }
2235     }
2236   }
2237   else {
2238     for(auto it = model->firstVertex(); it != model->lastVertex(); ++it)
2239       if((*it)->geomType() != GEntity::PartitionPoint &&
2240          (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0)))
2241         vertices.insert(*it);
2242     for(auto it = model->firstEdge(); it != model->lastEdge(); ++it)
2243       if((*it)->geomType() != GEntity::PartitionCurve &&
2244          (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0) ||
2245           (*it)->geomType() == GEntity::GhostCurve))
2246         edges.insert(*it);
2247     for(auto it = model->firstFace(); it != model->lastFace(); ++it)
2248       if((*it)->geomType() != GEntity::PartitionSurface &&
2249          (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0) ||
2250           (*it)->geomType() == GEntity::GhostSurface))
2251         faces.insert(*it);
2252     for(auto it = model->firstRegion(); it != model->lastRegion(); ++it)
2253       if((*it)->geomType() != GEntity::PartitionVolume &&
2254          (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0) ||
2255           (*it)->geomType() == GEntity::GhostVolume))
2256         regions.insert(*it);
2257   }
2258 }
2259 
writeMSH4Nodes(GModel * const model,FILE * fp,bool partitioned,int partitionToSave,bool binary,int saveParametric,double scalingFactor,bool saveAll,double version)2260 static void writeMSH4Nodes(GModel *const model, FILE *fp, bool partitioned,
2261                            int partitionToSave, bool binary, int saveParametric,
2262                            double scalingFactor, bool saveAll, double version)
2263 {
2264   std::set<GRegion *, GEntityPtrLessThan> regions;
2265   std::set<GFace *, GEntityPtrLessThan> faces;
2266   std::set<GEdge *, GEntityPtrLessThan> edges;
2267   std::set<GVertex *, GEntityPtrLessThan> vertices;
2268   getEntitiesToSave(model, partitioned, partitionToSave, saveAll, regions,
2269                     faces, edges, vertices);
2270 
2271   std::size_t numNodes = (saveAll && !partitioned) ?
2272     model->getNumMeshVertices() :
2273     getAdditionalEntities(regions, faces, edges, vertices);
2274 
2275   if(!numNodes) return;
2276 
2277   fprintf(fp, "$Nodes\n");
2278 
2279   std::size_t minTag = std::numeric_limits<std::size_t>::max(), maxTag = 0;
2280   for(auto it = vertices.begin(); it != vertices.end(); ++it) {
2281     for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) {
2282       minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum());
2283       maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum());
2284     }
2285   }
2286   for(auto it = edges.begin(); it != edges.end(); ++it) {
2287     for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) {
2288       minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum());
2289       maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum());
2290     }
2291   }
2292   for(auto it = faces.begin(); it != faces.end(); ++it) {
2293     for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) {
2294       minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum());
2295       maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum());
2296     }
2297   }
2298   for(auto it = regions.begin(); it != regions.end(); ++it) {
2299     for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) {
2300       minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum());
2301       maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum());
2302     }
2303   }
2304 
2305   if(binary) {
2306     std::size_t numSection =
2307       vertices.size() + edges.size() + faces.size() + regions.size();
2308     fwrite(&numSection, sizeof(std::size_t), 1, fp);
2309     fwrite(&numNodes, sizeof(std::size_t), 1, fp);
2310     fwrite(&minTag, sizeof(std::size_t), 1, fp);
2311     fwrite(&maxTag, sizeof(std::size_t), 1, fp);
2312   }
2313   else {
2314     if(version >= 4.1) {
2315       fprintf(fp, "%lu %lu %lu %lu\n",
2316               vertices.size() + edges.size() + faces.size() + regions.size(),
2317               numNodes, minTag, maxTag);
2318     }
2319     else {
2320       fprintf(fp, "%lu %lu\n",
2321               vertices.size() + edges.size() + faces.size() + regions.size(),
2322               numNodes);
2323     }
2324   }
2325 
2326   for(auto it = vertices.begin(); it != vertices.end(); ++it) {
2327     writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor,
2328                          version);
2329   }
2330   for(auto it = edges.begin(); it != edges.end(); ++it) {
2331     writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor,
2332                          version);
2333   }
2334   for(auto it = faces.begin(); it != faces.end(); ++it) {
2335     writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor,
2336                          version);
2337   }
2338   for(auto it = regions.begin(); it != regions.end(); ++it) {
2339     writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor,
2340                          version);
2341   }
2342 
2343   if(binary) fprintf(fp, "\n");
2344 
2345   fprintf(fp, "$EndNodes\n");
2346 }
2347 
writeMSH4Elements(GModel * const model,FILE * fp,bool partitioned,int partitionToSave,bool binary,bool saveAll,double version)2348 static void writeMSH4Elements(GModel *const model, FILE *fp, bool partitioned,
2349                               int partitionToSave, bool binary, bool saveAll,
2350                               double version)
2351 {
2352   std::set<GRegion *, GEntityPtrLessThan> regions;
2353   std::set<GFace *, GEntityPtrLessThan> faces;
2354   std::set<GEdge *, GEntityPtrLessThan> edges;
2355   std::set<GVertex *, GEntityPtrLessThan> vertices;
2356   getEntitiesToSave(model, partitioned, partitionToSave, saveAll, regions,
2357                     faces, edges, vertices);
2358 
2359   std::map<std::pair<int, int>, std::vector<MElement *> > elementsByType[4];
2360   std::size_t numElements = 0;
2361 
2362   for(auto it = vertices.begin(); it != vertices.end(); ++it) {
2363     if(!saveAll && (*it)->physicals.size() == 0) continue;
2364 
2365     numElements += (*it)->points.size();
2366     for(std::size_t i = 0; i < (*it)->points.size(); i++) {
2367       std::pair<int, int> p((*it)->tag(), (*it)->points[i]->getTypeForMSH());
2368       elementsByType[0][p].push_back((*it)->points[i]);
2369     }
2370   }
2371 
2372   for(auto it = edges.begin(); it != edges.end(); ++it) {
2373     if(!saveAll && (*it)->physicals.size() == 0 &&
2374        (*it)->geomType() != GEntity::GhostCurve)
2375       continue;
2376 
2377     numElements += (*it)->lines.size();
2378     for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
2379       std::pair<int, int> p((*it)->tag(), (*it)->lines[i]->getTypeForMSH());
2380       elementsByType[1][p].push_back((*it)->lines[i]);
2381     }
2382   }
2383 
2384   for(auto it = faces.begin(); it != faces.end(); ++it) {
2385     if(!saveAll && (*it)->physicals.size() == 0 &&
2386        (*it)->geomType() != GEntity::GhostSurface)
2387       continue;
2388 
2389     numElements += (*it)->triangles.size();
2390     for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
2391       std::pair<int, int> p((*it)->tag(), (*it)->triangles[i]->getTypeForMSH());
2392       elementsByType[2][p].push_back((*it)->triangles[i]);
2393     }
2394     numElements += (*it)->quadrangles.size();
2395     for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++) {
2396       std::pair<int, int> p((*it)->tag(),
2397                             (*it)->quadrangles[i]->getTypeForMSH());
2398       elementsByType[2][p].push_back((*it)->quadrangles[i]);
2399     }
2400   }
2401 
2402   for(auto it = regions.begin(); it != regions.end(); ++it) {
2403     if(!saveAll && (*it)->physicals.size() == 0 &&
2404        (*it)->geomType() != GEntity::GhostVolume)
2405       continue;
2406 
2407     numElements += (*it)->tetrahedra.size();
2408     for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++) {
2409       std::pair<int, int> p((*it)->tag(),
2410                             (*it)->tetrahedra[i]->getTypeForMSH());
2411       elementsByType[3][p].push_back((*it)->tetrahedra[i]);
2412     }
2413     numElements += (*it)->hexahedra.size();
2414     for(std::size_t i = 0; i < (*it)->hexahedra.size(); i++) {
2415       std::pair<int, int> p((*it)->tag(), (*it)->hexahedra[i]->getTypeForMSH());
2416       elementsByType[3][p].push_back((*it)->hexahedra[i]);
2417     }
2418     numElements += (*it)->prisms.size();
2419     for(std::size_t i = 0; i < (*it)->prisms.size(); i++) {
2420       std::pair<int, int> p((*it)->tag(), (*it)->prisms[i]->getTypeForMSH());
2421       elementsByType[3][p].push_back((*it)->prisms[i]);
2422     }
2423     numElements += (*it)->pyramids.size();
2424     for(std::size_t i = 0; i < (*it)->pyramids.size(); i++) {
2425       std::pair<int, int> p((*it)->tag(), (*it)->pyramids[i]->getTypeForMSH());
2426       elementsByType[3][p].push_back((*it)->pyramids[i]);
2427     }
2428     numElements += (*it)->trihedra.size();
2429     for(std::size_t i = 0; i < (*it)->trihedra.size(); i++) {
2430       std::pair<int, int> p((*it)->tag(), (*it)->trihedra[i]->getTypeForMSH());
2431       elementsByType[3][p].push_back((*it)->trihedra[i]);
2432     }
2433   }
2434 
2435   if(!numElements) return;
2436 
2437   fprintf(fp, "$Elements\n");
2438 
2439   std::size_t numSection = 0;
2440   for(int dim = 0; dim <= 3; dim++) numSection += elementsByType[dim].size();
2441 
2442   std::size_t minTag = std::numeric_limits<std::size_t>::max(), maxTag = 0;
2443   for(int dim = 0; dim <= 3; dim++) {
2444     for(auto it = elementsByType[dim].begin(); it != elementsByType[dim].end();
2445         ++it) {
2446       for(std::size_t i = 0; i < it->second.size(); i++) {
2447         minTag = std::min(minTag, it->second[i]->getNum());
2448         maxTag = std::max(maxTag, it->second[i]->getNum());
2449       }
2450     }
2451   }
2452 
2453   if(binary) {
2454     fwrite(&numSection, sizeof(std::size_t), 1, fp);
2455     fwrite(&numElements, sizeof(std::size_t), 1, fp);
2456     fwrite(&minTag, sizeof(std::size_t), 1, fp);
2457     fwrite(&maxTag, sizeof(std::size_t), 1, fp);
2458   }
2459   else {
2460     if(version >= 4.1)
2461       fprintf(fp, "%lu %lu %lu %lu\n", numSection, numElements, minTag, maxTag);
2462     else
2463       fprintf(fp, "%lu %lu\n", numSection, numElements);
2464   }
2465 
2466   for(int dim = 0; dim <= 3; dim++) {
2467     for(auto it = elementsByType[dim].begin(); it != elementsByType[dim].end();
2468         ++it) {
2469       int entityTag = it->first.first;
2470       int elmType = it->first.second;
2471       std::size_t numElm = it->second.size();
2472       if(binary) {
2473         fwrite(&dim, sizeof(int), 1, fp);
2474         fwrite(&entityTag, sizeof(int), 1, fp);
2475         fwrite(&elmType, sizeof(int), 1, fp);
2476         fwrite(&numElm, sizeof(std::size_t), 1, fp);
2477       }
2478       else {
2479         fprintf(fp, "%d %d %d %lu\n", (version >= 4.1) ? dim : entityTag,
2480                 (version >= 4.1) ? entityTag : dim, elmType, numElm);
2481       }
2482 
2483       std::size_t N = it->second.size();
2484       if(binary) {
2485         const int numVertPerElm = MElement::getInfoMSH(elmType);
2486         std::size_t n = 1 + numVertPerElm;
2487         std::vector<std::size_t> tags(N * n);
2488         std::size_t k = 0;
2489         for(std::size_t i = 0; i < N; i++) {
2490           MElement *e = it->second[i];
2491           tags[k] = e->getNum();
2492           for(int j = 0; j < numVertPerElm; j++) {
2493             tags[k + 1 + j] = e->getVertex(j)->getNum();
2494           }
2495           k += n;
2496         }
2497         fwrite(&tags[0], sizeof(std::size_t), N * n, fp);
2498       }
2499       else {
2500         for(std::size_t i = 0; i < N; i++) {
2501           MElement *e = it->second[i];
2502           fprintf(fp, "%lu ", e->getNum());
2503           for(std::size_t i = 0; i < e->getNumVertices(); i++) {
2504             fprintf(fp, "%lu ", e->getVertex(i)->getNum());
2505           }
2506           fprintf(fp, "\n");
2507         }
2508       }
2509     }
2510   }
2511 
2512   if(binary) fprintf(fp, "\n");
2513 
2514   fprintf(fp, "$EndElements\n");
2515 }
2516 
writeMSH4PeriodicNodes(GModel * const model,FILE * fp,bool binary,double version)2517 static void writeMSH4PeriodicNodes(GModel *const model, FILE *fp,
2518                                    bool binary, double version)
2519 {
2520   // To avoid saving correspondances bwteen nodes that are not saved (either in
2521   // the same file or not at all, e.g. in the partitioned case, or simply if
2522   // some physical entities are not defined), we could only apply the code below
2523   // to the entities returned by getEntitiesForNodes().
2524 
2525   // The current choice is to save everything, and not complain if a node is not
2526   // found when reading the file. This should be reevaluated at some point.
2527 
2528   std::size_t count = 0;
2529   std::vector<GEntity *> entities;
2530   model->getEntities(entities);
2531   for(std::size_t i = 0; i < entities.size(); i++)
2532     if(entities[i]->getMeshMaster() != entities[i]) count++;
2533 
2534   if(!count) return;
2535 
2536   fprintf(fp, "$Periodic\n");
2537 
2538   if(binary) { fwrite(&count, sizeof(std::size_t), 1, fp); }
2539   else {
2540     fprintf(fp, "%lu\n", count);
2541   }
2542 
2543   for(std::size_t i = 0; i < entities.size(); i++) {
2544     GEntity *g_slave = entities[i];
2545     GEntity *g_master = g_slave->getMeshMaster();
2546 
2547     if(g_slave != g_master) {
2548       std::map<MVertex *, MVertex *, MVertexPtrLessThan> corrVert = g_slave->correspondingVertices;
2549       if(CTX::instance()->mesh.hoSavePeriodic)
2550         corrVert.insert(g_slave->correspondingHighOrderVertices.begin(),
2551                         g_slave->correspondingHighOrderVertices.end());
2552 
2553       if(binary) {
2554         int gSlaveDim = g_slave->dim();
2555         int gSlaveTag = g_slave->tag();
2556         int gMasterTag = g_master->tag();
2557         fwrite(&gSlaveDim, sizeof(int), 1, fp);
2558         fwrite(&gSlaveTag, sizeof(int), 1, fp);
2559         fwrite(&gMasterTag, sizeof(int), 1, fp);
2560 
2561         if(g_slave->affineTransform.size() == 16) {
2562           std::size_t numAffine = 16;
2563           fwrite(&numAffine, sizeof(std::size_t), 1, fp);
2564           for(int j = 0; j < 16; j++) {
2565             double value = g_slave->affineTransform[j];
2566             fwrite(&value, sizeof(double), 1, fp);
2567           }
2568         }
2569         else {
2570           std::size_t numAffine = 0;
2571           fwrite(&numAffine, sizeof(std::size_t), 1, fp);
2572         }
2573 
2574         std::size_t corrVertSize = corrVert.size();
2575         fwrite(&corrVertSize, sizeof(std::size_t), 1, fp);
2576 
2577         for(auto it = corrVert.begin(); it != corrVert.end(); ++it) {
2578           std::size_t numFirst = it->first->getNum();
2579           std::size_t numSecond = it->second->getNum();
2580           fwrite(&numFirst, sizeof(std::size_t), 1, fp);
2581           fwrite(&numSecond, sizeof(std::size_t), 1, fp);
2582         }
2583       }
2584       else {
2585         fprintf(fp, "%d %d %d\n", g_slave->dim(), g_slave->tag(),
2586                 g_master->tag());
2587 
2588         if(version >= 4.1) {
2589           if(g_slave->affineTransform.size() == 16) {
2590             fprintf(fp, "16");
2591             for(int j = 0; j < 16; j++)
2592               fprintf(fp, " %.16g", g_slave->affineTransform[j]);
2593             fprintf(fp, "\n");
2594           }
2595           else {
2596             fprintf(fp, "0\n");
2597           }
2598         }
2599         else {
2600           if(g_slave->affineTransform.size() == 16) {
2601             fprintf(fp, "Affine");
2602             for(int j = 0; j < 16; j++)
2603               fprintf(fp, " %.16g", g_slave->affineTransform[j]);
2604             fprintf(fp, "\n");
2605           }
2606         }
2607 
2608         fprintf(fp, "%lu\n", corrVert.size());
2609 
2610         for(auto it = corrVert.begin(); it != corrVert.end(); ++it) {
2611           fprintf(fp, "%lu %lu\n", it->first->getNum(), it->second->getNum());
2612         }
2613       }
2614     }
2615   }
2616 
2617   if(binary) fprintf(fp, "\n");
2618   fprintf(fp, "$EndPeriodic\n");
2619 }
2620 
writeMSH4GhostCells(GModel * const model,FILE * fp,int partitionToSave,bool binary)2621 static void writeMSH4GhostCells(GModel *const model, FILE *fp,
2622                                 int partitionToSave, bool binary)
2623 {
2624   std::vector<GEntity *> entities;
2625   model->getEntities(entities);
2626   std::map<MElement *, std::vector<int> > ghostCells;
2627 
2628   for(std::size_t i = 0; i < entities.size(); i++) {
2629     std::map<MElement *, int> ghostElements;
2630     int partition;
2631 
2632     if(entities[i]->geomType() == GEntity::GhostCurve) {
2633       ghostElements = static_cast<ghostEdge *>(entities[i])->getGhostCells();
2634       partition = static_cast<ghostEdge *>(entities[i])->getPartition();
2635     }
2636     else if(entities[i]->geomType() == GEntity::GhostSurface) {
2637       ghostElements = static_cast<ghostFace *>(entities[i])->getGhostCells();
2638       partition = static_cast<ghostFace *>(entities[i])->getPartition();
2639     }
2640     else if(entities[i]->geomType() == GEntity::GhostVolume) {
2641       ghostElements = static_cast<ghostRegion *>(entities[i])->getGhostCells();
2642       partition = static_cast<ghostRegion *>(entities[i])->getPartition();
2643     }
2644 
2645     if(!partitionToSave || partitionToSave == partition) {
2646       for(auto it = ghostElements.begin(); it != ghostElements.end(); ++it) {
2647         if(ghostCells[it->first].size() == 0)
2648           ghostCells[it->first].push_back(it->second);
2649         ghostCells[it->first].push_back(partition);
2650       }
2651     }
2652   }
2653 
2654   if(ghostCells.size() != 0) {
2655     fprintf(fp, "$GhostElements\n");
2656     if(binary) {
2657       std::size_t ghostCellsSize = ghostCells.size();
2658       fwrite(&ghostCellsSize, sizeof(std::size_t), 1, fp);
2659 
2660       for(auto it = ghostCells.begin(); it != ghostCells.end(); ++it) {
2661         std::size_t elmTag = it->first->getNum();
2662         int partNum = it->second[0];
2663         std::size_t numGhostPartitions = it->second.size() - 1;
2664         fwrite(&elmTag, sizeof(std::size_t), 1, fp);
2665         fwrite(&partNum, sizeof(int), 1, fp);
2666         fwrite(&numGhostPartitions, sizeof(std::size_t), 1, fp);
2667         for(std::size_t i = 1; i < it->second.size(); i++) {
2668           fwrite(&it->second[i], sizeof(int), 1, fp);
2669         }
2670       }
2671       fprintf(fp, "\n");
2672     }
2673     else {
2674       fprintf(fp, "%ld\n", ghostCells.size());
2675 
2676       for(auto it = ghostCells.begin(); it != ghostCells.end(); ++it) {
2677         fprintf(fp, "%lu %d %ld", it->first->getNum(), it->second[0],
2678                 it->second.size() - 1);
2679         for(std::size_t i = 1; i < it->second.size(); i++) {
2680           fprintf(fp, " %d", it->second[i]);
2681         }
2682         fprintf(fp, "\n");
2683       }
2684     }
2685     fprintf(fp, "$EndGhostElements\n");
2686   }
2687 }
2688 
writeMSH4Parametrizations(GModel * const model,FILE * fp,bool binary)2689 static void writeMSH4Parametrizations(GModel *const model, FILE *fp,
2690                                       bool binary)
2691 {
2692   std::size_t nParamE = 0, nParamF = 0;
2693 
2694   for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
2695     discreteEdge *de = dynamic_cast<discreteEdge *>(*it);
2696     if(de && de->haveParametrization()) { nParamE++; }
2697   }
2698   for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
2699     discreteFace *df = dynamic_cast<discreteFace *>(*it);
2700     if(df && df->haveParametrization()) { nParamF++; }
2701   }
2702 
2703   if(!nParamE && !nParamF) return;
2704 
2705   fprintf(fp, "$Parametrizations\n");
2706 
2707   if(binary) {
2708     fwrite(&nParamE, sizeof(std::size_t), 1, fp);
2709     fwrite(&nParamF, sizeof(std::size_t), 1, fp);
2710   }
2711   else {
2712     fprintf(fp, "%lu %lu\n", nParamE, nParamF);
2713   }
2714 
2715   for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
2716     discreteEdge *de = dynamic_cast<discreteEdge *>(*it);
2717     if(de && de->haveParametrization()) {
2718       int t = de->tag();
2719       if(binary)
2720         fwrite(&t, sizeof(int), 1, fp);
2721       else
2722         fprintf(fp, "%d\n", t);
2723       de->writeParametrization(fp, binary);
2724     }
2725   }
2726   for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
2727     discreteFace *df = dynamic_cast<discreteFace *>(*it);
2728     if(df && df->haveParametrization()) {
2729       int t = df->tag();
2730       if(binary)
2731         fwrite(&t, sizeof(int), 1, fp);
2732       else
2733         fprintf(fp, "%d\n", t);
2734       df->writeParametrization(fp, binary);
2735     }
2736   }
2737 
2738   if(binary) fprintf(fp, "\n");
2739 
2740   fprintf(fp, "$EndParametrizations\n");
2741 }
2742 
_writeMSH4(const std::string & name,double version,bool binary,bool saveAll,bool saveParametric,double scalingFactor,bool append,int partitionToSave)2743 int GModel::_writeMSH4(const std::string &name, double version, bool binary,
2744                        bool saveAll, bool saveParametric, double scalingFactor,
2745                        bool append, int partitionToSave)
2746 {
2747   FILE *fp = nullptr;
2748   if(append)
2749     fp = Fopen(name.c_str(), binary ? "ab" : "a");
2750   else
2751     fp = Fopen(name.c_str(), binary ? "wb" : "w");
2752 
2753   if(!fp) {
2754     Msg::Error("Unable to open file '%s'", name.c_str());
2755     return 0;
2756   }
2757 
2758   if(version < 4.1 && binary) {
2759     Msg::Error("Can only write MSH 4.0 format in ASCII mode");
2760     return 0;
2761   }
2762 
2763   // if there are no physicals we save all the elements
2764   if(noPhysicalGroups()) saveAll = true;
2765 
2766   // header
2767   fprintf(fp, "$MeshFormat\n");
2768   fprintf(fp, "%g %d %lu\n", version, (binary ? 1 : 0), sizeof(std::size_t));
2769   if(binary) {
2770     int one = 1;
2771     fwrite(&one, sizeof(int), 1, fp); // swapping byte
2772     fprintf(fp, "\n");
2773   }
2774   fprintf(fp, "$EndMeshFormat\n");
2775 
2776   // physicals
2777   if(numPhysicalNames() > 0) {
2778     fprintf(fp, "$PhysicalNames\n");
2779     fprintf(fp, "%d\n", numPhysicalNames());
2780     for(auto it = firstPhysicalName(); it != lastPhysicalName(); ++it) {
2781       std::string name = it->second;
2782       if(name.size() > 128) name.resize(128);
2783       fprintf(fp, "%d %d \"%s\"\n", it->first.first, it->first.second,
2784               name.c_str());
2785     }
2786     fprintf(fp, "$EndPhysicalNames\n");
2787   }
2788 
2789   // entities
2790   writeMSH4Entities(this, fp, false, binary, scalingFactor, version);
2791 
2792   // check if the mesh is partitioned... and if we actually have elements in the
2793   // partitioned entities
2794   bool partitioned = getNumPartitions() > 0;
2795   if(partitioned) {
2796     std::vector<GEntity *> entities;
2797     getEntities(entities);
2798     std::size_t partEnt = 0;
2799     for(auto &ge : entities) {
2800       if(ge->geomType() == GEntity::PartitionPoint ||
2801          ge->geomType() == GEntity::PartitionCurve ||
2802          ge->geomType() == GEntity::PartitionSurface ||
2803          ge->geomType() == GEntity::PartitionVolume ||
2804          ge->geomType() == GEntity::GhostCurve ||
2805          ge->geomType() == GEntity::GhostSurface ||
2806          ge->geomType() == GEntity::GhostVolume)
2807         partEnt++;
2808     }
2809     if(!partEnt) {
2810       // this can happen when e.g. loading an old MSH2 files with partition tags
2811       // stored in elements
2812       Msg::Warning("No partition entities found, saving mesh as unpartitioned");
2813       partitioned = false;
2814     }
2815   }
2816 
2817   // partitioned entities
2818   if(partitioned)
2819     writeMSH4Entities(this, fp, true, binary, scalingFactor, version);
2820 
2821   // nodes
2822   writeMSH4Nodes(this, fp, partitioned, partitionToSave, binary,
2823                  saveParametric ? 1 : 0, scalingFactor, saveAll, version);
2824 
2825   // elements
2826   writeMSH4Elements(this, fp, partitioned, partitionToSave, binary, saveAll,
2827                     version);
2828 
2829   // periodic
2830   writeMSH4PeriodicNodes(this, fp, binary, version);
2831 
2832   // ghostCells
2833   writeMSH4GhostCells(this, fp, partitionToSave, binary);
2834 
2835   // parametrizations
2836   writeMSH4Parametrizations(this, fp, binary);
2837 
2838   fclose(fp);
2839 
2840   return 1;
2841 }
2842 
_writePartitionedMSH4(const std::string & baseName,double version,bool binary,bool saveAll,bool saveParametric,double scalingFactor)2843 int GModel::_writePartitionedMSH4(const std::string &baseName, double version,
2844                                   bool binary, bool saveAll,
2845                                   bool saveParametric, double scalingFactor)
2846 {
2847 #pragma omp parallel for
2848   for(std::size_t part = 1; part <= getNumPartitions(); part++) {
2849     std::ostringstream sstream;
2850     sstream << baseName << "_" << part << ".msh";
2851     if(getNumPartitions() > 100) {
2852       if(part % 100 == 1) {
2853         Msg::Info("Writing partition %d/%d in file '%s'", part, getNumPartitions(),
2854                   sstream.str().c_str());
2855       }
2856     }
2857     else {
2858       Msg::Info("Writing partition %d in file '%s'", part, sstream.str().c_str());
2859     }
2860     _writeMSH4(sstream.str(), version, binary, saveAll, saveParametric,
2861                scalingFactor, false, part);
2862   }
2863 
2864   return 1;
2865 }
2866 
getPhyscialNameInfo(const std::string & name,int & parentPhysicalTag,std::vector<int> & partitions)2867 static bool getPhyscialNameInfo(const std::string &name, int &parentPhysicalTag,
2868                                 std::vector<int> &partitions)
2869 {
2870   if(name[0] != '_') return false;
2871 
2872   const std::string part = "_part{";
2873   const std::string physical = "_physical{";
2874 
2875   size_t firstPart = name.find(part) + part.size();
2876   size_t lastPart = name.find_first_of('}', firstPart);
2877   const std::string partString = name.substr(firstPart, lastPart - firstPart);
2878 
2879   size_t firstPhysical = name.find(physical) + physical.size();
2880   size_t lastPhysical = name.find_first_of('}', firstPhysical);
2881   const std::string physicalString =
2882     name.substr(firstPhysical, lastPhysical - firstPhysical);
2883 
2884   std::string number;
2885   for(size_t i = 0; i < partString.size(); ++i) {
2886     if(partString[i] == ',') {
2887       partitions.push_back(atoi(number.c_str()));
2888       number.clear();
2889     }
2890     else {
2891       number += partString[i];
2892     }
2893   }
2894   partitions.push_back(atoi(number.c_str()));
2895 
2896   parentPhysicalTag = atoi(physicalString.c_str());
2897 
2898   return true;
2899 }
2900 
writePartitionedTopology(std::string & name)2901 int GModel::writePartitionedTopology(std::string &name)
2902 {
2903   Msg::Info("Writing '%s'", name.c_str());
2904 
2905   std::vector<std::map<int, std::pair<int, std::vector<int> > > > allParts(4);
2906   std::vector<GEntity *> entities;
2907   getEntities(entities);
2908   for(size_t i = 0; i < entities.size(); i++) {
2909     std::vector<int> physicals = entities[i]->getPhysicalEntities();
2910     for(size_t j = 0; j < physicals.size(); ++j) {
2911       const std::string phyName =
2912         this->getPhysicalName(entities[i]->dim(), physicals[j]);
2913       int parentPhysicalTag;
2914       std::vector<int> partitions;
2915       if(getPhyscialNameInfo(phyName, parentPhysicalTag, partitions)) {
2916         allParts[entities[i]->dim()].insert
2917           (std::make_pair(physicals[j],
2918                           std::make_pair(parentPhysicalTag, partitions)));
2919       }
2920     }
2921   }
2922 
2923   FILE *fp = Fopen(name.c_str(), "w");
2924   if(!fp) {
2925     Msg::Error("Could not open file '%s'", name.c_str());
2926     return 0;
2927   }
2928 
2929   fprintf(fp, "Group{\n");
2930   fprintf(fp, "  // Part~{dim}~{parentPhysicalTag}~{part1}~{part2}~...\n\n");
2931   std::vector<std::map<int, std::string> > tagToString(4);
2932   for(size_t i = 4; i > 0; --i) {
2933     fprintf(fp, "  // Dim %lu\n", i - 1);
2934     for(auto it = allParts[i - 1].begin(); it != allParts[i - 1].end(); ++it) {
2935       std::string partName = "Part~{" + std::to_string(i - 1) + "}~{" +
2936                              std::to_string(it->second.first) + "}";
2937       fprintf(fp, "  Part~{%lu}~{%d}", i - 1, it->second.first);
2938       for(size_t j = 0; j < it->second.second.size(); ++j) {
2939         partName += "~{" + std::to_string(it->second.second[j]) + "}";
2940         fprintf(fp, "~{%d}", it->second.second[j]);
2941       }
2942       tagToString[i - 1].insert(std::make_pair(it->first, partName));
2943       fprintf(fp, " = Region[{%d}];\n", it->first);
2944     }
2945     fprintf(fp, "\n");
2946   }
2947 
2948   fprintf(fp, "  // Global names\n\n");
2949   std::map<int, std::vector<int> > omegas;
2950   std::map<std::pair<int, int>, std::vector<int> > sigmasij;
2951   std::map<int, std::vector<int> > sigmas;
2952   std::map<int, std::set<int> > neighbors;
2953   std::size_t omegaDim = 0;
2954   for(size_t i = 4; i > 0; --i) {
2955     if(allParts[i - 1].size() != 0) {
2956       omegaDim = i - 1;
2957       break;
2958     }
2959   }
2960 
2961   // omega
2962   for(auto it = allParts[omegaDim].begin(); it != allParts[omegaDim].end();
2963       ++it) {
2964     if(it->second.second.size() == 1) {
2965       omegas[it->second.second[0]].push_back(it->first);
2966     }
2967   }
2968   fprintf(fp, "  // Omega\n");
2969   for(auto it = omegas.begin(); it != omegas.end(); ++it) {
2970     fprintf(fp, "  Omega~{%d} = Region[{", it->first);
2971     for(size_t j = 0; j < it->second.size(); ++j) {
2972       if(j == 0)
2973         fprintf(fp, "%s", tagToString[omegaDim][it->second[j]].c_str());
2974       else
2975         fprintf(fp, ", %s", tagToString[omegaDim][it->second[j]].c_str());
2976     }
2977     fprintf(fp, "}];\n");
2978   }
2979   fprintf(fp, "\n");
2980 
2981   if(omegaDim > 0) {
2982     // sigma
2983     for(auto it = allParts[omegaDim - 1].begin();
2984         it != allParts[omegaDim - 1].end(); ++it) {
2985       if(it->second.second.size() == 2) {
2986         sigmasij[std::make_pair(it->second.second[0],
2987                                 it->second.second[1])]
2988           .push_back(it->first);
2989         sigmasij[std::make_pair(it->second.second[1],
2990                                 it->second.second[0])]
2991           .push_back(it->first);
2992         sigmas[it->second.second[0]].push_back(it->first);
2993         sigmas[it->second.second[1]].push_back(it->first);
2994       }
2995     }
2996     fprintf(fp, "  // Sigma\n");
2997     for(auto it = sigmasij.begin(); it != sigmasij.end(); ++it) {
2998       fprintf(fp, "  Sigma~{%d}~{%d} = Region[{", it->first.first,
2999               it->first.second);
3000       for(size_t j = 0; j < it->second.size(); ++j) {
3001         if(j == 0)
3002           fprintf(fp, "%s", tagToString[omegaDim - 1][it->second[j]].c_str());
3003         else
3004           fprintf(fp, ", %s", tagToString[omegaDim - 1][it->second[j]].c_str());
3005       }
3006       fprintf(fp, "}];\n");
3007     }
3008     fprintf(fp, "\n");
3009 
3010     for(auto it = sigmas.begin(); it != sigmas.end(); ++it) {
3011       fprintf(fp, "  Sigma~{%d} = Region[{", it->first);
3012       for(size_t j = 0; j < it->second.size(); ++j) {
3013         if(j == 0)
3014           fprintf(fp, "%s", tagToString[omegaDim - 1][it->second[j]].c_str());
3015         else
3016           fprintf(fp, ", %s", tagToString[omegaDim - 1][it->second[j]].c_str());
3017       }
3018       fprintf(fp, "}];\n");
3019     }
3020     fprintf(fp, "\n");
3021   }
3022 
3023   // D
3024   fprintf(fp, "  D() = {");
3025   for(size_t i = 1; i <= getNumPartitions(); ++i) {
3026     if(i != 1) fprintf(fp, ", ");
3027     fprintf(fp, "%lu", i);
3028   }
3029   fprintf(fp, "};\n");
3030 
3031   if(omegaDim > 0) {
3032     // D~{i}
3033     for(auto it = allParts[omegaDim - 1].begin();
3034         it != allParts[omegaDim - 1].end(); ++it) {
3035       if(it->second.second.size() == 2) {
3036         neighbors[it->second.second[0]].insert(it->second.second[1]);
3037         neighbors[it->second.second[1]].insert(it->second.second[0]);
3038       }
3039     }
3040     for(size_t i = 1; i <= getNumPartitions(); ++i) {
3041       fprintf(fp, "  D~{%lu}() = {", i);
3042       for(auto it = neighbors[i].begin(); it != neighbors[i].end(); ++it) {
3043         if(it != neighbors[i].begin()) fprintf(fp, ", ");
3044         fprintf(fp, "%d", *it);
3045       }
3046       fprintf(fp, "};\n");
3047     }
3048   }
3049 
3050   fprintf(fp, "}\n\n");
3051 
3052   fclose(fp);
3053 
3054   Msg::Info("Done writing '%s'", name.c_str());
3055 
3056   return 1;
3057 }
3058