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