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 "GmshConfig.h"
7 #include "GmshMessage.h"
8 #include "PViewDataGModel.h"
9 #include "MVertex.h"
10 #include "MElement.h"
11 #include "Numeric.h"
12 #include "StringUtils.h"
13 #include "OS.h"
14 #include "Context.h"
15 #include "fullMatrix.h"
16 #include "CGNSCommon.h"
17 #include "CGNSConventions.h"
18
addData(GModel * model,const std::map<int,std::vector<double>> & data,int step,double time,int partition,int numComp)19 bool PViewDataGModel::addData(GModel *model,
20 const std::map<int, std::vector<double> > &data,
21 int step, double time, int partition, int numComp)
22 {
23 if(data.empty()) return false;
24 if(step < 0) return false;
25
26 if(numComp < 0) {
27 numComp = 9;
28 for(auto it = data.begin(); it != data.end(); it++)
29 numComp = std::min(numComp, (int)it->second.size());
30 }
31
32 while(step >= (int)_steps.size())
33 _steps.push_back(new stepData<double>(model, numComp));
34 _steps[step]->fillEntities();
35 _steps[step]->computeBoundingBox();
36 _steps[step]->setTime(time);
37
38 int numEnt = (_type == NodeData) ? model->getNumMeshVertices() :
39 model->getNumMeshElements();
40 _steps[step]->resizeData(numEnt);
41
42 for(auto it = data.begin(); it != data.end(); it++) {
43 int mult = it->second.size() / numComp;
44 double *d = _steps[step]->getData(it->first, true, mult);
45 for(int j = 0; j < numComp * mult; j++) d[j] = it->second[j];
46 }
47 if(partition >= 0) _steps[step]->getPartitions().insert(partition);
48 finalize();
49 return true;
50 }
51
addData(GModel * model,const std::vector<std::size_t> & tags,const std::vector<std::vector<double>> & data,int step,double time,int partition,int numComp)52 bool PViewDataGModel::addData(GModel *model,
53 const std::vector<std::size_t> &tags,
54 const std::vector<std::vector<double> > &data,
55 int step, double time, int partition, int numComp)
56 {
57 if(data.empty() || tags.empty() || data.size() != tags.size()) return false;
58
59 if(numComp < 0) {
60 if(_type == ElementNodeData) {
61 numComp = 1; // cannot infer, as we can have different element types
62 }
63 else {
64 numComp = 9;
65 for(std::size_t i = 0; i < data.size(); i++)
66 numComp = std::min(numComp, (int)data[i].size());
67 }
68 }
69
70 while(step >= (int)_steps.size())
71 _steps.push_back(new stepData<double>(model, numComp));
72 _steps[step]->fillEntities();
73 _steps[step]->computeBoundingBox();
74 _steps[step]->setTime(time);
75
76 int numEnt = (_type == NodeData) ? model->getNumMeshVertices() :
77 model->getNumMeshElements();
78 _steps[step]->resizeData(numEnt);
79
80 for(std::size_t i = 0; i < data.size(); i++) {
81 int mult = data[i].size() / numComp;
82 double *d = _steps[step]->getData(tags[i], true, mult);
83 for(int j = 0; j < numComp * mult; j++) d[j] = data[i][j];
84 }
85 if(partition >= 0) _steps[step]->getPartitions().insert(partition);
86 finalize();
87 return true;
88 }
89
addData(GModel * model,const std::vector<std::size_t> & tags,const std::vector<double> & data,int step,double time,int partition,int numComp)90 bool PViewDataGModel::addData(GModel *model,
91 const std::vector<std::size_t> &tags,
92 const std::vector<double> &data, int step,
93 double time, int partition, int numComp)
94 {
95 if(data.empty() || tags.empty()) return false;
96
97 std::size_t stride = data.size() / tags.size();
98 if(stride < 1) return false;
99
100 if(numComp < 0) {
101 if(_type == ElementNodeData) {
102 numComp = 1; // cannot infer, as we can have different element types
103 }
104 else {
105 numComp = (int)stride;
106 }
107 }
108
109 while(step >= (int)_steps.size())
110 _steps.push_back(new stepData<double>(model, numComp));
111 _steps[step]->fillEntities();
112 _steps[step]->computeBoundingBox();
113 _steps[step]->setTime(time);
114
115 int numEnt = (_type == NodeData) ? model->getNumMeshVertices() :
116 model->getNumMeshElements();
117 _steps[step]->resizeData(numEnt);
118
119 int mult = stride / numComp;
120 for(std::size_t i = 0; i < tags.size(); i++) {
121 double *d = _steps[step]->getData(tags[i], true, mult);
122 int k = i * stride;
123 for(std::size_t j = 0; j < stride; j++) d[j] = data[k + j];
124 }
125 if(partition >= 0) _steps[step]->getPartitions().insert(partition);
126 finalize();
127 return true;
128 }
129
destroyData()130 void PViewDataGModel::destroyData()
131 {
132 for(std::size_t i = 0; i < _steps.size(); i++) _steps[i]->destroyData();
133 }
134
readMSH(const std::string & viewName,const std::string & fileName,int fileIndex,FILE * fp,bool binary,bool swap,int step,double time,int partition,int numComp,int numEnt,const std::string & interpolationScheme)135 bool PViewDataGModel::readMSH(const std::string &viewName,
136 const std::string &fileName, int fileIndex,
137 FILE *fp, bool binary, bool swap, int step,
138 double time, int partition, int numComp,
139 int numEnt,
140 const std::string &interpolationScheme)
141 {
142 Msg::Debug("Reading view `%s' step %d (time %g) partition %d: %d records",
143 viewName.c_str(), step, time, partition, numEnt);
144
145 while(step >= (int)_steps.size())
146 _steps.push_back(new stepData<double>(GModel::current(), numComp));
147 _steps[step]->fillEntities();
148 _steps[step]->computeBoundingBox();
149 _steps[step]->setFileName(fileName);
150 _steps[step]->setFileIndex(fileIndex);
151 _steps[step]->setTime(time);
152
153 /*
154 // if we already have maxSteps for this view, return
155 int numSteps = 0, maxSteps = 1000000000;
156 for(std::size_t i = 0; i < _steps.size(); i++)
157 numSteps += _steps[i]->getNumData() ? 1 : 0;
158 if(numSteps > maxSteps) return true;
159 */
160
161 _steps[step]->resizeData(numEnt);
162
163 Msg::StartProgressMeter(numEnt);
164 for(int i = 0; i < numEnt; i++) {
165 int num;
166 if(binary) {
167 if(fread(&num, sizeof(int), 1, fp) != 1) return false;
168 if(swap) SwapBytes((char *)&num, sizeof(int), 1);
169 }
170 else {
171 if(fscanf(fp, "%d", &num) != 1) return false;
172 }
173 if(num < 0) return false;
174 int mult = 1;
175 if(_type == ElementNodeData || _type == GaussPointData) {
176 if(binary) {
177 if(fread(&mult, sizeof(int), 1, fp) != 1) return false;
178 if(swap) SwapBytes((char *)&mult, sizeof(int), 1);
179 }
180 else {
181 if(fscanf(fp, "%d", &mult) != 1) return false;
182 }
183 }
184 double *d = _steps[step]->getData(num, true, mult);
185 if(binary) {
186 if((int)fread(d, sizeof(double), numComp * mult, fp) != numComp * mult)
187 return false;
188 if(swap) SwapBytes((char *)d, sizeof(double), numComp * mult);
189 }
190 else {
191 for(int j = 0; j < numComp * mult; j++)
192 if(fscanf(fp, "%lf", &d[j]) != 1) return false;
193 }
194 // compute min/max here to avoid calling finalize(true) later:
195 // this would be very slow for large multi-step, multi-partition
196 // datasets (since we would recompute the min/max for all the
197 // previously loaded steps/partitions, and thus loop over all the
198 // elements many times)
199 for(int j = 0; j < mult; j++) {
200 double val = ComputeScalarRep(numComp, &d[numComp * j]);
201 _steps[step]->setMin(std::min(_steps[step]->getMin(), val));
202 _steps[step]->setMax(std::max(_steps[step]->getMax(), val));
203 _min = std::min(_min, val);
204 _max = std::max(_max, val);
205 }
206 if(numEnt > 100000) Msg::ProgressMeter(i + 1, true, "Reading data");
207 }
208 Msg::StopProgressMeter();
209 if(partition >= 0) _steps[step]->getPartitions().insert(partition);
210
211 finalize(false, interpolationScheme);
212 return true;
213 }
214
writeMSH(const std::string & fileName,double version,bool binary,bool saveMesh,bool multipleView,int partitionNum,bool saveInterpolationMatrices,bool forceNodeData,bool forceElementData)215 bool PViewDataGModel::writeMSH(const std::string &fileName, double version,
216 bool binary, bool saveMesh, bool multipleView,
217 int partitionNum, bool saveInterpolationMatrices,
218 bool forceNodeData, bool forceElementData)
219 {
220 if(_steps.empty()) return true;
221
222 if(hasMultipleMeshes()) {
223 Msg::Info("Exporting multi-mesh view in separate files");
224 }
225
226 if(forceNodeData && _type != NodeData) {
227 Msg::Warning("Cannot force NodeData for this dataset: saving native data");
228 }
229
230 if(forceElementData && _type != ElementData) {
231 Msg::Warning(
232 "Cannot force ElementData for this dataset: saving native data");
233 }
234
235 FILE *fp = nullptr;
236 GModel *model0 = _steps[0]->getModel();
237 int numFile = 0;
238
239 for(std::size_t step = 0; step < _steps.size(); step++) {
240 int numEnt = 0, numComp = _steps[step]->getNumComponents();
241 for(std::size_t i = 0; i < _steps[step]->getNumData(); i++)
242 if(_steps[step]->getData(i)) numEnt++;
243 if(!numEnt) continue; // skip step
244
245 // open file, save mesh and save interpolation matrices
246 if(!fp || _steps[step]->getModel() != model0) {
247 if(fp) fclose(fp);
248 std::string stepFileName = fileName;
249 if(hasMultipleMeshes()) {
250 std::ostringstream sstream;
251 std::vector<std::string> n = SplitFileName(fileName);
252 sstream << n[0] << n[1] << "_" << numFile++ << n[2];
253 stepFileName = sstream.str();
254 model0 = _steps[step]->getModel();
255 }
256 if(saveMesh) {
257 if(!_steps[step]->getModel()->writeMSH(stepFileName, version, binary,
258 false, false, 1.0, 0, 0,
259 multipleView))
260 return false;
261 // append data
262 fp = Fopen(stepFileName.c_str(), binary ? "ab" : "a");
263 if(!fp) {
264 Msg::Error("Unable to open file '%s'", stepFileName.c_str());
265 return false;
266 }
267 }
268 else {
269 if(multipleView) {
270 fp = Fopen(stepFileName.c_str(), binary ? "ab" : "a");
271 if(!fp) {
272 Msg::Error("Unable to open file '%s'", stepFileName.c_str());
273 return false;
274 }
275 }
276 else {
277 fp = Fopen(stepFileName.c_str(), binary ? "wb" : "w");
278 if(!fp) {
279 Msg::Error("Unable to open file '%s'", stepFileName.c_str());
280 return false;
281 }
282 fprintf(fp, "$MeshFormat\n");
283 fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0,
284 (int)sizeof(double));
285 if(binary) {
286 int one = 1;
287 fwrite(&one, sizeof(int), 1, fp);
288 fprintf(fp, "\n");
289 }
290 fprintf(fp, "$EndMeshFormat\n");
291 }
292 }
293 // save the interpolation matrix?
294 if(saveInterpolationMatrices && haveInterpolationMatrices()) {
295 fprintf(fp, "$InterpolationScheme\n");
296 fprintf(fp, "\"INTERPOLATION_SCHEME\"\n");
297 fprintf(fp, "%d\n", (int)_interpolation.size());
298 for(auto it = _interpolation.begin(); it != _interpolation.end();
299 it++) {
300 if(it->second.size() >= 2) {
301 fprintf(fp, "%d\n2\n", it->first);
302 for(int mat = 0; mat < 2; mat++) {
303 int m = it->second[mat]->size1(), n = it->second[mat]->size2();
304 fprintf(fp, "%d %d\n", m, n);
305 for(int i = 0; i < m; i++) {
306 for(int j = 0; j < n; j++)
307 fprintf(fp, "%.16g ", it->second[mat]->get(i, j));
308 fprintf(fp, "\n");
309 }
310 }
311 }
312 }
313 fprintf(fp, "$EndInterpolationScheme\n");
314 }
315 }
316
317 // save the data
318 if(_type == NodeData) {
319 fprintf(fp, "$NodeData\n");
320 fprintf(fp, "1\n\"%s\"\n", getName().c_str());
321 fprintf(fp, "1\n%.16g\n", _steps[step]->getTime());
322 if(partitionNum > 0) {
323 fprintf(fp, "4\n%lu\n%d\n%d\n%d\n", step, numComp, numEnt,
324 partitionNum);
325 }
326 else if(_steps[step]->getPartitions().size() == 1) {
327 int p = *_steps[step]->getPartitions().begin();
328 fprintf(fp, "4\n%lu\n%d\n%d\n%d\n", step, numComp, numEnt, p);
329 }
330 else {
331 fprintf(fp, "3\n%lu\n%d\n%d\n", step, numComp, numEnt);
332 }
333 for(std::size_t i = 0; i < _steps[step]->getNumData(); i++) {
334 if(_steps[step]->getData(i)) {
335 MVertex *v = _steps[step]->getModel()->getMeshVertexByTag(i);
336 if(!v) {
337 Msg::Error("Unknown node %d in data", i);
338 fclose(fp);
339 return false;
340 }
341 int num = (version >= 3.0) ? v->getNum() : v->getIndex();
342 if(binary) {
343 fwrite(&num, sizeof(int), 1, fp);
344 fwrite(_steps[step]->getData(i), sizeof(double), numComp, fp);
345 }
346 else {
347 fprintf(fp, "%d", num);
348 for(int k = 0; k < numComp; k++)
349 fprintf(fp, " %.16g", _steps[step]->getData(i)[k]);
350 fprintf(fp, "\n");
351 }
352 }
353 }
354 if(binary) fprintf(fp, "\n");
355 fprintf(fp, "$EndNodeData\n");
356 }
357 else {
358 if(_type == ElementNodeData)
359 fprintf(fp, "$ElementNodeData\n");
360 else
361 fprintf(fp, "$ElementData\n");
362 if(saveInterpolationMatrices && haveInterpolationMatrices())
363 fprintf(fp, "2\n\"%s\"\n\"INTERPOLATION_SCHEME\"\n", getName().c_str());
364 else
365 fprintf(fp, "1\n\"%s\"\n", getName().c_str());
366
367 fprintf(fp, "1\n%.16g\n", _steps[step]->getTime());
368 if(partitionNum > 0) {
369 fprintf(fp, "4\n%lu\n%d\n%d\n%d\n", step, numComp, numEnt,
370 partitionNum);
371 }
372 else if(_steps[step]->getPartitions().size() == 1) {
373 int p = *_steps[step]->getPartitions().begin();
374 fprintf(fp, "4\n%lu\n%d\n%d\n%d\n", step, numComp, numEnt, p);
375 }
376 else {
377 fprintf(fp, "3\n%lu\n%d\n%d\n", step, numComp, numEnt);
378 }
379 for(std::size_t i = 0; i < _steps[step]->getNumData(); i++) {
380 if(_steps[step]->getData(i)) {
381 MElement *e = _steps[step]->getModel()->getMeshElementByTag(i);
382 if(!e) {
383 Msg::Error("Unknown element %d in data", i);
384 fclose(fp);
385 return false;
386 }
387 int mult = _steps[step]->getMult(i);
388 int num = (version >= 3.0) ?
389 e->getNum() :
390 _steps[step]->getModel()->getMeshElementIndex(e);
391 if(binary) {
392 fwrite(&num, sizeof(int), 1, fp);
393 if(_type == ElementNodeData) fwrite(&mult, sizeof(int), 1, fp);
394 fwrite(_steps[step]->getData(i), sizeof(double), numComp * mult,
395 fp);
396 }
397 else {
398 fprintf(fp, "%d", num);
399 if(_type == ElementNodeData) fprintf(fp, " %d", mult);
400 for(int k = 0; k < numComp * mult; k++)
401 fprintf(fp, " %.16g", _steps[step]->getData(i)[k]);
402 fprintf(fp, "\n");
403 }
404 }
405 }
406 if(binary) fprintf(fp, "\n");
407 if(_type == ElementNodeData)
408 fprintf(fp, "$EndElementNodeData\n");
409 else
410 fprintf(fp, "$EndElementData\n");
411 }
412 }
413
414 fclose(fp);
415 return true;
416 }
417
importLists(int N[24],std::vector<double> * V[24])418 void PViewDataGModel::importLists(int N[24], std::vector<double> *V[24])
419 {
420 for(int idxtype = 0; idxtype < 24; idxtype++) {
421 int nbe = N[idxtype];
422 if(!nbe) continue;
423 std::vector<double> *list = V[idxtype];
424 int nc = 0, nn = 0;
425 switch(idxtype) {
426 case 0:
427 nc = 1;
428 nn = 1;
429 break; // SP
430 case 1:
431 nc = 3;
432 nn = 1;
433 break; // VP
434 case 2:
435 nc = 9;
436 nn = 1;
437 break; // TP
438 case 3:
439 nc = 1;
440 nn = 2;
441 break; // SL
442 case 4:
443 nc = 3;
444 nn = 2;
445 break; // VL
446 case 5:
447 nc = 9;
448 nn = 2;
449 break; // TL
450 case 6:
451 nc = 1;
452 nn = 3;
453 break; // ST
454 case 7:
455 nc = 3;
456 nn = 3;
457 break; // VT
458 case 8:
459 nc = 9;
460 nn = 3;
461 break; // TT
462 case 9:
463 nc = 1;
464 nn = 4;
465 break; // SQ
466 case 10:
467 nc = 3;
468 nn = 4;
469 break; // VQ
470 case 11:
471 nc = 9;
472 nn = 4;
473 break; // TQ
474 case 12:
475 nc = 1;
476 nn = 4;
477 break; // SS
478 case 13:
479 nc = 3;
480 nn = 4;
481 break; // VS
482 case 14:
483 nc = 9;
484 nn = 4;
485 break; // TS
486 case 15:
487 nc = 1;
488 nn = 8;
489 break; // SH
490 case 16:
491 nc = 3;
492 nn = 8;
493 break; // VH
494 case 17:
495 nc = 9;
496 nn = 8;
497 break; // TH
498 case 18:
499 nc = 1;
500 nn = 6;
501 break; // SI
502 case 19:
503 nc = 3;
504 nn = 6;
505 break; // VI
506 case 20:
507 nc = 9;
508 nn = 6;
509 break; // TI
510 case 21:
511 nc = 1;
512 nn = 5;
513 break; // SY
514 case 22:
515 nc = 3;
516 nn = 5;
517 break; // VY
518 case 23:
519 nc = 9;
520 nn = 5;
521 break; // TY
522 }
523 int stride = list->size() / nbe;
524 int numSteps = (stride - 1) / nc / nn;
525 for(int step = 0; step < numSteps; step++) {
526 _steps.push_back(new stepData<double>(GModel::current(), nc));
527 _steps[step]->fillEntities();
528 _steps[step]->computeBoundingBox();
529 _steps[step]->setTime(step);
530 _steps[step]->resizeData(nbe);
531 for(std::size_t j = 0; j < list->size(); j += stride) {
532 double *tmp = &(*list)[j];
533 int num = (int)tmp[0];
534 double *d = _steps[step]->getData(num, true, nn);
535 for(int k = 0; k < nc * nn; k++) { d[k] = tmp[1 + nc * nn * step + k]; }
536 }
537 }
538 }
539
540 finalize();
541 }
542
543 #if defined(HAVE_MED)
544
545 #include <med.h>
546
547 #if(MED_MAJOR_NUM >= 3)
548 // To avoid too many ifdefs below we use defines for the bits of the
549 // API that did not change too much between MED2 and MED3. If we
550 // remove MED2 support at some point, please remove these defines and
551 // replace the symbols accordingly.
552 #define med_geometrie_element med_geometry_type
553 #define med_entite_maillage med_entity_type
554 #define med_type_champ med_field_type
555 #define MED_TAILLE_NOM MED_NAME_SIZE
556 #define MED_TAILLE_PNOM MED_SNAME_SIZE
557 #define MED_LECTURE MED_ACC_RDONLY
558 #define MED_LECTURE_AJOUT MED_ACC_RDEXT
559 #define MED_NOEUD MED_NODE
560 #define MED_MAILLE MED_CELL
561 #define MED_NOEUD_MAILLE MED_NODE_ELEMENT
562 #define MED_NOPFL MED_NO_PROFILE
563 #define MEDouvrir MEDfileOpen
564 #define MEDfermer MEDfileClose
565 #define MEDnChamp MEDfieldnComponent
566 #define MEDnValProfil MEDprofileSizeByName
567 #endif
568
569 extern int med2mshElementType(med_geometrie_element med);
570 extern int med2mshNodeIndex(med_geometrie_element med, int k);
571
medGetFieldNames(const std::string & fileName)572 std::vector<std::string> medGetFieldNames(const std::string &fileName)
573 {
574 std::vector<std::string> fieldNames;
575
576 #if(MED_MAJOR_NUM >= 3)
577 med_idt fid = MEDfileOpen(fileName.c_str(), MED_ACC_RDONLY);
578 #else
579 med_idt fid = MEDouvrir((char *)fileName.c_str(), MED_LECTURE);
580 #endif
581 if(fid < 0) {
582 Msg::Error("Unable to open file '%s'", fileName.c_str());
583 return fieldNames;
584 }
585
586 #if(MED_MAJOR_NUM >= 3)
587 med_int numFields = MEDnField(fid);
588 #else
589 med_int numFields = MEDnChamp(fid, 0);
590 #endif
591
592 for(int index = 0; index < numFields; index++) {
593 med_int numComp = MEDnChamp(fid, index + 1);
594 if(numComp <= 0) {
595 Msg::Error("Could not get number of components for MED field");
596 return fieldNames;
597 }
598
599 char name[MED_TAILLE_NOM + 1], meshName[MED_TAILLE_NOM + 1];
600 char dtUnit[MED_TAILLE_PNOM + 1];
601 std::vector<char> compName(numComp * MED_TAILLE_PNOM + 1);
602 std::vector<char> compUnit(numComp * MED_TAILLE_PNOM + 1);
603 med_int numSteps = 0;
604 med_type_champ type;
605 #if(MED_MAJOR_NUM >= 3)
606 med_bool localMesh;
607 if(MEDfieldInfo(fid, index + 1, name, meshName, &localMesh, &type,
608 &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
609 #else
610 if(MEDchampInfo(fid, index + 1, name, &type, &compName[0], &compUnit[0],
611 numComp) < 0) {
612 #endif
613 Msg::Error("Could not get MED field info");
614 return fieldNames;
615 }
616 fieldNames.push_back(name);
617 }
618
619 #if(MED_MAJOR_NUM >= 3)
620 if(MEDfileClose(fid) < 0) {
621 #else
622 if(MEDfermer(fid) < 0) {
623 #endif
624 Msg::Error("Unable to close file '%s'", fileName.c_str());
625 }
626 return fieldNames;
627 }
628
629 bool PViewDataGModel::readMED(const std::string &fileName, int fileIndex)
630 {
631 med_idt fid = MEDouvrir((char *)fileName.c_str(), MED_LECTURE);
632 if(fid < 0) {
633 Msg::Error("Unable to open file '%s'", fileName.c_str());
634 return false;
635 }
636
637 med_int numComp = MEDnChamp(fid, fileIndex + 1);
638 if(numComp <= 0) {
639 Msg::Error("Could not get number of components for MED field");
640 return false;
641 }
642
643 char name[MED_TAILLE_NOM + 1], meshName[MED_TAILLE_NOM + 1];
644 char dtUnit[MED_TAILLE_PNOM + 1];
645 std::vector<char> compName(numComp * MED_TAILLE_PNOM + 1);
646 std::vector<char> compUnit(numComp * MED_TAILLE_PNOM + 1);
647 med_int numSteps = 0;
648 med_type_champ type;
649 #if(MED_MAJOR_NUM >= 3)
650 med_bool localMesh;
651 if(MEDfieldInfo(fid, fileIndex + 1, name, meshName, &localMesh, &type,
652 &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
653 #else
654 if(MEDchampInfo(fid, fileIndex + 1, name, &type, &compName[0], &compUnit[0],
655 numComp) < 0) {
656 #endif
657 Msg::Error("Could not get MED field info");
658 return false;
659 }
660
661 Msg::Info("Reading %d-component field '%s'", numComp, name);
662 setName(name);
663 setFileName(fileName);
664 setFileIndex(fileIndex);
665
666 int numCompMsh = (numComp <= 1) ? 1 :
667 (numComp <= 3) ? 3 :
668 (numComp <= 9) ? 9 :
669 numComp;
670
671 if(numCompMsh > 9)
672 Msg::Info(
673 "More than 9 components in field: you will probably want to force "
674 "the field type to scalar, vector or tensor in the options");
675
676 // Warning! The ordering of the elements in the last two lists is
677 // important: it should match the ordering of the MSH element types
678 // (when elements are saved without tags, this governs the order
679 // with which we implicitly index them in GModel::readMED)
680 const med_entite_maillage entType[] = {MED_NOEUD, MED_MAILLE,
681 MED_NOEUD_MAILLE};
682 #if(MED_MAJOR_NUM >= 3)
683 const med_geometrie_element eleType[] = {
684 MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
685 MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_QUAD9, MED_TETRA10,
686 MED_HEXA27, MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
687 const int nodesPerEle[] = {0, 2, 3, 4, 4, 8, 6, 5, 3,
688 6, 9, 10, 27, 1, 8, 20, 15, 13};
689 #else
690 const med_geometrie_element eleType[] = {
691 MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
692 MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_TETRA10, MED_POINT1,
693 MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
694 const int nodesPerEle[] = {0, 2, 3, 4, 4, 8, 6, 5,
695 3, 6, 10, 1, 8, 20, 15, 13};
696 #endif
697
698 std::vector<std::pair<int, int> > pairs;
699 for(std::size_t i = 0; i < sizeof(entType) / sizeof(entType[0]); i++) {
700 for(std::size_t j = 0; j < sizeof(eleType) / sizeof(eleType[0]); j++) {
701 if((!i && !j) || j) {
702 #if(MED_MAJOR_NUM >= 3)
703 med_int n = numSteps;
704 #else
705 med_int n = MEDnPasdetemps(fid, name, entType[i], eleType[j]);
706 #endif
707 if(n > 0) {
708 pairs.push_back(std::make_pair(i, j));
709 numSteps = std::max(numSteps, n);
710 }
711 if(!i && !j) break; // MED_NOEUD does not care about eleType
712 }
713 }
714 }
715
716 if(numSteps < 1 || pairs.empty()) {
717 Msg::Error("Nothing to import from MED file");
718 return false;
719 }
720
721 for(int step = 0; step < numSteps; step++) {
722 // FIXME: in MED3 we might want to loop over all profiles instead
723 // of relying of the default one
724
725 // FIXME: MED3 allows to store multi-step meshes; we should
726 // interface this with our own gmodel-per-step structure
727
728 for(std::size_t pair = 0; pair < pairs.size(); pair++) {
729 // get step info
730 med_entite_maillage ent = entType[pairs[pair].first];
731 med_geometrie_element ele = eleType[pairs[pair].second];
732 med_int numdt, numit, ngauss;
733 med_float dt;
734 #if(MED_MAJOR_NUM >= 3)
735 if(MEDfieldComputingStepInfo(fid, name, step + 1, &numdt, &numit, &dt) <
736 0) {
737 #else
738 char dtunit[MED_TAILLE_PNOM + 1];
739 med_booleen local;
740 med_int numMeshes;
741 if(MEDpasdetempsInfo(fid, name, ent, ele, step + 1, &ngauss, &numdt,
742 &numit, dtunit, &dt, meshName, &local,
743 &numMeshes) < 0) {
744 #endif
745 Msg::Error("Could not read step info");
746 return false;
747 }
748 // create step data
749 if(!pair) {
750 GModel *m = GModel::findByName(meshName);
751 if(!m) {
752 Msg::Error("Could not find mesh '%s'", meshName);
753 return false;
754 }
755 while(step >= (int)_steps.size())
756 _steps.push_back(new stepData<double>(m, numCompMsh));
757 _steps[step]->fillEntities();
758 _steps[step]->computeBoundingBox();
759 _steps[step]->setFileName(fileName);
760 _steps[step]->setFileIndex(fileIndex);
761 _steps[step]->setTime(dt);
762 }
763
764 char locName[MED_TAILLE_NOM + 1], profileName[MED_TAILLE_NOM + 1];
765
766 // get number of values in the field (numVal takes the number of
767 // Gauss points or the number of nodes per element into account,
768 // but not the number of components)
769 #if(MED_MAJOR_NUM >= 3)
770 med_int profileSize;
771 med_int numVal = MEDfieldnValueWithProfile(
772 fid, name, numdt, numit, ent, ele, 1, MED_COMPACT_STMODE, profileName,
773 &profileSize, locName, &ngauss);
774 numVal *= ngauss;
775 #else
776 med_int numVal =
777 MEDnVal(fid, name, ent, ele, numdt, numit, meshName, MED_COMPACT);
778 #endif
779 if(numVal <= 0) continue;
780
781 _type = (ent == MED_NOEUD) ? NodeData :
782 (ent == MED_MAILLE) ? ElementData :
783 ElementNodeData;
784 int mult = 1;
785 if(ent == MED_NOEUD_MAILLE) { mult = nodesPerEle[pairs[pair].second]; }
786 else if(ngauss != 1) {
787 mult = ngauss;
788 _type = GaussPointData;
789 }
790 _steps[step]->resizeData(numVal / mult);
791
792 // read field data
793 std::vector<double> val(numVal * numComp);
794 #if(MED_MAJOR_NUM >= 3)
795 if(MEDfieldValueWithProfileRd(fid, name, numdt, numit, ent, ele,
796 MED_COMPACT_STMODE, profileName,
797 MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
798 (unsigned char *)&val[0]) < 0) {
799 #else
800 if(MEDchampLire(fid, meshName, name, (unsigned char *)&val[0],
801 MED_FULL_INTERLACE, MED_ALL, locName, profileName,
802 MED_COMPACT, ent, ele, numdt, numit) < 0) {
803 #endif
804 Msg::Error("Could not read field values");
805 return false;
806 }
807
808 Msg::Debug(
809 "MED: eletyp=%d entity=%d (0:cell, 3:node, 4:elenode) ngauss=%d "
810 "localizationName=%s profileName=%s -- stepDataType=%d",
811 ele, ent, ngauss, locName, profileName, _type);
812
813 // read Gauss point data
814 if(_type == GaussPointData) {
815 std::vector<double> &p(
816 _steps[step]->getGaussPoints(med2mshElementType(ele)));
817 if(std::string(locName) == MED_GAUSS_ELNO) {
818 // special case: the gauss points are the vertices of the
819 // element; in this case no explicit localization has to be
820 // created in MED
821 p.resize(ngauss * 3, 1.e22);
822 }
823 else {
824 int dim = ele / 100;
825 std::vector<med_float> refcoo((ele % 100) * dim);
826 std::vector<med_float> gscoo(ngauss * dim);
827 std::vector<med_float> wg(ngauss);
828 #if(MED_MAJOR_NUM >= 3)
829 if(MEDlocalizationRd(fid, locName, MED_FULL_INTERLACE, &refcoo[0],
830 &gscoo[0], &wg[0]) < 0) {
831 #else
832 if(MEDgaussLire(fid, &refcoo[0], &gscoo[0], &wg[0],
833 MED_FULL_INTERLACE, locName) < 0) {
834 #endif
835 Msg::Error("Could not read Gauss points");
836 return false;
837 }
838 // FIXME: we should check that refcoo corresponds to our
839 // internal reference element
840 for(int i = 0; i < (int)gscoo.size(); i++) {
841 p.push_back(gscoo[i]);
842 if(i % dim == dim - 1)
843 for(int j = 0; j < 3 - dim; j++) p.push_back(0.);
844 }
845 }
846 }
847
848 // compute profile (indices in full array of entities of given type)
849 std::vector<med_int> profile;
850 if(std::string(profileName) != MED_NOPFL) {
851 med_int n = MEDnValProfil(fid, profileName);
852 if(n > 0) {
853 Msg::Debug("MED has full profile");
854 profile.resize(n);
855 #if(MED_MAJOR_NUM >= 3)
856 if(MEDprofileRd(fid, profileName, &profile[0]) < 0) {
857 #else
858 if(MEDprofilLire(fid, &profile[0], profileName) < 0) {
859 #endif
860 Msg::Error("Could not read profile");
861 return false;
862 }
863 }
864 }
865 if(profile.empty()) {
866 Msg::Debug("MED profile is empty -- using continuous sequence");
867 profile.resize(numVal / mult);
868 for(std::size_t i = 0; i < profile.size(); i++) profile[i] = i + 1;
869 }
870
871 // get size of full array and tags (if any) of entities
872 bool nodal = (ent == MED_NOEUD);
873 #if(MED_MAJOR_NUM >= 3)
874 med_bool changeOfCoord;
875 med_bool geoTransform;
876 med_int numEnt = MEDmeshnEntity(
877 fid, meshName, MED_NO_DT, MED_NO_IT, nodal ? MED_NODE : MED_CELL,
878 nodal ? MED_NO_GEOTYPE : ele, nodal ? MED_COORDINATE : MED_CONNECTIVITY,
879 nodal ? MED_NO_CMODE : MED_NODAL, &changeOfCoord, &geoTransform);
880 #else
881 med_int numEnt =
882 MEDnEntMaa(fid, meshName, nodal ? MED_COOR : MED_CONN,
883 nodal ? MED_NOEUD : MED_MAILLE, nodal ? MED_NONE : ele,
884 nodal ? (med_connectivite)0 : MED_NOD);
885 #endif
886 std::vector<med_int> tags(numEnt);
887 #if(MED_MAJOR_NUM >= 3)
888 if(MEDmeshEntityNumberRd(fid, meshName, MED_NO_DT, MED_NO_IT,
889 nodal ? MED_NODE : MED_CELL,
890 nodal ? MED_NO_GEOTYPE : ele, &tags[0]) < 0)
891 #else
892 if(MEDnumLire(fid, meshName, &tags[0], numEnt,
893 nodal ? MED_NOEUD : MED_MAILLE, nodal ? MED_NONE : ele) < 0)
894 #endif
895 tags.clear();
896
897 // if we don't have tags, compute the starting index (i.e., how many
898 // elements of different type are in the mesh before these ones)
899 std::size_t startIndex = 0;
900 if(tags.empty()) {
901 std::size_t maxv, maxe;
902 _steps[step]->getModel()->getCheckPointedMaxNumbers(maxv, maxe);
903 if(nodal) { startIndex += maxv; }
904 else {
905 for(int i = 1; i < pairs[pair].second; i++) {
906 #if(MED_MAJOR_NUM >= 3)
907 med_int n = MEDmeshnEntity(
908 fid, meshName, MED_NO_DT, MED_NO_IT, MED_CELL, eleType[i],
909 MED_CONNECTIVITY, MED_NODAL, &changeOfCoord, &geoTransform);
910 #else
911 med_int n = MEDnEntMaa(fid, meshName, MED_CONN, MED_MAILLE,
912 eleType[i], MED_NOD);
913 #endif
914 if(n > 0) startIndex += n;
915 }
916 startIndex += maxe;
917 }
918 Msg::Debug("MED has no tags -- assuming starting index %lu",
919 startIndex);
920 }
921
922 // compute entity numbers using profile, then fill step data
923 for(std::size_t i = 0; i < profile.size(); i++) {
924 std::size_t num;
925 if(tags.empty()) { num = startIndex + profile[i]; }
926 else {
927 if(profile[i] == 0 || profile[i] > (int)tags.size()) {
928 Msg::Error("Wrong index in profile");
929 return false;
930 }
931 num = tags[profile[i] - 1];
932 }
933
934 double *d = _steps[step]->getData(num, true, mult);
935 for(int j = 0; j < mult; j++) {
936 // reorder nodes if we have ElementNode data
937 int j2 = (ent == MED_NOEUD_MAILLE) ? med2mshNodeIndex(ele, j) : j;
938 for(int k = 0; k < numComp; k++)
939 d[numCompMsh * j + k] = val[numComp * mult * i + numComp * j2 + k];
940 }
941 }
942 }
943 }
944
945 finalize();
946
947 if(MEDfermer(fid) < 0) {
948 Msg::Error("Unable to close file '%s'", (char *)fileName.c_str());
949 return false;
950 }
951 return true;
952 }
953
954 bool PViewDataGModel::writeMED(const std::string &fileName)
955 {
956 if(_steps.empty()) return true;
957
958 if(hasMultipleMeshes()) {
959 Msg::Error("Export not done for multi-mesh views");
960 return false;
961 }
962
963 if(_type != NodeData) {
964 Msg::Error("Can only export node-based datasets for now");
965 return false;
966 }
967
968 GModel *model = _steps[0]->getModel();
969
970 // save the mesh
971 if(!model->writeMED(fileName, true)) return false;
972
973 std::string meshName(model->getName());
974 std::string fieldName(getName());
975
976 #if (MED_MAJOR_NUM >= 4) || ((MED_MAJOR_NUM >= 3) && (MED_MINOR_NUM >= 3))
977 // MEDfileVersionOpen actually appeared in MED 3.2.1
978 med_int major = MED_MAJOR_NUM, minor = MED_MINOR_NUM,
979 release = MED_RELEASE_NUM;
980 if(CTX::instance()->mesh.medFileMinorVersion >= 0) {
981 minor = (int)CTX::instance()->mesh.medFileMinorVersion;
982 Msg::Info("Forcing MED file version to %d.%d", major, minor);
983 }
984 med_idt fid = MEDfileVersionOpen((char *)fileName.c_str(), MED_ACC_RDEXT,
985 major, minor, release);
986 #else
987 med_idt fid = MEDouvrir((char *)fileName.c_str(), MED_LECTURE_AJOUT);
988 #endif
989
990 if(fid < 0) {
991 Msg::Error("Unable to open file '%s'", fileName.c_str());
992 return false;
993 }
994
995 // compute profile
996 char *profileName = (char *)"nodeProfile";
997 std::vector<med_int> profile, indices;
998 for(std::size_t i = 0; i < _steps[0]->getNumData(); i++) {
999 if(_steps[0]->getData(i)) {
1000 MVertex *v = _steps[0]->getModel()->getMeshVertexByTag(i);
1001 if(!v) {
1002 Msg::Error("Unknown node %d in data", i);
1003 return false;
1004 }
1005 profile.push_back(v->getIndex());
1006 indices.push_back(i);
1007 }
1008 }
1009
1010 if(profile.empty()) {
1011 Msg::Error("Nothing to save");
1012 return false;
1013 }
1014
1015 #if(MED_MAJOR_NUM >= 3)
1016 if(MEDprofileWr(fid, profileName, (med_int)profile.size(), &profile[0]) < 0) {
1017 #else
1018 if(MEDprofilEcr(fid, &profile[0], (med_int)profile.size(), profileName) < 0) {
1019 #endif
1020 Msg::Error("Could not create MED profile");
1021 return false;
1022 }
1023
1024 int numComp = _steps[0]->getNumComponents();
1025 #if(MED_MAJOR_NUM >= 3)
1026 if(MEDfieldCr(fid, (char *)fieldName.c_str(), MED_FLOAT64, (med_int)numComp,
1027 "unknown", "unknown", "unknown",
1028 (char *)meshName.c_str()) < 0) {
1029 #else
1030 if(MEDchampCr(fid, (char *)fieldName.c_str(), MED_FLOAT64, (char *)"unknown",
1031 (char *)"unknown", (med_int)numComp) < 0) {
1032 #endif
1033 Msg::Error("Could not create MED field");
1034 return false;
1035 }
1036
1037 #if(MED_MAJOR_NUM >= 3)
1038 med_bool changeOfCoord, geoTransform;
1039 med_int numNodes =
1040 MEDmeshnEntity(fid, (char *)meshName.c_str(), MED_NO_DT, MED_NO_IT,
1041 MED_NODE, MED_NO_GEOTYPE, MED_COORDINATE, MED_NO_CMODE,
1042 &changeOfCoord, &geoTransform);
1043 #else
1044 med_int numNodes = MEDnEntMaa(fid, (char *)meshName.c_str(), MED_COOR,
1045 MED_NOEUD, MED_NONE, (med_connectivite)0);
1046 #endif
1047 if(numNodes <= 0) {
1048 Msg::Error("Could not get valid number of nodes in mesh");
1049 return false;
1050 }
1051 for(std::size_t step = 0; step < _steps.size(); step++) {
1052 std::size_t n = 0;
1053 for(std::size_t i = 0; i < _steps[step]->getNumData(); i++)
1054 if(_steps[step]->getData(i)) n++;
1055 if(n != profile.size() || numComp != _steps[step]->getNumComponents()) {
1056 Msg::Error("Skipping incompatible step");
1057 continue;
1058 }
1059 double time = _steps[step]->getTime();
1060 std::vector<double> val(profile.size() * numComp);
1061 for(std::size_t i = 0; i < profile.size(); i++)
1062 for(int k = 0; k < numComp; k++)
1063 val[i * numComp + k] = _steps[step]->getData(indices[i])[k];
1064 #if(MED_MAJOR_NUM >= 3)
1065 if(MEDfieldValueWithProfileWr(
1066 fid, (char *)fieldName.c_str(), (med_int)(step + 1), MED_NO_IT, time,
1067 MED_NODE, MED_NO_GEOTYPE, MED_COMPACT_STMODE, profileName, "",
1068 MED_FULL_INTERLACE, MED_ALL_CONSTITUENT, numNodes,
1069 (unsigned char *)&val[0]) < 0) {
1070 #else
1071 if(MEDchampEcr(fid, (char *)meshName.c_str(), (char *)fieldName.c_str(),
1072 (unsigned char *)&val[0], MED_FULL_INTERLACE, numNodes,
1073 (char *)MED_NOGAUSS, MED_ALL, profileName, MED_COMPACT,
1074 MED_NOEUD, MED_NONE, (med_int)step, (char *)"unknown", time,
1075 MED_NONOR) < 0) {
1076 #endif
1077 Msg::Error("Could not write MED field");
1078 return false;
1079 }
1080 }
1081
1082 if(MEDfermer(fid) < 0) {
1083 Msg::Error("Unable to close file '%s'", (char *)fileName.c_str());
1084 return false;
1085 }
1086 return true;
1087 }
1088
1089 #else
1090
1091 bool PViewDataGModel::readMED(const std::string &fileName, int fileIndex)
1092 {
1093 Msg::Error("Gmsh must be compiled with MED support to read '%s'",
1094 fileName.c_str());
1095 return false;
1096 }
1097
1098 bool PViewDataGModel::writeMED(const std::string &fileName)
1099 {
1100 Msg::Error("Gmsh must be compiled with MED support to write '%s'",
1101 fileName.c_str());
1102 return false;
1103 }
1104
1105 #endif
1106
1107 bool PViewDataGModel::readPCH(const std::string &fileName, int fileIndex)
1108 {
1109 Msg::Info("Placeholder for reading punch file '%s'", fileName.c_str());
1110
1111 std::map<int, std::vector<double> > data;
1112 for(int i = 1; i < 200; i++) data[i].push_back(1.234);
1113 addData(GModel::current(), data, 0, 0.0, 1, 1);
1114
1115 return true;
1116 }
1117
1118 void PViewDataGModel::sendToServer(const std::string &name)
1119 {
1120 if(_steps.empty()) return;
1121
1122 if(_type != NodeData) {
1123 Msg::Error("sendToServer currently only implemented for nodal datasets");
1124 return;
1125 }
1126
1127 int numEnt = 0, numComp = 0;
1128 for(std::size_t step = 0; step < _steps.size(); step++) {
1129 int nc = _steps[step]->getNumComponents();
1130 int ne = 0;
1131 for(std::size_t i = 0; i < _steps[step]->getNumData(); i++)
1132 if(_steps[step]->getData(i)) ne++;
1133 if(!step) {
1134 numEnt = ne;
1135 numComp = nc;
1136 }
1137 else {
1138 if(ne != numEnt || nc != numComp) {
1139 Msg::Error("Can not send heterogeneous view to server");
1140 return;
1141 }
1142 }
1143 }
1144
1145 std::vector<double> exp;
1146 exp.push_back(numEnt);
1147
1148 for(std::size_t i = 0; i < _steps[0]->getNumData(); i++) {
1149 if(_steps[0]->getData(i)) {
1150 MVertex *v = _steps[0]->getModel()->getMeshVertexByTag(i);
1151 if(!v) {
1152 Msg::Error("Unknown node %d in data", i);
1153 return;
1154 }
1155 int num = v->getNum();
1156 exp.push_back(num);
1157 for(std::size_t step = 0; step < _steps.size(); step++) {
1158 for(int k = 0; k < numComp; k++) {
1159 double data = _steps[step]->getData(i)[k];
1160 exp.push_back(data);
1161 }
1162 }
1163 }
1164 }
1165
1166 Msg::SetOnelabNumber(name, exp, false);
1167 }
1168