1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5
6 #include <sstream>
7 #include <map>
8 #include <vector>
9 #include <list>
10 #include <string.h>
11 #include <math.h>
12 #include <set>
13
14 #include "GetDPVersion.h"
15 #include "GetDPConfig.h"
16 #include "ProData.h"
17 #include "ProDefine.h"
18 #include "ProParser.h"
19 #include "GeoData.h"
20 #include "DofData.h"
21 #include "Pos_Iso.h"
22 #include "Pos_Format.h"
23 #include "Pos_Element.h"
24 #include "Pos_Formulation.h"
25 #include "F.h"
26 #include "Cal_Value.h"
27 #include "Cal_Quantity.h"
28 #include "MallocUtils.h"
29 #include "Message.h"
30 #include "kissfft.hh"
31
32 #if defined(HAVE_GMSH)
33 #include <gmsh/PView.h>
34 #include <gmsh/PViewDataList.h>
35 #include <gmsh/PViewDataGModel.h>
36 #endif
37
38 #define TWO_PI 6.2831853071795865
39
40 #define NBR_MAX_ISO 200
41
42 #define SQU(a) ((a)*(a))
43
44 extern struct Problem Problem_S ;
45 extern struct CurrentData Current ;
46
47 extern int Flag_BIN, Flag_GMSH_VERSION;
48
49 extern FILE *PostStream ;
50
51 static List_T *PostElement_L = NULL ;
52 static List_T *TimeValue_L = NULL ;
53
54 /* ------------------------------------------------------------------------ */
55 /* Gmsh formats */
56 /* ------------------------------------------------------------------------ */
57
58 // global static lists for new-style Gmsh output (cannot be saved incrementally
59 // for each element)
60 static int Gmsh_StartNewView = 0 ;
61 static int NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT;
62 static int NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH;
63 static int NbSI, NbVI, NbTI, NbSY, NbVY, NbTY;
64 static int NbT2;
65 static std::vector<double> SP, VP, TP, SL, VL, TL, ST, VT, TT, SQ, VQ, TQ;
66 static std::vector<double> SS, VS, TS1 /* for petsc */, SH, VH, TH;
67 static std::vector<double> SI, VI, TI, SY, VY, TY, T2D;
68 static std::vector<char> T2C;
69 static char CurrentName[256] = "";
70 static int CurrentPartitionNumber = 0;
71
Gmsh_ResetStaticLists()72 static void Gmsh_ResetStaticLists()
73 {
74 NbSP = NbVP = NbTP = NbSL = NbVL = NbTL = 0;
75 NbST = NbVT = NbTT = NbSQ = NbVQ = NbTQ = 0;
76 NbSS = NbVS = NbTS = NbSH = NbVH = NbTH = 0;
77 NbSI = NbVI = NbTI = NbSY = NbVY = NbTY = 0;
78 NbT2 = 0;
79 SP.clear(); VP.clear(); TP.clear(); SL.clear(); VL.clear(); TL.clear();
80 ST.clear(); VT.clear(); TT.clear(); SQ.clear(); VQ.clear(); TQ.clear();
81 SS.clear(); VS.clear(); TS1.clear(); SH.clear(); VH.clear(); TH.clear();
82 SI.clear(); VI.clear(); TI.clear(); SY.clear(); VY.clear(); TY.clear();
83 T2D.clear(); T2C.clear();
84 if(!TimeValue_L) TimeValue_L = List_Create(100,1000000,sizeof(double));
85 else List_Reset(TimeValue_L);
86 }
87
Gmsh_StringStart(int Format,double x,double y,double style)88 static void Gmsh_StringStart(int Format, double x, double y, double style)
89 {
90 if(Flag_BIN){ /* bricolage: should use Format instead */
91 T2D.push_back(x);
92 T2D.push_back(y);
93 T2D.push_back(style);
94 T2D.push_back(T2C.size());
95 NbT2++;
96 }
97 else if(PostStream){
98 fprintf(PostStream, "T2(%g,%g,%g){", x, y, style);
99 }
100 }
101
Gmsh_StringAdd(int Format,int first,char * text)102 static void Gmsh_StringAdd(int Format, int first, char *text)
103 {
104 int i;
105 if(Flag_BIN){ /* bricolage: should use Format instead */
106 for(i = 0; i < (int)strlen(text)+1; i++)
107 T2C.push_back(text[i]);
108 }
109 else if(PostStream){
110 if(!first)
111 fprintf(PostStream, ",");
112 fprintf(PostStream, "\"%s\"", text);
113 }
114 }
115
Gmsh_StringEnd(int Format)116 static void Gmsh_StringEnd(int Format)
117 {
118 if(Flag_BIN){ /* bricolage: should use Format instead */
119 }
120 else if(PostStream){
121 fprintf(PostStream, "};\n") ;
122 }
123 }
124
Gmsh_PrintElementNodeData(struct PostSubOperation * PSO_P,int numTimeStep,int numComp,int Nb[8],std::vector<double> * L[8])125 static void Gmsh_PrintElementNodeData(struct PostSubOperation *PSO_P, int numTimeStep,
126 int numComp, int Nb[8], std::vector<double> *L[8])
127 {
128 if(!PostStream) return;
129 int N = 0;
130 for(int i = 0; i < 8; i++) N += Nb[i];
131 if(!N) return;
132 int step = 0;
133 for (int ts = 0; ts < numTimeStep; ts++) {
134 Pos_InitAllSolutions(PSO_P->TimeStep_L, ts);
135 for(int har = 0; har < Current.NbrHar; har++){
136 fprintf(PostStream, "$ElementNodeData\n");
137 fprintf(PostStream, "1\n");
138 fprintf(PostStream, "\"%s\"\n", CurrentName);
139 fprintf(PostStream, "1\n");
140 fprintf(PostStream, "%.16g\n", Current.Time);
141 fprintf(PostStream, "4\n");
142 fprintf(PostStream, "%d\n", (PSO_P->OverrideTimeStepValue >= 0) ?
143 PSO_P->OverrideTimeStepValue : (Current.NbrHar > 1 ? step :
144 (int)Current.TimeStep));
145 fprintf(PostStream, "%d\n", numComp);
146 fprintf(PostStream, "%d\n", N);
147 fprintf(PostStream, "%d\n", CurrentPartitionNumber);
148 for(int i = 0; i < 8; i++){
149 if(!Nb[i]) continue;
150 int stride = (*L[i]).size() / Nb[i];
151 for(unsigned int j = 0; j < (*L[i]).size(); j += stride){
152 double *tmp = &(*L[i])[j];
153 int num = (int)tmp[0];
154 int mult = (stride - 1) / numTimeStep / Current.NbrHar / numComp;
155 if(Flag_BIN){
156 fwrite(&num, sizeof(int), 1, PostStream);
157 fwrite(&mult, sizeof(int), 1, PostStream);
158 fwrite(&tmp[1 + step * mult * numComp], sizeof(double), mult * numComp,
159 PostStream);
160 }
161 else{
162 fprintf(PostStream, "%d %d", num, mult);
163 for(int k = 0; k < mult * numComp; k++)
164 fprintf(PostStream, " %.16g", tmp[1 + step * mult * numComp + k]);
165 fprintf(PostStream, "\n");
166 }
167 }
168 }
169 fprintf(PostStream, "$EndElementNodeData\n");
170 step++;
171 }
172 }
173 }
174
GmshParsed_PrintElement(double Time,int TimeStep,int NbTimeStep,int NbHarmonic,int HarmonicToTime,int Type,int NbrNodes,double * x,double * y,double * z,struct Value * Value)175 static void GmshParsed_PrintElement(double Time, int TimeStep, int NbTimeStep, int NbHarmonic,
176 int HarmonicToTime, int Type, int NbrNodes,
177 double *x, double *y, double *z, struct Value *Value)
178 {
179 int i,j,k,jj ;
180 double TimeMH ;
181 struct Value TmpValue ;
182 int symIndex[9] = {0, 1, 2, 1, 3, 4, 2, 4, 5} ;
183 int diagIndex[9] = {0, -1, -1, -1, 1, -1, -1, -1, 2} ;
184
185 if(Gmsh_StartNewView){
186 Gmsh_StartNewView = 0 ;
187 Gmsh_ResetStaticLists();
188 }
189
190 if (HarmonicToTime == 1){
191 if(NbHarmonic == 2 && NbTimeStep == 1){ // classical complex case
192 double zero = 0., one = 1.;
193 List_Put(TimeValue_L, 0, &zero);
194 List_Put(TimeValue_L, 1, &one);
195 }
196 else{
197 for(k = 0 ; k < NbHarmonic ; k++)
198 List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
199 }
200 }
201 else
202 for(k = 0 ; k < HarmonicToTime ; k++)
203 List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);
204
205 if(!PostStream) return;
206
207 switch (Value[0].Type) {
208
209 case SCALAR :
210
211 if(TimeStep == 0){
212 switch(Type){
213 case POINT_ELEMENT : fprintf(PostStream, "SP("); break;
214 case LINE : fprintf(PostStream, "SL("); break;
215 case TRIANGLE : fprintf(PostStream, "ST("); break;
216 case QUADRANGLE : fprintf(PostStream, "SQ("); break;
217 case TETRAHEDRON : fprintf(PostStream, "SS("); break;
218 case HEXAHEDRON : fprintf(PostStream, "SH("); break;
219 case PRISM : fprintf(PostStream, "SI("); break;
220 case PYRAMID : fprintf(PostStream, "SY("); break;
221 }
222 for(i = 0 ; i < NbrNodes ; i++){
223 if(i) fprintf(PostStream, ",");
224 fprintf(PostStream, "%.16g,%.16g,%.16g", x[i], y[i], z[i]);
225 }
226 fprintf(PostStream, "){");
227 }
228
229 if (HarmonicToTime == 1) {
230 for(k = 0 ; k < NbHarmonic ; k++) {
231 if(k || TimeStep) fprintf(PostStream, ",");
232 for(i = 0 ; i < NbrNodes ; i++){
233 if(i) fprintf(PostStream, ",");
234 fprintf(PostStream, "%.16g", Value[i].Val[MAX_DIM*k]);
235 }
236 }
237 }
238 else {
239 for(k = 0 ; k < HarmonicToTime ; k++){
240 if(k || TimeStep) fprintf(PostStream, ",");
241 for(i = 0 ; i < NbrNodes ; i++){
242 if(i) fprintf(PostStream, ",");
243 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
244 fprintf(PostStream, "%.16g", TmpValue.Val[0]);
245 }
246 }
247 }
248
249 if(TimeStep == NbTimeStep-1){
250 fprintf(PostStream, "};\n") ;
251 }
252 break ;
253
254 case VECTOR :
255
256 if(TimeStep == 0){
257 switch(Type){
258 case POINT_ELEMENT : fprintf(PostStream, "VP("); break;
259 case LINE : fprintf(PostStream, "VL("); break;
260 case TRIANGLE : fprintf(PostStream, "VT("); break;
261 case QUADRANGLE : fprintf(PostStream, "VQ("); break;
262 case TETRAHEDRON : fprintf(PostStream, "VS("); break;
263 case HEXAHEDRON : fprintf(PostStream, "VH("); break;
264 case PRISM : fprintf(PostStream, "VI("); break;
265 case PYRAMID : fprintf(PostStream, "VY("); break;
266 }
267 for(i = 0 ; i < NbrNodes ; i++){
268 if(i) fprintf(PostStream, ",");
269 fprintf(PostStream, "%.16g,%.16g,%.16g", x[i], y[i], z[i]);
270 }
271 fprintf(PostStream, "){");
272 }
273
274 if (HarmonicToTime == 1) {
275 for(k = 0 ; k < NbHarmonic ; k++) {
276 if(k || TimeStep) fprintf(PostStream, ",");
277 for(i = 0 ; i < NbrNodes ; i++){
278 if(i) fprintf(PostStream, ",");
279 for(j = 0 ; j < 3 ; j++){
280 if(j) fprintf(PostStream, ",");
281 fprintf(PostStream, "%.16g", Value[i].Val[MAX_DIM*k+j]);
282 }
283 }
284 }
285 }
286 else {
287 for(k = 0 ; k < HarmonicToTime ; k++){
288 if(k || TimeStep) fprintf(PostStream, ",");
289 for(i = 0 ; i < NbrNodes ; i++){
290 if(i) fprintf(PostStream, ",");
291 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
292 for(j = 0 ; j < 3 ; j++){
293 if(j) fprintf(PostStream, ",");
294 fprintf(PostStream, "%.16g", TmpValue.Val[j]);
295 }
296 }
297 }
298 }
299
300 if(TimeStep == NbTimeStep-1){
301 fprintf(PostStream, "};\n") ;
302 }
303 break ;
304
305 case TENSOR_DIAG :
306 case TENSOR_SYM :
307 case TENSOR :
308
309 if(TimeStep == 0){
310 switch(Type){
311 case POINT_ELEMENT : fprintf(PostStream, "TP("); break;
312 case LINE : fprintf(PostStream, "TL("); break;
313 case TRIANGLE : fprintf(PostStream, "TT("); break;
314 case QUADRANGLE : fprintf(PostStream, "TQ("); break;
315 case TETRAHEDRON : fprintf(PostStream, "TS("); break;
316 case HEXAHEDRON : fprintf(PostStream, "TH("); break;
317 case PRISM : fprintf(PostStream, "TI("); break;
318 case PYRAMID : fprintf(PostStream, "TY("); break;
319 }
320 for(i = 0 ; i < NbrNodes ; i++){
321 if(i) fprintf(PostStream, ",");
322 fprintf(PostStream, "%.16g,%.16g,%.16g", x[i], y[i], z[i]);
323 }
324 fprintf(PostStream, "){");
325 }
326
327 if (HarmonicToTime == 1) {
328 for(k = 0 ; k < NbHarmonic ; k++) {
329 if(k || TimeStep) fprintf(PostStream, ",");
330 for(i = 0 ; i < NbrNodes ; i++){
331 if(i) fprintf(PostStream, ",");
332 for(j = 0 ; j < 9 ; j++){
333 if(j) fprintf(PostStream, ",");
334 if(Value[0].Type != TENSOR_DIAG) {
335 if(Value[0].Type == TENSOR_SYM) jj = symIndex[j];
336 else jj = j;
337 fprintf(PostStream, "%.16g", Value[i].Val[MAX_DIM*k+jj]);
338 }
339 else {
340 jj = diagIndex[j];
341 if(jj == -1) fprintf(PostStream, "%.16g", 0.);
342 else fprintf(PostStream, "%.16g", Value[i].Val[MAX_DIM*k+jj]);
343 }
344 }
345 }
346 }
347 }
348 else {
349 for(k = 0 ; k < HarmonicToTime ; k++){
350 if(k || TimeStep) fprintf(PostStream, ",");
351 for(i = 0 ; i < NbrNodes ; i++){
352 if(i) fprintf(PostStream, ",");
353 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
354 for(j = 0 ; j < 9 ; j++){
355 if(j) fprintf(PostStream, ",");
356 if(Value[0].Type != TENSOR_DIAG) {
357 if(Value[0].Type == TENSOR_SYM) jj = symIndex[j];
358 else jj = j;
359 jj = symIndex[j];
360 fprintf(PostStream, "%.16g", TmpValue.Val[jj]);
361 }
362 else {
363 jj = diagIndex[j];
364 if(jj == -1) fprintf(PostStream, "%.16g", 0.);
365 else fprintf(PostStream, "%.16g", TmpValue.Val[jj]);
366 }
367 }
368 }
369 }
370 }
371
372 if(TimeStep == NbTimeStep-1){
373 fprintf(PostStream, "};\n") ;
374 }
375 break ;
376
377 }
378
379 }
380
Gmsh_PrintElement(double Time,int TimeStep,int NbTimeStep,int NbHarmonic,int HarmonicToTime,int Type,int ElementNum,int NbrNodes,double * x,double * y,double * z,struct Value * Value,struct PostSubOperation * PSO_P,int Store)381 static void Gmsh_PrintElement(double Time, int TimeStep, int NbTimeStep, int NbHarmonic,
382 int HarmonicToTime, int Type, int ElementNum, int NbrNodes,
383 double *x, double *y, double *z, struct Value *Value,
384 struct PostSubOperation *PSO_P, int Store)
385 {
386 int i,j,k,jj ;
387 double TimeMH ;
388 struct Value TmpValue ;
389 static std::vector<double> *Current_L ;
390 int symIndex[9] = {0, 1, 2, 1, 3, 4, 2, 4, 5} ;
391 int diagIndex[9] = {0, -1, -1, -1, 1, -1, -1, -1, 2} ;
392
393 if(Gmsh_StartNewView){
394 Gmsh_StartNewView = 0 ;
395 Gmsh_ResetStaticLists();
396 }
397
398 switch (Value[0].Type) {
399
400 case SCALAR :
401
402 if(TimeStep == 0){
403 switch(Type){
404 case POINT_ELEMENT : Current_L = &SP ; NbSP++ ; break ;
405 case LINE : Current_L = &SL ; NbSL++ ; break ;
406 case TRIANGLE : Current_L = &ST ; NbST++ ; break ;
407 case QUADRANGLE : Current_L = &SQ ; NbSQ++ ; break ;
408 case TETRAHEDRON : Current_L = &SS ; NbSS++ ; break ;
409 case HEXAHEDRON : Current_L = &SH ; NbSH++ ; break ;
410 case PRISM : Current_L = &SI ; NbSI++ ; break ;
411 case PYRAMID : Current_L = &SY ; NbSY++ ; break ;
412
413 case LINE_2 : Current_L = &SL ; NbSL++ ; break ;
414 case TRIANGLE_2 : Current_L = &ST ; NbST++ ; break ;
415 case QUADRANGLE_2 : Current_L = &SQ ; NbSQ++ ; break ;
416 case QUADRANGLE_2_8N: Current_L = &SQ ; NbSQ++ ; break ;
417 case TETRAHEDRON_2 : Current_L = &SS ; NbSS++ ; break ;
418 case HEXAHEDRON_2 : Current_L = &SH ; NbSH++ ; break ;
419 case HEXAHEDRON_2_20N: Current_L = &SH ; NbSH++ ; break ;
420 case PRISM_2 : Current_L = &SI ; NbSI++ ; break ;
421 case PRISM_2_15N : Current_L = &SI ; NbSI++ ; break ;
422 case PYRAMID_2 : Current_L = &SY ; NbSY++ ; break ;
423 case PYRAMID_2_13N : Current_L = &SY ; NbSY++ ; break ;
424
425 case LINE_3 : Current_L = &SL ; NbSL++ ; break ;
426 case TRIANGLE_3 : Current_L = &ST ; NbST++ ; break ;
427 case QUADRANGLE_3 : Current_L = &SQ ; NbSQ++ ; break ;
428 case TETRAHEDRON_3 : Current_L = &SS ; NbSS++ ; break ;
429 case HEXAHEDRON_3 : Current_L = &SH ; NbSH++ ; break ;
430 case PRISM_3 : Current_L = &SI ; NbSI++ ; break ;
431 case PYRAMID_3 : Current_L = &SY ; NbSY++ ; break ;
432
433 case LINE_4 : Current_L = &SL ; NbSL++ ; break ;
434 case TRIANGLE_4 : Current_L = &ST ; NbST++ ; break ;
435 case QUADRANGLE_4 : Current_L = &SQ ; NbSQ++ ; break ;
436 case TETRAHEDRON_4 : Current_L = &SS ; NbSS++ ; break ;
437 case HEXAHEDRON_4 : Current_L = &SH ; NbSH++ ; break ;
438 case PRISM_4 : Current_L = &SI ; NbSI++ ; break ;
439 //case PYRAMID_4 : Current_L = &SY ; NbSY++ ; break ;
440 }
441 if(Flag_GMSH_VERSION != 2){
442 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(x[i]);
443 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(y[i]);
444 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(z[i]);
445 }
446 else{
447 double tmp = ElementNum;
448 Current_L->push_back(tmp);
449 }
450 }
451 if (HarmonicToTime == 1)
452 for(k = 0 ; k < NbHarmonic ; k++){
453 List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
454 for(i = 0 ; i < NbrNodes ; i++)
455 Current_L->push_back(Value[i].Val[MAX_DIM*k]);
456 }
457 else
458 for(k = 0 ; k < HarmonicToTime ; k++){
459 List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);
460 for(i = 0 ; i < NbrNodes ; i++){
461 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
462 Current_L->push_back(TmpValue.Val[0]);
463 }
464 }
465 break ;
466
467 case VECTOR :
468
469 if(TimeStep == 0){
470 switch(Type){
471 case POINT_ELEMENT : Current_L = &VP ; NbVP++ ; break ;
472
473 case LINE : Current_L = &VL ; NbVL++ ; break ;
474 case TRIANGLE : Current_L = &VT ; NbVT++ ; break ;
475 case QUADRANGLE : Current_L = &VQ ; NbVQ++ ; break ;
476 case TETRAHEDRON : Current_L = &VS ; NbVS++ ; break ;
477 case HEXAHEDRON : Current_L = &VH ; NbVH++ ; break ;
478 case PRISM : Current_L = &VI ; NbVI++ ; break ;
479 case PYRAMID : Current_L = &VY ; NbVY++ ; break ;
480
481 case LINE_2 : Current_L = &VL ; NbVL++ ; break ;
482 case TRIANGLE_2 : Current_L = &VT ; NbVT++ ; break ;
483 case QUADRANGLE_2 : Current_L = &VQ ; NbVQ++ ; break ;
484 case QUADRANGLE_2_8N: Current_L = &VQ ; NbVQ++ ; break ;
485 case TETRAHEDRON_2 : Current_L = &VS ; NbVS++ ; break ;
486 case HEXAHEDRON_2 : Current_L = &VH ; NbVH++ ; break ;
487 case HEXAHEDRON_2_20N: Current_L = &VH ; NbVH++ ; break ;
488 case PRISM_2 : Current_L = &VI ; NbVI++ ; break ;
489 case PRISM_2_15N : Current_L = &VI ; NbVI++ ; break ;
490 case PYRAMID_2 : Current_L = &VY ; NbVY++ ; break ;
491 case PYRAMID_2_13N : Current_L = &VY ; NbVY++ ; break ;
492
493 case LINE_3 : Current_L = &VL ; NbVL++ ; break ;
494 case TRIANGLE_3 : Current_L = &VT ; NbVT++ ; break ;
495 case QUADRANGLE_3 : Current_L = &VQ ; NbVQ++ ; break ;
496 case TETRAHEDRON_3 : Current_L = &VS ; NbVS++ ; break ;
497 case HEXAHEDRON_3 : Current_L = &VH ; NbVH++ ; break ;
498 case PRISM_3 : Current_L = &VI ; NbVI++ ; break ;
499 case PYRAMID_3 : Current_L = &VY ; NbVY++ ; break ;
500
501 case LINE_4 : Current_L = &VL ; NbVL++ ; break ;
502 case TRIANGLE_4 : Current_L = &VT ; NbVT++ ; break ;
503 case QUADRANGLE_4 : Current_L = &VQ ; NbVQ++ ; break ;
504 case TETRAHEDRON_4 : Current_L = &VS ; NbVS++ ; break ;
505 case HEXAHEDRON_4 : Current_L = &VH ; NbVH++ ; break ;
506 case PRISM_4 : Current_L = &VI ; NbVI++ ; break ;
507 //case PYRAMID_4 : Current_L = &VY ; NbVY++ ; break ;
508 }
509 if(Flag_GMSH_VERSION != 2){
510 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(x[i]);
511 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(y[i]);
512 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(z[i]);
513 }
514 else{
515 double tmp = ElementNum;
516 Current_L->push_back(tmp);
517 }
518 }
519 if (HarmonicToTime == 1)
520 for(k = 0 ; k < NbHarmonic ; k++){
521 List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
522 for(i = 0 ; i < NbrNodes ; i++)
523 for(j = 0 ; j < 3 ; j++)
524 Current_L->push_back(Value[i].Val[MAX_DIM*k+j]);
525 }
526 else
527 for(k = 0 ; k < HarmonicToTime ; k++){
528 List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);
529 for(i = 0 ; i < NbrNodes ; i++){
530 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
531 for(j = 0 ; j < 3 ; j++)
532 Current_L->push_back(TmpValue.Val[j]);
533 }
534 }
535 break ;
536
537 case TENSOR_DIAG :
538 case TENSOR_SYM :
539 case TENSOR :
540
541 if(TimeStep == 0){
542 switch(Type){
543 case POINT_ELEMENT : Current_L = &TP ; NbTP++ ; break ;
544
545 case LINE : Current_L = &TL ; NbTL++ ; break ;
546 case TRIANGLE : Current_L = &TT ; NbTT++ ; break ;
547 case QUADRANGLE : Current_L = &TQ ; NbTQ++ ; break ;
548 case TETRAHEDRON : Current_L = &TS1 ; NbTS++ ; break ;
549 case HEXAHEDRON : Current_L = &TH ; NbTH++ ; break ;
550 case PRISM : Current_L = &TI ; NbTI++ ; break ;
551 case PYRAMID : Current_L = &TY ; NbTY++ ; break ;
552
553 case LINE_2 : Current_L = &TL ; NbTL++ ; break ;
554 case TRIANGLE_2 : Current_L = &TT ; NbTT++ ; break ;
555 case QUADRANGLE_2 : Current_L = &TQ ; NbTQ++ ; break ;
556 case QUADRANGLE_2_8N: Current_L = &TQ ; NbTQ++ ; break ;
557 case TETRAHEDRON_2 : Current_L = &TS1 ; NbTS++ ; break ;
558 case HEXAHEDRON_2 : Current_L = &TH ; NbTH++ ; break ;
559 case HEXAHEDRON_2_20N: Current_L = &TH ; NbTH++ ; break ;
560 case PRISM_2 : Current_L = &TI ; NbTI++ ; break ;
561 case PRISM_2_15N : Current_L = &TI ; NbTI++ ; break ;
562 case PYRAMID_2 : Current_L = &TY ; NbTY++ ; break ;
563 case PYRAMID_2_13N : Current_L = &TY ; NbTY++ ; break ;
564
565 case LINE_3 : Current_L = &TL ; NbTL++ ; break ;
566 case TRIANGLE_3 : Current_L = &TT ; NbTT++ ; break ;
567 case QUADRANGLE_3 : Current_L = &TQ ; NbTQ++ ; break ;
568 case TETRAHEDRON_3 : Current_L = &TS1 ; NbTS++ ; break ;
569 case HEXAHEDRON_3 : Current_L = &TH ; NbTH++ ; break ;
570 case PRISM_3 : Current_L = &TI ; NbTI++ ; break ;
571 case PYRAMID_3 : Current_L = &TY ; NbTY++ ; break ;
572
573 case LINE_4 : Current_L = &TL ; NbTL++ ; break ;
574 case TRIANGLE_4 : Current_L = &TT ; NbTT++ ; break ;
575 case QUADRANGLE_4 : Current_L = &TQ ; NbTQ++ ; break ;
576 case TETRAHEDRON_4 : Current_L = &TS1 ; NbTS++ ; break ;
577 case HEXAHEDRON_4 : Current_L = &TH ; NbTH++ ; break ;
578 case PRISM_4 : Current_L = &TI ; NbTI++ ; break ;
579 //case PYRAMID_4 : Current_L = &TY ; NbTY++ ; break ;
580 }
581 if(Flag_GMSH_VERSION != 2){
582 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(x[i]);
583 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(y[i]);
584 for(i = 0 ; i < NbrNodes ; i++) Current_L->push_back(z[i]);
585 }
586 else{
587 double tmp = ElementNum;
588 Current_L->push_back(tmp);
589 }
590 }
591 if (HarmonicToTime == 1)
592 for(k = 0 ; k < NbHarmonic ; k++){
593 List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
594 for(i = 0 ; i < NbrNodes ; i++){
595 for(j = 0 ; j < 9 ; j++){
596 if(Value[0].Type != TENSOR_DIAG) {
597 if(Value[0].Type == TENSOR_SYM) jj = symIndex[j];
598 else jj = j;
599 Current_L->push_back(Value[i].Val[MAX_DIM*k+jj]);
600 }
601 else {
602 jj = diagIndex[j];
603 if(jj == -1) Current_L->push_back(0.);
604 else Current_L->push_back(Value[i].Val[MAX_DIM*k+jj]);
605 }
606 }
607 }
608 }
609 else
610 for(k = 0 ; k < HarmonicToTime ; k++){
611 List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);
612 for(i = 0 ; i < NbrNodes ; i++){
613 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
614 for(j = 0 ; j < 9 ; j++) {
615 if(Value[0].Type != TENSOR_DIAG) {
616 if(Value[0].Type == TENSOR_SYM) jj = symIndex[j];
617 else jj = j;
618 Current_L->push_back(TmpValue.Val[jj]);
619 }
620 else {
621 jj = diagIndex[j];
622 if(jj == -1) Current_L->push_back(0.);
623 else Current_L->push_back(TmpValue.Val[jj]);
624 }
625 }
626 }
627 }
628
629 }
630
631 // reduce memory requirements by automatically partitioning large
632 // output views into chunks not larger than 1Gb
633 if(Flag_GMSH_VERSION == 2 && TimeStep == NbTimeStep - 1 &&
634 Current_L->size() > (int)(1024 * 1024 * 1024 / sizeof(double))){
635 Format_PostFooter(PSO_P, Store);
636 CurrentPartitionNumber++;
637 Gmsh_StartNewView = 1;
638 }
639
640 }
641
dVecWrite(std::vector<double> & v,FILE * fp,bool binary)642 static void dVecWrite(std::vector<double> &v, FILE *fp, bool binary)
643 {
644 if(v.empty()) return;
645 if(binary)
646 fwrite(&v[0], sizeof(double), v.size(), fp);
647 else
648 for(unsigned i = 0; i < v.size(); i++)
649 fprintf(fp, " %.16g", v[i]);
650 }
651
cVecWrite(std::vector<char> & v,FILE * fp,bool binary)652 static void cVecWrite(std::vector<char> &v, FILE *fp, bool binary)
653 {
654 if(v.empty()) return;
655 if(binary)
656 fwrite(&v[0], sizeof(char), v.size(), fp);
657 else
658 for(unsigned i = 0; i < v.size(); i++)
659 fputc(v[i], fp);
660 }
661
662 /* ------------------------------------------------------------------------ */
663 /* Gnuplot format */
664 /* ------------------------------------------------------------------------ */
665
Gnuplot_PrintElement(int Format,double Time,int TimeStep,int NbrTimeSteps,int NbrHarmonics,int HarmonicToTime,int ElementType,int NumElement,int NbrNodes,double * x,double * y,double * z,double * Dummy,struct Value * Value)666 static void Gnuplot_PrintElement(int Format, double Time, int TimeStep, int NbrTimeSteps,
667 int NbrHarmonics, int HarmonicToTime,
668 int ElementType, int NumElement, int NbrNodes,
669 double *x, double *y, double *z, double *Dummy,
670 struct Value *Value)
671 {
672 static int Size, TmpIndex ;
673 static double * TmpValues ;
674 int i, j, k, t, i2, k2 ;
675 double TimeMH ;
676 struct Value TmpValue ;
677
678 if(!PostStream) return;
679
680 if(TimeStep == 0){
681 switch(Value->Type){
682 case SCALAR : Size = 1 ; break ;
683 case VECTOR : Size = 3 ; break ;
684 case TENSOR_DIAG : Size = 3 ; break ;
685 case TENSOR_SYM : Size = 6 ; break ;
686 case TENSOR : Size = 9 ; break ;
687 }
688 TmpValues = (double*) Malloc(NbrTimeSteps*NbrNodes*NbrHarmonics*Size*sizeof(double));
689 TmpIndex = 0;
690 }
691
692 for(i = 0 ; i < NbrNodes ; i++)
693 for(k = 0 ; k < NbrHarmonics ; k++)
694 for(j = 0 ; j < Size ; j++)
695 TmpValues[TmpIndex++] = Value[i].Val[MAX_DIM*k+j];
696
697 if(TimeStep == NbrTimeSteps-1){
698
699 for(i = 0 ; i <= NbrNodes ; i++){ /* New line for each node, closed loop for tri/qua */
700
701 if(i != NbrNodes)
702 i2 = i ;
703 else{
704 if(NbrNodes < 3) break ;
705 else i2 = 0 ;
706 }
707
708 fprintf(PostStream, "%d %d ", GetDP2Gmsh(ElementType), NumElement);
709 fprintf(PostStream, " %.16g %.16g %.16g ", x[i2], y[i2], z[i2]);
710 if(Dummy){
711 if(Dummy[3]<0){
712 if(!i)
713 fprintf(PostStream, " %.16g %.16g 0 ", Dummy[0], Dummy[2]);
714 else
715 fprintf(PostStream, " %.16g %.16g 0 ", Dummy[1], Dummy[2]);
716 }
717 else
718 fprintf(PostStream, " %.16g %.16g %.16g ", Dummy[0], Dummy[1], Dummy[2]);
719 }
720 else
721 fprintf(PostStream, " 0 0 0 ");
722
723 for(t = 0 ; t < NbrTimeSteps ; t++){
724
725 if (HarmonicToTime == 1) {
726 for(k = 0 ; k < NbrHarmonics ; k++) {
727 for(j = 0 ; j < Size ; j++){
728 fprintf(PostStream, " %.16g",
729 TmpValues[ t*NbrNodes*NbrHarmonics*Size
730 + i2*NbrHarmonics*Size
731 + k*Size
732 + j ]);
733 }
734 fprintf(PostStream, " ");
735 }
736 }
737 else {
738 TmpValue.Type = Value->Type ;
739 for(k = 0 ; k < HarmonicToTime ; k++){
740
741 for(k2 = 0 ; k2 < NbrHarmonics ; k2++)
742 for(j = 0 ; j < Size ; j++)
743 TmpValue.Val[MAX_DIM*k2+j] =
744 TmpValues[ t*NbrNodes*NbrHarmonics*Size
745 + i2*NbrHarmonics*Size
746 + k2*Size
747 + j ] ;
748
749 F_MHToTime0(k, &TmpValue, &TmpValue, k, HarmonicToTime, &TimeMH) ;
750 for(j = 0 ; j < Size ; j++)
751 fprintf(PostStream, "%.16g", TmpValue.Val[0]);
752 fprintf(PostStream, " ");
753 }
754 }
755 fprintf(PostStream, " ");
756
757 } /* for t */
758 fprintf(PostStream, "\n");
759
760 } /* for i */
761 if(NbrNodes > 1) fprintf(PostStream, "\n");
762
763 Free(TmpValues);
764 }
765 }
766
767 /* ------------------------------------------------------------------------ */
768 /* Tabular format */
769 /* ------------------------------------------------------------------------ */
770
771 // global static list for tables output
772 static int TableList_StartNew = 0;
773 static std::list<double> TableList;
774
Tabular_PrintElement(struct PostSubOperation * PSO_P,int Format,double Time,int TimeStep,int NbrTimeSteps,int NbrHarmonics,int HarmonicToTime,int ElementType,int NumElement,int NbrNodes,double * x,double * y,double * z,double * Dummy,struct Value * Value)775 static void Tabular_PrintElement(struct PostSubOperation *PSO_P,
776 int Format, double Time, int TimeStep, int NbrTimeSteps,
777 int NbrHarmonics, int HarmonicToTime,
778 int ElementType, int NumElement, int NbrNodes,
779 double *x, double *y, double *z, double *Dummy,
780 struct Value *Value)
781 {
782 static int Size ;
783 int i,j,k ;
784 double TimeMH ;
785 struct Value TmpValue ;
786
787 if(TableList_StartNew){
788 TableList_StartNew = 0 ;
789 TableList.clear();
790 }
791
792 if(!PostStream) return;
793
794 if(TimeStep == 0){
795 switch(Value->Type){
796 case SCALAR : Size = 1 ; break ;
797 case VECTOR : Size = 3 ; break ;
798 case TENSOR_DIAG : Size = 3 ; break ;
799 case TENSOR_SYM : Size = 6 ; break ;
800 case TENSOR : Size = 9 ; break ;
801 }
802 }
803
804 if(Format == FORMAT_SPACE_TABLE || Format == FORMAT_SIMPLE_SPACE_TABLE
805 || Format == FORMAT_VALUE_ONLY){
806 if(TimeStep == 0){
807 if(Format != FORMAT_SIMPLE_SPACE_TABLE && Format != FORMAT_VALUE_ONLY){
808 fprintf(PostStream, "%d %d ", GetDP2Gmsh(ElementType), NumElement);
809 TableList.push_back(GetDP2Gmsh(ElementType));
810 TableList.push_back(NumElement);
811 }
812 if(Format != FORMAT_VALUE_ONLY)
813 for(i=0 ; i<NbrNodes ; i++){
814 fprintf(PostStream, "%.16g %.16g %.16g ", x[i], y[i], z[i]);
815 TableList.push_back(x[i]);
816 TableList.push_back(y[i]);
817 TableList.push_back(z[i]);
818 }
819 if(Format != FORMAT_SIMPLE_SPACE_TABLE && Format != FORMAT_VALUE_ONLY){
820 if(Dummy){
821 fprintf(PostStream, "%.16g %.16g %.16g ", Dummy[0], Dummy[1], Dummy[2]);
822 TableList.push_back(Dummy[0]);
823 TableList.push_back(Dummy[1]);
824 TableList.push_back(Dummy[2]);
825 }
826 else{
827 fprintf(PostStream, "0 0 0 ");
828 TableList.push_back(0);
829 TableList.push_back(0);
830 TableList.push_back(0);
831 }
832 }
833 }
834 }
835
836 if (HarmonicToTime == 1) {
837 if(Format == FORMAT_TIME_TABLE){
838 fprintf(PostStream, "%d %.16g ", TimeStep, Time);
839 TableList.push_back(TimeStep);
840 TableList.push_back(Time);
841 for(i=0 ; i<NbrNodes ; i++){
842 fprintf(PostStream, " %.16g %.16g %.16g ", x[i], y[i], z[i]);
843 TableList.push_back(x[i]);
844 TableList.push_back(y[i]);
845 TableList.push_back(z[i]);
846 }
847 }
848 for(k = 0 ; k < NbrHarmonics ; k++) {
849 for(i = 0 ; i < NbrNodes ; i++){
850 for(j = 0 ; j < Size ; j++){
851 if (Format != FORMAT_VALUE_ONLY){
852 fprintf(PostStream, " %.16g", Value[i].Val[MAX_DIM*k+j]);
853 TableList.push_back(Value[i].Val[MAX_DIM*k+j]);
854 }
855 else{
856 fprintf(PostStream, " %s_%d_%d = %.16g",
857 PSO_P->ValueName, j, PSO_P->ValueIndex, Value[i].Val[MAX_DIM*k+j]);
858 TableList.push_back(Value[i].Val[MAX_DIM*k+j]);
859 if (j<Size-1)
860 fprintf(PostStream, "\n");
861 }
862 }
863 fprintf(PostStream, " ");
864 }
865 fprintf(PostStream, " ");
866 }
867 }
868 else {
869 for(k = 0 ; k < HarmonicToTime ; k++){
870 for(i = 0 ; i < NbrNodes ; i++){
871 F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
872 if (!i && Format == FORMAT_TIME_TABLE) {
873 fprintf(PostStream, "%d %.16g ", TimeStep, TimeMH);
874 TableList.push_back(TimeStep);
875 TableList.push_back(TimeMH);
876 for(i=0 ; i<NbrNodes ; i++){
877 fprintf(PostStream, " %.16g %.16g %.16g ", x[i], y[i], z[i]);
878 TableList.push_back(x[i]);
879 TableList.push_back(y[i]);
880 TableList.push_back(z[i]);
881 }
882 }
883 for(j = 0 ; j < Size ; j++){
884 fprintf(PostStream, " %.16g", TmpValue.Val[j]) ;
885 TableList.push_back(TmpValue.Val[j]);
886 }
887 fprintf(PostStream, " ");
888 }
889 if(Format == FORMAT_TIME_TABLE)
890 fprintf(PostStream, "\n");
891 }
892 }
893
894 if(TimeStep == NbrTimeSteps-1 || Format == FORMAT_TIME_TABLE)
895 fprintf(PostStream, "\n");
896 }
897
898 /* ------------------------------------------------------------------------ */
899 /* NodeTable format */
900 /* ------------------------------------------------------------------------ */
901
902 // global static map for node table output (cannot be saved incrementally for
903 // each element)
904 static int NodeTable_StartNew = 0;
905 static std::map<int, std::vector<double> > NodeTable;
906
NodeTable_PrintElement(int TimeStep,int NbTimeStep,int NbrHarmonics,struct PostElement * PE)907 static void NodeTable_PrintElement(int TimeStep, int NbTimeStep, int NbrHarmonics,
908 struct PostElement *PE)
909 {
910 if(NodeTable_StartNew){
911 NodeTable_StartNew = 0 ;
912 NodeTable.clear();
913 }
914 for(int i = 0 ; i < PE->NbrNodes ; i++){
915 int n = PE->NumNodes[i];
916 int Size = 0;
917 switch(PE->Value[0].Type){
918 case SCALAR : Size = 1 ; break ;
919 case VECTOR : Size = 3 ; break ;
920 case TENSOR_DIAG : Size = 3 ; break ;
921 case TENSOR_SYM : Size = 6 ; break ;
922 case TENSOR : Size = 9 ; break ;
923 }
924 if(n > 0 && Size){ // we have data on an actual node
925 NodeTable[n].resize(NbTimeStep * NbrHarmonics * Size, 0.);
926 for(int k = 0 ; k < NbrHarmonics ; k++){
927 for(int j = 0 ; j < Size ; j++){
928 double val = PE->Value[i].Val[MAX_DIM * k + j];
929 int idx = NbrHarmonics * Size * TimeStep + k * Size + j;
930 NodeTable[n][idx] = val;
931 }
932 }
933 }
934 }
935
936 }
937
938 /* ------------------------------------------------------------------------ */
939 /* ElementTable format */
940 /* ------------------------------------------------------------------------ */
941
942 // global static map for element table output (cannot be saved incrementally for
943 // each element)
944 static int ElementTable_StartNew = 0;
945 static std::map<int, std::vector<double> > ElementTable;
946
ElementTable_PrintElement(int TimeStep,int NbTimeStep,int NbrHarmonics,struct PostElement * PE)947 static void ElementTable_PrintElement(int TimeStep, int NbTimeStep, int NbrHarmonics,
948 struct PostElement *PE)
949 {
950 if(ElementTable_StartNew){
951 ElementTable_StartNew = 0 ;
952 ElementTable.clear();
953 }
954 int numEle = -1;
955 if(PE->Index >= 0 && PE->Index < Geo_GetNbrGeoElements()){
956 numEle = Geo_GetGeoElement(PE->Index)->Num;
957 int Size = 0;
958 switch(PE->Value[0].Type){
959 case SCALAR : Size = 1 ; break ;
960 case VECTOR : Size = 3 ; break ;
961 case TENSOR_DIAG : Size = 3 ; break ;
962 case TENSOR_SYM : Size = 6 ; break ;
963 case TENSOR : Size = 9 ; break ;
964 }
965 if(Size){
966 ElementTable[numEle].resize
967 (NbTimeStep * PE->NbrNodes * NbrHarmonics * Size, 0.);
968 for(int i = 0 ; i < PE->NbrNodes ; i++){
969 for(int k = 0 ; k < NbrHarmonics ; k++){
970 for(int j = 0 ; j < Size ; j++){
971 double val = PE->Value[i].Val[MAX_DIM * k + j];
972 int idx = TimeStep * PE->NbrNodes * NbrHarmonics * Size +
973 i * NbrHarmonics * Size + k * Size + j;
974 ElementTable[numEle][idx] = val;
975 }
976 }
977 }
978 }
979 }
980 }
981
982 /* ------------------------------------------------------------------------ */
983 /* S t o r e P o s t O p R e s u l t */
984 /* ------------------------------------------------------------------------ */
985
986 static List_T *PostOpResults_L=NULL;
987
StorePostOpResult(int NbrHarmonics,struct PostElement * PE)988 static void StorePostOpResult(int NbrHarmonics, struct PostElement *PE)
989 {
990 int Size;
991 double val;
992
993 if(!PostOpResults_L)
994 PostOpResults_L = List_Create(1000,1000,sizeof(double));
995
996 for(int i = 0 ; i < PE->NbrNodes ; i++){
997 Size = 0;
998 switch(PE->Value[0].Type){
999 case SCALAR : Size = 1 ; break ;
1000 case VECTOR : Size = 3 ; break ;
1001 case TENSOR_DIAG : Size = 3 ; break ;
1002 case TENSOR_SYM : Size = 6 ; break ;
1003 case TENSOR : Size = 9 ; break ;
1004 }
1005 if(Size){ // we have data
1006 for(int k = 0 ; k < NbrHarmonics ; k++){
1007 for(int j = 0 ; j < Size ; j++){
1008 val = PE->Value[i].Val[MAX_DIM * k + j];
1009 List_Add(PostOpResults_L, &val);
1010 }
1011 }
1012 }
1013 }
1014 }
1015
StorePostOpResult(int NbrHarmonics,struct Value * Value)1016 static void StorePostOpResult(int NbrHarmonics, struct Value *Value)
1017 {
1018 int Size;
1019 double val;
1020
1021 if(!PostOpResults_L)
1022 PostOpResults_L = List_Create(1000,1000,sizeof(double));
1023
1024 Size = 0;
1025 switch(Value[0].Type){
1026 case SCALAR : Size = 1 ; break ;
1027 case VECTOR : Size = 3 ; break ;
1028 case TENSOR_DIAG : Size = 3 ; break ;
1029 case TENSOR_SYM : Size = 6 ; break ;
1030 case TENSOR : Size = 9 ; break ;
1031 }
1032 if(Size){ // we have data
1033 for(int k = 0 ; k < NbrHarmonics ; k++){
1034 for(int j = 0 ; j < Size ; j++){
1035 val = Value[0].Val[MAX_DIM * k + j];
1036 List_Add(PostOpResults_L, &val);
1037 }
1038 }
1039 }
1040 }
1041
1042 /* ------------------------------------------------------------------------ */
1043 /* UNV format */
1044 /* ------------------------------------------------------------------------ */
1045
1046 #if !defined(HAVE_NX)
1047 #define NX { Message::Error("UNV export not available in this version"); }
1048 #else
1049 #define NX ;
1050 #endif
1051
1052 double NXUnv_UnitFactor = 1;
1053 int NXUnv_DatasetLocation = 3; //1: Data at nodes, 2: Data on elements, 3: Data at nodes on elements
1054
Unv_PrintHeader(FILE * PostStream,char * name,double Time,int TimeStep,double & NXUnv_UnitFactor,int & NXUnv_DatasetLocation)1055 void Unv_PrintHeader(FILE *PostStream, char *name, double Time, int TimeStep, double& NXUnv_UnitFactor, int& NXUnv_DatasetLocation) NX
1056 void Unv_PrintFooter(FILE *PostStream) NX
1057 void Unv_PrintElement(FILE *PostStream, int Num_Element, int NbrNodes, struct Value *Value, int NbrHarmonics, int& NXUnv_DatasetLocation, double& NXUnv_UnitFactor) NX
1058 void Unv_PrintNodeTable(FILE *PostStream, std::map<int, std::vector<double>> &NodeTable, double& NXUnv_UnitFactor) NX
1059 void Unv_PrintRegion(FILE *PostStream, int Flag_Comma, int numRegion, int NbrHarmonics, int Size, struct Value *Value, double& NXUnv_UnitFactor) NX
1060 #undef NX
1061
1062 /* ------------------------------------------------------------------------ */
1063 /* F o r m a t _ P o s t F o r m a t */
1064 /* ------------------------------------------------------------------------ */
1065
1066 void Format_PostFormat(struct PostSubOperation *PSO_P)
1067 {
1068 if(!PostStream || PSO_P->Type == POP_EXPRESSION) return;
1069
1070 int Format = PSO_P->Format;
1071 int NoMesh = PSO_P->NoMesh;
1072
1073 switch(Format){
1074 case FORMAT_GMSH :
1075 if((PSO_P->StoreInField >= 0 || PSO_P->StoreInMeshBasedField >= 0) &&
1076 !PSO_P->FileOut) break;
1077 if(Flag_GMSH_VERSION == 2){
1078 fprintf(PostStream, "$MeshFormat\n") ;
1079 fprintf(PostStream, "2.2 %d %d\n", Flag_BIN, (int)sizeof(double)) ;
1080 if(Flag_BIN){
1081 int one=1;
1082 fwrite(&one, sizeof(int), 1, PostStream);
1083 fprintf(PostStream, "\n");
1084 }
1085 fprintf(PostStream, "$EndMeshFormat\n") ;
1086 if(!NoMesh){
1087 std::vector<Geo_Element*> elements;
1088 std::set<int> nodes;
1089 if(PSO_P->SubType == PRINT_ONELEMENTSOF){
1090 List_T *Region_L =
1091 ((struct Group *)List_Pointer(Problem_S.Group,
1092 PSO_P->Case.OnRegion.RegionIndex))->InitialList ;
1093 for(int i = 0 ; i < Geo_GetNbrGeoElements() ; i++) {
1094 Geo_Element *Geo_Element = Geo_GetGeoElement(i) ;
1095 if(List_Search(Region_L, &Geo_Element->Region, fcmp_int)){
1096 elements.push_back(Geo_Element);
1097 for (int j = 0 ; j < Geo_Element->NbrNodes ; j++)
1098 nodes.insert(Geo_Element->NumNodes[j]) ;
1099 }
1100 }
1101 }
1102 else{
1103 for(int i = 0 ; i < Geo_GetNbrGeoElements() ; i++) {
1104 Geo_Element *Geo_Element = Geo_GetGeoElement(i) ;
1105 elements.push_back(Geo_Element);
1106 for (int j = 0 ; j < Geo_Element->NbrNodes ; j++)
1107 nodes.insert(Geo_Element->NumNodes[j]) ;
1108 }
1109 }
1110 fprintf(PostStream, "$Nodes\n%d\n", (int)nodes.size());
1111 for (int i = 0 ; i < Geo_GetNbrGeoNodes() ; i++) {
1112 Geo_Node *Geo_Node = Geo_GetGeoNode(i);
1113 if(nodes.find(Geo_Node->Num) != nodes.end()){
1114 if(Flag_BIN){
1115 fwrite(&Geo_Node->Num, sizeof(int), 1, PostStream);
1116 double data[3] = {Geo_Node->x, Geo_Node->y, Geo_Node->z};
1117 fwrite(data, sizeof(double), 3, PostStream);
1118 }
1119 else{
1120 fprintf(PostStream, "%d %.16g %.16g %.16g\n",
1121 Geo_Node->Num, Geo_Node->x, Geo_Node->y, Geo_Node->z) ;
1122 }
1123 }
1124 }
1125 fprintf(PostStream, "$EndNodes\n$Elements\n%d\n", (int)elements.size());
1126 for (unsigned int i = 0 ; i < elements.size() ; i++) {
1127 Geo_Element *Geo_Element = elements[i];
1128 int Type = GetDP2Gmsh(Geo_Element->Type) ;
1129 if(Flag_BIN){
1130 int blob[6] = {Type, 1, 2, Geo_Element->Num, Geo_Element->Region,
1131 Geo_Element->ElementaryRegion};
1132 fwrite(blob, sizeof(int), 6, PostStream);
1133 std::vector<int> verts(Geo_Element->NbrNodes);
1134 for (int j = 0 ; j < Geo_Element->NbrNodes ; j++)
1135 verts[j] = Geo_Element->NumNodes[j] ;
1136 fwrite(&verts[0], sizeof(int), Geo_Element->NbrNodes, PostStream);
1137 }
1138 else{
1139 fprintf(PostStream, "%d %d 2 %d %d ", Geo_Element->Num,
1140 Type, Geo_Element->Region, Geo_Element->ElementaryRegion) ;
1141 for (int j = 0 ; j < Geo_Element->NbrNodes ; j++)
1142 fprintf(PostStream, "%d ", Geo_Element->NumNodes[j]) ;
1143 fprintf(PostStream, "\n") ;
1144 }
1145 }
1146 fprintf(PostStream, "$EndElements\n");
1147 }
1148 }
1149 else if(PostStream && Flag_BIN){/* bricolage */
1150 fprintf(PostStream, "$PostFormat /* Gmsh 1.2, %s */\n",
1151 Flag_BIN ? "binary" : "ascii") ;
1152 fprintf(PostStream, "1.2 %d %d\n", Flag_BIN, (int)sizeof(double)) ;
1153 fprintf(PostStream, "$EndPostFormat\n") ;
1154 }
1155 break ;
1156 case FORMAT_GNUPLOT :
1157 fprintf(PostStream, "# GetDP %s, %s\n", GETDP_VERSION,
1158 Flag_BIN ? "binary" : "ascii") ;
1159 break ;
1160 }
1161 }
1162
1163 /* ------------------------------------------------------------------------ */
1164 /* F o r m a t _ P o s t H e a d e r */
1165 /* ------------------------------------------------------------------------ */
1166
Format_PostHeader(struct PostSubOperation * PSO_P,int NbTimeStep,int Order,char * Name1,char * Name2)1167 void Format_PostHeader(struct PostSubOperation *PSO_P, int NbTimeStep,
1168 int Order, char *Name1, char *Name2)
1169 {
1170 int Format = PSO_P->Format;
1171 double Time = Current.Time;
1172 int TimeStep = Current.TimeStep;
1173 int Contour = PSO_P->Iso;
1174 int Type = PSO_P->CombinationType;
1175
1176 char name[256] ;
1177
1178 CurrentPartitionNumber = 0;
1179
1180 if(Contour){
1181 if(!PostElement_L)
1182 PostElement_L = List_Create(20, 20, sizeof(struct PostElement*));
1183 else
1184 List_Reset(PostElement_L);
1185 }
1186
1187 if(Name1 && Name2) {
1188 strcpy(name, Order ? Name1 : Name2) ;
1189 strcat(name, Get_StringForDefine(PostSubOperation_CombinationType, Type)) ;
1190 strcat(name, Order ? Name2 : Name1) ;
1191 }
1192 else if(Name1)
1193 strcpy(name, Name1) ;
1194 else if(Name2)
1195 strcpy(name, Name2) ;
1196 else
1197 strcpy(name, "unnamed");
1198
1199 strcpy(CurrentName, name);
1200
1201 switch(Format){
1202 case FORMAT_GMSH_PARSED :
1203 if(PostStream) fprintf(PostStream, "View \"%s\" {\n", name) ;
1204 Gmsh_StartNewView = 1 ;
1205 break ;
1206 case FORMAT_GMSH :
1207 Gmsh_StartNewView = 1 ;
1208 if((PSO_P->StoreInField >= 0 || PSO_P->StoreInMeshBasedField >= 0) &&
1209 !PSO_P->FileOut) break;
1210 if(PostStream && Flag_GMSH_VERSION != 2){
1211 if(Flag_BIN){ /* bricolage */
1212 fprintf(PostStream, "$View /* %s */\n", name);
1213 fprintf(PostStream, "%s ", name);
1214 }
1215 else {
1216 fprintf(PostStream, "View \"%s\" {\n", name) ;
1217 }
1218 }
1219 break ;
1220 case FORMAT_NXUNV :
1221 if(PostStream){
1222 NodeTable_StartNew = 1 ;
1223 Unv_PrintHeader(PostStream, name, Time, TimeStep, NXUnv_UnitFactor, NXUnv_DatasetLocation);
1224 }
1225 break ;
1226 case FORMAT_GNUPLOT :
1227 if(PostStream){
1228 fprintf(PostStream, "# PostData '%s'\n", name);
1229 fprintf(PostStream, "# Type Num X Y Z N1 N2 N3 Values <Values>...\n");
1230 }
1231 break ;
1232 case FORMAT_NODE_TABLE :
1233 NodeTable_StartNew = 1 ;
1234 break ;
1235 case FORMAT_ELEMENT_TABLE :
1236 ElementTable_StartNew = 1 ;
1237 break ;
1238 case FORMAT_SPACE_TABLE :
1239 case FORMAT_TIME_TABLE :
1240 case FORMAT_SIMPLE_SPACE_TABLE :
1241 case FORMAT_VALUE_ONLY :
1242 TableList_StartNew = 1 ;
1243 break ;
1244 case FORMAT_ADAPT :
1245 if(PostStream) fprintf(PostStream, "$Adapt /* %s */\n", name) ;
1246 break ;
1247 }
1248 }
1249
1250 /* ------------------------------------------------------------------------ */
1251 /* F o r m a t _ P o s t F o o t e r */
1252 /* ------------------------------------------------------------------------ */
1253
Format_PostFooter(struct PostSubOperation * PSO_P,int Store,bool SendToServer)1254 void Format_PostFooter(struct PostSubOperation *PSO_P, int Store,
1255 bool SendToServer)
1256 {
1257 List_T *Iso_L[NBR_MAX_ISO], *Solutions_L;
1258 double IsoMin = 1.e200, IsoMax = -1.e200, IsoVal = 0.0, freq, valr, vali ;
1259 int NbrIso = 0 ;
1260 int iPost, iNode, iIso, iTime, One=1, i, j, NbTimeStep ;
1261 char tmp[1024];
1262 bool PostOpSolutionGenerated;
1263 struct PostElement *PE ;
1264 struct Solution *Solution_P=NULL, Solution_S;
1265
1266
1267 if( !(NbTimeStep = List_Nbr(PSO_P->TimeStep_L)) )
1268 NbTimeStep = List_Nbr(Current.DofData->Solutions);
1269
1270 if ( (PSO_P->Format == FORMAT_GMSH || PSO_P->Format == FORMAT_GMSH_PARSED) &&
1271 Flag_GMSH_VERSION != 2 ){
1272
1273 switch(PSO_P->Legend){
1274
1275 case LEGEND_TIME:
1276 Gmsh_StringStart(PSO_P->Format, PSO_P->LegendPosition[0],
1277 PSO_P->LegendPosition[1], PSO_P->LegendPosition[2]);
1278 for (i = 0 ; i < NbTimeStep ; i++) {
1279 Pos_InitAllSolutions(PSO_P->TimeStep_L, i) ;
1280 valr = Current.DofData->CurrentSolution->Time ;
1281 for (j = 0 ; j < Current.NbrHar ; j++){
1282 sprintf(tmp, "Step %d/%d: Time = %g", i+1, NbTimeStep, valr);
1283 Gmsh_StringAdd(PSO_P->Format, (!i && !j), tmp);
1284 }
1285 }
1286 Gmsh_StringEnd(PSO_P->Format);
1287 break;
1288
1289 case LEGEND_FREQUENCY:
1290 if(Current.NbrHar > 1) {
1291 Gmsh_StringStart(PSO_P->Format, PSO_P->LegendPosition[0],
1292 PSO_P->LegendPosition[1], PSO_P->LegendPosition[2]);
1293 for (i = 0 ; i < NbTimeStep ; i++) {
1294 Pos_InitAllSolutions(PSO_P->TimeStep_L, i) ;
1295 for (j = 0 ; j < Current.NbrHar ; j+=2) {
1296 freq = 0.5*Current.DofData->Val_Pulsation[j/2]/M_PI ;
1297 sprintf(tmp, "%g Hz (real part: COSINUS)", freq);
1298 Gmsh_StringAdd(PSO_P->Format, (!i && !j), tmp);
1299 sprintf(tmp, "%g Hz (imag part: -SINUS)", freq);
1300 Gmsh_StringAdd(PSO_P->Format, 0, tmp);
1301 }
1302 }
1303 Gmsh_StringEnd(PSO_P->Format);
1304 }
1305 break;
1306
1307 case LEGEND_EIGENVALUES:
1308 Gmsh_StringStart(PSO_P->Format, PSO_P->LegendPosition[0],
1309 PSO_P->LegendPosition[1], PSO_P->LegendPosition[2]);
1310 for (i = 0 ; i < NbTimeStep ; i++) {
1311 Pos_InitAllSolutions(PSO_P->TimeStep_L, i) ;
1312 valr = Current.DofData->CurrentSolution->Time ;
1313 vali = Current.DofData->CurrentSolution->TimeImag ;
1314 if(Current.NbrHar == 1){
1315 sprintf(tmp, "Eigenvalue %d/%d: %g",
1316 i+1, NbTimeStep, valr);
1317 Gmsh_StringAdd(PSO_P->Format, !i, tmp);
1318 }
1319 else{
1320 for (j = 0 ; j < Current.NbrHar ; j++) {
1321 if(!(j % 2)){
1322 sprintf(tmp, "Eigenvalue %d/%d: %g %s i * %g (Real Part)",
1323 i+1, NbTimeStep, valr, (vali > 0) ? "+" : "-",
1324 (vali > 0) ? vali : -vali);
1325 Gmsh_StringAdd(PSO_P->Format, (!i && !j), tmp);
1326 }
1327 else{
1328 sprintf(tmp, "Eigenvalue %d/%d: %g %s i * %g (Imaginary Part)",
1329 i+1, NbTimeStep, valr, (vali > 0) ? "+" : "-",
1330 (vali > 0) ? vali : -vali);
1331 Gmsh_StringAdd(PSO_P->Format, 0, tmp);
1332 }
1333 }
1334 }
1335 }
1336 Gmsh_StringEnd(PSO_P->Format);
1337 break;
1338 }
1339 }
1340
1341 if(PSO_P->Iso){
1342
1343 for(iPost = 0 ; iPost < List_Nbr(PostElement_L) ; iPost++){
1344 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
1345 for (iNode = 0 ; iNode < PE->NbrNodes ; iNode++ ){
1346 IsoMin = std::min(IsoMin, PE->Value[iNode].Val[0]) ;
1347 IsoMax = std::max(IsoMax, PE->Value[iNode].Val[0]) ;
1348 }
1349 }
1350
1351 if((NbrIso = PSO_P->Iso) < 0)
1352 NbrIso = List_Nbr(PSO_P->Iso_L) ;
1353 if(NbrIso > NBR_MAX_ISO){
1354 Message::Error("Too many Iso values");
1355 NbrIso = NBR_MAX_ISO;
1356 }
1357
1358 if(PostStream && PSO_P->Format == FORMAT_GNUPLOT)
1359 fprintf(PostStream, "# NbIso = %d, Min = %g, Max = %g\n",
1360 NbrIso, IsoMin, IsoMax) ;
1361
1362 for(iIso = 0 ; iIso < NbrIso ; iIso++)
1363 Iso_L[iIso] = List_Create(10, 10, sizeof(struct PostElement*)) ;
1364
1365 for(iPost = 0 ; iPost < List_Nbr(PostElement_L) ; iPost++){
1366 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
1367 for(iIso = 0 ; iIso < NbrIso ; iIso++){
1368 if(PSO_P->Iso > 0){
1369 Cal_Iso(PE, Iso_L[iIso], IsoMin+iIso*(IsoMax-IsoMin)/(double)(NbrIso-1),
1370 IsoMin, IsoMax, PSO_P->DecomposeInSimplex) ;
1371 }
1372 else{
1373 List_Read(PSO_P->Iso_L, iIso, &IsoVal) ;
1374 Cal_Iso(PE, Iso_L[iIso], IsoVal, IsoMin, IsoMax, PSO_P->DecomposeInSimplex) ;
1375 }
1376 }
1377 if(!Store) Destroy_PostElement(PE);
1378 }
1379
1380 for(iIso = 0 ; iIso < NbrIso ; iIso++){
1381 for(iPost = 0 ; iPost < List_Nbr(Iso_L[iIso]) ; iPost++){
1382 PE = *(struct PostElement**)List_Pointer(Iso_L[iIso], iPost) ;
1383 Format_PostElement(PSO_P, 0, 0,
1384 Current.Time, 0, 1,
1385 Current.NbrHar, PSO_P->HarmonicToTime,
1386 NULL, PE);
1387 Destroy_PostElement(PE) ;
1388 }
1389 List_Delete(Iso_L[iIso]) ;
1390 if(PostStream && PSO_P->Format == FORMAT_GNUPLOT) fprintf(PostStream, "\n") ;
1391 }
1392 }
1393
1394 switch(PSO_P->Format){
1395 case FORMAT_GMSH_PARSED :
1396 if(PostStream && List_Nbr(TimeValue_L) > 1){
1397 fprintf(PostStream, "TIME{");
1398 for(iTime = 0; iTime < List_Nbr(TimeValue_L); iTime++){
1399 if(iTime) fprintf(PostStream, ",");
1400 fprintf(PostStream, "%.16g", *(double*)List_Pointer(TimeValue_L, iTime));
1401 }
1402 fprintf(PostStream, "};\n");
1403 }
1404 fprintf(PostStream, "};\n") ;
1405 break ;
1406 case FORMAT_GMSH :
1407 if(Gmsh_StartNewView) Gmsh_ResetStaticLists(); // nothing to print!
1408 if(PSO_P->StoreInField >= 0 || PSO_P->StoreInMeshBasedField >= 0){
1409 #if defined(HAVE_GMSH)
1410 int field = (PSO_P->StoreInField >= 0) ? PSO_P->StoreInField :
1411 PSO_P->StoreInMeshBasedField;
1412 Message::Info("Storing data in field %d (%s)", field,
1413 PSO_P->StoreInField >= 0 ? "list-based" : "mesh-based");
1414 int NS[24] = {NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT,
1415 NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH,
1416 NbSI, NbVI, NbTI, NbSY, NbVY, NbTY};
1417 std::vector<double> *LS[24] = {&SP, &VP, &TP, &SL, &VL, &TL, &ST, &VT, &TT,
1418 &SQ, &VQ, &TQ, &SS, &VS, &TS1, &SH, &VH, &TH,
1419 &SI, &VI, &TI, &SY, &VY, &TY};
1420 PViewData *data;
1421 if(PSO_P->StoreInField >= 0)
1422 data = new PViewDataList();
1423 else
1424 data = new PViewDataGModel(PViewDataGModel::ElementNodeData);
1425 data->importLists(NS, LS);
1426 new PView(data, field);
1427 #else
1428 Message::Error("GetDP must be compiled with Gmsh support to store data as field");
1429 #endif
1430 if(!PSO_P->FileOut) break;
1431 }
1432 if(Flag_GMSH_VERSION == 2){
1433 int NS[8] = {NbSP, NbSL, NbST, NbSQ, NbSS, NbSH, NbSI, NbSY};
1434 std::vector<double> *LS[8] = {&SP, &SL, &ST, &SQ, &SS, &SH, &SI, &SY};
1435 Gmsh_PrintElementNodeData(PSO_P, NbTimeStep, 1, NS, LS);
1436
1437 int NV[8] = {NbVP, NbVL, NbVT, NbVQ, NbVS, NbVH, NbVI, NbVY};
1438 std::vector<double> *LV[8] = {&VP, &VL, &VT, &VQ, &VS, &VH, &VI, &VY};
1439 Gmsh_PrintElementNodeData(PSO_P, NbTimeStep, 3, NV, LV);
1440
1441 int NT[8] = {NbTP, NbTL, NbTT, NbTQ, NbTS, NbTH, NbTI, NbTY};
1442 std::vector<double> *LT[8] = {&TP, &TL, &TT, &TQ, &TS1, &TH, &TI, &TY};
1443 Gmsh_PrintElementNodeData(PSO_P, NbTimeStep, 9, NT, LT);
1444 }
1445 else if(PostStream && Flag_BIN){ /* bricolage */
1446 fprintf(PostStream, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
1447 "%d %d %d %d %d %d %d %d %d %d %d 0 0\n",
1448 List_Nbr(TimeValue_L),
1449 NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT,
1450 NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH,
1451 NbSI, NbVI, NbTI, NbSY, NbVY, NbTY, NbT2, (int)T2C.size());
1452 fwrite(&One, sizeof(int), 1, PostStream);
1453 List_WriteToFile(TimeValue_L, PostStream, LIST_FORMAT_BINARY);
1454 bool f = true;
1455 dVecWrite(SP, PostStream, f); dVecWrite(VP, PostStream, f);
1456 dVecWrite(TP, PostStream, f);
1457 dVecWrite(SL, PostStream, f); dVecWrite(VL, PostStream, f);
1458 dVecWrite(TL, PostStream, f);
1459 dVecWrite(ST, PostStream, f); dVecWrite(VT, PostStream, f);
1460 dVecWrite(TT, PostStream, f);
1461 dVecWrite(SQ, PostStream, f); dVecWrite(VQ, PostStream, f);
1462 dVecWrite(TQ, PostStream, f);
1463 dVecWrite(SS, PostStream, f); dVecWrite(VS, PostStream, f);
1464 dVecWrite(TS1, PostStream, f);
1465 dVecWrite(SH, PostStream, f); dVecWrite(VH, PostStream, f);
1466 dVecWrite(TH, PostStream, f);
1467 dVecWrite(SI, PostStream, f); dVecWrite(VI, PostStream, f);
1468 dVecWrite(TI, PostStream, f);
1469 dVecWrite(SY, PostStream, f); dVecWrite(VY, PostStream, f);
1470 dVecWrite(TY, PostStream, f);
1471 dVecWrite(T2D, PostStream, f); cVecWrite(T2C, PostStream, f);
1472 fprintf(PostStream, "\n");
1473 fprintf(PostStream, "$EndView\n");
1474 }
1475 else if(PostStream){
1476 if(List_Nbr(TimeValue_L) > 1){
1477 fprintf(PostStream, "TIME{");
1478 for(iTime = 0; iTime < List_Nbr(TimeValue_L); iTime++){
1479 if(iTime) fprintf(PostStream, ",");
1480 fprintf(PostStream, "%.16g", *(double*)List_Pointer(TimeValue_L, iTime));
1481 }
1482 fprintf(PostStream, "};\n");
1483 }
1484 fprintf(PostStream, "};\n") ;
1485 }
1486 break ;
1487 case FORMAT_ADAPT :
1488 if(PostStream) fprintf(PostStream, "$EndAdapt\n");
1489 break ;
1490 case FORMAT_NXUNV :
1491 if(PostStream){
1492 if(NXUnv_DatasetLocation == 1) //Data at nodes
1493 Unv_PrintNodeTable(PostStream, NodeTable, NXUnv_UnitFactor);
1494 Unv_PrintFooter(PostStream);
1495 }
1496 break ;
1497 case FORMAT_NODE_TABLE :
1498 if(PostStream && NodeTable.size()){
1499 fprintf(PostStream, "%d\n", (int)NodeTable.size());
1500 for(std::map<int, std::vector<double> >::iterator it = NodeTable.begin();
1501 it != NodeTable.end(); it++){
1502 fprintf(PostStream, "%d", it->first);
1503 for(unsigned int i = 0; i < it->second.size(); i++)
1504 fprintf(PostStream, " %.16g", it->second[i]);
1505 fprintf(PostStream, "\n");
1506 }
1507 }
1508 {
1509 std::vector<double> exp;
1510 exp.push_back(NodeTable.size());
1511 for(std::map<int, std::vector<double> >::iterator it = NodeTable.begin();
1512 it != NodeTable.end(); it++){
1513 exp.push_back(it->first);
1514 for(unsigned int i = 0; i < it->second.size(); i++)
1515 exp.push_back(it->second[i]);
1516 }
1517 GetDPNumbers[CurrentName] = exp;
1518 GetDPNumbersMap[CurrentName] = NodeTable;
1519 if(SendToServer && PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No"))
1520 Message::AddOnelabNumberChoice(PSO_P->SendToServer, exp, PSO_P->Color,
1521 PSO_P->Units, PSO_P->Label, PSO_P->Visible,
1522 PSO_P->Closed);
1523 }
1524 break;
1525 case FORMAT_ELEMENT_TABLE :
1526 if(PostStream){
1527 fprintf(PostStream, "%d\n", (int)ElementTable.size());
1528 for(std::map<int, std::vector<double> >::iterator it = ElementTable.begin();
1529 it != ElementTable.end(); it++){
1530 fprintf(PostStream, "%d", it->first);
1531 for(unsigned int i = 0; i < it->second.size(); i++)
1532 fprintf(PostStream, " %.16g", it->second[i]);
1533 fprintf(PostStream, "\n");
1534 }
1535 }
1536 {
1537 std::vector<double> exp;
1538 exp.push_back(ElementTable.size());
1539 for(std::map<int, std::vector<double> >::iterator it = ElementTable.begin();
1540 it != ElementTable.end(); it++){
1541 exp.push_back(it->first);
1542 for(unsigned int i = 0; i < it->second.size(); i++)
1543 exp.push_back(it->second[i]);
1544 }
1545 GetDPNumbers[CurrentName] = exp;
1546 GetDPNumbersMap[CurrentName] = ElementTable;
1547 if(SendToServer && PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No"))
1548 Message::AddOnelabNumberChoice(PSO_P->SendToServer, exp, PSO_P->Color,
1549 PSO_P->Units, PSO_P->Label, PSO_P->Visible,
1550 PSO_P->Closed);
1551 }
1552 break;
1553 case FORMAT_SPACE_TABLE :
1554 case FORMAT_TIME_TABLE :
1555 case FORMAT_SIMPLE_SPACE_TABLE :
1556 case FORMAT_VALUE_ONLY :
1557 {
1558 if(TableList.size()) {
1559 std::vector<double> v(TableList.begin(), TableList.end());
1560 GetDPNumbers[CurrentName] = v;
1561 if(SendToServer && PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No"))
1562 Message::AddOnelabNumberChoice(PSO_P->SendToServer, v, PSO_P->Color,
1563 PSO_P->Units, PSO_P->Label, PSO_P->Visible,
1564 PSO_P->Closed);
1565 }
1566 }
1567 break;
1568 case FORMAT_LOOP_ERROR :
1569 Solutions_L = ((struct PostOpSolutions*)
1570 List_Pointer(Current.PostOpData_L, Current.PostOpDataIndex))->Solutions_L;
1571 PostOpSolutionGenerated = false;
1572 if(List_Nbr(Solutions_L)>0) {
1573 Solution_P = (struct Solution*)List_Pointer(Solutions_L, List_Nbr(Solutions_L)-1);
1574 PostOpSolutionGenerated = (Solution_P->TimeStep == (int)Current.TimeStep);
1575 }
1576 if (!PostOpSolutionGenerated) {
1577 Solution_S.Time = Current.Time;
1578 Solution_S.TimeImag = Current.TimeImag;
1579 Solution_S.TimeStep = Current.TimeStep;
1580 Solution_S.SolutionExist = 1;
1581 Solution_S.TimeFunctionValues = NULL;
1582 LinAlg_CreateVector(&Solution_S.x, &Current.DofData->Solver,
1583 List_Nbr(PostOpResults_L));
1584 for(int i=0; i<List_Nbr(PostOpResults_L); i++){
1585 List_Read(PostOpResults_L, i, &valr);
1586 LinAlg_SetDoubleInVector(valr, &Solution_S.x, i);
1587 }
1588 List_Add(Solutions_L, &Solution_S);
1589 }
1590 else
1591 for(int i=0; i<List_Nbr(PostOpResults_L); i++){
1592 List_Read(PostOpResults_L, i, &valr);
1593 LinAlg_SetDoubleInVector(valr, &Solution_P->x, i);
1594 }
1595
1596 List_Delete(PostOpResults_L);
1597 PostOpResults_L = NULL;
1598 break;
1599 }
1600 }
1601
1602 /* ------------------------------------------------------------------------ */
1603 /* F o r m a t _ P o s t E l e m e n t */
1604 /* ------------------------------------------------------------------------ */
1605
Format_PostElement(struct PostSubOperation * PSO_P,int Contour,int Store,double Time,int TimeStep,int NbTimeStep,int NbrHarmonics,int HarmonicToTime,double * Dummy,struct PostElement * PE)1606 void Format_PostElement(struct PostSubOperation *PSO_P, int Contour, int Store,
1607 double Time, int TimeStep, int NbTimeStep,
1608 int NbrHarmonics, int HarmonicToTime, double *Dummy,
1609 struct PostElement * PE)
1610 {
1611 int i, j, k, l, Num_Element ;
1612 struct PostElement * PE2 ;
1613 struct Value Value ;
1614
1615 static int Warning_FirstHarmonic = 0 ;
1616
1617 /* TODO
1618
1619 static int Size = 0 ;
1620
1621 int flag_storeAllTimeResults, indexInTmpValues;
1622 static struct Value TmpValue, *TmpValues ;
1623 static double *Times ;
1624 struct Value *FourierValues;
1625
1626 flag_storeAllTimeResults = PSO_P->TimeToHarmonic ;
1627 indexInTmpValues = flag_storeAllTimeResults? iTime * NbrRegion : 0 ;
1628
1629 if(1){
1630 switch(PE->Value[0].Type){
1631 case SCALAR : Size = 1 ; break ;
1632 case VECTOR : Size = 3 ; break ;
1633 case TENSOR_DIAG : Size = 3 ; break ;
1634 case TENSOR_SYM : Size = 6 ; break ;
1635 case TENSOR : Size = 9 ; break ;
1636 default : Size = 9 ; break ;
1637 }
1638 }
1639 */
1640
1641 if(PE->Index != NO_ELEMENT)
1642 Num_Element = Geo_GetGeoElement(PE->Index)->Num ;
1643 else
1644 Num_Element = 0 ;
1645
1646 if(Contour){
1647 if(PE->Value[0].Type != SCALAR){
1648 Message::Error("Non scalar Element %d in contour creation", Num_Element);
1649 return;
1650 }
1651 if(NbTimeStep != 1){
1652 Message::Error("Contour creation not allowed for multiple time steps");
1653 return;
1654 }
1655 if(Current.NbrHar != 1 && !Warning_FirstHarmonic){
1656 Message::Warning("Contour creation done only for first harmonic (use Re[] or Im[])");
1657 Warning_FirstHarmonic = 1 ;
1658 }
1659 if(Store)
1660 List_Add(PostElement_L, &PE) ;
1661 else{
1662 PE2 = PartialCopy_PostElement(PE) ;
1663 List_Add(PostElement_L, &PE2) ;
1664 }
1665 return ;
1666 }
1667
1668 if(PSO_P->ChangeOfCoordinates[0] >= 0){
1669 for(i=0 ; i<PE->NbrNodes ; i++){
1670 Current.x = PE->x[i];
1671 Current.y = PE->y[i];
1672 Current.z = PE->z[i];
1673 for(j = 0; j<9 ; j++) Current.Val[j] = PE->Value[i].Val[j];
1674 Get_ValueOfExpressionByIndex(PSO_P->ChangeOfCoordinates[0], NULL, 0., 0., 0., &Value) ;
1675 PE->x[i] = Value.Val[0];
1676 Get_ValueOfExpressionByIndex(PSO_P->ChangeOfCoordinates[1], NULL, 0., 0., 0., &Value) ;
1677 PE->y[i] = Value.Val[0];
1678 Get_ValueOfExpressionByIndex(PSO_P->ChangeOfCoordinates[2], NULL, 0., 0., 0., &Value) ;
1679 PE->z[i] = Value.Val[0];
1680 }
1681 }
1682
1683 if(PSO_P->ChangeOfValues && List_Nbr(PSO_P->ChangeOfValues) > 0){
1684 for(i=0 ; i < PE->NbrNodes ; i++){
1685 Current.x = PE->x[i];
1686 Current.y = PE->y[i];
1687 Current.z = PE->z[i];
1688 for(k=0 ; k<Current.NbrHar ; k++){
1689 for(j = 0; j<9 ; j++) Current.Val[j] = PE->Value[i].Val[MAX_DIM*k+j];
1690 for(l=0 ; l<List_Nbr(PSO_P->ChangeOfValues) ; l++){
1691 Get_ValueOfExpressionByIndex(*(int*)List_Pointer(PSO_P->ChangeOfValues, l),
1692 NULL, 0., 0., 0., &Value) ;
1693 PE->Value[i].Val[MAX_DIM*k+l] = Value.Val[0];
1694 }
1695 }
1696 }
1697 }
1698
1699 switch(PSO_P->Format){
1700 case FORMAT_GMSH_PARSED :
1701 GmshParsed_PrintElement(Time, TimeStep, NbTimeStep, NbrHarmonics, HarmonicToTime,
1702 PE->Type, PE->NbrNodes, PE->x, PE->y, PE->z,
1703 PE->Value) ;
1704 break ;
1705 case FORMAT_NXUNV :
1706 if(PostStream){
1707 if(NXUnv_DatasetLocation == 1) //Data at nodes
1708 NodeTable_PrintElement(TimeStep, NbTimeStep, NbrHarmonics, PE);
1709 else if(NXUnv_DatasetLocation == 2 || NXUnv_DatasetLocation == 3) //Data at elements or nodes on elements
1710 Unv_PrintElement(PostStream, Num_Element, PE->NbrNodes, PE->Value, NbrHarmonics, NXUnv_DatasetLocation, NXUnv_UnitFactor) ;
1711 }
1712 break ;
1713 case FORMAT_GMSH :
1714 if(PSO_P->StoreInField >= 0 || PSO_P->StoreInMeshBasedField >= 0){
1715 Gmsh_PrintElement(Time, TimeStep, NbTimeStep, NbrHarmonics, HarmonicToTime,
1716 PE->Type, Num_Element, PE->NbrNodes, PE->x, PE->y, PE->z,
1717 PE->Value, PSO_P, Store) ;
1718 if(!PSO_P->FileOut || Flag_GMSH_VERSION == 2 || Flag_BIN) break;
1719 }
1720 if(Flag_GMSH_VERSION == 2 || Flag_BIN){ /* bricolage */
1721 Gmsh_PrintElement(Time, TimeStep, NbTimeStep, NbrHarmonics, HarmonicToTime,
1722 PE->Type, Num_Element, PE->NbrNodes, PE->x, PE->y, PE->z,
1723 PE->Value, PSO_P, Store) ;
1724 }
1725 else{
1726 GmshParsed_PrintElement(Time, TimeStep, NbTimeStep, NbrHarmonics, HarmonicToTime,
1727 PE->Type, PE->NbrNodes, PE->x, PE->y, PE->z,
1728 PE->Value) ;
1729 }
1730 break ;
1731 case FORMAT_GNUPLOT :
1732 Gnuplot_PrintElement(PSO_P->Format, Time, TimeStep, NbTimeStep, NbrHarmonics,
1733 HarmonicToTime, PE->Type, Num_Element, PE->NbrNodes,
1734 PE->x, PE->y, PE->z, Dummy, PE->Value) ;
1735 break ;
1736 case FORMAT_SPACE_TABLE :
1737 case FORMAT_TIME_TABLE :
1738 case FORMAT_SIMPLE_SPACE_TABLE :
1739 case FORMAT_VALUE_ONLY :
1740 Tabular_PrintElement(PSO_P, PSO_P->Format, Time, TimeStep, NbTimeStep,
1741 NbrHarmonics, HarmonicToTime, PE->Type, Num_Element,
1742 PE->NbrNodes, PE->x, PE->y, PE->z, Dummy, PE->Value) ;
1743 break ;
1744 case FORMAT_NODE_TABLE :
1745 NodeTable_PrintElement(TimeStep, NbTimeStep, NbrHarmonics, PE);
1746 break;
1747 case FORMAT_ELEMENT_TABLE :
1748 ElementTable_PrintElement(TimeStep, NbTimeStep, NbrHarmonics, PE);
1749 break;
1750 case FORMAT_LOOP_ERROR :
1751 StorePostOpResult(NbrHarmonics, PE);
1752 break;
1753 case FORMAT_ADAPT:
1754 if(PostStream){
1755 if(Dummy[4]) fprintf(PostStream, "%d\n", (int)Dummy[4]) ;
1756 fprintf(PostStream, "%d %g %g %g\n",
1757 (int)Dummy[0], Dummy[1], Dummy[2], Dummy[3]);
1758 }
1759 break ;
1760 default :
1761 Message::Error("Unknown format in Format_PostElement");
1762 }
1763
1764 if (PE->NbrNodes == 1 &&
1765 PSO_P->Format != FORMAT_NODE_TABLE &&
1766 PSO_P->Format != FORMAT_ELEMENT_TABLE){
1767 if(PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No")){
1768 std::vector<double> v;
1769 Export_Value(&PE->Value[0], v, PSO_P->SendToServerList);
1770 Message::AddOnelabNumberChoice(PSO_P->SendToServer, v, PSO_P->Color,
1771 PSO_P->Units, PSO_P->Label, PSO_P->Visible,
1772 PSO_P->Closed);
1773 }
1774 }
1775 }
1776
1777 /* ------------------------------------------------------------------------ */
1778 /* F o r m a t _ P o s t V a l u e */
1779 /* ------------------------------------------------------------------------ */
1780
Format_PostValue(struct PostQuantity * PQ_P,struct PostSubOperation * PSO_P,int Format,char * Comma,int Group_FunctionType,int iTime,double Time,int NbrTimeStep,int iRegion,int numRegion,int NbrRegion,int NbrHarmonics,int HarmonicToTime,int FourierTransform,int Flag_NoNewLine,struct Value * Value)1781 void Format_PostValue(struct PostQuantity *PQ_P,
1782 struct PostSubOperation *PSO_P,
1783 int Format, char *Comma, int Group_FunctionType,
1784 int iTime, double Time, int NbrTimeStep,
1785 int iRegion, int numRegion, int NbrRegion,
1786 int NbrHarmonics, int HarmonicToTime, int FourierTransform,
1787 int Flag_NoNewLine,
1788 struct Value * Value)
1789 {
1790 static int Size ;
1791 int j, k ;
1792 double TimeMH, Freq ;
1793 double x, y, z ;
1794 int flag_storeAllTimeResults, indexInTmpValues ;
1795 static struct Value TmpValue, *TmpValues ;
1796 static double *Times ;
1797
1798 flag_storeAllTimeResults = FourierTransform || PSO_P->TimeToHarmonic ;
1799 indexInTmpValues = flag_storeAllTimeResults? iTime * NbrRegion : 0 ;
1800
1801 if(iRegion == 0){
1802 switch(Value->Type){
1803 case SCALAR : Size = 1 ; break ;
1804 case VECTOR : Size = 3 ; break ;
1805 case TENSOR_DIAG : Size = 3 ; break ;
1806 case TENSOR_SYM : Size = 6 ; break ;
1807 case TENSOR : Size = 9 ; break ;
1808 }
1809 }
1810
1811 if (Format == FORMAT_REGION_TABLE) {
1812 if(iRegion == 0){
1813 if(PostStream == stdout || PostStream == stderr) {
1814 Message::Direct("%d%s", NbrRegion, Comma ? Comma : "");
1815 }
1816 else if(PostStream) {
1817 fprintf(PostStream, "%d%s\n", NbrRegion, Comma ? Comma : "") ;
1818 }
1819 }
1820 std::ostringstream sstream;
1821 sstream.precision(16);
1822 sstream << numRegion;
1823 for (k = 0 ; k < NbrHarmonics ; k++) {
1824 for(j = 0 ; j < Size ; j++) {
1825 sstream << " " << Value->Val[MAX_DIM*k+j] ;
1826 if (Comma) sstream << Comma;
1827 }
1828 }
1829 if(PostStream == stdout || PostStream == stderr)
1830 Message::Direct(sstream.str().c_str());
1831 else if(PostStream)
1832 fprintf(PostStream, "%s\n", sstream.str().c_str()) ;
1833 }
1834 else if (Format == FORMAT_GETDP) {
1835 std::ostringstream sstream;
1836 sstream.precision(16);
1837 sstream << (PSO_P->ValueName? PSO_P->ValueName : PQ_P->Name);
1838 if (numRegion != NO_REGION) sstream << "~{" << numRegion << "}";
1839
1840 if (NbrHarmonics >= 2 || Size > 1) sstream << "() = {";
1841 else sstream << " =";
1842
1843 for (k = 0 ; k < NbrHarmonics ; k++) {
1844 for(j = 0 ; j < Size ; j++) {
1845 if (k || j) sstream << ",";
1846 sstream << " " << Value->Val[MAX_DIM*k+j] ;
1847 }
1848 }
1849
1850 if (NbrHarmonics >= 2 || Size > 1) sstream << " };";
1851
1852 if(PostStream == stdout || PostStream == stderr)
1853 Message::Direct(sstream.str().c_str());
1854 else if(PostStream)
1855 fprintf(PostStream, "%s\n", sstream.str().c_str()) ;
1856 }
1857 else if (Format == FORMAT_GMSH && Flag_GMSH_VERSION != 2) {
1858 if (Group_FunctionType == NODESOF)
1859 Geo_GetNodesCoordinates(1, &numRegion, &x, &y, &z) ;
1860 else {
1861 x = y = z = 0.;
1862 Message::Warning("Post Format \'Gmsh\' not adapted for global quantities supported"
1863 " by Regions. Zero coordinates are considered.") ;
1864 }
1865 GmshParsed_PrintElement(Time, 0, 1, NbrHarmonics, HarmonicToTime,
1866 POINT_ELEMENT, 1, &x, &y, &z,
1867 Value) ;
1868 }
1869 else if (Format == FORMAT_NXUNV) {
1870 if(PostStream)
1871 Unv_PrintRegion(PostStream, Comma ? 1 : 0, numRegion, NbrHarmonics,
1872 Size, Value, NXUnv_UnitFactor);
1873 }
1874 else if (Format == FORMAT_LOOP_ERROR) {
1875 StorePostOpResult(NbrHarmonics, Value);
1876 }
1877 else if (Format == FORMAT_NODE_TABLE) {
1878 // FIXME: this leads to output files without the total number of nodes at
1879 // the beginning (i.e. not compatible with NodeTable obtained e.g. for
1880 // OnElementsOf)
1881 fprintf(PostStream, "%d", numRegion) ;
1882 Geo_GetNodesCoordinates(1, &numRegion, &x, &y, &z) ;
1883 fprintf(PostStream, " %.16g %.16g %.16g", x,y,z) ;
1884 for (k = 0 ; k < NbrHarmonics ; k++) {
1885 for(j = 0 ; j < Size ; j++) {
1886 fprintf(PostStream, " %.16g", Value->Val[MAX_DIM*k+j]) ;
1887 }
1888 }
1889 fprintf(PostStream, "\n") ;
1890 }
1891 // else, for other FORMATs, e.g., FORMAT_FREQUENCY_TABLE
1892 else {
1893 if(iRegion == 0){
1894 if (!flag_storeAllTimeResults)
1895 TmpValues = (struct Value*) Malloc(NbrRegion*sizeof(struct Value)) ;
1896 else{
1897 if (iTime == 0){
1898 TmpValues = (struct Value*) Malloc(NbrTimeStep*NbrRegion*sizeof(struct Value)) ;
1899 Times = (double*) Malloc(NbrTimeStep*sizeof(double)) ;
1900 }
1901 Times[iTime] = Time ;
1902 }
1903 }
1904
1905 Cal_CopyValue(Value, &TmpValues[indexInTmpValues+iRegion]) ;
1906
1907 if (!flag_storeAllTimeResults && iRegion == NbrRegion-1) {
1908
1909 if (PostStream && HarmonicToTime == 1) {
1910 switch (Format) {
1911 case FORMAT_FREQUENCY_TABLE :
1912 case FORMAT_FREQUENCY_REGION_VALUE :
1913 if (NbrHarmonics == 1){
1914 Message::Error("FrequencyTable format not allowed (only one harmonic)") ;
1915 return;
1916 }
1917 break ;
1918 case FORMAT_VALUE_ONLY :
1919 break;
1920 default :
1921 fprintf(PostStream, " %.16g", Time) ;
1922 if (Comma) fprintf(PostStream, "%s", Comma);
1923 break ;
1924 }
1925 for (iRegion = 0 ; iRegion < NbrRegion ; iRegion++) {
1926 for (k = 0 ; k < NbrHarmonics ; k++) {
1927 if ((Format == FORMAT_FREQUENCY_TABLE ||
1928 Format == FORMAT_FREQUENCY_REGION_VALUE)
1929 && !(k % 2) && iRegion == 0) {
1930 Freq = Current.DofData->Val_Pulsation[0] / TWO_PI ;
1931 fprintf(PostStream, " %.16g", Freq) ;
1932 if (Comma) fprintf(PostStream, "%s", Comma);
1933 }
1934 for(j = 0 ; j < Size ; j++) {
1935 if (Format != FORMAT_REGION_VALUE &&
1936 Format != FORMAT_FREQUENCY_REGION_VALUE) {
1937 fprintf(PostStream, " %.16g",
1938 TmpValues[indexInTmpValues+iRegion].Val[MAX_DIM*k+j]) ;
1939 if (Comma) fprintf(PostStream, "%s", Comma);
1940 }
1941 }
1942 }
1943 }
1944 if (Flag_NoNewLine ||
1945 Format == FORMAT_REGION_VALUE || Format == FORMAT_FREQUENCY_REGION_VALUE)
1946 fprintf(PostStream, " ") ;
1947 else
1948 fprintf(PostStream, "\n") ;
1949 }
1950 else if(PostStream){
1951 for(k = 0 ; k < HarmonicToTime ; k++) {
1952 for (iRegion = 0 ; iRegion < NbrRegion ; iRegion++) {
1953 F_MHToTime0(k+iRegion, &TmpValues[indexInTmpValues+iRegion], &TmpValue,
1954 k, HarmonicToTime, &TimeMH) ;
1955 if (iRegion == 0) {
1956 fprintf(PostStream, " %.16g", TimeMH) ;
1957 if (Comma) fprintf(PostStream, "%s", Comma);
1958 }
1959 for(j = 0 ; j < Size ; j++) {
1960 fprintf(PostStream, " %.16g", TmpValue.Val[j]) ;
1961 if (Comma) fprintf(PostStream, "%s", Comma);
1962 }
1963 }
1964 fprintf(PostStream, "\n") ;
1965 }
1966 }
1967
1968 if (flag_storeAllTimeResults) Free(Times) ;
1969 Free(TmpValues) ;
1970 }
1971
1972 else if (flag_storeAllTimeResults &&
1973 iTime == NbrTimeStep-1 && iRegion == NbrRegion-1) {
1974
1975 Pos_FourierTransform(NbrTimeStep, NbrRegion, Times, TmpValues, Size,
1976 1, PSO_P->TimeToHarmonic, NULL, NULL, NULL);
1977
1978 Free(Times);
1979 Free(TmpValues) ;
1980 }
1981
1982 }
1983 }
1984
1985
1986 /* ------------------------------------------------------------------------ */
1987 /* P o s _ F o u r i e r T r a n s f o r m */
1988 /* ------------------------------------------------------------------------ */
1989
Pos_FourierTransform(int NbrTimeStep,int NbrRegion,double * Times,struct Value * TmpValues,int Size,int TypeOutput,int Nb_Freq_Select_0,int * NbrFreq,double ** Frequencies,struct Value ** OutValues)1990 void Pos_FourierTransform(int NbrTimeStep, int NbrRegion,
1991 double *Times, struct Value *TmpValues, int Size,
1992 int TypeOutput, int Nb_Freq_Select_0,
1993 int *NbrFreq, double **Frequencies, struct Value **OutValues)
1994 {
1995 #if NEW_CODE
1996 *NbrFreq = (NbrTimeStep-1)/2+1;
1997 *Frequencies = (double *)Malloc(*NbrFreq*sizeof(double));
1998 *OutValues = (struct Value *)Malloc(*NbrFreq*sizeof(struct Value));
1999
2000 int nfft = *NbrFreq;
2001 kissfft<double> fft(nfft, false);
2002 std::vector<std::complex<double> > inbuf(nfft);
2003 std::vector<std::complex<double> > outbuf(nfft);
2004 for (int k = 0; k < nfft; ++k)
2005 inbuf[k]= std::complex<double>(rand()/(double)RAND_MAX - .5,
2006 rand()/(double)RAND_MAX - .5);
2007 fft.transform(&inbuf[0], &outbuf[0]);
2008 #else
2009
2010 int iTime, iRegion, k_fc, i_k, j, k;
2011 int N, Nhalf, NbrFourierComps;
2012 double *val_FourierComps;
2013 double val, val_r, val_i, norm, Period, w, v_cos, v_sin;
2014
2015 N = NbrTimeStep-1;
2016 Nhalf = N/2;
2017 // Nhalf = 2;
2018 NbrFourierComps = Nhalf*2;
2019
2020 Period = Times[NbrTimeStep-1] - Times[0];
2021 w = TWO_PI/Period;
2022
2023 val_FourierComps = (double*) Malloc(NbrFourierComps*MAX_DIM*2*sizeof(double)) ;
2024
2025 for (k_fc=-Nhalf; k_fc<Nhalf; k_fc++){
2026 i_k = Nhalf+k_fc;
2027 for (k=0; k<2; k++){
2028 for (j = 0 ; j < Size ; j++){
2029 val_FourierComps[(2*i_k+k)*MAX_DIM + j] = 0.;
2030 }
2031 }
2032 }
2033
2034 for (iTime=0; iTime<N; iTime++){
2035 iRegion = 0; // only for 1 region now!
2036
2037 for (k_fc=-Nhalf; k_fc<Nhalf; k_fc++){
2038 i_k = Nhalf+k_fc;
2039
2040 v_cos = cos( k_fc*w*Times[iTime]);
2041 v_sin = sin(-k_fc*w*Times[iTime]);
2042
2043 for (j = 0 ; j < Size ; j++){
2044 val = TmpValues[iTime*NbrRegion+iRegion].Val[j];
2045
2046 val_FourierComps[(2*i_k+0)*MAX_DIM + j] += val * v_cos;
2047 val_FourierComps[(2*i_k+1)*MAX_DIM + j] += val * v_sin;
2048 }
2049 }
2050
2051 }
2052
2053 for (k_fc=-Nhalf; k_fc<Nhalf; k_fc++){
2054 i_k = Nhalf+k_fc;
2055 for (k=0; k<2; k++){
2056 for (j = 0 ; j < Size ; j++){
2057 val_FourierComps[(2*i_k+k)*MAX_DIM + j] /= N;
2058 }
2059 }
2060 }
2061
2062 if(!PostStream) TypeOutput = 2;
2063
2064 int Nb_Freq_Select = (Nb_Freq_Select_0>0)? Nb_Freq_Select_0 : Nhalf-1;
2065
2066 // Limited to cosine transform now
2067 if (TypeOutput == 2){
2068 *NbrFreq = Nhalf+1;
2069 *Frequencies = (double *)Malloc(*NbrFreq*sizeof(double));
2070 *OutValues = (struct Value *)Malloc(*NbrFreq*sizeof(struct Value));
2071 }
2072
2073 // for (k_fc=-Nhalf; k_fc<Nhalf; k_fc++){
2074 for (k_fc=0; k_fc<=Nb_Freq_Select; k_fc++){
2075 i_k = Nhalf+k_fc;
2076
2077 if (TypeOutput == 1) {
2078 fprintf(PostStream, "%.16g", k_fc*w/TWO_PI) ;
2079 }
2080 else if (TypeOutput == 2) {
2081 (*Frequencies)[k_fc] = k_fc*w/TWO_PI;
2082 (*OutValues)[k_fc].Type = TmpValues[0].Type;
2083 }
2084
2085 for (k=0; k<2; k++){
2086 for(j = 0 ; j < Size ; j++){
2087 /*
2088 if (k_fc != 0)
2089 val = ((k==0)?1:-1) * 2 * val_FourierComps[(2*i_k+k)*MAX_DIM + j];
2090 else
2091 val = val_FourierComps[(2*i_k+k)*MAX_DIM + j];
2092 */
2093
2094 if (k_fc != 0) {
2095 val_r = 2 * val_FourierComps[(2*i_k+0)*MAX_DIM + j];
2096 val_i = -2 * val_FourierComps[(2*i_k+1)*MAX_DIM + j];
2097 norm = sqrt(SQU(val_r) + SQU(val_i));
2098 // val = (k==0)? norm : -asin(val_i/norm); // Phase for CosineTransform
2099 val = (k==0)? norm : atan2(val_i,val_r); // Phase for FourierTransform
2100 }
2101 else {
2102 val = (k==0)? val_FourierComps[(2*i_k+k)*MAX_DIM + j] : 0.;
2103 }
2104
2105 if (TypeOutput == 1) {
2106 fprintf(PostStream, " %.16g", val);
2107 }
2108 if (k == 0 && TypeOutput == 2) {
2109 (*OutValues)[k_fc].Val[j] = val;
2110 }
2111
2112 }
2113 }
2114 if (TypeOutput == 1) {
2115 fprintf(PostStream, "\n") ;
2116 }
2117 }
2118
2119 Free(val_FourierComps);
2120 #endif
2121 }
2122