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 <string.h>
8 #include <math.h>
9 #include "ProData.h"
10 #include "ProParser.h"
11 #include "GeoData.h"
12 #include "DofData.h"
13 #include "Cal_PostQuantity.h"
14 #include "Get_Geometry.h"
15 #include "Get_DofOfElement.h"
16 #include "Get_FunctionValue.h"
17 #include "ExtendedGroup.h"
18 #include "Cal_Quantity.h"
19 #include "Cal_Value.h"
20 #include "Pos_Formulation.h"
21 #include "Pos_Element.h"
22 #include "Pos_Search.h"
23 #include "Pos_Print.h"
24 #include "Pos_Format.h"
25 #include "Adapt.h"
26 #include "MallocUtils.h"
27 #include "Message.h"
28
29 #define SQU(a) ((a)*(a))
30
31 extern struct Problem Problem_S ;
32 extern struct CurrentData Current ;
33
34 extern int Flag_BIN, Flag_GMSH_VERSION;
35
36 extern FILE *PostStream ;
37
38 /*
39 Print OnElementsOf
40 ------------------
41 expl: plot on elements, belonging to the current mesh, where
42 the solution was computed during the processing stage
43 args: list of groups of region type
44
45 Print OnSection
46 ---------------
47 expl: plot an a 'real' cut of the mesh, i.e. computation on the
48 intersections of the mesh with a cutting entity (plane, line)
49 args: 2 (not done) or 3 points, specifying the cutting line or the cutting plane
50
51 Print OnGrid
52 ------------
53 expl: reinterpolate the solution on a grid
54 args: - a list of groups of region type (belonging to a mesh, where the
55 solution will be reinterpolated)
56 - 3 expressions (using $S and $T) and 2 intervals for the parametric
57 grid definition
58
59 Print OnPoint, OnLine, OnPlane, OnBox
60 -------------------------------------
61 expl: reinterpolate the solution on a grid (particular cases)
62 args: 1, 2, 3 or 4 points (0d, 1d, 2d or 3d grid) and the associated
63 number of divisions
64
65 Print OnRegion
66 --------------
67 expl: print Global Quantities associated with Regions
68 args: list of groups of region type
69
70 */
71
72 /* ------------------------------------------------------------------------ */
73 /* P o s _ P r i n t O n E l e m e n t s O f */
74 /* ------------------------------------------------------------------------ */
75
76 struct CutEdge {
77 int nbc ;
78 double x[2],y[2],z[2] ;
79 double xc,yc,zc ;
80 double uc,vc,wc ;
81 struct Value *Value ;
82 } ;
83
84 struct xyzv {
85 double x,y,z;
86 struct Value v;
87 /*int nbvals; for time domain -> malloc Value *v... */
88 int nboccurences;
89 };
90
91 struct ValMinMax {
92 struct Value Val, ValX, ValY, ValZ;
93 };
94
CompareValue(const Value * valA_P,const Value * valB_P)95 int CompareValue(const Value * valA_P, const Value * valB_P)
96 {
97 double cmp=0, VecLengthSquA, VecLengthSquB;
98
99 //if (Current.NbrHar != 1)
100 // Message::Error("Cannot compare multi-harmonic values");
101 // -> we compare the real part in this case
102
103 switch (valA_P->Type) {
104 case SCALAR:
105 cmp = valA_P->Val[0] - valB_P->Val[0];
106 break;
107 case VECTOR:
108 VecLengthSquA = valA_P->Val[0] * valA_P->Val[0] +
109 valA_P->Val[1] * valA_P->Val[1] +
110 valA_P->Val[2] * valA_P->Val[2];
111 VecLengthSquB = valB_P->Val[0] * valB_P->Val[0] +
112 valB_P->Val[1] * valB_P->Val[1] +
113 valB_P->Val[2] * valB_P->Val[2];
114 cmp = VecLengthSquA - VecLengthSquB;
115 break;
116 default:
117 Message::Error("Cannot compare values other than SCALAR and VECTOR");
118 }
119 if(cmp > 1.e-16)
120 return 1;
121 else if(cmp < -1.e-16)
122 return -1;
123 else
124 return 0;
125 }
126
SetValMinMax(struct PostElement * PE_P,int iNode,struct ValMinMax * ValueMinMax_P)127 void SetValMinMax(struct PostElement *PE_P,
128 int iNode,
129 struct ValMinMax *ValueMinMax_P)
130 {
131 Cal_CopyValue(&PE_P->Value[iNode], &ValueMinMax_P->Val);
132 ValueMinMax_P->ValX.Val[0] = PE_P->x[iNode];
133 ValueMinMax_P->ValY.Val[0] = PE_P->y[iNode];
134 ValueMinMax_P->ValZ.Val[0] = PE_P->z[iNode];
135 }
136
InitValMinMax(struct ValMinMax * ValueMinMax_P,struct PostElement * PE_P)137 void InitValMinMax(struct ValMinMax *ValueMinMax_P,
138 struct PostElement *PE_P)
139 {
140 // Init ValueMin and ValueMax
141 ValueMinMax_P->ValX.Type = SCALAR;
142 ValueMinMax_P->ValY.Type = SCALAR;
143 ValueMinMax_P->ValZ.Type = SCALAR;
144 Cal_ZeroValue(&ValueMinMax_P->ValX);
145 Cal_ZeroValue(&ValueMinMax_P->ValY);
146 Cal_ZeroValue(&ValueMinMax_P->ValZ);
147 SetValMinMax(PE_P, 0, ValueMinMax_P);
148 }
149
EvalMinMax(struct PostElement * PE_P,struct ValMinMax * ValueMin_P,struct ValMinMax * ValueMax_P)150 void EvalMinMax(struct PostElement *PE_P,
151 struct ValMinMax *ValueMin_P,
152 struct ValMinMax *ValueMax_P)
153 {
154 for(int iNode = 0 ; iNode < PE_P->NbrNodes ; iNode++) {
155 if (CompareValue(&PE_P->Value[iNode], &ValueMin_P->Val) < 0)
156 SetValMinMax(PE_P, iNode, ValueMin_P);
157 if (CompareValue(&PE_P->Value[iNode], &ValueMax_P->Val) > 0)
158 SetValMinMax(PE_P, iNode, ValueMax_P);
159 }
160 }
161
fcmp_xyzv(const void * a,const void * b)162 static int fcmp_xyzv(const void * a, const void * b)
163 {
164 struct xyzv *p1, *p2;
165 double TOL = Current.GeoData->CharacteristicLength * 1.e-8;
166 p1 = (struct xyzv*)a;
167 p2 = (struct xyzv*)b;
168 if(p1->x - p2->x > TOL) return 1;
169 if(p1->x - p2->x <-TOL) return -1;
170 if(p1->y - p2->y > TOL) return 1;
171 if(p1->y - p2->y <-TOL) return -1;
172 if(p1->z - p2->z > TOL) return 1;
173 if(p1->z - p2->z <-TOL) return -1;
174 return 0;
175 }
176
177 static List_T * SkinPostElement_L ;
178 static int SkinDepth ;
179
Cut_SkinPostElement(void * a,void * b)180 static void Cut_SkinPostElement(void *a, void *b)
181 {
182 struct PostElement * PE ;
183
184 PE = *(struct PostElement**)a ;
185
186 Cut_PostElement(PE, Geo_GetGeoElement(PE->Index), SkinPostElement_L,
187 PE->Index, SkinDepth, 0, 1) ;
188 }
189
Decompose_SkinPostElement(void * a,void * b)190 static void Decompose_SkinPostElement(void *a, void *b)
191 {
192 struct PostElement * PE, * PE2 ;
193
194 PE = *(struct PostElement**)a ;
195
196 if(PE->Type != QUADRANGLE) return;
197 /* change the quad to a tri */
198 PE->Type = TRIANGLE;
199 PE->NbrNodes = 3;
200 /* create a second tri */
201 PE2 = NodeCopy_PostElement(PE) ;
202 PE2->NumNodes[1] = PE->NumNodes[2];
203 PE2->u[1] = PE->u[2]; PE2->x[1] = PE->x[2];
204 PE2->v[1] = PE->v[2]; PE2->y[1] = PE->y[2];
205 PE2->w[1] = PE->w[2]; PE2->z[1] = PE->z[2];
206 PE2->NumNodes[2] = PE->NumNodes[3];
207 PE2->u[2] = PE->u[3]; PE2->x[2] = PE->x[3];
208 PE2->v[2] = PE->v[3]; PE2->y[2] = PE->y[3];
209 PE2->w[2] = PE->w[3]; PE2->z[2] = PE->z[3];
210 List_Add(SkinPostElement_L, &PE2);
211 }
212
Pos_PrintOnElementsOf(struct PostQuantity * NCPQ_P,struct PostQuantity * CPQ_P,int Order,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostSubOperation * PSO_P)213 void Pos_PrintOnElementsOf(struct PostQuantity *NCPQ_P,
214 struct PostQuantity *CPQ_P,
215 int Order,
216 struct DefineQuantity *DefineQuantity_P0,
217 struct QuantityStorage *QuantityStorage_P0,
218 struct PostSubOperation *PSO_P)
219 {
220 Tree_T * PostElement_T ;
221 List_T * PostElement_L, * Region_L ;
222 struct Group * Group_P ;
223
224 struct Element Element ;
225 struct PostElement * PE ;
226 struct Value * CumulativeValues ;
227 struct xyzv xyzv, *xyzv_P ;
228 struct ValMinMax ValueMin, ValueMax ;
229 Tree_T * xyzv_T ;
230 double * Error=NULL, Dummy[5], d, x1, x2 ;
231 int jj, NbrGeo, iGeo, incGeo, NbrPost=0, iPost ;
232 int NbrTimeStep, iTime, iNode ;
233 int Store = 0, DecomposeInSimplex = 0, Depth ;
234 bool StoreMinMax, ValueMinMaxInitialized;
235
236 /* Do we have to store min. and max. values? */
237 if (PSO_P->StoreMinInRegister >= 0 || PSO_P->StoreMaxInRegister >= 0 ||
238 PSO_P->StoreMinXinRegister >= 0 || PSO_P->StoreMaxXinRegister >= 0 ||
239 PSO_P->StoreMinYinRegister >= 0 || PSO_P->StoreMaxYinRegister >= 0 ||
240 PSO_P->StoreMinZinRegister >= 0 || PSO_P->StoreMaxZinRegister >= 0) {
241 StoreMinMax = true;
242 }
243 else
244 StoreMinMax = false;
245
246 /* Select the TimeSteps */
247
248 NbrTimeStep = Pos_InitTimeSteps(PSO_P);
249
250 /* Print the header */
251
252 NbrGeo = Geo_GetNbrGeoElements() ;
253
254 Format_PostHeader(PSO_P, NbrTimeStep, Order,
255 PSO_P->Label ? PSO_P->Label :
256 (NCPQ_P ? NCPQ_P->Name : NULL),
257 PSO_P->Label ? NULL :
258 (CPQ_P ? CPQ_P->Name : NULL));
259
260 /* Get the region */
261
262 Group_P = (struct Group *)List_Pointer(Problem_S.Group,
263 PSO_P->Case.OnRegion.RegionIndex);
264 Region_L = Group_P->InitialList ;
265 Get_InitDofOfElement(&Element) ;
266
267 /* Compute the Cumulative quantity, if any */
268
269 if(CPQ_P){
270 Cal_PostCumulativeQuantity(Region_L,
271 PSO_P->PostQuantitySupport[Order],
272 PSO_P->TimeStep_L,
273 CPQ_P, DefineQuantity_P0,
274 QuantityStorage_P0, &CumulativeValues);
275 }
276
277 /* If we compute a skin, apply smoothing, sort the results, or
278 perform adaption, we'll need to store all the PostElements */
279
280 if(PSO_P->Smoothing || PSO_P->Skin ||
281 PSO_P->Adapt || PSO_P->Sort)
282 Store = 1 ;
283
284 /* Check if everything is OK for Adaption */
285
286 if(PSO_P->Adapt){
287 if(PSO_P->Dimension == DIM_ALL){
288 Message::Error("You have to specify a Dimension for the adaptation (2 or 3)");
289 return;
290 }
291 if(PSO_P->Target < 0.){
292 Message::Error("You have to specify a Target for the adaptation (e.g. 0.01)");
293 return;
294 }
295 if(NbrTimeStep > 1){
296 Message::Error("Adaption not ready with more than one time step");
297 return;
298 }
299 }
300
301 /* Check if we should decompose all PostElements to simplices */
302
303 if(!PSO_P->Skin && PSO_P->DecomposeInSimplex)
304 DecomposeInSimplex = 1 ;
305
306 /* Check for de-refinement */
307
308 if(PSO_P->Depth < 0)
309 incGeo = - PSO_P->Depth ;
310 else
311 incGeo = 1 ;
312
313 /* Create the list of PostElements */
314
315 PostElement_L = List_Create(Store ? NbrGeo/10 : 10, Store ? NbrGeo/10 : 10,
316 sizeof(struct PostElement *)) ;
317
318 if(Store){
319
320 /* If we have a Skin, we will divide after the skin extraction */
321
322 if(PSO_P->Skin && PSO_P->Depth > 1)
323 Depth = 1;
324 else
325 Depth = PSO_P->Depth;
326
327 /* Generate all PostElements */
328
329 Message::ResetProgressMeter();
330 for(iGeo = 0 ; iGeo < NbrGeo ; iGeo += incGeo) {
331 Element.GeoElement = Geo_GetGeoElement(iGeo) ;
332 if ((Group_P->Type != ELEMENTLIST &&
333 List_Search(Region_L, &Element.GeoElement->Region, fcmp_int))
334 ||
335 (Group_P->Type == ELEMENTLIST &&
336 Check_IsEntityInExtendedGroup(Group_P, Element.GeoElement->Num, 0))
337 ) {
338 Fill_PostElement(Element.GeoElement, PostElement_L, iGeo,
339 Depth, PSO_P->Skin,
340 DecomposeInSimplex, 0, PSO_P->Gauss) ;
341 }
342 Message::ProgressMeter(iGeo + 1, NbrGeo, "Post-processing (Generate)");
343 if(Message::GetErrorCount()) break;
344 }
345
346 /* Compute the skin */
347
348 if(PSO_P->Skin){
349 PostElement_T = Tree_Create(sizeof(struct PostElement *), fcmp_PostElement);
350
351 Message::ResetProgressMeter();
352 for(iPost = 0 ; iPost < List_Nbr(PostElement_L) ; iPost++){
353 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
354 if(Tree_PQuery(PostElement_T, &PE)){
355 Tree_Suppress(PostElement_T, &PE);
356 Destroy_PostElement(PE) ;
357 }
358 else
359 Tree_Add(PostElement_T, &PE);
360 Message::ProgressMeter(iPost + 1, List_Nbr(PostElement_L), "Post-processing (Skin)");
361 if(Message::GetErrorCount()) break;
362 }
363
364 /* only decompose in simplices (triangles!) now */
365 if(PSO_P->DecomposeInSimplex){
366 List_Reset(PostElement_L);
367 SkinPostElement_L = PostElement_L ;
368 Tree_Action(PostElement_T, Decompose_SkinPostElement);
369 for(iPost = 0 ; iPost < List_Nbr(SkinPostElement_L) ; iPost++)
370 Tree_Add(PostElement_T,
371 (struct PostElement**)List_Pointer(SkinPostElement_L, iPost)) ;
372 }
373
374 if(PSO_P->Depth > 1){
375 List_Reset(PostElement_L);
376 SkinPostElement_L = PostElement_L ;
377 SkinDepth = PSO_P->Depth ;
378 Tree_Action(PostElement_T, Cut_SkinPostElement) ;
379 }
380 else{
381 List_Delete(PostElement_L);
382 PostElement_L = Tree2List(PostElement_T);
383 }
384
385 Tree_Delete(PostElement_T);
386 }
387
388 } /* if Store */
389
390 /* Loop on GeoElements */
391
392 Message::ResetProgressMeter();
393 ValueMinMaxInitialized = false;
394 for(iGeo = 0 ; iGeo < NbrGeo ; iGeo += incGeo) {
395
396 if(Store){
397 if(iGeo) break ;
398 }
399 else{
400 List_Reset(PostElement_L) ;
401 Element.GeoElement = Geo_GetGeoElement(iGeo) ;
402 if ((Group_P->Type != ELEMENTLIST &&
403 List_Search(Region_L, &Element.GeoElement->Region, fcmp_int))
404 ||
405 (Group_P->Type == ELEMENTLIST &&
406 Check_IsEntityInExtendedGroup(Group_P, Element.GeoElement->Num, 0))
407 ) {
408 int HighOrder =
409 (PSO_P->Format == FORMAT_GMSH && (PSO_P->StoreInField >= 0 ||
410 PSO_P->StoreInMeshBasedField >= 0 ||
411 Flag_GMSH_VERSION == 2)) ? 1 : 0;
412 Fill_PostElement(Element.GeoElement, PostElement_L, iGeo,
413 PSO_P->Depth, PSO_P->Skin,
414 DecomposeInSimplex, HighOrder, PSO_P->Gauss) ;
415 }
416 }
417
418 NbrPost = List_Nbr(PostElement_L) ;
419
420 /* Loop on PostElements */
421
422 for(iPost = 0 ; iPost < NbrPost ; iPost++) {
423
424 PE = *(struct PostElement **)List_Pointer(PostElement_L, iPost) ;
425
426 if(!NCPQ_P){ /* Only one Cumulative */
427 for (iTime = 0 ; iTime < NbrTimeStep ; iTime++){
428 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++)
429 Cal_CopyValue(&CumulativeValues[iTime], &PE->Value[iNode]);
430 if(!Store)
431 Format_PostElement(PSO_P, PSO_P->Iso, 0,
432 Current.Time, iTime, NbrTimeStep,
433 Current.NbrHar, PSO_P->HarmonicToTime,
434 NULL, PE);
435 }
436 }
437 else{ /* There is one non-cumulative */
438
439 if(PSO_P->SubType == PRINT_ONGRID){ /* We re-interpolate */
440 for (iTime = 0 ; iTime < NbrTimeStep ; iTime++) {
441 Pos_InitAllSolutions(PSO_P->TimeStep_L, iTime) ;
442 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
443 InWhichElement(&Current.GeoData->Grid, Region_L, &Element,
444 PSO_P->Dimension,
445 PE->x[iNode], PE->y[iNode], PE->z[iNode],
446 &PE->u[iNode], &PE->v[iNode], &PE->w[iNode]) ;
447 Current.Region = Element.Region ;
448 Current.x = PE->x[iNode] ;
449 Current.y = PE->y[iNode] ;
450 Current.z = PE->z[iNode] ;
451 Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,
452 NULL, &Element,
453 PE->u[iNode], PE->v[iNode], PE->w[iNode], &PE->Value[iNode]);
454 if(CPQ_P)
455 Combine_PostQuantity(PSO_P->CombinationType, Order,
456 &PE->Value[iNode], &CumulativeValues[iNode]) ;
457 }
458 if (StoreMinMax) {
459 if (!ValueMinMaxInitialized){
460 // Init ValueMin and ValueMax
461 InitValMinMax(&ValueMin, PE);
462 InitValMinMax(&ValueMax, PE);
463 ValueMinMaxInitialized = true;
464 }
465 EvalMinMax(PE, &ValueMin, &ValueMax);
466 }
467
468 if(!Store)
469 Format_PostElement(PSO_P, PSO_P->Iso, 0,
470 Current.Time, iTime, NbrTimeStep,
471 Current.NbrHar, PSO_P->HarmonicToTime,
472 NULL, PE);
473 }
474 }
475 else{ /* PRINT_ONREGION: We work on the real mesh */
476 Element.GeoElement = Geo_GetGeoElement(PE->Index) ;
477 Element.Num = Element.GeoElement->Num ;
478 Element.Type = Element.GeoElement->Type ;
479 Current.Region = Element.Region = Element.GeoElement->Region ;
480 Get_NodesCoordinatesOfElement(&Element) ;
481
482 for (iTime = 0 ; iTime < NbrTimeStep ; iTime++) {
483
484 Pos_InitAllSolutions(PSO_P->TimeStep_L, iTime) ;
485
486 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
487 Current.x = PE->x[iNode] ;
488 Current.y = PE->y[iNode] ;
489 Current.z = PE->z[iNode] ;
490 Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,
491 NULL, &Element,
492 PE->u[iNode], PE->v[iNode], PE->w[iNode], &PE->Value[iNode]);
493 if(CPQ_P)
494 Combine_PostQuantity(PSO_P->CombinationType, Order,
495 &PE->Value[iNode], &CumulativeValues[iTime]) ;
496 }
497 if (StoreMinMax) {
498 if (!ValueMinMaxInitialized){
499 // Init ValueMin and ValueMax
500 InitValMinMax(&ValueMin, PE);
501 InitValMinMax(&ValueMax, PE);
502 ValueMinMaxInitialized = true;
503 }
504 EvalMinMax(PE, &ValueMin, &ValueMax);
505 }
506
507 if(!Store)
508 Format_PostElement(PSO_P, PSO_P->Iso, 0,
509 Current.Time, iTime, NbrTimeStep,
510 Current.NbrHar, PSO_P->HarmonicToTime,
511 NULL, PE);
512 }
513 }
514 }
515
516 if(!Store) Destroy_PostElement(PE) ;
517
518 }
519 Message::ProgressMeter(iGeo + 1, NbrGeo, "Post-processing (Compute)");
520 if(Message::GetErrorCount()) break;
521 } /* for iGeo */
522
523 /* Store minimum or maximum value in register */
524 if (StoreMinMax) {
525 if (PSO_P->StoreMinInRegister >= 0)
526 Cal_StoreInRegister(&ValueMin.Val, PSO_P->StoreMinInRegister) ;
527 if (PSO_P->StoreMinXinRegister >= 0)
528 Cal_StoreInRegister(&ValueMin.ValX, PSO_P->StoreMinXinRegister) ;
529 if (PSO_P->StoreMinYinRegister >= 0)
530 Cal_StoreInRegister(&ValueMin.ValY, PSO_P->StoreMinYinRegister) ;
531 if (PSO_P->StoreMinZinRegister >= 0)
532 Cal_StoreInRegister(&ValueMin.ValZ, PSO_P->StoreMinZinRegister) ;
533 if (PSO_P->StoreMaxInRegister >= 0)
534 Cal_StoreInRegister(&ValueMax.Val, PSO_P->StoreMaxInRegister) ;
535 if (PSO_P->StoreMaxXinRegister >= 0)
536 Cal_StoreInRegister(&ValueMax.ValX, PSO_P->StoreMaxXinRegister) ;
537 if (PSO_P->StoreMaxYinRegister >= 0)
538 Cal_StoreInRegister(&ValueMax.ValY, PSO_P->StoreMaxYinRegister) ;
539 if (PSO_P->StoreMaxZinRegister >= 0)
540 Cal_StoreInRegister(&ValueMax.ValZ, PSO_P->StoreMaxZinRegister) ;
541 }
542
543 /* Perform Smoothing */
544
545 if(PSO_P->Smoothing){
546
547 Message::Info("Smoothing (phase 1)");
548
549 xyzv_T = Tree_Create(sizeof(struct xyzv), fcmp_xyzv);
550
551 for (iPost = 0 ; iPost < NbrPost ; iPost++){
552 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
553 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++) {
554 xyzv.x = PE->x[iNode];
555 xyzv.y = PE->y[iNode];
556 xyzv.z = PE->z[iNode];
557 if((xyzv_P = (struct xyzv*)Tree_PQuery(xyzv_T, &xyzv))){
558 x1 = (double)(xyzv_P->nboccurences)/ (double)(xyzv_P->nboccurences + 1.);
559 x2 = 1./(double)(xyzv_P->nboccurences + 1);
560 Cal_AddMultValue2(&xyzv_P->v, x1, &PE->Value[iNode], x2);
561 xyzv_P->nboccurences++;
562 }
563 else{
564 Cal_CopyValue(&PE->Value[iNode],&xyzv.v);
565 xyzv.nboccurences = 1;
566 Tree_Add(xyzv_T, &xyzv);
567 }
568 }
569 }
570
571 Message::Info("Smoothing (phase 2)");
572
573 for(iPost = 0 ; iPost < NbrPost ; iPost++){
574 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
575 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++) {
576 xyzv.x = PE->x[iNode];
577 xyzv.y = PE->y[iNode];
578 xyzv.z = PE->z[iNode];
579 if((xyzv_P = (struct xyzv*)Tree_PQuery(xyzv_T, &xyzv))){
580 Cal_CopyValue(&xyzv_P->v, &PE->Value[iNode]);
581 }
582 else
583 Message::Warning("Node (%g,%g,%g) not found", xyzv.x, xyzv.y, xyzv.z);
584 }
585 }
586
587 Tree_Delete(xyzv_T);
588
589 } /* if Smoothing */
590
591 /* Perform Adaption */
592
593 if(PSO_P->Adapt){
594
595 if(!Current.GeoData->H)
596 Current.GeoData->H = (double*)Malloc((NbrGeo+2)*sizeof(double)) ;
597
598 if(!Current.GeoData->P)
599 Current.GeoData->P = (double*)Malloc((NbrGeo+2)*sizeof(double)) ;
600
601 Error = (double*)Malloc((NbrGeo+1)*sizeof(double)) ;
602
603 /* All elements are perfect... */
604 for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++){
605 Element.GeoElement = Geo_GetGeoElement(iGeo) ;
606 Element.Num = Element.GeoElement->Num ;
607 Element.Type = Element.GeoElement->Type ;
608 Element.Region = Element.GeoElement->Region ;
609 Get_NodesCoordinatesOfElement(&Element) ;
610
611 Current.GeoData->H[iGeo+1] = Cal_MaxEdgeLength(&Element) ;
612 Current.GeoData->P[iGeo+1] = 1. ;
613 Error[iGeo+1] = PSO_P->Target ;
614 }
615
616 /* ...except those we want to optimize */
617 for(iPost = 0 ; iPost < NbrPost ; iPost++){
618 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
619 Error[PE->Index+1] = 0. ;
620 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++)
621 Error[PE->Index+1] += PE->Value[iNode].Val[0] ;
622 Error[PE->Index+1] /= (double)PE->NbrNodes ;
623 }
624
625 Adapt (NbrGeo, PSO_P->Adapt,
626 PSO_P->Dimension, Error,
627 Current.GeoData->H, Current.GeoData->P,
628 PSO_P->Target);
629
630 /* Clean up the interpolation orders to fit to what's available */
631 if(List_Nbr(PSO_P->Value_L)){
632 for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++){
633 for(jj = List_Nbr(PSO_P->Value_L)-1 ; jj >= 0 ; jj--){
634 d = *(double*)List_Pointer(PSO_P->Value_L, jj);
635 if(Current.GeoData->P[iGeo+1] > d || jj == 0){
636 Current.GeoData->P[iGeo+1] = d ;
637 break ;
638 }
639 }
640 }
641 }
642 } /* if Adapt */
643
644 /* Print everything if we are in Store mode */
645
646 if(Store){
647
648 /* Sort the elements */
649
650 switch(PSO_P->Sort){
651 case SORT_BY_POSITION : List_Sort(PostElement_L, fcmp_PostElement) ; break ;
652 case SORT_BY_CONNECTIVITY : Sort_PostElement_Connectivity(PostElement_L) ; break ;
653 }
654
655 Dummy[0] = Dummy[1] = Dummy[2] = Dummy[3] = Dummy[4] = 0. ;
656
657 for(iPost = 0 ; iPost < NbrPost ; iPost++){
658 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
659
660 /* Get the values from adaption */
661 if(PSO_P->Adapt){
662 Element.GeoElement = Geo_GetGeoElement(PE->Index) ;
663
664 Dummy[0] = Element.GeoElement->Num ;
665 Dummy[1] = Error[PE->Index+1] ;
666 Dummy[2] = Current.GeoData->H[PE->Index+1] ;
667 Dummy[3] = Current.GeoData->P[PE->Index+1] ;
668 Dummy[4] = iPost ? 0 : NbrPost ;
669
670 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
671 PE->Value[iNode].Type = SCALAR ;
672 if(PSO_P->Adapt == ADAPT_H1 ||
673 PSO_P->Adapt == ADAPT_H2)
674 PE->Value[iNode].Val[0] = Dummy[2] ;
675 else
676 PE->Value[iNode].Val[0] = Dummy[3] ;
677 }
678 }
679
680 /* Compute curvilinear coord if connection sort */
681 if(PSO_P->Sort == SORT_BY_CONNECTIVITY){
682 Dummy[0] = Dummy[1] ;
683 Dummy[1] = Dummy[0] + sqrt(SQU(PE->x[1]-PE->x[0])+
684 SQU(PE->y[1]-PE->y[0])+
685 SQU(PE->z[1]-PE->z[0])) ;
686 Dummy[2] = PE->v[0] ;
687 Dummy[3] = -1. ;
688 }
689
690 Format_PostElement(PSO_P, PSO_P->Iso, 1,
691 Current.Time, 0, 1,
692 Current.NbrHar, PSO_P->HarmonicToTime,
693 Dummy, PE);
694 }
695 }
696
697 Format_PostFooter(PSO_P, Store);
698
699 if(Store)
700 for(iPost = 0 ; iPost < NbrPost ; iPost++){
701 PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
702 Destroy_PostElement(PE) ;
703 }
704
705 List_Delete(PostElement_L);
706
707 if(CPQ_P) Free(CumulativeValues);
708 if(PSO_P->Adapt) Free(Error) ;
709 }
710
711 /* ------------------------------------------------------------------------ */
712 /* P o s _ P r i n t O n S e c t i o n */
713 /* ------------------------------------------------------------------------ */
714
Plane(double a,double b,double c,double d,double x,double y,double z)715 double Plane(double a, double b, double c, double d,
716 double x, double y, double z)
717 {
718 return (a*x+b*y+c*z+d);
719 }
720
721 static double DIRX[3], DIRY[3], DIRZ[3], XCP, YCP ;
722
fcmp_Angle(const void * a,const void * b)723 int fcmp_Angle (const void *a, const void *b)
724 {
725 struct CutEdge *q , *w;
726 double x1,y1,x2,y2,ang1,ang2;
727
728 q = (struct CutEdge*)a;
729 w = (struct CutEdge*)b;
730
731 x1 = q->xc*DIRX[0] + q->yc*DIRX[1] + q->zc*DIRX[2];
732 y1 = q->xc*DIRY[0] + q->yc*DIRY[1] + q->zc*DIRY[2];
733 x2 = w->xc*DIRX[0] + w->yc*DIRX[1] + w->zc*DIRX[2];
734 y2 = w->xc*DIRY[0] + w->yc*DIRY[1] + w->zc*DIRY[2];
735
736 ang1 = atan2(y1-YCP,x1-XCP);
737 ang2 = atan2(y2-YCP,x2-XCP);
738
739 if(ang1>ang2)return 1;
740 return -1;
741 }
742
prodvec(double * a,double * b,double * c)743 void prodvec (double *a , double *b , double *c)
744 {
745 c[0] = a[1]*b[2]-a[2]*b[1];
746 c[1] = a[2]*b[0]-a[0]*b[2];
747 c[2] = a[0]*b[1]-a[1]*b[0];
748 }
749
normvec(double * a)750 void normvec(double *a)
751 {
752 double mod;
753 mod = sqrt(SQU(a[0])+SQU(a[1])+SQU(a[2]));
754 a[0]/=mod;
755 a[1]/=mod;
756 a[2]/=mod;
757 }
758
759 #define NBR_MAX_CUT 10
760
761 #define LETS_PRINT_THE_RESULT \
762 List_Reset(PE_L); \
763 if(PSO_P->Depth < 2) \
764 List_Add(PE_L, &PE) ; \
765 else \
766 Cut_PostElement(PE, Element.GeoElement, PE_L, PE->Index, \
767 PSO_P->Depth, 0, 1) ; \
768 for(iPost = 0 ; iPost < List_Nbr(PE_L) ; iPost++){ \
769 PE = *(struct PostElement **)List_Pointer(PE_L, iPost) ; \
770 for(iTime = 0 ; iTime < NbTimeStep ; iTime++){ \
771 Pos_InitAllSolutions(PSO_P->TimeStep_L, iTime) ; \
772 for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){ \
773 if(NCPQ_P){ \
774 Current.x = PE->x[iNode] ; \
775 Current.y = PE->y[iNode] ; \
776 Current.z = PE->z[iNode] ; \
777 Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0, \
778 NULL, &Element, PE->u[iNode], PE->v[iNode], PE->w[iNode], \
779 &PE->Value[iNode]); \
780 if(CPQ_P) \
781 Combine_PostQuantity(PSO_P->CombinationType, Order, \
782 &PE->Value[iNode], &CumulativeValues[iTime]) ; \
783 } \
784 else \
785 Cal_CopyValue(&CumulativeValues[iTime],&PE->Value[iNode]); \
786 } \
787 Format_PostElement(PSO_P, PSO_P->Iso,0, \
788 Current.Time, iTime, NbTimeStep, \
789 Current.NbrHar, PSO_P->HarmonicToTime, \
790 NULL, PE); \
791 } \
792 } \
793 for(iPost = 0 ; iPost < List_Nbr(PE_L) ; iPost++) \
794 Destroy_PostElement(*(struct PostElement **)List_Pointer(PE_L, iPost));
795
796
Pos_PrintOnSection(struct PostQuantity * NCPQ_P,struct PostQuantity * CPQ_P,int Order,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostSubOperation * PSO_P)797 void Pos_PrintOnSection(struct PostQuantity *NCPQ_P,
798 struct PostQuantity *CPQ_P,
799 int Order,
800 struct DefineQuantity *DefineQuantity_P0,
801 struct QuantityStorage *QuantityStorage_P0,
802 struct PostSubOperation *PSO_P)
803 {
804 struct CutEdge e[NBR_MAX_CUT];
805 struct Element Element ;
806 struct PostElement * PE ;
807 struct Value * CumulativeValues ;
808 List_T * PE_L ;
809
810 int NbGeoElement, NbTimeStep, NbCut, * NumNodes ;
811 int iPost, iNode, iGeo, iCut, iEdge, iTime ;
812 double A, B, C, D, d1, d2, u, xcg, ycg, zcg ;
813 double x[3], y[3], z[3] ;
814
815 NbTimeStep = Pos_InitTimeSteps(PSO_P);
816
817 PE_L = List_Create(10, 10, sizeof(struct PostElement *)) ;
818
819 for(iCut = 0 ; iCut < NBR_MAX_CUT ; iCut++)
820 e[iCut].Value = (struct Value*) Malloc(NbTimeStep*sizeof(struct Value)) ;
821
822 Format_PostHeader(PSO_P, NbTimeStep, Order,
823 PSO_P->Label ? PSO_P->Label :
824 (NCPQ_P ? NCPQ_P->Name : NULL),
825 PSO_P->Label ? NULL :
826 (CPQ_P ? CPQ_P->Name : NULL));
827
828 if(CPQ_P){
829 Cal_PostCumulativeQuantity(NULL,
830 PSO_P->PostQuantitySupport[Order],
831 PSO_P->TimeStep_L,
832 CPQ_P, DefineQuantity_P0,
833 QuantityStorage_P0, &CumulativeValues);
834 }
835
836 switch(PSO_P->SubType) {
837
838 case PRINT_ONSECTION_1D :
839 Message::Error("Print on 1D cuts not done (use Print OnLine instead)");
840 break;
841
842 case PRINT_ONSECTION_2D :
843
844 /* Ax+By+Cz+D=0 from (x0,y0,z0),(x1,y1,z1),(x2,y2,z2) */
845
846 x[0] = PSO_P->Case.OnSection.x[0] ;
847 y[0] = PSO_P->Case.OnSection.y[0] ;
848 z[0] = PSO_P->Case.OnSection.z[0] ;
849 x[1] = PSO_P->Case.OnSection.x[1] ;
850 y[1] = PSO_P->Case.OnSection.y[1] ;
851 z[1] = PSO_P->Case.OnSection.z[1] ;
852 x[2] = PSO_P->Case.OnSection.x[2] ;
853 y[2] = PSO_P->Case.OnSection.y[2] ;
854 z[2] = PSO_P->Case.OnSection.z[2] ;
855 A = (y[1]-y[0])*(z[2]-z[0]) - (z[1]-z[0])*(y[2]-y[0]) ;
856 B = -(x[1]-x[0])*(z[2]-z[0]) + (z[1]-z[0])*(x[2]-x[0]) ;
857 C = (x[1]-x[0])*(y[2]-y[0]) - (y[1]-y[0])*(x[2]-x[0]) ;
858 D = -A*x[0]-B*y[0]-C*z[0] ;
859
860 /* Cut each element */
861
862 NbGeoElement = Geo_GetNbrGeoElements() ;
863
864 Message::ResetProgressMeter();
865 for(iGeo = 0 ; iGeo < NbGeoElement ; iGeo++) {
866 Element.GeoElement = Geo_GetGeoElement(iGeo) ;
867 Element.Num = Element.GeoElement->Num ;
868 Element.Type = Element.GeoElement->Type ;
869 Current.Region = Element.Region = Element.GeoElement->Region ;
870
871 if((PSO_P->Dimension == DIM_ALL &&
872 (Element.GeoElement->Type != POINT_ELEMENT)) ||
873 (PSO_P->Dimension == DIM_3D &&
874 (Element.GeoElement->Type & (TETRAHEDRON|HEXAHEDRON|PRISM|PYRAMID))) ||
875 (PSO_P->Dimension == DIM_2D &&
876 (Element.GeoElement->Type & (TRIANGLE|QUADRANGLE))) ||
877 (PSO_P->Dimension == DIM_1D &&
878 (Element.GeoElement->Type & LINE))){
879
880 Get_NodesCoordinatesOfElement(&Element) ;
881
882 if(Element.GeoElement->NbrEdges == 0)
883 Geo_CreateEdgesOfElement(Element.GeoElement) ;
884
885 NbCut = 0;
886
887 for(iEdge = 0 ; iEdge < Element.GeoElement->NbrEdges ; iEdge++){
888 NumNodes = Geo_GetNodesOfEdgeInElement(Element.GeoElement, iEdge) ;
889 e[NbCut].x[0] = Element.x[abs(NumNodes[0])-1] ;
890 e[NbCut].y[0] = Element.y[abs(NumNodes[0])-1] ;
891 e[NbCut].z[0] = Element.z[abs(NumNodes[0])-1] ;
892 e[NbCut].x[1] = Element.x[abs(NumNodes[1])-1] ;
893 e[NbCut].y[1] = Element.y[abs(NumNodes[1])-1] ;
894 e[NbCut].z[1] = Element.z[abs(NumNodes[1])-1] ;
895 d1 = Plane(A,B,C,D,e[NbCut].x[0],e[NbCut].y[0],e[NbCut].z[0]);
896 d2 = Plane(A,B,C,D,e[NbCut].x[1],e[NbCut].y[1],e[NbCut].z[1]);
897
898 if(d1*d2 <= 0) {
899 if(d1*d2 < 0.) u = -d2/(d1-d2) ;
900 else if(d1 == 0.) u = 1. ;
901 else u = 0. ;
902 e[NbCut].xc = u*e[NbCut].x[0] + (1.-u)*e[NbCut].x[1];
903 e[NbCut].yc = u*e[NbCut].y[0] + (1.-u)*e[NbCut].y[1];
904 e[NbCut].zc = u*e[NbCut].z[0] + (1.-u)*e[NbCut].z[1];
905
906 if(NCPQ_P)
907 xyz2uvwInAnElement(&Element, e[NbCut].xc, e[NbCut].yc, e[NbCut].zc,
908 &e[NbCut].uc, &e[NbCut].vc, &e[NbCut].wc);
909 NbCut++;
910 }
911 }
912
913 if(NbCut > 3){
914 xcg = ycg = zcg = 0.;
915 for(iCut = 0 ; iCut < NbCut ; iCut++){
916 xcg += e[iCut].xc; ycg += e[iCut].yc; zcg += e[iCut].zc;
917 }
918 xcg /= (double)NbCut; ycg /= (double)NbCut; zcg /= (double)NbCut;
919 DIRZ[0] = A; DIRY[0] = xcg-e[0].xc;
920 DIRZ[1] = B; DIRY[1] = ycg-e[0].yc;
921 DIRZ[2] = C; DIRY[2] = zcg-e[0].zc;
922 normvec(DIRZ); normvec(DIRY); prodvec(DIRY,DIRZ,DIRX); normvec(DIRX);
923 XCP = xcg*DIRX[0] + ycg*DIRX[1] + zcg*DIRX[2];
924 YCP = xcg*DIRY[0] + ycg*DIRY[1] + zcg*DIRY[2];
925 qsort(e,NbCut,sizeof(struct CutEdge), fcmp_Angle);
926 }
927
928 if(NbCut > 2){
929 iCut = 2;
930 while(iCut < NbCut){
931 if(PSO_P->Depth > 0){
932 PE = Create_PostElement(iGeo, TRIANGLE, 3, 1) ;
933 PE->x[0] = e[0].xc; PE->x[1] = e[iCut-1].xc; PE->x[2] = e[iCut].xc;
934 PE->y[0] = e[0].yc; PE->y[1] = e[iCut-1].yc; PE->y[2] = e[iCut].yc;
935 PE->z[0] = e[0].zc; PE->z[1] = e[iCut-1].zc; PE->z[2] = e[iCut].zc;
936 PE->u[0] = e[0].uc; PE->u[1] = e[iCut-1].uc; PE->u[2] = e[iCut].uc;
937 PE->v[0] = e[0].vc; PE->v[1] = e[iCut-1].vc; PE->v[2] = e[iCut].vc;
938 PE->w[0] = e[0].wc; PE->w[1] = e[iCut-1].wc; PE->w[2] = e[iCut].wc;
939 LETS_PRINT_THE_RESULT ;
940 }
941 else{
942 PE = Create_PostElement(iGeo, POINT_ELEMENT, 1, 0) ;
943 PE->x[0] = (e[0].xc + e[iCut-1].xc + e[iCut].xc) / 3. ;
944 PE->y[0] = (e[0].yc + e[iCut-1].yc + e[iCut].yc) / 3. ;
945 PE->z[0] = (e[0].zc + e[iCut-1].zc + e[iCut].zc) / 3. ;
946 PE->u[0] = (e[0].uc + e[iCut-1].uc + e[iCut].uc) / 3. ;
947 PE->v[0] = (e[0].vc + e[iCut-1].vc + e[iCut].vc) / 3. ;
948 PE->w[0] = (e[0].wc + e[iCut-1].wc + e[iCut].wc) / 3. ;
949 LETS_PRINT_THE_RESULT ;
950 }
951 iCut++;
952 }
953 }
954
955 if(NbCut == 2){
956 if(PSO_P->Depth > 0){
957 PE = Create_PostElement(iGeo, LINE, 2, 1) ;
958 PE->x[0] = e[0].xc; PE->x[1] = e[1].xc;
959 PE->y[0] = e[0].yc; PE->y[1] = e[1].yc;
960 PE->z[0] = e[0].zc; PE->z[1] = e[1].zc;
961 PE->u[0] = e[0].uc; PE->u[1] = e[1].uc;
962 PE->v[0] = e[0].vc; PE->v[1] = e[1].vc;
963 PE->w[0] = e[0].wc; PE->w[1] = e[1].wc;
964 LETS_PRINT_THE_RESULT ;
965 }
966 else{
967 PE = Create_PostElement(iGeo, POINT_ELEMENT, 1, 0) ;
968 PE->x[0] = (e[0].xc + e[1].xc) / 2. ;
969 PE->y[0] = (e[0].yc + e[1].yc) / 2. ;
970 PE->z[0] = (e[0].zc + e[1].zc) / 2. ;
971 PE->u[0] = (e[0].uc + e[1].uc) / 2. ;
972 PE->v[0] = (e[0].vc + e[1].vc) / 2. ;
973 PE->w[0] = (e[0].wc + e[1].wc) / 2. ;
974 LETS_PRINT_THE_RESULT ;
975 }
976 }
977
978 }
979 Message::ProgressMeter(iGeo + 1, NbGeoElement, "Post-processing (Cut)") ;
980 if(Message::GetErrorCount()) break;
981 }
982 Format_PostFooter(PSO_P, 0);
983 break;
984
985 default :
986 Message::Error("Unknown operation in Print OnSection");
987 break;
988 }
989
990 List_Delete(PE_L) ;
991 if(CPQ_P) Free(CumulativeValues);
992 for(iCut = 0 ; iCut < NBR_MAX_CUT ; iCut++) Free(e[iCut].Value) ;
993 }
994
995 #undef NBR_MAX_CUT
996 #undef LETS_PRINT_THE_RESULT
997
998 /* ------------------------------------------------------------------------ */
999 /* P o s _ P r i n t O n G r i d */
1000 /* ------------------------------------------------------------------------ */
1001
1002
1003 #define LETS_PRINT_THE_RESULT \
1004 PE->x[0] = Current.xp = Current.x ; \
1005 PE->y[0] = Current.yp = Current.y ; \
1006 PE->z[0] = Current.zp = Current.z ; \
1007 if(!NCPQ_P){ \
1008 for (ts = 0 ; ts < NbTimeStep ; ts++){ \
1009 PE->Value[0] = CumulativeValues[ts] ; \
1010 Format_PostElement(PSO_P, PSO_P->Iso, 0, \
1011 Current.Time, ts, NbTimeStep, \
1012 Current.NbrHar, PSO_P->HarmonicToTime, \
1013 Normal, PE); \
1014 } \
1015 } \
1016 else{ \
1017 InWhichElement(&Current.GeoData->Grid, NULL, &Element, PSO_P->Dimension, \
1018 Current.x, Current.y, Current.z, &u, &v, &w) ; \
1019 Current.Region = Element.Region ; \
1020 if(Element.Num != NO_ELEMENT) \
1021 PE->Index = Geo_GetGeoElementIndex(Element.GeoElement) ; \
1022 else \
1023 PE->Index = NO_ELEMENT ; \
1024 for (ts = 0 ; ts < NbTimeStep ; ts++) { \
1025 Pos_InitAllSolutions(PSO_P->TimeStep_L, ts) ; \
1026 Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0, \
1027 NULL, &Element, u, v, w, &PE->Value[0]); \
1028 if(CPQ_P) \
1029 Combine_PostQuantity(PSO_P->CombinationType, Order, \
1030 &PE->Value[0], &CumulativeValues[ts]) ; \
1031 Format_PostElement(PSO_P, PSO_P->Iso, 0, \
1032 Current.Time, ts, NbTimeStep, \
1033 Current.NbrHar, PSO_P->HarmonicToTime, \
1034 Normal, PE); \
1035 } \
1036 }
1037
1038 #define ARRAY(i,j,k,t) \
1039 Array[ (t) * Current.NbrHar * ((int)N[0]+1) * ((int)N[1]+1) + \
1040 (k) * ((int)N[0]+1) * ((int)N[1]+1) + \
1041 (j) * ((int)N[0]+1) + \
1042 (i) ]
1043
1044 #define LETS_STORE_THE_RESULT \
1045 if(!NCPQ_P){ \
1046 if(CumulativeValues[0].Type != SCALAR) \
1047 Message::Error("Print OnPlane not designed for non scalars with Depth > 1"); \
1048 else \
1049 for (ts = 0 ; ts < NbTimeStep ; ts++) \
1050 for(k = 0 ; k < Current.NbrHar ; k++) \
1051 ARRAY(i1,i2,k,ts) = (float)CumulativeValues[ts].Val[MAX_DIM*k] ; \
1052 } \
1053 else{ \
1054 InWhichElement(&Current.GeoData->Grid, NULL, &Element, PSO_P->Dimension, \
1055 Current.x, Current.y, Current.z, &u, &v, &w) ; \
1056 Current.Region = Element.Region ; \
1057 for (ts = 0 ; ts < NbTimeStep ; ts++) { \
1058 Pos_InitAllSolutions(PSO_P->TimeStep_L, ts) ; \
1059 Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0, \
1060 NULL, &Element, u, v, w, &PE->Value[0]); \
1061 if(PE->Value[0].Type != SCALAR) \
1062 Message::Error("Print OnPlane not designed for non scalars with Depth > 1"); \
1063 if(CPQ_P) \
1064 Combine_PostQuantity(PSO_P->CombinationType, Order, \
1065 &PE->Value[0], &CumulativeValues[ts]) ; \
1066 for(k = 0 ; k < Current.NbrHar ; k++) \
1067 ARRAY(i1,i2,k,ts) = (float)PE->Value[0].Val[MAX_DIM*k] ; \
1068 } \
1069 }
1070
Pos_PrintOnGrid(struct PostQuantity * NCPQ_P,struct PostQuantity * CPQ_P,int Order,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostSubOperation * PSO_P)1071 void Pos_PrintOnGrid(struct PostQuantity *NCPQ_P,
1072 struct PostQuantity *CPQ_P,
1073 int Order,
1074 struct DefineQuantity *DefineQuantity_P0,
1075 struct QuantityStorage *QuantityStorage_P0,
1076 struct PostSubOperation *PSO_P)
1077 {
1078 struct Element Element ;
1079 struct Value * CumulativeValues, Value ;
1080 struct PostElement * PE , * PE2 ;
1081
1082 int i1, i2, i3, k, NbTimeStep, ts ;
1083 float *Array = NULL ;
1084 double u, v, w, Length, Normal[4] = {0., 0., 0., 0.} ;
1085 double X[4], Y[4], Z[4], S[4], N[4], tmp[3];
1086
1087 Get_InitDofOfElement(&Element) ;
1088
1089 NbTimeStep = Pos_InitTimeSteps(PSO_P);
1090
1091 if(CPQ_P){
1092 Cal_PostCumulativeQuantity(NULL, PSO_P->PostQuantitySupport[Order],
1093 PSO_P->TimeStep_L,
1094 CPQ_P, DefineQuantity_P0,
1095 QuantityStorage_P0, &CumulativeValues);
1096 }
1097
1098 Format_PostHeader(PSO_P, NbTimeStep, Order,
1099 PSO_P->Label ? PSO_P->Label :
1100 (NCPQ_P ? NCPQ_P->Name : NULL),
1101 PSO_P->Label ? NULL :
1102 (CPQ_P ? CPQ_P->Name : NULL));
1103
1104 PE = Create_PostElement(0, POINT_ELEMENT, 1, 0) ;
1105
1106 switch(PSO_P->SubType) {
1107
1108 case PRINT_ONGRID_0D :
1109 Current.x = PSO_P->Case.OnGrid.x[0] ;
1110 Current.y = PSO_P->Case.OnGrid.y[0] ;
1111 Current.z = PSO_P->Case.OnGrid.z[0] ;
1112 Normal[0] = Normal[1] = Normal[2] = 0.0 ;
1113 LETS_PRINT_THE_RESULT ;
1114
1115 if (PSO_P->StoreInRegister >= 0)
1116 Cal_StoreInRegister(&PE->Value[0], PSO_P->StoreInRegister) ;
1117 if (PSO_P->StoreInVariable)
1118 Cal_StoreInVariable(&PE->Value[0], PSO_P->StoreInVariable) ;
1119 break;
1120
1121 case PRINT_ONGRID_1D :
1122 X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
1123 Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
1124 Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
1125 N[0] = PSO_P->Case.OnGrid.n[0] ;
1126 Normal[1] = Normal[2] = 0.0 ;
1127 Length = sqrt(SQU(X[1]-X[0]) + SQU(Y[1]-Y[0]) + SQU(Z[1]-Z[0])) ;
1128 for (i1 = 0 ; i1 <= N[0] ; i1++) {
1129 S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
1130 Normal[0] = Length * S[0] ;
1131 Current.x = X[0] + (X[1] - X[0]) * S[0] ;
1132 Current.y = Y[0] + (Y[1] - Y[0]) * S[0] ;
1133 Current.z = Z[0] + (Z[1] - Z[0]) * S[0] ;
1134 LETS_PRINT_THE_RESULT ;
1135 }
1136 break;
1137
1138 case PRINT_ONGRID_2D :
1139 X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
1140 Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
1141 Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
1142 X[2] = PSO_P->Case.OnGrid.x[2] ;
1143 Y[2] = PSO_P->Case.OnGrid.y[2] ;
1144 Z[2] = PSO_P->Case.OnGrid.z[2] ;
1145
1146 S[0] = X[1]-X[0]; S[1] = Y[1]-Y[0]; S[2] = Z[1]-Z[0];
1147 N[0] = X[2]-X[0]; N[1] = Y[2]-Y[0]; N[2] = Z[2]-Z[0];
1148 prodvec(S,N,Normal);
1149 Length = sqrt(SQU(Normal[0])+SQU(Normal[1])+SQU(Normal[2]));
1150 if(!Length){
1151 Message::Warning("Bad plane (null normal)");
1152 return ;
1153 }
1154 Normal[0]/=Length ; Normal[1]/=Length ; Normal[2]/=Length ;
1155
1156 N[0] = PSO_P->Case.OnGrid.n[0] ; N[1] = PSO_P->Case.OnGrid.n[1] ;
1157
1158 if(PSO_P->Depth > 1)
1159 Array = (float*)
1160 Malloc(NbTimeStep*Current.NbrHar*(int)((N[0]+1)*(N[1]+1))*sizeof(float)) ;
1161
1162 for (i1 = 0 ; i1 <= N[0] ; i1++) {
1163 S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
1164 for (i2 = 0 ; i2 <= N[1] ; i2++) {
1165 S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
1166 Current.x = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[1] ;
1167 Current.y = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[1] ;
1168 Current.z = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[1] ;
1169 if(PSO_P->Depth > 1){
1170 LETS_STORE_THE_RESULT ;
1171 }
1172 else{
1173 LETS_PRINT_THE_RESULT ;
1174 }
1175 }
1176 if(PostStream && PSO_P->Depth < 2 && !Flag_BIN) fprintf(PostStream, "\n");
1177 }
1178
1179 if(PSO_P->Depth > 1){
1180 PE2 = Create_PostElement(0, TRIANGLE, 3, 0);
1181 PE2->Value[0].Type = PE2->Value[1].Type = PE2->Value[2].Type = SCALAR ;
1182 for (i1 = 0 ; i1 < N[0] ; i1++) {
1183 S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
1184 S[2] = (double)(i1+1) / (double)(N[0] ? N[0] : 1) ;
1185 for (i2 = 0 ; i2 < N[1] ; i2++) {
1186 S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
1187 S[3] = (double)(i2+1) / (double)(N[1] ? N[1] : 1) ;
1188 PE2->x[0] = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[1] ;
1189 PE2->y[0] = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[1] ;
1190 PE2->z[0] = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[1] ;
1191 PE2->x[1] = X[0] + (X[1] - X[0]) * S[2] + (X[2] - X[0]) * S[1] ;
1192 PE2->y[1] = Y[0] + (Y[1] - Y[0]) * S[2] + (Y[2] - Y[0]) * S[1] ;
1193 PE2->z[1] = Z[0] + (Z[1] - Z[0]) * S[2] + (Z[2] - Z[0]) * S[1] ;
1194 PE2->x[2] = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[3] ;
1195 PE2->y[2] = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[3] ;
1196 PE2->z[2] = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[3] ;
1197 for (ts = 0 ; ts < NbTimeStep ; ts++){
1198 for(k = 0 ; k < Current.NbrHar ; k++){
1199 PE2->Value[0].Val[MAX_DIM*k] = ARRAY(i1,i2,k,ts) ;
1200 PE2->Value[1].Val[MAX_DIM*k] = ARRAY(i1+1,i2,k,ts) ;
1201 PE2->Value[2].Val[MAX_DIM*k] = ARRAY(i1,i2+1,k,ts) ;
1202 }
1203 Format_PostElement(PSO_P, PSO_P->Iso, 0,
1204 Current.Time, ts, NbTimeStep,
1205 Current.NbrHar, PSO_P->HarmonicToTime,
1206 Normal, PE2);
1207 }
1208 PE2->x[0] = X[0] + (X[1] - X[0]) * S[2] + (X[2] - X[0]) * S[3] ;
1209 PE2->y[0] = Y[0] + (Y[1] - Y[0]) * S[2] + (Y[2] - Y[0]) * S[3] ;
1210 PE2->z[0] = Z[0] + (Z[1] - Z[0]) * S[2] + (Z[2] - Z[0]) * S[3] ;
1211 tmp[0] = PE2->x[1]; PE2->x[1] = PE2->x[2]; PE2->x[2] = tmp[0];
1212 tmp[1] = PE2->y[1]; PE2->y[1] = PE2->y[2]; PE2->y[2] = tmp[1];
1213 tmp[2] = PE2->z[1]; PE2->z[1] = PE2->z[2]; PE2->z[2] = tmp[2];
1214 for (ts = 0 ; ts < NbTimeStep ; ts++){
1215 for(k = 0 ; k < Current.NbrHar ; k++){
1216 PE2->Value[0].Val[MAX_DIM*k] = ARRAY(i1+1,i2+1,k,ts) ;
1217 PE2->Value[1].Val[MAX_DIM*k] = ARRAY(i1,i2+1,k,ts) ;
1218 PE2->Value[2].Val[MAX_DIM*k] = ARRAY(i1+1,i2,k,ts) ;
1219 }
1220 Format_PostElement(PSO_P, PSO_P->Iso, 0,
1221 Current.Time, ts, NbTimeStep,
1222 Current.NbrHar, PSO_P->HarmonicToTime,
1223 Normal, PE2);
1224 }
1225 }
1226 }
1227 Destroy_PostElement(PE2) ;
1228 Free(Array) ;
1229 }
1230 break;
1231
1232 case PRINT_ONGRID_3D :
1233 X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
1234 Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
1235 Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
1236 X[2] = PSO_P->Case.OnGrid.x[2] ; X[3] = PSO_P->Case.OnGrid.x[3] ;
1237 Y[2] = PSO_P->Case.OnGrid.y[2] ; Y[3] = PSO_P->Case.OnGrid.y[3] ;
1238 Z[2] = PSO_P->Case.OnGrid.z[2] ; Z[3] = PSO_P->Case.OnGrid.z[3] ;
1239 N[0] = PSO_P->Case.OnGrid.n[0] ;
1240 N[1] = PSO_P->Case.OnGrid.n[1] ;
1241 N[2] = PSO_P->Case.OnGrid.n[2] ;
1242 Normal[0] = Normal[1] = Normal[2] = 0.0 ;
1243 for (i1 = 0 ; i1 <= N[0] ; i1++) {
1244 S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
1245 for (i2 = 0 ; i2 <= N[1] ; i2++) {
1246 S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
1247 for (i3 = 0 ; i3 <= N[2] ; i3++) {
1248 S[2] = (double)i3 / (double)(N[2] ? N[2] : 1) ;
1249 Current.x = X[0] + (X[1]-X[0])*S[0] + (X[2]-X[0])*S[1] + (X[3]-X[0])*S[2] ;
1250 Current.y = Y[0] + (Y[1]-Y[0])*S[0] + (Y[2]-Y[0])*S[1] + (Y[3]-Y[0])*S[2] ;
1251 Current.z = Z[0] + (Z[1]-Z[0])*S[0] + (Z[2]-Z[0])*S[1] + (Z[3]-Z[0])*S[2] ;
1252 LETS_PRINT_THE_RESULT ;
1253 }
1254 if(PostStream && !Flag_BIN) fprintf(PostStream, "\n");
1255 }
1256 if(PostStream && !Flag_BIN) fprintf(PostStream, "\n\n");
1257 /* two blanks lines for -index in gnuplot */
1258 }
1259 break;
1260
1261 case PRINT_ONGRID_PARAM :
1262 for (i1 = 0 ; i1 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[0]) ; i1++) {
1263 List_Read(PSO_P->Case.OnParamGrid.ParameterValue[0], i1, &Current.a) ;
1264 for (i2 = 0 ; i2 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[1]) ; i2++) {
1265 List_Read(PSO_P->Case.OnParamGrid.ParameterValue[1], i2, &Current.b) ;
1266 for (i3 = 0 ; i3 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[2]) ; i3++) {
1267 List_Read(PSO_P->Case.OnParamGrid.ParameterValue[2], i3, &Current.c) ;
1268 Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[0],
1269 NULL, 0., 0., 0., &Value) ;
1270 Current.x = Value.Val[0];
1271 Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[1],
1272 NULL, 0., 0., 0., &Value) ;
1273 Current.y = Value.Val[0];
1274 Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[2],
1275 NULL, 0., 0., 0., &Value) ;
1276 Current.z = Value.Val[0];
1277 Normal[0] = Current.a ;
1278 Normal[1] = Current.b ;
1279 Normal[2] = Current.c ;
1280 LETS_PRINT_THE_RESULT ;
1281 }
1282 if(PostStream && List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[2])>1 && !Flag_BIN)
1283 fprintf(PostStream, "\n");
1284 }
1285 if(PostStream && List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[1])>1 && !Flag_BIN)
1286 fprintf(PostStream, "\n\n");
1287 /* two blanks lines for -index in gnuplot */
1288 }
1289 break;
1290 }
1291
1292 Destroy_PostElement(PE) ;
1293
1294 Format_PostFooter(PSO_P, 0);
1295
1296 if(CPQ_P) Free(CumulativeValues);
1297 }
1298
1299 #undef LETS_PRINT_THE_RESULT
1300 #undef LETS_STORE_THE_RESULT
1301 #undef ARRAY
1302
1303
1304 /* ------------------------------------------------------------------------ */
1305 /* P o s _ P r i n t O n R e g i o n */
1306 /* ------------------------------------------------------------------------ */
1307
Pos_PrintOnRegion(struct PostQuantity * NCPQ_P,struct PostQuantity * CPQ_P,int Order,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostSubOperation * PSO_P)1308 void Pos_PrintOnRegion(struct PostQuantity *NCPQ_P,
1309 struct PostQuantity *CPQ_P,
1310 int Order,
1311 struct DefineQuantity *DefineQuantity_P0,
1312 struct QuantityStorage *QuantityStorage_P0,
1313 struct PostSubOperation *PSO_P)
1314 {
1315 struct Element Element ;
1316 struct Value Value, ValueSummed ;
1317 struct PostQuantity *PQ_P ;
1318 struct Group * Group_P ;
1319
1320 List_T *Region_L, *Support_L ;
1321 int i, iTime, NbrTimeStep ;
1322 int Nbr_Region=0, Num_Region, Group_FunctionType ;
1323 int Type_Evaluation=0;
1324 double u, v, w;
1325
1326 NbrTimeStep = Pos_InitTimeSteps(PSO_P);
1327
1328 if (CPQ_P && NCPQ_P){
1329 Message::Error("Only one PostProcessing Quantity allowed in PostOperation") ;
1330 return;
1331 }
1332
1333 if (CPQ_P) {
1334 PQ_P = CPQ_P ;
1335 /*
1336 If the PostQuantityTerm is an Integral quantity to be integrated
1337 over a support in this 'Print' PostOperation
1338 (i.e. syntax: PQ[ Support ]), the InitialList (list of regions) of the
1339 group 'Support' is affected to Support_L.
1340 If however the group 'Support' is of type ELEMENTLIST,
1341 the latter is substituted to the ->InIndex of the PostQuantityTerm here,
1342 i.e., just before evaluation.
1343
1344 For consistency, it should be checked that all ellements
1345 in the ELEMENTLIST Support are indeed in the region PostQuantityTerm->InIndex.
1346 This is not done.
1347 */
1348
1349 /* code original - enlever apres vérification FH
1350 Support_L = // for e.g. PQ[ Support ] ...
1351 ((struct Group *)
1352 List_Pointer(Problem_S.Group,
1353 PSO_P->PostQuantitySupport[Order]))->InitialList ;
1354 */
1355
1356 // FIXME reuse Group_P instead of definig a specific variable ?
1357 struct Group * SupportGroup_P = (struct Group *)
1358 List_Pointer(Problem_S.Group, PSO_P->PostQuantitySupport[Order]);
1359
1360 Support_L = SupportGroup_P->InitialList; // for e.g. PQ[ Support ] ...
1361
1362 if( SupportGroup_P->Type == ELEMENTLIST ) {
1363 ((struct PostQuantityTerm *)
1364 List_Pointer(PQ_P->PostQuantityTerm, 0))->InIndex = SupportGroup_P->Num;
1365 // FIXME What if PostQuantity has several PostQuantityTerm's ?
1366 // Here only the first term (index 0) is considered.
1367 }
1368
1369
1370 }
1371 else {
1372 PQ_P = NCPQ_P ; Support_L = NULL ;
1373 }
1374
1375 Format_PostHeader(PSO_P, NbrTimeStep, Order,
1376 PSO_P->Label ? PSO_P->Label :
1377 (NCPQ_P ? NCPQ_P->Name : NULL),
1378 PSO_P->Label ? NULL :
1379 (CPQ_P ? CPQ_P->Name : NULL));
1380
1381 Group_P = (PSO_P->Case.OnRegion.RegionIndex < 0)? NULL :
1382 (struct Group *)
1383 List_Pointer(Problem_S.Group,
1384 PSO_P->Case.OnRegion.RegionIndex);
1385
1386 Region_L = Group_P? Group_P->InitialList : NULL ;
1387 Group_FunctionType = Group_P? Group_P->FunctionType : REGION;
1388
1389 if (!Support_L &&
1390 List_Nbr(NCPQ_P->PostQuantityTerm) &&
1391 (
1392 ((struct PostQuantityTerm *)List_Pointer(NCPQ_P->PostQuantityTerm, 0))
1393 ->Type == LOCALQUANTITY &&
1394 ((struct PostQuantityTerm *)List_Pointer(NCPQ_P->PostQuantityTerm, 0))
1395 ->EvaluationType == LOCAL)
1396 ) {
1397 if (Group_FunctionType == REGION){
1398 Message::Error("Print OnRegion not valid for PostProcessing Quantity '%s'",
1399 NCPQ_P->Name);
1400 return;
1401 }
1402 else
1403 Type_Evaluation = LOCAL;
1404 }
1405 else
1406 Type_Evaluation = GLOBAL;
1407
1408 if (Region_L) {
1409 if (Group_P->FunctionType == REGION) {
1410 List_Sort(Region_L, fcmp_int) ;
1411 Nbr_Region = List_Nbr(Region_L) ;
1412
1413 if (!PSO_P->NoTitle &&
1414 PSO_P->Format != FORMAT_SPACE_TABLE &&
1415 PSO_P->Format != FORMAT_VALUE_ONLY &&
1416 PSO_P->Format != FORMAT_GETDP) {
1417 std::ostringstream sstream;
1418 if (PSO_P->Format == FORMAT_GMSH)
1419 sstream << "// ";
1420 else
1421 sstream << "# ";
1422 sstream << PQ_P->Name << " on";
1423 for(i = 0 ; i < Nbr_Region ; i++) {
1424 List_Read(Region_L, i, &Num_Region) ;
1425 sstream << " " << Num_Region;
1426 }
1427 if(PostStream == stdout || PostStream == stderr)
1428 Message::Direct(sstream.str().c_str());
1429 else if(PostStream)
1430 fprintf(PostStream, "%s\n", sstream.str().c_str()) ;
1431 }
1432 }
1433 else if (Group_P->FunctionType == NODESOF) {
1434 if (!Group_P->ExtendedList) Generate_ExtendedGroup(Group_P) ;
1435 Region_L = Group_P->ExtendedList ; /* Attention: new Region_L */
1436 Nbr_Region = List_Nbr(Region_L) ;
1437 }
1438 else {
1439 Message::Error("Function type (%d) not allowed for PrintOnRegion",
1440 Group_P->FunctionType) ;
1441 return;
1442 }
1443 }
1444 else
1445 Nbr_Region = 1 ;
1446
1447 for (iTime = 0 ; iTime < NbrTimeStep ; iTime++) {
1448
1449 Pos_InitAllSolutions(PSO_P->TimeStep_L, iTime) ;
1450
1451 if (PSO_P->Format == FORMAT_REGION_VALUE ||
1452 PSO_P->Format == FORMAT_FREQUENCY_REGION_VALUE) {
1453 Cal_ZeroValue(&ValueSummed) ;
1454 }
1455
1456 if(Nbr_Region > 1)
1457 Message::ResetProgressMeter();
1458
1459 for(i = 0 ; i < Nbr_Region ; i++) {
1460
1461 if (Region_L)
1462 List_Read(Region_L, i, &Num_Region) ;
1463 else
1464 Num_Region = NO_REGION ;
1465 Current.SubRegion = Num_Region ; /* Region being a GlobalQuantity Entity no */
1466
1467 Current.NumEntity = Num_Region ; /* for OnRegion NodesOf */
1468
1469 Element.GeoElement = NULL ;
1470 Element.Num = NO_ELEMENT ;
1471 Element.Type = -1 ;
1472 Current.Region = Element.Region = Num_Region ;
1473 Current.x = Current.y = Current.z = 0. ;
1474
1475 if (Type_Evaluation == GLOBAL) {
1476 Cal_PostQuantity(PQ_P, DefineQuantity_P0, QuantityStorage_P0,
1477 Support_L, &Element, 0., 0., 0., &Value) ;
1478 }
1479 else {
1480 if (Group_FunctionType == NODESOF)
1481 Geo_GetNodesCoordinates(1, &Num_Region,
1482 &Current.x, &Current.y, &Current.z) ;
1483 InWhichElement(&Current.GeoData->Grid, NULL, &Element,
1484 PSO_P->Dimension,
1485 Current.x, Current.y, Current.z, &u, &v, &w) ;
1486
1487 Cal_PostQuantity(PQ_P, DefineQuantity_P0, QuantityStorage_P0,
1488 Support_L, &Element, u, v, w, &Value) ;
1489 }
1490
1491 if (PSO_P->Format != FORMAT_REGION_VALUE &&
1492 PSO_P->Format != FORMAT_FREQUENCY_REGION_VALUE) {
1493 if (PSO_P->StoreInRegister >= 0)
1494 Cal_StoreInRegister(&Value, PSO_P->StoreInRegister) ;
1495 if (PSO_P->StoreInVariable)
1496 Cal_StoreInVariable(&Value, PSO_P->StoreInVariable) ;
1497 if (PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No")){
1498 std::vector<double> v;
1499 Export_Value(&Value, v, PSO_P->SendToServerList);
1500 Message::AddOnelabNumberChoice(PSO_P->SendToServer, v, PSO_P->Color,
1501 PSO_P->Units, PSO_P->Label, PSO_P->Visible,
1502 PSO_P->Closed);
1503 }
1504 }
1505
1506 Format_PostValue(PQ_P, PSO_P, PSO_P->Format, PSO_P->Comma,
1507 Group_FunctionType,
1508 iTime, Current.Time, NbrTimeStep,
1509 i, Current.NumEntity, Nbr_Region,
1510 Current.NbrHar, PSO_P->HarmonicToTime,
1511 PSO_P->FourierTransform,
1512 PSO_P->NoNewLine,
1513 &Value) ;
1514
1515 if (PSO_P->Format == FORMAT_REGION_VALUE ||
1516 PSO_P->Format == FORMAT_FREQUENCY_REGION_VALUE) {
1517 ValueSummed.Type = Value.Type ;
1518 Cal_AddValue(&ValueSummed, &Value, &ValueSummed);
1519 }
1520
1521 if(Nbr_Region > 1)
1522 Message::ProgressMeter(i + 1, Nbr_Region, "Post-processing (OnRegion)");
1523 }
1524
1525 if (PostStream && (PSO_P->Format == FORMAT_REGION_VALUE ||
1526 PSO_P->Format == FORMAT_FREQUENCY_REGION_VALUE)) {
1527 fprintf(PostStream, "%s", Print_Value_ToString(&ValueSummed).c_str());
1528 if (PSO_P->StoreInRegister >= 0)
1529 Cal_StoreInRegister(&ValueSummed, PSO_P->StoreInRegister) ;
1530 if (PSO_P->StoreInVariable)
1531 Cal_StoreInVariable(&ValueSummed, PSO_P->StoreInVariable) ;
1532 if (PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No")){
1533 std::vector<double> v;
1534 Export_Value(&ValueSummed, v, PSO_P->SendToServerList);
1535 Message::AddOnelabNumberChoice(PSO_P->SendToServer, v, PSO_P->Color,
1536 PSO_P->Units, PSO_P->Label, PSO_P->Visible,
1537 PSO_P->Closed);
1538 }
1539 }
1540 }
1541
1542 // prevent SendToServer here, as we have alredy done it
1543 Format_PostFooter(PSO_P, 0, false);
1544 }
1545
1546 /* ------------------------------------------------------------------------ */
1547 /* P o s _ P r i n t W i t h A r g u m e n t */
1548 /* ------------------------------------------------------------------------ */
1549
Pos_PrintWithArgument(struct PostQuantity * NCPQ_P,struct PostQuantity * CPQ_P,int Order,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostSubOperation * PSO_P)1550 void Pos_PrintWithArgument(struct PostQuantity *NCPQ_P,
1551 struct PostQuantity *CPQ_P,
1552 int Order,
1553 struct DefineQuantity *DefineQuantity_P0,
1554 struct QuantityStorage *QuantityStorage_P0,
1555 struct PostSubOperation *PSO_P)
1556 {
1557 struct Element Element ;
1558 struct Value Value ;
1559
1560 struct Expression * Expression_P ;
1561 List_T *Region_L ;
1562 int i, N, Num_Region ;
1563 double X[2], S, x ;
1564
1565 if(CPQ_P){
1566 Message::Error("Cumulative PostProcessing Quantity in PrintWithArgument not done") ;
1567 return;
1568 }
1569
1570 X[0] = PSO_P->Case.WithArgument.x[0] ;
1571 X[1] = PSO_P->Case.WithArgument.x[1] ;
1572 N = PSO_P->Case.WithArgument.n ;
1573
1574 Expression_P = (struct Expression *)
1575 List_Pointer(Problem_S.Expression,
1576 PSO_P->Case.WithArgument.ArgumentIndex) ;
1577
1578 Region_L = ((struct Group *)
1579 List_Pointer(Problem_S.Group,
1580 PSO_P->Case.WithArgument.RegionIndex))
1581 ->InitialList ;
1582
1583 if (List_Nbr(Region_L))
1584 List_Read(Region_L, 0, &Num_Region) ;
1585 else
1586 Num_Region = NO_REGION ;
1587
1588 for (i = 0 ; i <= N ; i++) {
1589 S = (double)i / (double)(N ? N : 1) ;
1590 x = X[0] + (X[1] - X[0]) * S ;
1591 Expression_P->Case.Constant = x ;
1592
1593 Element.GeoElement = NULL ;
1594 Element.Num = NO_ELEMENT ;
1595 Element.Type = -1 ;
1596 Current.Region = Element.Region = Num_Region ;
1597 Current.x = Current.y = Current.z = 0. ;
1598
1599 Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,
1600 NULL, &Element, 0., 0., 0., &Value) ;
1601
1602 Format_PostValue(NCPQ_P, PSO_P, PSO_P->Format, PSO_P->Comma,
1603 REGION,
1604 0, x, 1, 0, 0, 1,
1605 Current.NbrHar, PSO_P->HarmonicToTime,
1606 PSO_P->FourierTransform,
1607 PSO_P->NoNewLine,
1608 &Value) ;
1609 }
1610 }
1611
1612 /* ------------------------------------------------------------------------ */
1613 /* P o s _ P r i n t E x p r e s s i o n */
1614 /* ------------------------------------------------------------------------ */
1615
Pos_PrintExpression(struct PostSubOperation * PSO_P)1616 void Pos_PrintExpression(struct PostSubOperation *PSO_P)
1617 {
1618 int NbrTimeStep, iTime;
1619 struct Value Value;
1620 char *str = PSO_P->Case.Expression.String;
1621 char *str2 = PSO_P->Case.Expression.String2;
1622 List_T *expr = PSO_P->Case.Expression.Expressions;
1623
1624 if((!str || !strlen(str)) && (!str2 || !strlen(str2)) && !List_Nbr(expr))
1625 return; // nothing to print; useful to request merging an existing file
1626
1627 NbrTimeStep = Pos_InitTimeSteps(PSO_P);
1628
1629 for(iTime = 0; iTime < NbrTimeStep; iTime++){
1630 Pos_InitAllSolutions(PSO_P->TimeStep_L, iTime) ;
1631 if(List_Nbr(expr) && str2){ // old style, unformatted single expession output
1632 int j; List_Read(expr, 0, &j);
1633 Get_ValueOfExpressionByIndex(j, NULL, 0., 0., 0., &Value) ;
1634 if(PostStream){
1635 if(str) fprintf(PostStream, "%s", str);
1636 fprintf(PostStream, "%s", Print_Value_ToString(&Value).c_str());
1637 }
1638 }
1639 else if(List_Nbr(expr) && str){ // new style, same syntax as in resolution operations
1640 List_T *list = List_Create(10, 10, sizeof(double));
1641 for(int i = 0; i < List_Nbr(expr); i++){
1642 int j; List_Read(expr, i, &j) ;
1643 Get_ValueOfExpressionByIndex(j, NULL, 0., 0., 0., &Value) ;
1644 List_Add(list, &Value.Val[0]);
1645 }
1646 char buffer[1024];
1647 Print_ListOfDouble(str, list, buffer);
1648 if(PostStream)
1649 fprintf(PostStream, "%s", buffer);
1650 List_Delete(list);
1651 }
1652 else if(str2){
1653 if(PostStream){
1654 if(str) fprintf(PostStream, "%s", str);
1655 fprintf(PostStream, "%s", str2);
1656 }
1657 }
1658 else if(str){
1659 if(PostStream)
1660 fprintf(PostStream, "%s", str);
1661 }
1662
1663 if(PostStream){
1664 if(PSO_P->NoNewLine)
1665 fprintf(PostStream, " ") ;
1666 else
1667 fprintf(PostStream, "\n") ;
1668 }
1669 }
1670 }
1671
1672 /* ------------------------------------------------------------------------ */
1673 /* P o s _ P r i n t G r o u p */
1674 /* ------------------------------------------------------------------------ */
1675
Pos_PrintGroup(struct PostSubOperation * PSO_P)1676 void Pos_PrintGroup(struct PostSubOperation *PSO_P)
1677 {
1678 struct Group *Group_P;
1679 struct Geo_Element *GeoElement;
1680 struct PostElement *SL, *ST, *SQ;
1681 List_T *Region_L;
1682 int i, NbrGeo, iGeo, *NumNodes;
1683 double x [NBR_MAX_NODES_IN_ELEMENT] ;
1684 double y [NBR_MAX_NODES_IN_ELEMENT] ;
1685 double z [NBR_MAX_NODES_IN_ELEMENT] ;
1686
1687 NbrGeo = Geo_GetNbrGeoElements() ;
1688
1689 Format_PostHeader(PSO_P, 1, 0, PSO_P->Label, NULL);
1690
1691 Region_L = ((struct Group *)
1692 List_Pointer(Problem_S.Group,
1693 PSO_P->Case.Group.GroupIndex))->InitialList ;
1694 Group_P = (struct Group *)
1695 List_Pointer(Problem_S.Group,
1696 PSO_P->Case.Group.ExtendedGroupIndex);
1697
1698 SL = Create_PostElement(0, LINE, 2, 1) ;
1699 ST = Create_PostElement(0, TRIANGLE, 3, 1) ;
1700 SQ = Create_PostElement(0, QUADRANGLE, 4, 1) ;
1701
1702 if(!Group_P->ExtendedList) Generate_ExtendedGroup(Group_P) ;
1703
1704 Message::ResetProgressMeter();
1705 for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++) {
1706
1707 GeoElement = Geo_GetGeoElement(iGeo) ;
1708 if(List_Search(Region_L, &GeoElement->Region, fcmp_int)){
1709
1710 Geo_GetNodesCoordinates
1711 (GeoElement->NbrNodes, GeoElement->NumNodes, x, y, z) ;
1712
1713 switch (Group_P->FunctionType) {
1714
1715 case EDGESOF : case EDGESOFTREEIN :
1716 if(!GeoElement->NbrEdges) Geo_CreateEdgesOfElement(GeoElement) ;
1717 for(i=0 ; i<GeoElement->NbrEdges ; i++){
1718 if(List_Search(Group_P->ExtendedList, &GeoElement->NumEdges[i], fcmp_absint)){
1719 NumNodes = Geo_GetNodesOfEdgeInElement(GeoElement, i) ;
1720 SL->Index = iGeo;
1721 SL->x[0] = x[abs(NumNodes[0])-1]; SL->x[1] = x[abs(NumNodes[1])-1];
1722 SL->y[0] = y[abs(NumNodes[0])-1]; SL->y[1] = y[abs(NumNodes[1])-1];
1723 SL->z[0] = z[abs(NumNodes[0])-1]; SL->z[1] = z[abs(NumNodes[1])-1];
1724 SL->Value[0].Type = SL->Value[1].Type = SCALAR ;
1725 SL->Value[0].Val[0] = SL->Value[1].Val[0] = GeoElement->NumEdges[i];
1726 Format_PostElement(PSO_P, PSO_P->Iso, 0,
1727 0, 0, 1, 1, 1,
1728 NULL, SL);
1729 }
1730 }
1731 break ;
1732
1733 case FACETSOFTREEIN :
1734 if(!GeoElement->NbrFacets) Geo_CreateFacetsOfElement(GeoElement) ;
1735 for(i=0 ; i<GeoElement->NbrFacets ; i++){
1736 if(List_Search(Group_P->ExtendedList, &GeoElement->NumFacets[i], fcmp_absint)){
1737 NumNodes = Geo_GetNodesOfFacetInElement(GeoElement, i) ;
1738 if(!NumNodes[3]){ // we have triangle
1739 ST->Index = iGeo;
1740 ST->x[0] = x[abs(NumNodes[0])-1]; ST->x[1] = x[abs(NumNodes[1])-1];
1741 ST->y[0] = y[abs(NumNodes[0])-1]; ST->y[1] = y[abs(NumNodes[1])-1];
1742 ST->z[0] = z[abs(NumNodes[0])-1]; ST->z[1] = z[abs(NumNodes[1])-1];
1743 ST->x[2] = x[abs(NumNodes[2])-1];
1744 ST->y[2] = y[abs(NumNodes[2])-1];
1745 ST->z[2] = z[abs(NumNodes[2])-1];
1746 ST->Value[0].Type = ST->Value[1].Type = ST->Value[2].Type = SCALAR ;
1747 ST->Value[0].Val[0] = ST->Value[1].Val[0] = ST->Value[2].Val[0] =
1748 GeoElement->NumFacets[i];
1749 Format_PostElement(PSO_P, PSO_P->Iso, 0,
1750 0, 0, 1, 1, 1,
1751 NULL, ST);
1752 }
1753 else{ // we have a quad
1754 SQ->Index = iGeo;
1755 SQ->x[0] = x[abs(NumNodes[0])-1]; SQ->x[1] = x[abs(NumNodes[1])-1];
1756 SQ->y[0] = y[abs(NumNodes[0])-1]; SQ->y[1] = y[abs(NumNodes[1])-1];
1757 SQ->z[0] = z[abs(NumNodes[0])-1]; SQ->z[1] = z[abs(NumNodes[1])-1];
1758 SQ->x[2] = x[abs(NumNodes[2])-1]; SQ->x[3] = x[abs(NumNodes[3])-1];
1759 SQ->y[2] = y[abs(NumNodes[2])-1]; SQ->y[3] = y[abs(NumNodes[3])-1];
1760 SQ->z[2] = z[abs(NumNodes[2])-1]; SQ->z[3] = z[abs(NumNodes[3])-1];
1761 SQ->Value[0].Type = SQ->Value[1].Type = SQ->Value[2].Type =
1762 SQ->Value[3].Type = SCALAR ;
1763 SQ->Value[0].Val[0] = SQ->Value[1].Val[0] = SQ->Value[2].Val[0] =
1764 SQ->Value[3].Val[0] = GeoElement->NumFacets[i];
1765 Format_PostElement(PSO_P, PSO_P->Iso, 0,
1766 0, 0, 1, 1, 1,
1767 NULL, SQ);
1768 }
1769 }
1770 }
1771 break ;
1772
1773 default :
1774 Message::Error("Print function not implemented for this kind of Group");
1775 break ;
1776
1777 }
1778
1779 }
1780 Message::ProgressMeter(iGeo + 1, NbrGeo, "Post-processing (Compute)") ;
1781 if(Message::GetErrorCount()) break;
1782 }
1783
1784 Destroy_PostElement(SL) ;
1785
1786 Format_PostFooter(PSO_P, 0);
1787 }
1788