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