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 <string.h>
7 #include <set>
8 #include "PViewDataList.h"
9 #include "MElement.h"
10 #include "Numeric.h"
11 #include "StringUtils.h"
12 #include "GmshMessage.h"
13 #include "GmshDefines.h"
14 #include "MVertexRTree.h"
15 #include "Context.h"
16 #include "adaptiveData.h"
17 #include "OS.h"
18 
dVecRead(std::vector<double> & v,int n,FILE * fp,bool binary,int swap)19 static void dVecRead(std::vector<double> &v, int n, FILE *fp, bool binary,
20                      int swap)
21 {
22   if(n <= 0) return;
23   v.resize(n);
24   if(binary) {
25     if(!fread(&v[0], sizeof(double), n, fp)) Msg::Error("Read error");
26     if(swap) SwapBytes((char *)&v[0], sizeof(double), n);
27   }
28   else {
29     for(int i = 0; i < n; i++) {
30       if(fscanf(fp, "%lf", &v[i]) != 1) {
31         Msg::Error("Read error");
32         break;
33       }
34     }
35   }
36 }
37 
cVecRead(std::vector<char> & v,int n,FILE * fp,bool binary,int swap,bool oldStyle)38 static void cVecRead(std::vector<char> &v, int n, FILE *fp, bool binary,
39                      int swap, bool oldStyle)
40 {
41   if(n <= 0) return;
42   v.resize(n);
43   if(binary) {
44     if(!fread(&v[0], sizeof(char), n, fp)) Msg::Error("Read error");
45     if(swap) SwapBytes((char *)&v[0], sizeof(char), n);
46   }
47   else {
48     if(oldStyle) {
49       for(int i = 0; i < n; i++) {
50         if(fscanf(fp, "%c", &v[i]) != 1) {
51           Msg::Error("Read error");
52           break;
53         }
54         if(v[i] == '^') v[i] = '\0';
55       }
56     }
57     else {
58       for(int i = 0; i < n; i++) {
59         char c = (char)fgetc(fp);
60         if(c == EOF) {
61           Msg::Error("Read error");
62           break;
63         }
64         else {
65           v[i] = c;
66         }
67       }
68     }
69   }
70 }
71 
dVecWrite(std::vector<double> & v,FILE * fp,bool binary)72 static void dVecWrite(std::vector<double> &v, FILE *fp, bool binary)
73 {
74   if(v.empty()) return;
75   if(binary)
76     fwrite(&v[0], sizeof(double), v.size(), fp);
77   else
78     for(unsigned i = 0; i < v.size(); i++) fprintf(fp, " %.16g", v[i]);
79 }
80 
cVecWrite(std::vector<char> & v,FILE * fp,bool binary)81 static void cVecWrite(std::vector<char> &v, FILE *fp, bool binary)
82 {
83   if(v.empty()) return;
84   if(binary)
85     fwrite(&v[0], sizeof(char), v.size(), fp);
86   else
87     for(unsigned i = 0; i < v.size(); i++) fputc(v[i], fp);
88 }
89 
readPOS(FILE * fp,double version,bool binary)90 bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
91 {
92   char name[256];
93   int t2l, t3l;
94 
95   int NbSL2 = 0, NbVL2 = 0, NbTL2 = 0, NbST2 = 0, NbVT2 = 0, NbTT2 = 0;
96   int NbSQ2 = 0, NbVQ2 = 0, NbTQ2 = 0, NbSS2 = 0, NbVS2 = 0, NbTS2 = 0;
97   int NbSH2 = 0, NbVH2 = 0, NbTH2 = 0, NbSI2 = 0, NbVI2 = 0, NbTI2 = 0;
98   int NbSY2 = 0, NbVY2 = 0, NbTY2 = 0;
99   std::vector<double> SL2, VL2, TL2, ST2, VT2, TT2;
100   std::vector<double> SQ2, VQ2, TQ2, SS2, VS2, TS2;
101   std::vector<double> SH2, VH2, TH2, SI2, VI2, TI2;
102   std::vector<double> SY2, VY2, TY2;
103 
104   if(version <= 1.0) {
105     Msg::Debug("Detected post-processing view format <= 1.0");
106     if(fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d\n", name,
107               &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL, &NbST,
108               &NbVT, &NbTT, &NbSS, &NbVS, &NbTS) != 14) {
109       Msg::Error("Read error");
110       return false;
111     }
112     NbT2 = t2l = NbT3 = t3l = 0;
113   }
114   else if(version == 1.1) {
115     Msg::Debug("Detected post-processing view format 1.1");
116     if(fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
117               name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
118               &NbST, &NbVT, &NbTT, &NbSS, &NbVS, &NbTS, &NbT2, &t2l, &NbT3,
119               &t3l) != 18) {
120       Msg::Error("Read error");
121       return false;
122     }
123   }
124   else if(version == 1.2 || version == 1.3) {
125     Msg::Debug("Detected post-processing view format %g", version);
126     if(fscanf(fp,
127               "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
128               "%d %d %d %d %d %d %d %d %d %d %d %d %d\n",
129               name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
130               &NbST, &NbVT, &NbTT, &NbSQ, &NbVQ, &NbTQ, &NbSS, &NbVS, &NbTS,
131               &NbSH, &NbVH, &NbTH, &NbSI, &NbVI, &NbTI, &NbSY, &NbVY, &NbTY,
132               &NbT2, &t2l, &NbT3, &t3l) != 30) {
133       Msg::Error("Read error");
134       return false;
135     }
136   }
137   else if(version == 1.4) {
138     Msg::Debug("Detected post-processing view format 1.4");
139     if(fscanf(fp,
140               "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
141               "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
142               "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
143               name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
144               &NbST, &NbVT, &NbTT, &NbSQ, &NbVQ, &NbTQ, &NbSS, &NbVS, &NbTS,
145               &NbSH, &NbVH, &NbTH, &NbSI, &NbVI, &NbTI, &NbSY, &NbVY, &NbTY,
146               &NbSL2, &NbVL2, &NbTL2, &NbST2, &NbVT2, &NbTT2, &NbSQ2, &NbVQ2,
147               &NbTQ2, &NbSS2, &NbVS2, &NbTS2, &NbSH2, &NbVH2, &NbTH2, &NbSI2,
148               &NbVI2, &NbTI2, &NbSY2, &NbVY2, &NbTY2, &NbT2, &t2l, &NbT3,
149               &t3l) != 51) {
150       Msg::Error("Read error");
151       return false;
152     }
153   }
154   else {
155     Msg::Error("Unknown post-processing file format (version %g)", version);
156     return false;
157   }
158 
159   for(int i = 0; i < (int)strlen(name); i++)
160     if(name[i] == '^') name[i] = ' ';
161 
162   int swap = 0;
163   if(binary) {
164     int testone;
165     if(!fread(&testone, sizeof(int), 1, fp)) {
166       Msg::Error("Read error");
167       return false;
168     }
169     if(testone != 1) {
170       Msg::Info("Swapping bytes from binary file");
171       swap = 1;
172     }
173   }
174 
175   dVecRead(Time, NbTimeStep, fp, binary, swap);
176   dVecRead(SP, NbSP * (NbTimeStep * 1 + 3), fp, binary, swap);
177   dVecRead(VP, NbVP * (NbTimeStep * 3 + 3), fp, binary, swap);
178   dVecRead(TP, NbTP * (NbTimeStep * 9 + 3), fp, binary, swap);
179   dVecRead(SL, NbSL * (NbTimeStep * 2 * 1 + 6), fp, binary, swap);
180   dVecRead(VL, NbVL * (NbTimeStep * 2 * 3 + 6), fp, binary, swap);
181   dVecRead(TL, NbTL * (NbTimeStep * 2 * 9 + 6), fp, binary, swap);
182   dVecRead(ST, NbST * (NbTimeStep * 3 * 1 + 9), fp, binary, swap);
183   dVecRead(VT, NbVT * (NbTimeStep * 3 * 3 + 9), fp, binary, swap);
184   dVecRead(TT, NbTT * (NbTimeStep * 3 * 9 + 9), fp, binary, swap);
185   dVecRead(SQ, NbSQ * (NbTimeStep * 4 * 1 + 12), fp, binary, swap);
186   dVecRead(VQ, NbVQ * (NbTimeStep * 4 * 3 + 12), fp, binary, swap);
187   dVecRead(TQ, NbTQ * (NbTimeStep * 4 * 9 + 12), fp, binary, swap);
188   dVecRead(SS, NbSS * (NbTimeStep * 4 * 1 + 12), fp, binary, swap);
189   dVecRead(VS, NbVS * (NbTimeStep * 4 * 3 + 12), fp, binary, swap);
190   dVecRead(TS, NbTS * (NbTimeStep * 4 * 9 + 12), fp, binary, swap);
191   dVecRead(SH, NbSH * (NbTimeStep * 8 * 1 + 24), fp, binary, swap);
192   dVecRead(VH, NbVH * (NbTimeStep * 8 * 3 + 24), fp, binary, swap);
193   dVecRead(TH, NbTH * (NbTimeStep * 8 * 9 + 24), fp, binary, swap);
194   dVecRead(SI, NbSI * (NbTimeStep * 6 * 1 + 18), fp, binary, swap);
195   dVecRead(VI, NbVI * (NbTimeStep * 6 * 3 + 18), fp, binary, swap);
196   dVecRead(TI, NbTI * (NbTimeStep * 6 * 9 + 18), fp, binary, swap);
197   dVecRead(SY, NbSY * (NbTimeStep * 5 * 1 + 15), fp, binary, swap);
198   dVecRead(VY, NbVY * (NbTimeStep * 5 * 3 + 15), fp, binary, swap);
199   dVecRead(TY, NbTY * (NbTimeStep * 5 * 9 + 15), fp, binary, swap);
200 
201   // overwrite first order data with second order data (if any)
202   dVecRead(SL, NbSL2 * (NbTimeStep * 3 * 1 + 9), fp, binary, swap);
203   dVecRead(VL, NbVL2 * (NbTimeStep * 3 * 3 + 9), fp, binary, swap);
204   dVecRead(TL, NbTL2 * (NbTimeStep * 3 * 9 + 9), fp, binary, swap);
205   dVecRead(ST, NbST2 * (NbTimeStep * 6 * 1 + 18), fp, binary, swap);
206   dVecRead(VT, NbVT2 * (NbTimeStep * 6 * 3 + 18), fp, binary, swap);
207   dVecRead(TT, NbTT2 * (NbTimeStep * 6 * 9 + 18), fp, binary, swap);
208   dVecRead(SQ, NbSQ2 * (NbTimeStep * 9 * 1 + 27), fp, binary, swap);
209   dVecRead(VQ, NbVQ2 * (NbTimeStep * 9 * 3 + 27), fp, binary, swap);
210   dVecRead(TQ, NbTQ2 * (NbTimeStep * 9 * 9 + 27), fp, binary, swap);
211   dVecRead(SS, NbSS2 * (NbTimeStep * 10 * 1 + 30), fp, binary, swap);
212   dVecRead(VS, NbVS2 * (NbTimeStep * 10 * 3 + 30), fp, binary, swap);
213   dVecRead(TS, NbTS2 * (NbTimeStep * 10 * 9 + 30), fp, binary, swap);
214   dVecRead(SH, NbSH2 * (NbTimeStep * 27 * 1 + 81), fp, binary, swap);
215   dVecRead(VH, NbVH2 * (NbTimeStep * 27 * 3 + 81), fp, binary, swap);
216   dVecRead(TH, NbTH2 * (NbTimeStep * 27 * 9 + 81), fp, binary, swap);
217   dVecRead(SI, NbSI2 * (NbTimeStep * 18 * 1 + 54), fp, binary, swap);
218   dVecRead(VI, NbVI2 * (NbTimeStep * 18 * 3 + 54), fp, binary, swap);
219   dVecRead(TI, NbTI2 * (NbTimeStep * 18 * 9 + 54), fp, binary, swap);
220   dVecRead(SY, NbSY2 * (NbTimeStep * 14 * 1 + 42), fp, binary, swap);
221   dVecRead(VY, NbVY2 * (NbTimeStep * 14 * 3 + 42), fp, binary, swap);
222   dVecRead(TY, NbTY2 * (NbTimeStep * 14 * 9 + 42), fp, binary, swap);
223   if(NbSL2) {
224     NbSL = NbSL2;
225     setOrder2(TYPE_LIN);
226   }
227   if(NbVL2) {
228     NbVL = NbVL2;
229     setOrder2(TYPE_LIN);
230   }
231   if(NbTL2) {
232     NbTL = NbTL2;
233     setOrder2(TYPE_LIN);
234   }
235   if(NbST2) {
236     NbST = NbST2;
237     setOrder2(TYPE_TRI);
238   }
239   if(NbVT2) {
240     NbVT = NbVT2;
241     setOrder2(TYPE_TRI);
242   }
243   if(NbTT2) {
244     NbTT = NbTT2;
245     setOrder2(TYPE_TRI);
246   }
247   if(NbSQ2) {
248     NbSQ = NbSQ2;
249     setOrder2(TYPE_QUA);
250   }
251   if(NbVQ2) {
252     NbVQ = NbVQ2;
253     setOrder2(TYPE_QUA);
254   }
255   if(NbTQ2) {
256     NbTQ = NbTQ2;
257     setOrder2(TYPE_QUA);
258   }
259   if(NbSS2) {
260     NbSS = NbSS2;
261     setOrder2(TYPE_TET);
262   }
263   if(NbVS2) {
264     NbVS = NbVS2;
265     setOrder2(TYPE_TET);
266   }
267   if(NbTS2) {
268     NbTS = NbTS2;
269     setOrder2(TYPE_TET);
270   }
271   if(NbSH2) {
272     NbSH = NbSH2;
273     setOrder2(TYPE_HEX);
274   }
275   if(NbVH2) {
276     NbVH = NbVH2;
277     setOrder2(TYPE_HEX);
278   }
279   if(NbTH2) {
280     NbTH = NbTH2;
281     setOrder2(TYPE_HEX);
282   }
283   if(NbSI2) {
284     NbSI = NbSI2;
285     setOrder2(TYPE_PRI);
286   }
287   if(NbVI2) {
288     NbVI = NbVI2;
289     setOrder2(TYPE_PRI);
290   }
291   if(NbTI2) {
292     NbTI = NbTI2;
293     setOrder2(TYPE_PRI);
294   }
295   if(NbSY2) {
296     NbSY = NbSY2;
297     setOrder2(TYPE_PYR);
298   }
299   if(NbVY2) {
300     NbVY = NbVY2;
301     setOrder2(TYPE_PYR);
302   }
303   if(NbTY2) {
304     NbTY = NbTY2;
305     setOrder2(TYPE_PYR);
306   }
307 
308   dVecRead(T2D, NbT2 * 4, fp, binary, swap);
309   cVecRead(T2C, t2l, fp, binary, swap, (version <= 1.2));
310   dVecRead(T3D, NbT3 * 5, fp, binary, swap);
311   cVecRead(T3C, t3l, fp, binary, swap, (version <= 1.2));
312 
313   Msg::Debug("Read View '%s' (%d TimeSteps): "
314              "SP(%d/%d) VP(%d/%d) TP(%d/%d) "
315              "SL(%d/%d) VL(%d/%d) TL(%d/%d) "
316              "ST(%d/%d) VT(%d/%d) TT(%d/%d) "
317              "SQ(%d/%d) VQ(%d/%d) TQ(%d/%d) "
318              "SS(%d/%d) VS(%d/%d) TS(%d/%d) "
319              "SH(%d/%d) VH(%d/%d) TH(%d/%d) "
320              "SI(%d/%d) VI(%d/%d) TI(%d/%d) "
321              "SY(%d/%d) VY(%d/%d) TY(%d/%d) "
322              "T2(%d/%d/%d) T3(%d/%d/%d) ",
323              name, NbTimeStep, NbSP, SP.size(), NbVP, VP.size(), NbTP,
324              TP.size(), NbSL, SL.size(), NbVL, VL.size(), NbTL, TL.size(), NbST,
325              ST.size(), NbVT, VT.size(), NbTT, TT.size(), NbSQ, SQ.size(), NbVQ,
326              VQ.size(), NbTQ, TQ.size(), NbSS, SS.size(), NbVS, VS.size(), NbTS,
327              TS.size(), NbSH, SH.size(), NbVH, VH.size(), NbTH, TH.size(), NbSI,
328              SI.size(), NbVI, VI.size(), NbTI, TI.size(), NbSY, SY.size(), NbVY,
329              VY.size(), NbTY, TY.size(), NbT2, T2D.size(), T2C.size(), NbT3,
330              T3D.size(), T3C.size());
331 
332   setName(name);
333   finalize();
334   return true;
335 }
336 
writeTimePOS(FILE * fp,std::vector<double> & list)337 static void writeTimePOS(FILE *fp, std::vector<double> &list)
338 {
339   if(list.size() > 1) {
340     fprintf(fp, "TIME{");
341     for(std::size_t i = 0; i < list.size(); i++) {
342       if(i) fprintf(fp, ",");
343       fprintf(fp, "%.16g", list[i]);
344     }
345     fprintf(fp, "};\n");
346   }
347 }
348 
writeElementPOS(FILE * fp,const char * str,int nbnod,int nb,std::vector<double> & list)349 static void writeElementPOS(FILE *fp, const char *str, int nbnod, int nb,
350                             std::vector<double> &list)
351 {
352   if(nb) {
353     int n = list.size() / nb;
354     for(std::size_t i = 0; i < list.size(); i += n) {
355       double *x = &list[i];
356       double *y = &list[i + nbnod];
357       double *z = &list[i + 2 * nbnod];
358       fprintf(fp, "%s(", str);
359       for(int j = 0; j < nbnod; j++) {
360         if(j) fprintf(fp, ",");
361         fprintf(fp, "%.16g,%.16g,%.16g", x[j], y[j], z[j]);
362       }
363       fprintf(fp, "){");
364       for(int j = 3 * nbnod; j < n; j++) {
365         if(j - 3 * nbnod) fprintf(fp, ",");
366         fprintf(fp, "%.16g", list[i + j]);
367       }
368       fprintf(fp, "};\n");
369     }
370   }
371 }
372 
writeTextPOS(FILE * fp,int nbc,int nb,std::vector<double> & TD,std::vector<char> & TC)373 static void writeTextPOS(FILE *fp, int nbc, int nb, std::vector<double> &TD,
374                          std::vector<char> &TC)
375 {
376   if(!nb || (nbc != 4 && nbc != 5)) return;
377   for(std::size_t j = 0; j < TD.size(); j += nbc) {
378     double x = TD[j];
379     double y = TD[j + 1];
380     double z = (nbc == 5) ? TD[j + 2] : 0.;
381     double style = TD[j + nbc - 2];
382     if(nbc == 4)
383       fprintf(fp, "T2(%g,%g,%g){", x, y, style);
384     else
385       fprintf(fp, "T3(%g,%g,%g,%g){", x, y, z, style);
386     double start = TD[j + nbc - 1];
387     double end;
388     if(j + nbc * 2 - 1 < TD.size())
389       end = TD[j + nbc * 2 - 1];
390     else
391       end = TC.size();
392     int l = 0;
393     while(l < end - start) {
394       char *str = &TC[(int)start + l];
395       if(l) fprintf(fp, ",");
396       fprintf(fp, "\"%s\"", str);
397       l += strlen(str) + 1;
398     }
399     fprintf(fp, "};\n");
400   }
401 }
402 
writePOS(const std::string & fileName,bool binary,bool parsed,bool append)403 bool PViewDataList::writePOS(const std::string &fileName, bool binary,
404                              bool parsed, bool append)
405 {
406   if(_adaptive) {
407     Msg::Warning(
408       "Writing adapted dataset (will only export current time step)");
409     return _adaptive->getData()->writePOS(fileName, binary, parsed, append);
410   }
411 
412   if(haveInterpolationMatrices()) {
413     Msg::Error(
414       "Cannot export datasets with interpolation matrices in old POS format: "
415       "consider using the new mesh-based format instead, or select 'Adapt "
416       "post-processing data' before exporting");
417     return false;
418   }
419 
420   FILE *fp = Fopen(fileName.c_str(),
421                    append ? (binary ? "ab" : "a") : (binary ? "wb" : "w"));
422   if(!fp) {
423     Msg::Error("Unable to open file '%s'", fileName.c_str());
424     return false;
425   }
426 
427   if(!parsed && !append) {
428     fprintf(fp, "$PostFormat /* Gmsh 1.3, %s */\n",
429             binary ? "binary" : "ascii");
430     fprintf(fp, "1.3 %d %d\n", binary, (int)sizeof(double));
431     fprintf(fp, "$EndPostFormat\n");
432   }
433 
434   std::string str = getName();
435   for(std::size_t i = 0; i < str.size(); i++)
436     if(str[i] == ' ') str[i] = '^';
437 
438   if(!parsed) {
439     fprintf(fp, "$View /* %s */\n", getName().c_str());
440     if(str.empty())
441       fprintf(fp, "noname ");
442     else
443       fprintf(fp, "%s ", str.c_str());
444     fprintf(fp,
445             "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
446             "%d %d %d %d %d %d %d %d %d %d %d %d\n",
447             (int)Time.size(), NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT,
448             NbTT, NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH, NbSI,
449             NbVI, NbTI, NbSY, NbVY, NbTY, NbT2, (int)T2C.size(), NbT3,
450             (int)T3C.size());
451     if(binary) {
452       int one = 1;
453       if(!fwrite(&one, sizeof(int), 1, fp)) {
454         Msg::Error("Write error");
455         fclose(fp);
456         return false;
457       }
458     }
459     dVecWrite(Time, fp, binary);
460     dVecWrite(SP, fp, binary);
461     dVecWrite(VP, fp, binary);
462     dVecWrite(TP, fp, binary);
463     dVecWrite(SL, fp, binary);
464     dVecWrite(VL, fp, binary);
465     dVecWrite(TL, fp, binary);
466     dVecWrite(ST, fp, binary);
467     dVecWrite(VT, fp, binary);
468     dVecWrite(TT, fp, binary);
469     dVecWrite(SQ, fp, binary);
470     dVecWrite(VQ, fp, binary);
471     dVecWrite(TQ, fp, binary);
472     dVecWrite(SS, fp, binary);
473     dVecWrite(VS, fp, binary);
474     dVecWrite(TS, fp, binary);
475     dVecWrite(SH, fp, binary);
476     dVecWrite(VH, fp, binary);
477     dVecWrite(TH, fp, binary);
478     dVecWrite(SI, fp, binary);
479     dVecWrite(VI, fp, binary);
480     dVecWrite(TI, fp, binary);
481     dVecWrite(SY, fp, binary);
482     dVecWrite(VY, fp, binary);
483     dVecWrite(TY, fp, binary);
484     dVecWrite(T2D, fp, binary);
485     cVecWrite(T2C, fp, binary);
486     dVecWrite(T3D, fp, binary);
487     cVecWrite(T3C, fp, binary);
488     fprintf(fp, "\n");
489     fprintf(fp, "$EndView\n");
490   }
491   else {
492     fprintf(fp, "View \"%s\" {\n", getName().c_str());
493     writeTimePOS(fp, Time);
494     writeElementPOS(fp, "SP", 1, NbSP, SP);
495     writeElementPOS(fp, "VP", 1, NbVP, VP);
496     writeElementPOS(fp, "TP", 1, NbTP, TP);
497     writeElementPOS(fp, "SL", 2, NbSL, SL);
498     writeElementPOS(fp, "VL", 2, NbVL, VL);
499     writeElementPOS(fp, "TL", 2, NbTL, TL);
500     writeElementPOS(fp, "ST", 3, NbST, ST);
501     writeElementPOS(fp, "VT", 3, NbVT, VT);
502     writeElementPOS(fp, "TT", 3, NbTT, TT);
503     writeElementPOS(fp, "SQ", 4, NbSQ, SQ);
504     writeElementPOS(fp, "VQ", 4, NbVQ, VQ);
505     writeElementPOS(fp, "TQ", 4, NbTQ, TQ);
506     writeElementPOS(fp, "SS", 4, NbSS, SS);
507     writeElementPOS(fp, "VS", 4, NbVS, VS);
508     writeElementPOS(fp, "TS", 4, NbTS, TS);
509     writeElementPOS(fp, "SH", 8, NbSH, SH);
510     writeElementPOS(fp, "VH", 8, NbVH, VH);
511     writeElementPOS(fp, "TH", 8, NbTH, TH);
512     writeElementPOS(fp, "SI", 6, NbSI, SI);
513     writeElementPOS(fp, "VI", 6, NbVI, VI);
514     writeElementPOS(fp, "TI", 6, NbTI, TI);
515     writeElementPOS(fp, "SY", 5, NbSY, SY);
516     writeElementPOS(fp, "VY", 5, NbVY, VY);
517     writeElementPOS(fp, "TY", 5, NbTY, TY);
518     writeTextPOS(fp, 4, NbT2, T2D, T2C);
519     writeTextPOS(fp, 5, NbT3, T3D, T3C);
520     fprintf(fp, "};\n");
521   }
522 
523   fclose(fp);
524   return true;
525 }
526 
createVertices(std::vector<double> & list,int nbelm,int nbnod,std::vector<MVertex * > & nodes)527 static void createVertices(std::vector<double> &list, int nbelm, int nbnod,
528                            std::vector<MVertex *> &nodes)
529 {
530   if(!nbelm) return;
531   int nb = list.size() / nbelm;
532   for(std::size_t i = 0; i < list.size(); i += nb) {
533     double *x = &list[i];
534     double *y = &list[i + nbnod];
535     double *z = &list[i + 2 * nbnod];
536     for(int j = 0; j < nbnod; j++)
537       nodes.push_back(new MVertex(x[j], y[j], z[j]));
538   }
539 }
540 
541 class nodeData {
542 public:
543   int nbnod;
544   int nod;
545   double *data;
nodeData()546   nodeData() : nbnod(0), nod(0), data(nullptr) {}
nodeData(int _nbnod,int _nod,double * _data)547   nodeData(int _nbnod, int _nod, double *_data)
548     : nbnod(_nbnod), nod(_nod), data(_data)
549   {
550   }
551 };
552 
createElements(std::vector<double> & list,int nbelm,int nbnod,MVertexRTree & pos,std::vector<MElement * > & elements,int type,std::map<MVertex *,nodeData> * vertexData)553 static void createElements(std::vector<double> &list, int nbelm, int nbnod,
554                            MVertexRTree &pos, std::vector<MElement *> &elements,
555                            int type, std::map<MVertex *, nodeData> *vertexData)
556 {
557   if(!nbelm) return;
558   int t = 0;
559   // reverse-engineer geometrical element type according to the number
560   // of nodes (this should be completed, but is likely enough for most
561   // legacy .pos files out there...)
562   switch(type) {
563   case TYPE_PNT: t = MSH_PNT; break;
564   case TYPE_LIN:
565     switch(nbnod) {
566     case 2: t = MSH_LIN_2; break;
567     case 3: t = MSH_LIN_3; break;
568     }
569     break;
570   case TYPE_TRI:
571     switch(nbnod) {
572     case 3: t = MSH_TRI_3; break;
573     case 6: t = MSH_TRI_6; break;
574     }
575     break;
576   case TYPE_QUA:
577     switch(nbnod) {
578     case 4: t = MSH_QUA_4; break;
579     case 8: t = MSH_QUA_8; break;
580     case 9: t = MSH_QUA_9; break;
581     }
582     break;
583   case TYPE_TET:
584     switch(nbnod) {
585     case 4: t = MSH_TET_4; break;
586     case 10: t = MSH_TET_10; break;
587     }
588     break;
589   case TYPE_HEX:
590     switch(nbnod) {
591     case 8: t = MSH_HEX_8; break;
592     case 20: t = MSH_HEX_20; break;
593     case 27: t = MSH_HEX_27; break;
594     }
595     break;
596   case TYPE_PRI:
597     switch(nbnod) {
598     case 6: t = MSH_PRI_6; break;
599     case 15: t = MSH_PRI_15; break;
600     case 18: t = MSH_PRI_18; break;
601     }
602     break;
603   case TYPE_PYR:
604     switch(nbnod) {
605     case 5: t = MSH_PYR_5; break;
606     case 13: t = MSH_PYR_13; break;
607     case 14: t = MSH_PYR_14; break;
608     }
609     break;
610   }
611   if(!t) {
612     Msg::Warning("Discarding elements of type (%d nodes)", nbnod);
613     return;
614   }
615   MElementFactory factory;
616   int nb = list.size() / nbelm;
617   for(std::size_t i = 0; i < list.size(); i += nb) {
618     double *x = &list[i];
619     double *y = &list[i + nbnod];
620     double *z = &list[i + 2 * nbnod];
621     std::vector<MVertex *> verts(nbnod);
622     for(int j = 0; j < nbnod; j++) {
623       verts[j] = pos.find(x[j], y[j], z[j]);
624       if(vertexData)
625         (*vertexData)[verts[j]] = nodeData(nbnod, j, &list[i + 3 * nbnod]);
626     }
627     MElement *e = factory.create(t, verts);
628     elements.push_back(e);
629   }
630 }
631 
writeMSH(const std::string & fileName,double version,bool binary,bool saveMesh,bool multipleView,int partitionNum,bool saveInterpolationMatrices,bool forceNodeData,bool forceElementData)632 bool PViewDataList::writeMSH(const std::string &fileName, double version,
633                              bool binary, bool saveMesh, bool multipleView,
634                              int partitionNum, bool saveInterpolationMatrices,
635                              bool forceNodeData, bool forceElementData)
636 {
637   if(_adaptive) {
638     Msg::Warning(
639       "Writing adapted dataset (will only export current time step)");
640     return _adaptive->getData()->writeMSH(fileName, version, binary);
641   }
642 
643   FILE *fp = Fopen(fileName.c_str(), "w");
644   if(!fp) {
645     Msg::Error("Unable to open file '%s'", fileName.c_str());
646     return false;
647   }
648 
649   double tol = CTX::instance()->geom.tolerance;
650   double eps = norm(SVector3(BBox.max(), BBox.min())) * tol;
651 
652   std::vector<MVertex *> vertices;
653   std::vector<MElement *> elements;
654 
655   int numComponents = 9;
656   for(int i = 0; i < 24; i++) {
657     std::vector<double> *list = nullptr;
658     int *numEle = nullptr, numNodes, numComp;
659     _getRawData(i, &list, &numEle, &numComp, &numNodes);
660     if(*numEle) numComponents = std::min(numComponents, numComp);
661     createVertices(*list, *numEle, numNodes, vertices);
662   }
663   MVertexRTree pos(eps);
664   std::vector<MVertex *> unique;
665   for(std::size_t i = 0; i < vertices.size(); i++) {
666     if(!pos.insert(vertices[i])) unique.push_back(vertices[i]);
667   }
668   vertices.clear();
669 
670   std::map<MVertex *, nodeData> vertexData;
671 
672   for(int i = 0; i < 24; i++) {
673     std::vector<double> *list = nullptr;
674     int *numEle = nullptr, numComp, numNodes;
675     int typ = _getRawData(i, &list, &numEle, &numComp, &numNodes);
676     createElements(*list, *numEle, numNodes, pos, elements, typ,
677                    forceNodeData ? &vertexData : nullptr);
678   }
679 
680   int num = 0;
681   for(std::size_t i = 0; i < unique.size(); i++) unique[i]->setIndex(++num);
682 
683   if(version > 2.2)
684     Msg::Warning("Mesh-based export of list-based datasets not available with "
685                  "MSH %g: using MSH 2.2",
686                  version);
687 
688   fprintf(fp, "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n");
689 
690   if(saveMesh) {
691     fprintf(fp, "$Nodes\n");
692     fprintf(fp, "%d\n", (int)unique.size());
693     for(std::size_t i = 0; i < unique.size(); i++) {
694       MVertex *v = unique[i];
695       fprintf(fp, "%ld %.16g %.16g %.16g\n", v->getIndex(), v->x(), v->y(),
696               v->z());
697     }
698     fprintf(fp, "$EndNodes\n");
699 
700     fprintf(fp, "$Elements\n");
701     fprintf(fp, "%d\n", (int)elements.size());
702     for(std::size_t i = 0; i < elements.size(); i++) {
703       elements[i]->writeMSH2(fp, 2.2, false, i + 1);
704     }
705     fprintf(fp, "$EndElements\n");
706   }
707 
708   if(saveInterpolationMatrices && haveInterpolationMatrices() &&
709      !forceNodeData && !forceElementData) {
710     fprintf(fp, "$InterpolationScheme\n");
711     fprintf(fp, "\"INTERPOLATION_SCHEME\"\n");
712     fprintf(fp, "%d\n", (int)_interpolation.size());
713     for(auto it = _interpolation.begin(); it != _interpolation.end(); it++) {
714       if(it->second.size() >= 2) {
715         fprintf(fp, "%d\n2\n", it->first);
716         for(int mat = 0; mat < 2; mat++) {
717           int m = it->second[mat]->size1(), n = it->second[mat]->size2();
718           fprintf(fp, "%d %d\n", m, n);
719           for(int i = 0; i < m; i++) {
720             for(int j = 0; j < n; j++)
721               fprintf(fp, "%.16g ", it->second[mat]->get(i, j));
722             fprintf(fp, "\n");
723           }
724         }
725       }
726     }
727     fprintf(fp, "$EndInterpolationScheme\n");
728   }
729 
730   for(int ts = 0; ts < NbTimeStep; ts++) {
731     if(forceNodeData)
732       fprintf(fp, "$NodeData\n");
733     else if(forceElementData)
734       fprintf(fp, "$ElementData\n");
735     else
736       fprintf(fp, "$ElementNodeData\n");
737     if(saveInterpolationMatrices && haveInterpolationMatrices() &&
738        !forceNodeData && !forceElementData)
739       fprintf(fp, "2\n\"%s\"\n\"INTERPOLATION_SCHEME\"\n", getName().c_str());
740     else
741       fprintf(fp, "1\n\"%s\"\n", getName().c_str());
742     fprintf(fp, "1\n%.16g\n", getTime(ts));
743     int size = forceNodeData ? (int)unique.size() : (int)elements.size();
744     if(partitionNum > 0)
745       fprintf(fp, "4\n%d\n%d\n%d\n%d\n", ts, numComponents, size, partitionNum);
746     else
747       fprintf(fp, "3\n%d\n%d\n%d\n", ts, numComponents, size);
748 
749     if(forceNodeData) {
750       for(std::size_t i = 0; i < unique.size(); i++) {
751         MVertex *v = unique[i];
752         fprintf(fp, "%ld", v->getIndex());
753         int nbnod = vertexData[v].nbnod;
754         int nod = vertexData[v].nod;
755         double *d = vertexData[v].data;
756         for(int j = 0; j < numComponents; j++)
757           fprintf(fp, " %.16g",
758                   d[numComponents * nbnod * ts + numComponents * nod + j]);
759         fprintf(fp, "\n");
760       }
761       fprintf(fp, "$EndNodeData\n");
762     }
763     else {
764       int n = 0;
765       for(int i = 0; i < 24; i++) {
766         std::vector<double> *list = nullptr;
767         int *numEle = nullptr, numComp, numNodes;
768         int typ = _getRawData(i, &list, &numEle, &numComp, &numNodes);
769         if(*numEle) {
770           int mult = numNodes;
771           if(_interpolation.count(typ)) mult = _interpolation[typ][0]->size1();
772           int nb = list->size() / *numEle;
773           for(std::size_t i = 0; i < list->size(); i += nb) {
774             double *v = &(*list)[i + 3 * numNodes];
775             if(forceElementData) { // just keep first vertex value
776               fprintf(fp, "%d", ++n);
777               for(int j = 0; j < numComponents; j++)
778                 fprintf(fp, " %.16g", v[numComponents * mult * ts + j]);
779             }
780             else {
781               fprintf(fp, "%d %d", ++n, mult);
782               for(int j = 0; j < numComponents * mult; j++)
783                 fprintf(fp, " %.16g", v[numComponents * mult * ts + j]);
784             }
785             fprintf(fp, "\n");
786           }
787         }
788       }
789       if(forceElementData)
790         fprintf(fp, "$EndElementData\n");
791       else
792         fprintf(fp, "$EndElementNodeData\n");
793     }
794   }
795 
796   fclose(fp);
797   return true;
798 }
799 
importLists(int N[24],std::vector<double> * V[24])800 void PViewDataList::importLists(int N[24], std::vector<double> *V[24])
801 {
802   for(int i = 0; i < 24; i++) {
803     std::vector<double> *list = nullptr;
804     int *nbe = nullptr, nbc, nbn;
805     _getRawData(i, &list, &nbe, &nbc, &nbn);
806     *nbe = N[i];
807     *list = *V[i]; // deep copy
808   }
809   finalize();
810 }
811 
importList(int index,int n,const std::vector<double> & v,bool fin)812 void PViewDataList::importList(int index, int n, const std::vector<double> &v,
813                                bool fin)
814 {
815   if(index < 0 || index >= 24) {
816     Msg::Error("Wrong list index to import");
817     return;
818   }
819   std::vector<double> *list = nullptr;
820   int *nbe = nullptr, nbc, nbn;
821   _getRawData(index, &list, &nbe, &nbc, &nbn);
822   *nbe = n;
823   *list = v; // deep copy
824   if(fin) finalize();
825 }
826 
getListPointers(int N[24],std::vector<double> * V[24])827 void PViewDataList::getListPointers(int N[24], std::vector<double> *V[24])
828 {
829   for(int i = 0; i < 24; i++) {
830     std::vector<double> *list = nullptr;
831     int *nbe = nullptr, nbc, nbn;
832     _getRawData(i, &list, &nbe, &nbc, &nbn);
833     N[i] = *nbe;
834     V[i] = list; // copy pointer only
835   }
836 }
837