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