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 <string.h>
7 #include <math.h>
8 #include <map>
9 #include "GetDPConfig.h"
10 #include "ProData.h"
11 #include "DofData.h"
12 #include "F.h"
13 #include "Cal_Quantity.h"
14 #include "Cal_Value.h"
15 #include "MallocUtils.h"
16 #include "Message.h"
17 #include "ProParser.h"
18 
19 #if defined(HAVE_KERNEL)
20 #include "GeoData.h"
21 #include "Get_Geometry.h"
22 #include "Pos_FemInterpolation.h"
23 #include "Pos_Search.h"
24 #include "Get_FunctionValue.h"
25 #include "Pos_Format.h"
26 #endif
27 
28 extern struct Problem Problem_S ;
29 extern struct CurrentData Current ;
30 
31 #if defined(HAVE_KERNEL)
32 extern int TreatmentStatus ;
33 #else
34 const int TreatmentStatus = 0; // FIXME
35 #endif
36 
37 /* ------------------------------------------------------------------------ */
38 /*  G e t _ V a l u e O f E x p r e s s i o n                               */
39 /* ------------------------------------------------------------------------ */
40 
Get_ValueOfExpression(struct Expression * Expression_P,struct QuantityStorage * QuantityStorage_P0,double u,double v,double w,struct Value * Value,int NbrArguments,char * CallingExpressionName)41 void  Get_ValueOfExpression(struct Expression * Expression_P,
42                             struct QuantityStorage * QuantityStorage_P0,
43                             double u, double v, double w,
44                             struct Value * Value,
45                             int    NbrArguments,
46                             char   *CallingExpressionName)
47 {
48   static char *Flag_WarningUndefined = NULL ;
49 
50   switch (Expression_P->Type) {
51 
52   case CONSTANT :
53     if (Current.NbrHar == 1) {
54       Value->Val[0] = Expression_P->Case.Constant ;
55     }
56     else {
57       for (int k = 0 ; k < Current.NbrHar ; k += 2) {
58         Value->Val[MAX_DIM* k   ] = Expression_P->Case.Constant ;
59         Value->Val[MAX_DIM*(k+1)] = 0. ;
60       }
61     }
62     Value->Type = SCALAR ;
63     break ;
64 
65   case WHOLEQUANTITY :
66     Cal_WholeQuantity(Current.Element, QuantityStorage_P0,
67                       Expression_P->Case.WholeQuantity,
68                       u,v,w, -1, 0, Value, NbrArguments,
69                       CallingExpressionName ?
70                           CallingExpressionName : Expression_P->Name) ;
71     break ;
72 
73   case PIECEWISEFUNCTION :
74     struct Expression * NextExpression_P;
75     struct ExpressionPerRegion  * ExpressionPerRegion_P ;
76 
77     if (Current.Region == Expression_P->Case.PieceWiseFunction.NumLastRegion) {
78       NextExpression_P = Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion;
79     }
80     else {
81       if ((ExpressionPerRegion_P = (struct ExpressionPerRegion*)
82            List_PQuery(Expression_P->Case.PieceWiseFunction.ExpressionPerRegion,
83                        &Current.Region, fcmp_int))) {
84         Expression_P->Case.PieceWiseFunction.NumLastRegion = Current.Region ;
85         NextExpression_P =
86           Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion =
87           (struct Expression*)List_Pointer(Problem_S.Expression,
88                                            ExpressionPerRegion_P->ExpressionIndex) ;
89       }
90       else if (Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default >= 0) {
91         NextExpression_P =
92           (struct Expression*)
93           List_Pointer(Problem_S.Expression,
94                        Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default) ;
95       }
96       else {
97         NextExpression_P = NULL;
98         if(Current.Region == NO_REGION)
99           Message::Error("Function '%s' undefined in expressions without support",
100                          Expression_P->Name);
101         else
102           Message::Error("Function '%s' undefined in Region %d",
103                          Expression_P->Name, Current.Region);
104       }
105     }
106 
107     Get_ValueOfExpression
108       (NextExpression_P,
109        QuantityStorage_P0, u, v, w, Value, NbrArguments, Expression_P->Name) ;
110     break ;
111 
112   case PIECEWISEFUNCTION2 :
113     //    struct Expression * NextExpression_P;
114     struct ExpressionPerRegion2  * ExpressionPerRegion2_P ;
115     int twoRegions[2];
116 
117     if (Current.Region == Expression_P->Case.PieceWiseFunction2.NumLastRegion[0] &&
118         Current.SubRegion == Expression_P->Case.PieceWiseFunction2.NumLastRegion[1]) {
119       NextExpression_P = Expression_P->Case.PieceWiseFunction2.ExpressionForLastRegion;
120     }
121     else {
122       twoRegions[0] = Current.Region;
123       twoRegions[1] = Current.SubRegion;
124       if ((ExpressionPerRegion2_P = (struct ExpressionPerRegion2*)
125            List_PQuery(Expression_P->Case.PieceWiseFunction2.ExpressionPerRegion,
126                        twoRegions, fcmp_Integer2))) {
127         Expression_P->Case.PieceWiseFunction2.NumLastRegion[0] = Current.Region ;
128         Expression_P->Case.PieceWiseFunction2.NumLastRegion[1] = Current.SubRegion ;
129         NextExpression_P =
130           Expression_P->Case.PieceWiseFunction2.ExpressionForLastRegion =
131           (struct Expression*)List_Pointer(Problem_S.Expression,
132                                            ExpressionPerRegion2_P->ExpressionIndex) ;
133       }
134       else if (Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default >= 0) {
135         NextExpression_P =
136           (struct Expression*)
137           List_Pointer(Problem_S.Expression,
138                        Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default) ;
139       }
140       else {
141         NextExpression_P = NULL;
142         if(Current.Region == NO_REGION)
143           Message::Error("Function '%s' undefined in expressions without support",
144                          Expression_P->Name);
145         else
146           Message::Error("Function '%s' undefined in [ Region %d, SubRegion %d ]",
147                          Expression_P->Name, Current.Region, Current.SubRegion);
148       }
149     }
150 
151     Get_ValueOfExpression
152       (NextExpression_P,
153        QuantityStorage_P0, u, v, w, Value, NbrArguments, Expression_P->Name) ;
154     break ;
155 
156   case UNDEFINED_EXP :
157     if(!Flag_WarningUndefined || strcmp(Flag_WarningUndefined, Expression_P->Name)){
158       Message::Warning("Undefined expression '%s' (assuming zero)", Expression_P->Name) ;
159       Flag_WarningUndefined = Expression_P->Name;
160     }
161     Cal_ZeroValue(Value);
162     Value->Type = SCALAR ;
163     break;
164 
165   default :
166     Message::Error("Unknown type (%d) of Expression (%s)",
167                    Expression_P->Type, Expression_P->Name) ;
168     break;
169   }
170 }
171 
172 /* ------------------------------------------------------------------------ */
173 /*  G e t _ V a l u e O f E x p r e s s i o n B y I n d e x                 */
174 /* ------------------------------------------------------------------------ */
175 
Get_ValueOfExpressionByIndex(int Index_Expression,struct QuantityStorage * QuantityStorage_P0,double u,double v,double w,struct Value * Value)176 void Get_ValueOfExpressionByIndex(int Index_Expression,
177                                   struct QuantityStorage * QuantityStorage_P0,
178                                   double u, double v, double w,
179                                   struct Value * Value) {
180   Get_ValueOfExpression
181     ((struct Expression *)List_Pointer(Problem_S.Expression, Index_Expression),
182      QuantityStorage_P0,  u, v, w, Value) ;
183 }
184 
Is_ExpressionConstant(struct Expression * Expression_P)185 bool Is_ExpressionConstant(struct Expression * Expression_P)
186 {
187   if(Expression_P->Type == CONSTANT) return true;
188   if(Expression_P->Type == WHOLEQUANTITY){
189     for(int i = 0; i < List_Nbr(Expression_P->Case.WholeQuantity); i++){
190       struct WholeQuantity *WholeQuantity_P = (struct WholeQuantity*)
191         List_Pointer(Expression_P->Case.WholeQuantity, i);
192       if(WholeQuantity_P->Type != WQ_CONSTANT)
193         return false;
194     }
195     return true;
196   }
197   return false;
198 }
199 
200 /* ------------------------------------------------------------------------ */
201 /*  C a l _ S o l i d A n g l e                                             */
202 /* ------------------------------------------------------------------------ */
203 
Cal_SolidAngle(int Source,struct Element * Element,struct QuantityStorage * QuantityStorage,int Nbr_Dof,int Index,struct Value ** Stack)204 void Cal_SolidAngle(int Source, struct Element *Element,
205                     struct QuantityStorage *QuantityStorage,
206                     int Nbr_Dof, int Index,
207                     struct Value **Stack)
208 {
209 #if !defined(HAVE_KERNEL)
210   Message::Error("Cal_SolidAngle requires Kernel");
211 #else
212   struct Element     *Elt ;
213   struct Geo_Element *GeoElement ;
214   struct Geo_Node    *GeoNode1, *GeoNode2, *GeoNode3 ;
215 
216   double X, Y, Atan ;
217   int    In, TypEnt, NumNode, NbrElements, *NumElements ;
218   int    i, j ;
219 
220   if(Nbr_Dof != QuantityStorage->NbrElementaryBasisFunction){
221     Message::Error("Uncompatible Quantity (%s) in SolidAngle computation",
222                    QuantityStorage->DefineQuantity->Name);
223     return;
224   }
225 
226   if(Source){
227     Elt = Element->ElementSource ;
228     In  = Current.SourceIntegrationSupportIndex ;
229   }
230   else{
231     Elt = Element ;
232     In  = Current.IntegrationSupportIndex ;
233   }
234 
235   if (Elt->NumLastElementForSolidAngle == Elt->Num) {
236     for(i=0 ; i<Nbr_Dof ; i++) {
237       Stack[i][Index].Type = SCALAR ;
238       Stack[i][Index].Val[0] = Elt->angle[i] ;
239     }
240     return;
241   }
242 
243   for(i=0 ; i<Nbr_Dof ; i++){
244 
245     Stack[i][Index].Type = SCALAR ;
246 
247     TypEnt = ((struct Group*)
248               List_Pointer(Problem_S.Group,
249                            QuantityStorage->BasisFunction[i].
250                            BasisFunction->EntityIndex))->FunctionType ;
251 
252     if(TypEnt != NODESOF){
253 
254       if(Elt->Type == LINE){
255         Stack[i][Index].Val[0] = M_PI ;
256       }
257       else{
258         Stack[i][Index].Val[0] = 2.*M_PI ;
259       }
260 
261     }
262 
263     else{
264 
265       NumNode = Elt->GeoElement->
266         NumNodes[QuantityStorage->BasisFunction[i].NumEntityInElement] ;
267 
268       Geo_CreateNodesXElements(NumNode, In, &NbrElements, &NumElements) ;
269 
270       if(NbrElements != 2){
271         Message::Error("SolidAngle not done for incidence != 2 (%d)", NbrElements);
272         return;
273       }
274 
275       GeoNode2 = Geo_GetGeoNodeOfNum(NumNode) ;
276       GeoElement = Geo_GetGeoElementOfNum(abs(NumElements[0])) ;
277 
278       if(GeoElement->Type != LINE){
279         Message::Error("SolidAngle not done for Elements other than LINE");
280         return;
281       }
282 
283       if(NumElements[0]>0){
284         GeoNode1 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[0]) ;
285         GeoElement = Geo_GetGeoElementOfNum(abs(NumElements[1])) ;
286         GeoNode3 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[1]) ;
287 
288       }
289       else{
290         GeoNode3 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[1]) ;
291         GeoElement = Geo_GetGeoElementOfNum(NumElements[1]) ;
292         GeoNode1 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[0]) ;
293       }
294 
295       Y = (GeoNode1->y - GeoNode2->y) * (GeoNode3->x - GeoNode2->x) -
296           (GeoNode1->x - GeoNode2->x) * (GeoNode3->y - GeoNode2->y) ;
297       X = (GeoNode1->x - GeoNode2->x) * (GeoNode3->x - GeoNode2->x) +
298           (GeoNode1->y - GeoNode2->y) * (GeoNode3->y - GeoNode2->y) ;
299 
300       Atan = atan2(Y,X) ;
301 
302       Stack[i][Index].Val[0] = (Atan >= 0) ? Atan : (Atan+2.*M_PI) ;
303 
304       if(Message::GetVerbosity() > 5){
305         printf("Solid angle=%g node=%d, region=%s, elms=",
306                Stack[i][Index].Val[0] * 180/M_PI, NumNode,
307                ((struct Group*)List_Pointer(Problem_S.Group, In))->Name);
308         for(j=0 ; j<NbrElements ; j++) printf("%d ", NumElements[j]);
309         printf("\n");
310       }
311 
312     }
313 
314   }
315 
316   if (Elt->NumLastElementForSolidAngle != Elt->Num) {
317     Elt->NumLastElementForSolidAngle = Elt->Num ;
318     for(i=0 ; i<Nbr_Dof ; i++) Elt->angle[i] = Stack[i][Index].Val[0] ;
319   }
320 #endif
321 }
322 
323 /* ------------------------------------------------------------------------ */
324 /*  C a l _ W h o l e Q u a n t i t y                                       */
325 /* ------------------------------------------------------------------------ */
326 
327 #define CAST3V    void(*)(struct Value*, struct Value*, struct Value*)
328 #define CAST1V    void(*)(struct Value*)
329 #define CASTF2V   void(*)(struct Function*, struct Value*, struct Value*)
330 
331 // There can be at max one "Dof{op qty}" per WholeQuantity, but as
332 // many {op qty} as you want.
333 
334 static std::map<int, struct Value> ValueSaved ;
335 static std::map<std::string, struct Value> NamedValueSaved ;
336 
Cal_WholeQuantity(struct Element * Element,struct QuantityStorage * QuantityStorage_P0,List_T * WholeQuantity_L,double u,double v,double w,int DofIndexInWholeQuantity,int Nbr_Dof,struct Value DofValue[],int NbrArguments,char * ExpressionName)337 void Cal_WholeQuantity(struct Element * Element,
338                        struct QuantityStorage * QuantityStorage_P0,
339                        List_T * WholeQuantity_L,
340                        double u, double v, double w,
341                        int DofIndexInWholeQuantity,
342                        int Nbr_Dof, struct Value DofValue[],
343                        int NbrArguments, char *ExpressionName) {
344 
345   static int Flag_WarningMissSolForDt = 0 ;
346   static int Flag_WarningMissSolForTime_ntime = 0 ;
347   static int Flag_InfoForTime_ntime = 0 ;
348 
349   int     i_WQ, j, k, Flag_True, Index, DofIndex, Multi[MAX_STACK_SIZE] ;
350   int     Save_NbrHar, Save_Region, Type_Dimension, ntime ;
351   double  Save_Time, Save_TimeImag, Save_TimeStep, X, Y, Z, Order ;
352   double  Save_x, Save_y, Save_z ;
353 
354   struct WholeQuantity   *WholeQuantity_P0, *WholeQuantity_P ;
355   struct DofData         *Save_DofData ;
356   struct Solution        *Solution_P0, *Solution_PN ;
357 
358   struct Element* Save_CurrentElement ;
359 
360   // we could make this dynamic (with std::vector) to reduce stack usage, but
361   // the performance hit is important
362 
363   // ==> forcing a reduced size of stack for MH case...
364   // MAX_STACK_SIZE0 = 8 by default, 2 for MH
365   // segmentation violation and out of memory with high number of harmonics
366   // MAX_STACK_SIZE at least MAX_HAR_SIZE if harmonic function in formulation term
367   struct Value Stack[MAX_STACK_SIZE0][MAX_STACK_SIZE] ;
368 
369 
370   WholeQuantity_P0 = (struct WholeQuantity*)List_Pointer(WholeQuantity_L, 0) ;
371 
372   Index     = 0 ;
373   DofIndex  = -1 ;
374 
375   for (i_WQ = 0 ; i_WQ < List_Nbr(WholeQuantity_L) ; i_WQ++) {
376 
377     if(Index >= MAX_STACK_SIZE){
378       Message::Error("Stack size exceeded (%d>%d)", Index, MAX_STACK_SIZE);
379       return;
380     }
381 
382     WholeQuantity_P = WholeQuantity_P0 + i_WQ ;
383 
384     switch (WholeQuantity_P->Type) {
385 
386     case WQ_OPERATORANDQUANTITY : /* {op qty}  Dof{op qty}  BF{op qty} */
387       Save_Region = Current.Region ;
388       Save_CurrentElement = Current.Element ;
389       if (i_WQ != DofIndexInWholeQuantity){ /* Attention!!! || TreatmentStatus == STATUS_POS){  */
390 #if defined(HAVE_KERNEL)
391         Pos_FemInterpolation
392           (Element,
393            QuantityStorage_P0,
394            QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
395            WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
396            WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0,
397            u, v, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
398 #else
399         Message::Error("TODO Post_FemInterpolation");
400 #endif
401         Multi[Index] = 0 ;
402       }
403       else {
404         DofIndex = Index ;
405       }
406       Index++ ;
407       Current.Element = Save_CurrentElement ;
408       Current.Region = Save_Region ;
409       break ;
410 
411     case WQ_ORDER : /* Order[{qty}] */
412 #if defined(HAVE_KERNEL)
413       Order = Cal_InterpolationOrder
414         (Element, QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index) ;
415 #else
416         Message::Error("TODO Cal_InterpolationOrder");
417 #endif
418       for (k = 0 ; k < Current.NbrHar ; k += 2) {
419         Stack[0][Index].Val[MAX_DIM* k   ] = Order ;
420         Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ;
421       }
422       Stack[0][Index].Type = SCALAR ;
423       Multi[Index] = 0 ;
424       Index++ ;
425       break ;
426 
427     case WQ_OPERATORANDQUANTITYEVAL :
428       Save_Region = Current.Region ;
429       Save_CurrentElement = Current.Element ;
430       /* {op qty}[x,y,z], {op qty}[x,y,z,dimension]
431          or {op qty}[Vector[x,y,x],dimension]
432          or {op qty}[ntime] */
433       if (i_WQ != DofIndexInWholeQuantity || TreatmentStatus == STATUS_POS){
434         j = WholeQuantity_P->Case.OperatorAndQuantity.NbrArguments;
435         if (j == 2 || j == 3 || j == 4) {
436           if (j == 3 || j == 4) {
437             Index -= j ;
438             X = Stack[0][Index  ].Val[0] ;
439             Y = Stack[0][Index+1].Val[0] ;
440             Z = Stack[0][Index+2].Val[0] ;
441             if(j == 4)
442               Type_Dimension = (int)Stack[0][Index+3].Val[0] ;
443             else
444               Type_Dimension = -1 ;
445           }
446           else { /* j==2 */
447             Index -= j ;
448             X = Stack[0][Index  ].Val[0] ;
449             Y = Stack[0][Index  ].Val[1] ;
450             Z = Stack[0][Index  ].Val[2] ;
451             Type_Dimension = (int)Stack[0][Index+1].Val[0] ;
452           }
453 #if defined(HAVE_KERNEL)
454           Pos_FemInterpolation
455             (Element,
456              QuantityStorage_P0,
457              QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
458              WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
459              WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, Type_Dimension, 1,
460              u, v, w, X, Y, Z, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
461 #else
462           Message::Error("TODO Post_FemInterpolation");
463 #endif
464           Multi[Index] = 0 ;
465           Index++ ;
466         }
467         else if (j == 1) {
468           Index -= j ;
469           ntime = (int)Stack[0][Index].Val[0] ;
470 
471           for (k = 0 ; k < Current.NbrSystem ; k++){
472             if(!List_Nbr((Current.DofData_P0+k)->Solutions)) continue;
473             Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)->Solutions, 0);
474             if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) >= ntime){
475               ((Current.DofData_P0+k)->CurrentSolution) -= ntime ;
476               if (Flag_InfoForTime_ntime != List_Nbr((Current.DofData_P0+k)->Solutions)) {
477                 Message::Debug("Accessing solution from %d time steps ago", ntime);
478                 Message::Debug("  -> System %d/%d: TimeStep = %d, Time = %g + i * %g",
479                                k+1, Current.NbrSystem,
480                                (Current.DofData_P0+k)->CurrentSolution->TimeStep,
481                                (Current.DofData_P0+k)->CurrentSolution->Time,
482                                (Current.DofData_P0+k)->CurrentSolution->TimeImag);
483                 Flag_InfoForTime_ntime = List_Nbr((Current.DofData_P0+k)->Solutions);
484               }
485             }
486             else {
487               if (!Flag_WarningMissSolForTime_ntime) {
488                 Message::Warning("Missing solution for time step -%d computation (System #%d/%d)",
489                                  ntime, k+1, Current.NbrSystem);
490                 Flag_WarningMissSolForTime_ntime = 1 ;
491               }
492             }
493           }
494 
495 #if defined(HAVE_KERNEL)
496           Pos_FemInterpolation
497             (Element,
498              QuantityStorage_P0,
499              QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
500              WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
501              WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0,
502              u, v, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
503 #else
504           Message::Error("TODO Post_FemInterpolation");
505 #endif
506 
507           Multi[Index] = 0 ;
508           Index++ ;
509 
510           for (k = 0 ; k < Current.NbrSystem ; k++){
511             if(!List_Nbr((Current.DofData_P0+k)->Solutions)) continue;
512             Solution_PN = (struct Solution*)
513               List_Pointer((Current.DofData_P0+k)->Solutions,
514                            List_Nbr((Current.DofData_P0+k)->Solutions)-1);
515             if((Solution_PN - (Current.DofData_P0+k)->CurrentSolution) >= ntime)
516               ((Current.DofData_P0+k)->CurrentSolution) += ntime ;
517           }
518 
519         }
520         else
521           Message::Error("Explicit (x,y,z,time) evaluation not implemented");
522       }
523       else{
524         Message::Error("Explicit Dof{} evaluation out of context");
525       }
526       Current.Element = Save_CurrentElement ;
527       Current.Region = Save_Region ;
528       break ;
529 
530     case WQ_TRACE : /* Trace[WholeQuantity, Group] */
531       Save_Region = Current.Region ;
532 
533       if(!Element->ElementTrace){
534         Message::Error("Trace must act on discrete quantity (and not in post-processing)");
535         break;
536       }
537 
538       Current.Region = Element->ElementTrace->Region ;
539 
540       if(WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity >= 0){
541         Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0,
542                           WholeQuantity_P->Case.Trace.WholeQuantity,
543                           Current.ut, Current.vt, Current.wt,
544                           WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity,
545                           Nbr_Dof, DofValue, NbrArguments, ExpressionName) ;
546         DofIndexInWholeQuantity = DofIndex = Index ;
547       }
548       else{
549         Current.x = Current.y = Current.z = 0. ;
550         for (j = 0 ; j < Element->GeoElement->NbrNodes ; j++) {
551           Current.x += Element->x[j] * Element->n[j] ;
552           Current.y += Element->y[j] * Element->n[j] ;
553           Current.z += Element->z[j] * Element->n[j] ;
554         }
555 #if defined(HAVE_KERNEL)
556         xyz2uvwInAnElement(Element->ElementTrace, Current.x, Current.y, Current.z,
557                            &Current.ut, &Current.vt, &Current.wt) ;
558 #else
559         Message::Error("TODO xyz2uvwInAnElement");
560 #endif
561         for (j=0; j<NbrArguments; j++) {
562           Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
563           Multi[Index] = 0 ;
564           Index++ ;
565         }
566         Index -= NbrArguments ;
567         Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0,
568                           WholeQuantity_P->Case.Trace.WholeQuantity,
569                           Current.ut, Current.vt, Current.wt,
570                           -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
571       }
572 
573       Current.Region = Save_Region ;
574       Multi[Index] = 0 ;
575       Index++ ;
576       break ;
577 
578     case WQ_SOLIDANGLE : /* SolidAngle[{qty}] */
579       Cal_SolidAngle(0, Element, QuantityStorage_P0 +
580                      WholeQuantity_P->Case.OperatorAndQuantity.Index,
581                      Nbr_Dof, Index, (struct Value **)Stack);
582       Multi[Index] = 1 ;
583       Index++ ;
584       break ;
585 
586     case WQ_BINARYOPERATOR : /* + - * x / % ^ < > <= >= == != && || */
587       if (Index-2 != DofIndex && Index-1 != DofIndex){
588         if(!Multi[Index-1] && !Multi[Index-2])
589           ((CAST3V)WholeQuantity_P->Case.Operator.Function)
590             (&Stack[0][Index-2], &Stack[0][Index-1], &Stack[0][Index-2]) ;
591         else if(Multi[Index-1] && Multi[Index-2])
592           for(j=0 ; j<Nbr_Dof ; j++)
593             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
594               (&Stack[j][Index-2], &Stack[j][Index-1], &Stack[j][Index-2]) ;
595         else if(Multi[Index-2])
596           for(j=0 ; j<Nbr_Dof ; j++)
597             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
598               (&Stack[j][Index-2], &Stack[0][Index-1], &Stack[j][Index-2]) ;
599         else {
600           for(j=0 ; j<Nbr_Dof ; j++)
601             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
602               (&Stack[0][Index-2], &Stack[j][Index-1], &Stack[j][Index-2]) ;
603           Multi[Index-2] = 1 ;
604         }
605       }
606       else if (Index-1 == DofIndex) {
607         if(Multi[Index-2])
608           for (j = 0 ; j < Nbr_Dof ; j++)
609             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
610               (&Stack[j][Index-2], &DofValue[j], &DofValue[j]) ;
611         else
612           for (j = 0 ; j < Nbr_Dof ; j++)
613             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
614               (&Stack[0][Index-2], &DofValue[j], &DofValue[j]) ;
615         DofIndex-- ;
616       }
617       else { /* Index-2 == DofIndex */
618         if(Multi[Index-1])
619           for (j = 0 ; j < Nbr_Dof ; j++)
620             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
621               (&DofValue[j], &Stack[j][Index-1], &DofValue[j]) ;
622         else
623           for (j = 0 ; j < Nbr_Dof ; j++)
624             ((CAST3V)WholeQuantity_P->Case.Operator.Function)
625               (&DofValue[j], &Stack[0][Index-1], &DofValue[j]) ;
626       }
627       Index-- ;
628       break ;
629 
630     case WQ_UNARYOPERATOR : /* + - ! */
631       if (Index-1 == DofIndex)
632         for (j = 0 ; j < Nbr_Dof ; j++)
633           ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&DofValue[j]) ;
634       else if(Multi[Index-1])
635         for(j=0 ; j<Nbr_Dof ; j++)
636           ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[j][Index-1]) ;
637       else
638         ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[0][Index-1]) ;
639       break ;
640 
641       /* WARNING: all the rest assumes 0 multi status */
642 
643     case WQ_TEST :
644       Flag_True = (Stack[0][Index-1].Val[0] != 0.) ;
645       for (j=0; j<NbrArguments; j++) {
646         Cal_CopyValue(DofValue + j, &Stack[0][Index-1]) ;
647         Multi[Index-1] = 0 ;
648         Index++ ;
649       }
650       Index -= NbrArguments ;
651       Cal_WholeQuantity(Element, QuantityStorage_P0,
652                         (Flag_True) ? WholeQuantity_P->Case.Test.WholeQuantity_True :
653                                       WholeQuantity_P->Case.Test.WholeQuantity_False,
654                         u, v, w, -1, 0, &Stack[0][Index-1], NbrArguments, ExpressionName) ;
655       break ;
656 
657     case WQ_EXPRESSION :
658       Index -= WholeQuantity_P->Case.Expression.NbrArguments ;
659       Get_ValueOfExpression
660         ((struct Expression*)List_Pointer
661          (Problem_S.Expression, WholeQuantity_P->Case.Expression.Index),
662          QuantityStorage_P0, u, v, w, &Stack[0][Index],
663          WholeQuantity_P->Case.Expression.NbrArguments) ;
664       Multi[Index] = 0 ;
665       Index++ ;
666       break ;
667 
668     case WQ_BUILTINFUNCTION :
669       Index -= WholeQuantity_P->Case.Function.NbrArguments ;
670 
671       if (Index != DofIndex)
672         ((CASTF2V)WholeQuantity_P->Case.Function.Fct)
673           (&WholeQuantity_P->Case.Function, &Stack[0][Index], &Stack[0][Index]) ;
674       else /* Dof must be the only argument, only valid with linear functions */
675         for (j = 0 ; j < Nbr_Dof ; j++) {
676           Current.NumEntity = Current.NumEntities[j]; /* temp */
677           ((CASTF2V)WholeQuantity_P->Case.Function.Fct)
678             (&WholeQuantity_P->Case.Function, &DofValue[j], &DofValue[j]) ;
679         }
680 
681       Multi[Index] = 0 ;
682       Index++ ;
683       break ;
684 
685     case WQ_EXTERNBUILTINFUNCTION :
686       ((CASTF2V)WholeQuantity_P->Case.Function.Fct)
687         (&WholeQuantity_P->Case.Function, DofValue, &Stack[0][Index]) ;
688       Multi[Index] = 0 ;
689       Index++ ;
690       break ;
691 
692     case WQ_CONSTANT :
693       if (Current.NbrHar == 1) {
694         Stack[0][Index].Val[0] = WholeQuantity_P->Case.Constant ;
695       }
696       else {
697         for (k = 0 ; k < Current.NbrHar ; k += 2) {
698           Stack[0][Index].Val[MAX_DIM* k   ] = WholeQuantity_P->Case.Constant ;
699           Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ;
700         }
701       }
702       Stack[0][Index].Type = SCALAR ;
703       Multi[Index] = 0 ;
704       Index++ ;
705       break ;
706 
707     case WQ_CURRENTVALUE :
708       if (Current.NbrHar == 1) {
709         Stack[0][Index].Val[0] = *(WholeQuantity_P->Case.CurrentValue.Value) ;
710       }
711       else {
712         for (k = 0 ; k < Current.NbrHar ; k += 2) {
713           Stack[0][Index].Val[MAX_DIM* k   ] = *(WholeQuantity_P->Case.CurrentValue.Value) ;
714           Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ;
715         }
716       }
717       Stack[0][Index].Type = SCALAR ;
718       Multi[Index] = 0 ;
719       Index++ ;
720       break ;
721 
722     case WQ_ARGUMENT :
723       if (WholeQuantity_P->Case.Argument.Index > NbrArguments){
724         Message::Error("Function %s called with too few arguments.", ExpressionName);
725       }
726       Cal_CopyValue(DofValue + WholeQuantity_P->Case.Argument.Index - 1,
727                     &Stack[0][Index]) ;
728       Multi[Index] = 0 ;
729       Index++ ;
730       break ;
731 
732     case WQ_TIMEDERIVATIVE :
733       if (Current.TypeTime == TIME_GEAR) {
734         for (j=0; j<NbrArguments; j++) {
735           Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
736           Multi[Index] = 0 ;
737           Index++ ;
738         }
739         Index -= NbrArguments ;
740         Cal_WholeQuantity(Element, QuantityStorage_P0,
741                           WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
742                           u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName);
743 
744         for (k = 0 ; k < Current.NbrSystem ; k++)
745           (Current.DofData_P0+k)->Save_CurrentSolution =
746             (Current.DofData_P0+k)->CurrentSolution;
747         Save_TimeStep = Current.TimeStep ;
748         Save_Time = Current.Time ;
749         Save_TimeImag = Current.TimeImag ;
750 
751         for (int n = 0; n < Current.CorrOrder; n++) {
752 
753           for (k = 0 ; k < Current.NbrSystem ; k++){
754             Solution_P0 = (struct Solution*)List_Pointer
755               ((Current.DofData_P0+k)->Solutions, 0);
756             if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) > 0)
757               ((Current.DofData_P0+k)->CurrentSolution) -= 1 ;
758             else{
759               Message::Error("Too few solutions for Dt with Gear's method");
760               break;
761             }
762           }
763 
764           Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ;
765           Current.Time = Current.DofData->CurrentSolution->Time ;
766           Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ;
767 
768           for (j=0; j<NbrArguments; j++) {
769             Cal_CopyValue(DofValue + j, &Stack[0][Index+1]) ;
770             Multi[Index+1] = 0 ;
771             Index++ ;
772           }
773           Index -= NbrArguments ;
774           Cal_WholeQuantity(Element, QuantityStorage_P0,
775                             WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
776                             u, v, w, -1, 0, &Stack[0][Index+1], NbrArguments, ExpressionName);
777           Cal_AddMultValue(&Stack[0][Index], &Stack[0][Index+1],
778                            -Current.aCorrCoeff[n], &Stack[0][Index]);
779 
780         }
781 
782         Cal_MultValue(&Stack[0][Index], 1./(Current.bCorrCoeff*Current.DTime),
783                       &Stack[0][Index]);
784 
785         for (k = 0 ; k < Current.NbrSystem ; k++)
786           (Current.DofData_P0+k)->CurrentSolution =
787               (Current.DofData_P0+k)->Save_CurrentSolution;
788         Current.TimeStep = Save_TimeStep ;
789         Current.Time = Save_Time ;
790         Current.TimeImag = Save_TimeImag ;
791       }
792       else if (Current.NbrHar == 1) {
793         for (j=0; j<NbrArguments; j++) {
794           Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
795           Multi[Index] = 0 ;
796           Index++ ;
797         }
798         Index -= NbrArguments ;
799         Cal_WholeQuantity(Element, QuantityStorage_P0,
800                           WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
801                           u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
802 
803         for (k = 0 ; k < Current.NbrSystem ; k++){
804           (Current.DofData_P0+k)->Save_CurrentSolution =
805            (Current.DofData_P0+k)->CurrentSolution;
806 
807           if(List_Nbr((Current.DofData_P0+k)->Solutions) > 1){
808             Solution_P0 = (struct Solution*)List_Pointer
809               ((Current.DofData_P0+k)->Solutions, 0);
810             if ((Current.DofData_P0+k)->CurrentSolution != Solution_P0)
811               ((Current.DofData_P0+k)->CurrentSolution) -- ;
812           }
813           else {
814             if (!Flag_WarningMissSolForDt) {
815               Message::Warning("Missing solution for time derivative computation "
816                                "(Sys#%d/%d)", k+1, Current.NbrSystem);
817               Flag_WarningMissSolForDt = 1 ;
818             }
819           }
820         }
821 
822         Save_TimeStep = Current.TimeStep ;
823         Save_Time = Current.Time ;
824         Save_TimeImag = Current.TimeImag ;
825 
826         Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ;
827         Current.Time = Current.DofData->CurrentSolution->Time ;
828         Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ;
829         for (j=0; j<NbrArguments; j++) {
830           Cal_CopyValue(DofValue + j, &Stack[0][Index+1]) ;
831           Multi[Index+1] = 0 ;
832           Index++ ;
833         }
834         Index -= NbrArguments ;
835         Cal_WholeQuantity(Element, QuantityStorage_P0,
836                           WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
837                           u, v, w, -1, 0, &Stack[0][Index+1], NbrArguments, ExpressionName) ;
838         Cal_SubstractValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) ;
839         Stack[0][Index+1].Val[0] = Save_Time - Current.Time ;
840         Stack[0][Index+1].Type = SCALAR ;
841         if(Stack[0][Index+1].Val[0])
842           Cal_DivideValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) ;
843         else
844           Cal_ZeroValue(&Stack[0][Index]);
845 
846         for (k = 0 ; k < Current.NbrSystem ; k++)
847           (Current.DofData_P0+k)->CurrentSolution =
848             (Current.DofData_P0+k)->Save_CurrentSolution;
849         Current.TimeStep = Save_TimeStep ;
850         Current.Time = Save_Time ;
851         Current.TimeImag = Save_TimeImag ;
852       }
853       else {
854         for (j=0; j<NbrArguments; j++) {
855           Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
856           Multi[Index] = 0 ;
857           Index++ ;
858         }
859         Index -= NbrArguments ;
860         Cal_WholeQuantity(Element, QuantityStorage_P0,
861                           WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
862                           u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
863         for (k = 0 ; k < Current.NbrHar ; k += 2) {
864           Stack[0][Index+1].Val[MAX_DIM* k   ] = 0. ;
865           Stack[0][Index+1].Val[MAX_DIM*(k+1)] = Current.DofData->Val_Pulsation[k/2] ;
866         }
867         Stack[0][Index+1].Type = SCALAR ;
868         Cal_ProductValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) ;
869       }
870       Multi[Index] = 0 ;
871       Index++ ;
872       break ;
873 
874     case WQ_ATANTERIORTIMESTEP :
875       ntime = WholeQuantity_P->Case.AtAnteriorTimeStep.TimeStep ;
876 
877       for (k = 0 ; k < Current.NbrSystem ; k++){
878         Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)->Solutions, 0);
879         if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) >= ntime){
880           ((Current.DofData_P0+k)->CurrentSolution) -= ntime ;
881           if (Flag_InfoForTime_ntime != List_Nbr((Current.DofData_P0+k)->Solutions)) {
882             Message::Info("Accessing solution from %d time steps ago", ntime);
883             Message::Info("  -> System %d/%d: TimeStep = %d, Time = %g + i * %g",
884                           k+1, Current.NbrSystem,
885                           (Current.DofData_P0+k)->CurrentSolution->TimeStep,
886                           (Current.DofData_P0+k)->CurrentSolution->Time,
887                           (Current.DofData_P0+k)->CurrentSolution->TimeImag);
888             Flag_InfoForTime_ntime = List_Nbr((Current.DofData_P0+k)->Solutions);
889           }
890         }
891         else {
892           if (!Flag_WarningMissSolForTime_ntime) {
893             Message::Warning("Missing solution for time step -%d computation "
894                              "(System #%d/%d)", ntime, k+1, Current.NbrSystem);
895             Flag_WarningMissSolForTime_ntime = 1 ;
896           }
897         }
898       }
899 
900       Save_TimeStep = Current.TimeStep ;
901       Save_Time = Current.Time ;
902       Save_TimeImag = Current.TimeImag ;
903       Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ;
904       Current.Time = Current.DofData->CurrentSolution->Time ;
905       Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ;
906 
907       for (j=0; j<NbrArguments; j++) {
908         Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
909         Multi[Index] = 0 ;
910         Index++ ;
911       }
912       Index -= NbrArguments ;
913       Cal_WholeQuantity(Element, QuantityStorage_P0,
914                         WholeQuantity_P->Case.AtAnteriorTimeStep.WholeQuantity,
915                         u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
916 
917       Current.TimeStep = Save_TimeStep ;
918       Current.Time = Save_Time ;
919       Current.TimeImag = Save_TimeImag ;
920 
921       for (k = 0 ; k < Current.NbrSystem ; k++){
922         Solution_PN = (struct Solution*)
923           List_Pointer((Current.DofData_P0+k)->Solutions,
924                        List_Nbr((Current.DofData_P0+k)->Solutions)-1);
925         if((Solution_PN - (Current.DofData_P0+k)->CurrentSolution) >= ntime)
926           ((Current.DofData_P0+k)->CurrentSolution) += ntime ;
927       }
928 
929       Multi[Index] = 0 ;
930       Index++ ;
931       break ;
932 
933 
934     case WQ_MAXOVERTIME :
935       if (Current.NbrHar == 1) {
936         double time_init = WholeQuantity_P->Case.MaxOverTime.TimeInit;
937         double time_final = WholeQuantity_P->Case.MaxOverTime.TimeFinal;
938         /*
939         for (k = 0 ; k < Current.NbrSystem ; k++)
940           (Current.DofData_P0+k)->Save_CurrentSolution =
941            (Current.DofData_P0+k)->CurrentSolution;
942         */
943         Save_TimeStep = Current.TimeStep ;
944         Save_Time = Current.Time ;
945         Save_TimeImag = Current.TimeImag ;
946 
947         for (j=0; j<NbrArguments; j++) {
948           Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
949           Multi[Index] = 0 ;
950           Index++ ;
951         }
952         Index -= NbrArguments ;
953 
954         double val_maxInTime = -1.e99;
955 
956         for (j=1; j<List_Nbr((Current.DofData)->Solutions); j++) {
957 
958          Current.DofData->CurrentSolution =
959             (struct Solution*)List_Pointer((Current.DofData)->Solutions, j);
960 
961           //++++ Add: also for other systems!
962 
963           Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ;
964           Current.Time = Current.DofData->CurrentSolution->Time ;
965           Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ;
966 
967           //++++ test to do more accurately!
968           if (Current.Time >= time_init && Current.Time <= time_final) {
969             Cal_WholeQuantity(Element, QuantityStorage_P0,
970                               WholeQuantity_P->Case.MaxOverTime.WholeQuantity,
971                               u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
972 
973             if (Stack[0][Index].Type == SCALAR) {
974               if (Stack[0][Index].Val[0] > val_maxInTime){
975                 val_maxInTime = Stack[0][Index].Val[0];
976               }
977             }
978             else {
979               Message::Error("MaxOverTime can only be applied on scalar values") ;
980             }
981           }
982         }
983         Stack[0][Index].Val[0] = val_maxInTime;
984         /*
985         for (k = 0 ; k < Current.NbrSystem ; k++)
986           (Current.DofData_P0+k)->CurrentSolution =
987             (Current.DofData_P0+k)->Save_CurrentSolution;
988         */
989         Current.TimeStep = Save_TimeStep ;
990         Current.Time = Save_Time ;
991         Current.TimeImag = Save_TimeImag ;
992 
993         Multi[Index] = 0 ;
994         Index++ ;
995       }
996       else {
997         Message::Error("MaxOverTime can only be used in time domain") ;
998         break;
999       }
1000       break ;
1001 
1002     case WQ_FOURIERSTEINMETZ :
1003       if (Current.NbrHar == 1) {
1004         double time_init = WholeQuantity_P->Case.FourierSteinmetz.TimeInit;
1005         double time_final = WholeQuantity_P->Case.FourierSteinmetz.TimeFinal;
1006         int nbrFrequencyInFormula = WholeQuantity_P->Case.FourierSteinmetz.NbrFrequency;
1007         double exponent_f = WholeQuantity_P->Case.FourierSteinmetz.Exponent_f;
1008         double exponent_b = WholeQuantity_P->Case.FourierSteinmetz.Exponent_b;
1009 
1010         /*
1011         for (k = 0 ; k < Current.NbrSystem ; k++)
1012           (Current.DofData_P0+k)->Save_CurrentSolution =
1013            (Current.DofData_P0+k)->CurrentSolution;
1014         */
1015         Save_TimeStep = Current.TimeStep ;
1016         Save_Time = Current.Time ;
1017         Save_TimeImag = Current.TimeImag ;
1018 
1019         for (j=0; j<NbrArguments; j++) {
1020           Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
1021           Multi[Index] = 0 ;
1022           Index++ ;
1023         }
1024         Index -= NbrArguments ;
1025 
1026         int i_Solution_init = -1, i_Solution_final = -1;
1027         int NbrTimeStep, Size=-1;
1028 
1029         for (j=0; j<List_Nbr((Current.DofData)->Solutions); j++) {
1030           Current.DofData->CurrentSolution =
1031             (struct Solution*)List_Pointer((Current.DofData)->Solutions, j);
1032           Current.Time = Current.DofData->CurrentSolution->Time ;
1033           if (Current.Time >= time_init && i_Solution_init < 0)
1034             i_Solution_init = j;
1035           if (Current.Time <= time_final) {
1036             i_Solution_final = j;
1037           }
1038           if (Current.Time > time_final) {
1039             break;
1040           }
1041         }
1042         NbrTimeStep = i_Solution_final-i_Solution_init+1;
1043         if (NbrTimeStep < 2)
1044           Message::Error("Wrong time interval in Function FourierSteinmetz (%d,%d)",
1045                          i_Solution_init, i_Solution_final) ;
1046 
1047         double *Times = (double *)Malloc(NbrTimeStep*sizeof(double));
1048         struct Value *TmpValues = (struct Value *)Malloc(NbrTimeStep*sizeof(struct Value));
1049 
1050         for (j=0; j<NbrTimeStep; j++) {
1051           Current.DofData->CurrentSolution =
1052             (struct Solution*)List_Pointer((Current.DofData)->Solutions,
1053                                            i_Solution_init+j);
1054 
1055           //++++ Add: also for other systems!
1056 
1057           Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ;
1058           Current.Time = Current.DofData->CurrentSolution->Time ;
1059           Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ;
1060 
1061           Cal_WholeQuantity(Element, QuantityStorage_P0,
1062                             WholeQuantity_P->Case.MaxOverTime.WholeQuantity,
1063                             u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
1064           Times[j] = Current.Time ;
1065           Cal_CopyValue(&Stack[0][Index], &TmpValues[j]) ;
1066           if (Stack[0][Index].Type == SCALAR)
1067             Size = 1;
1068           else if (Stack[0][Index].Type == VECTOR)
1069             Size = 3;
1070           else
1071             Message::Error("FourierSteinmetz can only be applied on scalar or vector values") ;
1072         }
1073 
1074         // FourierTransform
1075         int NbrFreq ;
1076         double *Frequencies;
1077         struct Value *FourierValues;
1078 
1079 #if defined(HAVE_KERNEL)
1080         Pos_FourierTransform(NbrTimeStep, 1, Times, TmpValues, Size,
1081                              2, nbrFrequencyInFormula,
1082                              &NbrFreq, &Frequencies, &FourierValues);
1083 #else
1084         Message::Error("TODO Pos_FourierTransform");
1085 #endif
1086         /*
1087           we calculate the Sum for all frequencies of frequency_i^exponent_f *
1088           b_i^exponent_b
1089         */
1090 
1091         if (nbrFrequencyInFormula > NbrFreq-1)
1092           Message::Error("FourierSteinmetz: too many frequencies asked "
1093                          "(%d asked and only %d available)",
1094                          nbrFrequencyInFormula, NbrFreq-1) ;
1095 
1096         double val=0.;
1097         for (j=0; j<nbrFrequencyInFormula; j++) {
1098           if (Size==1) {
1099             val += pow(Frequencies[j+1], exponent_f) *
1100               pow(FourierValues[j+1].Val[0], exponent_b);
1101           }
1102           else {
1103             val +=
1104               pow(Frequencies[j+1], exponent_f) *
1105               pow(sqrt(FourierValues[j+1].Val[0]*FourierValues[j+1].Val[0]
1106                        + FourierValues[j+1].Val[1]*FourierValues[j+1].Val[1]
1107                        + FourierValues[j+1].Val[2]*FourierValues[j+1].Val[2]),
1108                   exponent_b);
1109           }
1110         }
1111         Stack[0][Index].Type = SCALAR;
1112         Stack[0][Index].Val[0] = val;
1113 
1114         Free(Frequencies);
1115         Free(FourierValues);
1116 
1117         Free(Times); Free(TmpValues) ;
1118 
1119         /*
1120         for (k = 0 ; k < Current.NbrSystem ; k++)
1121           (Current.DofData_P0+k)->CurrentSolution =
1122             (Current.DofData_P0+k)->Save_CurrentSolution;
1123         */
1124         Current.TimeStep = Save_TimeStep ;
1125         Current.Time = Save_Time ;
1126         Current.TimeImag = Save_TimeImag ;
1127 
1128         Multi[Index] = 0 ;
1129         Index++ ;
1130       }
1131       else {
1132         Message::Error("FourierSteinmetz can only be used in time domain") ;
1133         break;
1134       }
1135       break ;
1136 
1137     case WQ_MHTRANSFORM :
1138       if(Current.NbrHar == 1){
1139         Message::Error("MHTransform can only be used in complex (multi-harmonic)"
1140                        " calculations") ;
1141         break;
1142       }
1143 
1144       for (j=0; j<NbrArguments; j++) {
1145         Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
1146         Multi[Index] = 0 ;
1147         Index++ ;
1148       }
1149       Index -= NbrArguments ;
1150       {
1151         int N = List_Nbr(WholeQuantity_P->Case.MHTransform.WholeQuantity_L);
1152         std::vector<struct Value> MH_Values(N);
1153         for(int j = 0; j < N; j++){
1154           List_T *WQ; List_Read(WholeQuantity_P->Case.MHTransform.WholeQuantity_L, j, &WQ);
1155           Cal_WholeQuantity(Element, QuantityStorage_P0, WQ, u, v, w, -1, 0,
1156                             &MH_Values[j], NbrArguments, ExpressionName) ;
1157         }
1158         MHTransform(Element, QuantityStorage_P0, u, v, w, MH_Values,
1159                     (struct Expression *)List_Pointer(Problem_S.Expression,
1160                                                       WholeQuantity_P->Case.MHTransform.Index),
1161                     WholeQuantity_P->Case.MHTransform.NbrPoints, Stack[0][Index]) ;
1162       }
1163       Multi[Index] = 0 ;
1164       Index++ ;
1165       break ;
1166 
1167     case WQ_CAST :
1168       /* This should be changed... */
1169       Save_NbrHar = Current.NbrHar ;
1170       Save_DofData = Current.DofData ;
1171 
1172       if (!WholeQuantity_P->Case.Cast.NbrHar){
1173         Current.DofData =
1174           ((struct FunctionSpace *)
1175            List_Pointer(Problem_S.FunctionSpace,
1176                         WholeQuantity_P->Case.Cast.FunctionSpaceIndexForType))
1177           ->DofData ;
1178 
1179         Current.NbrHar = Current.DofData->NbrHar ;
1180       }
1181       else{
1182         Current.NbrHar = WholeQuantity_P->Case.Cast.NbrHar ;
1183       }
1184 
1185       for (j=0; j<NbrArguments; j++) {
1186         Cal_CopyValue(DofValue + j, &Stack[0][Index]) ;
1187         Multi[Index] = 0 ;
1188         Index++ ;
1189       }
1190       Index -= NbrArguments ;
1191       Cal_WholeQuantity(Element, QuantityStorage_P0,
1192                         WholeQuantity_P->Case.Cast.WholeQuantity,
1193                         u, v, w, -1, 0, &Stack[0][Index], NbrArguments, ExpressionName) ;
1194 
1195       if (Current.NbrHar < Save_NbrHar)  /* ne plus a completer ...?? */
1196         Cal_SetZeroHarmonicValue(&Stack[0][Index], Save_NbrHar) ;
1197 
1198       Current.NbrHar = Save_NbrHar ;
1199       Current.DofData = Save_DofData ;
1200       Multi[Index] = 0 ;
1201       Index++ ;
1202       break ;
1203 
1204     case WQ_CHANGECURRENTPOSITION :
1205       Save_x = Current.x ; Save_y = Current.y ; Save_z = Current.z ;
1206 
1207       Current.x = Stack[0][Index-1].Val[0] ;
1208       Current.y = Stack[0][Index-1].Val[1] ;
1209       Current.z = Stack[0][Index-1].Val[2] ;
1210 
1211       for (j=0; j<NbrArguments; j++) {
1212         Cal_CopyValue(DofValue + j, &Stack[0][Index-1]) ;
1213         Multi[Index-1] = 0 ;
1214         Index++ ;
1215       }
1216       Index -= NbrArguments ;
1217       Cal_WholeQuantity(Element, QuantityStorage_P0,
1218                         WholeQuantity_P->Case.ChangeCurrentPosition.WholeQuantity,
1219                         u, v, w, -1, 0, &Stack[0][Index-1], NbrArguments, ExpressionName) ;
1220 
1221       Current.x = Save_x ; Current.y = Save_y ; Current.z = Save_z ;
1222       break ;
1223 
1224     case WQ_SAVEVALUE :
1225       Cal_CopyValue(&Stack[0][Index-1],
1226                     &ValueSaved[WholeQuantity_P->Case.SaveValue.Index]) ;
1227       break ;
1228 
1229     case WQ_VALUESAVED :
1230       if(ValueSaved.count(WholeQuantity_P->Case.ValueSaved.Index))
1231         Cal_CopyValue(&ValueSaved[WholeQuantity_P->Case.ValueSaved.Index],
1232                       &Stack[0][Index]) ;
1233       else{
1234         if(TreatmentStatus != STATUS_PRE)
1235           Message::Warning("Empty register %d: assuming zero value",
1236                            WholeQuantity_P->Case.ValueSaved.Index + 1);
1237         Cal_ZeroValue(&Stack[0][Index]);
1238         Stack[0][Index].Type = SCALAR ;
1239       }
1240       Multi[Index] = 0 ;
1241       Index++ ;
1242       break ;
1243 
1244     case WQ_SAVENAMEDVALUE :
1245       Cal_CopyValue(&Stack[0][Index-1],
1246                     &NamedValueSaved[WholeQuantity_P->Case.NamedValue.Name]) ;
1247       break ;
1248 
1249     case WQ_NAMEDVALUESAVED :
1250       if(NamedValueSaved.count(WholeQuantity_P->Case.NamedValue.Name))
1251         Cal_CopyValue(&NamedValueSaved[WholeQuantity_P->Case.NamedValue.Name],
1252                       &Stack[0][Index]) ;
1253       else{
1254         if(TreatmentStatus != STATUS_PRE)
1255           Message::Warning("Unknown current value '$%s': assuming zero value",
1256                            WholeQuantity_P->Case.NamedValue.Name);
1257         Cal_ZeroValue(&Stack[0][Index]);
1258         Stack[0][Index].Type = SCALAR ;
1259       }
1260       Multi[Index] = 0 ;
1261       Index++ ;
1262       break ;
1263 
1264     case WQ_SHOWVALUE :
1265       if (Index-1 == DofIndex) {
1266         for(j=0 ; j<Nbr_Dof ; j++){
1267           fprintf(stderr, "##%d Dof %d ", WholeQuantity_P->Case.ShowValue.Index, j+1);
1268           Show_Value(&DofValue[j]);
1269         }
1270       }
1271       else {
1272         fprintf(stderr, "##%d ", WholeQuantity_P->Case.ShowValue.Index);
1273         Show_Value(&Stack[0][Index-1]);
1274       }
1275       break ;
1276 
1277     default :
1278       Message::Error("Unknown type of WholeQuantity (%d)", WholeQuantity_P->Type);
1279       break;
1280     }
1281   }
1282 
1283   if (DofIndexInWholeQuantity < 0) Cal_CopyValue(&Stack[0][0], &DofValue[0]) ;
1284 }
1285 
1286 /* ------------------------------------------------------------------------ */
1287 /*  C a l _ S t o r e I n R e g i s t e r                                   */
1288 /* ------------------------------------------------------------------------ */
1289 
Cal_StoreInRegister(struct Value * Value,int RegisterIndex)1290 void Cal_StoreInRegister(struct  Value  *Value, int RegisterIndex)
1291 {
1292   Cal_CopyValue(Value, &ValueSaved[RegisterIndex]) ;
1293 }
1294 
1295 /* ------------------------------------------------------------------------ */
1296 /*  C a l _ S t o r e I n V a r i a b l e                                   */
1297 /* ------------------------------------------------------------------------ */
1298 
Cal_StoreInVariable(struct Value * Value,const char * name)1299 void Cal_StoreInVariable(struct  Value  *Value, const char *name)
1300 {
1301   Cal_CopyValue(Value, &NamedValueSaved[name]) ;
1302   Export_Value(Value, GetDPNumbers[name], 0, false); // don't append
1303 }
1304 
Cal_GetValueSaved(struct Value * Value,const char * name)1305 void Cal_GetValueSaved(struct  Value  *Value, const char *name)
1306 {
1307   if(NamedValueSaved.count(name))
1308     Cal_CopyValue(&NamedValueSaved[name], Value) ;
1309   else{
1310     Message::Warning("Unknown current value '$%s': assuming zero value", name);
1311     Cal_ZeroValue(Value);
1312     Value->Type = SCALAR ;
1313   }
1314 }
1315 
Get_AllValueSaved()1316 std::map<std::string, struct Value> &Get_AllValueSaved()
1317 {
1318   return NamedValueSaved;
1319 }
1320