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, ¶metric,
527 &numNodes) != 4) {
528 delete[] vertexCache;
529 return nullptr;
530 }
531 }
532 else {
533 if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, ¶metric,
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(®ionsSize, 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(¶metric, 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> ®ions,
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> ®ions,
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