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 <stdlib.h>
7 #include <string.h>
8 #include <math.h>
9 #include "GetDPConfig.h"
10 #include "ProData.h"
11 #include "ProDefine.h"
12 #include "ProParser.h"
13 #include "F.h"
14 #include "Message.h"
15 #include "OS.h"
16 #include "Cal_Value.h"
17 #include "Cal_Quantity.h"
18 
19 extern struct CurrentData Current ;
20 
F_Printf(F_ARG)21 void F_Printf(F_ARG)
22 {
23   Print_Value(A) ;
24   printf("\n") ;
25 }
26 
F_Rand(F_ARG)27 void F_Rand(F_ARG)
28 {
29   int     k;
30 
31   if(A->Type != SCALAR)
32     Message::Error("Non scalar argument for function 'Rand");
33 
34   V->Val[0] = A->Val[0] * (double)rand() / (double)RAND_MAX;
35 
36   if (Current.NbrHar != 1){
37     V->Val[MAX_DIM] = 0. ;
38     for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
39       V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
40   }
41 
42   V->Type = SCALAR ;
43 }
44 
F_CompElementNum(F_ARG)45 void F_CompElementNum (F_ARG)
46 {
47   if(!Current.Element || !Current.ElementSource)
48     Message::Error("Uninitialized Element in 'F_CompElementNum'");
49 
50   V->Type = SCALAR ;
51   V->Val[0] = (Current.Element->Num == Current.ElementSource->Num) ;
52 }
53 
F_ElementNum(F_ARG)54 void F_ElementNum (F_ARG)
55 {
56   if(!Current.Element)
57     Message::Error("Uninitialized Element in 'F_ElementNum'");
58 
59   V->Type = SCALAR ;
60   V->Val[0] = Current.Element->Num ;
61 }
62 
F_QuadraturePointIndex(F_ARG)63 void F_QuadraturePointIndex (F_ARG)
64 {
65   V->Type = SCALAR ;
66   V->Val[0] = Current.QuadraturePointIndex ;
67 }
68 
F_GetCpuTime(F_ARG)69 void F_GetCpuTime (F_ARG)
70 {
71   double s = 0.;
72   long mem = 0;
73   GetResources(&s, &mem);
74   V->Type = SCALAR ;
75   V->Val[0] = s ;
76 }
77 
F_GetWallClockTime(F_ARG)78 void F_GetWallClockTime (F_ARG)
79 {
80   V->Type = SCALAR ;
81   V->Val[0] = Message::GetWallClockTime() ;
82 }
83 
F_GetMemory(F_ARG)84 void F_GetMemory (F_ARG)
85 {
86   double s = 0.;
87   long mem = 0;
88   GetResources(&s, &mem);
89   double val = mem / 1024. / 1024.;
90   V->Type = SCALAR ;
91   V->Val[0] = val ;
92 }
93 
F_SetNumberRunTime(F_ARG)94 void F_SetNumberRunTime (F_ARG)
95 {
96   double val = A->Val[0];
97   int type = A->Type;
98 
99   for (int k = 0; k < Current.NbrHar; k++)
100     V->Val[MAX_DIM * k] = 0. ;
101   V->Type = SCALAR;
102   if(type != SCALAR){
103     Message::Error("Non scalar argument for function 'SetNumberRunTime'");
104     return;
105   }
106   if(!Fct->String){
107     Message::Error("Missing ONELAB variable name: use SetNumberRunTime[value]{\"name\"}");
108     return;
109   }
110 
111   Message::SetOnelabNumber(Fct->String, val);
112   V->Val[0] = val ;
113 }
114 
F_SetNumberRunTimeWithChoices(F_ARG)115 void F_SetNumberRunTimeWithChoices (F_ARG)
116 {
117   double val = A->Val[0];
118   int type = A->Type;
119 
120   for (int k = 0; k < Current.NbrHar; k++)
121     V->Val[MAX_DIM * k] = 0. ;
122   V->Type = SCALAR;
123   if(type != SCALAR){
124     Message::Error("Non scalar argument for function 'SetNumberRunTime'");
125     return;
126   }
127   if(!Fct->String){
128     Message::Error("Missing ONELAB variable name: use SetNumberRunTime[value]{\"name\"}");
129     return;
130   }
131 
132   std::vector<double> v(1, val);
133   Message::AddOnelabNumberChoice(Fct->String, v);
134   V->Val[0] = val ;
135 }
136 
F_GetNumberRunTime(F_ARG)137 void F_GetNumberRunTime (F_ARG)
138 {
139   if(Fct->NbrArguments){
140     Cal_CopyValue(A, V);
141   }
142   else{
143     for (int k = 0; k < Current.NbrHar; k++)
144       V->Val[MAX_DIM * k] = 0. ;
145     V->Type = SCALAR;
146   }
147   if(!Fct->String){
148     Message::Error("Missing ONELAB variable name: use GetNumberRunTime[]{\"name\"}");
149     return;
150   }
151   if(Message::UseOnelab())
152     V->Val[0] = Message::GetOnelabNumber(Fct->String);
153 }
154 
F_SetVariable(F_ARG)155 void F_SetVariable (F_ARG)
156 {
157   if(!Fct->String){
158     Message::Error("Missing runtime variable name: use SetVariable[...]{$name}");
159     return;
160   }
161 
162   if(Fct->NbrArguments < 1){
163     Message::Error("Missing argument in SetVariable[...]{$name}");
164     return;
165   }
166 
167   Cal_CopyValue(A, V);
168 
169   char tmp[256];
170   strcpy(tmp, Fct->String);
171   for(int i = 1; i < Fct->NbrArguments; i++){
172     if((A+i)->Type != SCALAR){
173       Message::Error("Non-scalar argument in SetVariable");
174       return;
175     }
176     char tmp2[256];
177     sprintf(tmp2, "_%g", (A+i)->Val[0]);
178     strcat(tmp, tmp2);
179   }
180   Cal_StoreInVariable(V, tmp);
181 }
182 
F_SetCumulativeVariable(F_ARG)183 void F_SetCumulativeVariable (F_ARG)
184 {
185   if(Fct->NbrArguments < 1){
186     Message::Error("Missing argument in SetCumulativeVariable[...]{$name}");
187     return;
188   }
189 
190   Cal_CopyValue(A, V);
191 
192   char tmp[256];
193   strcpy(tmp, Fct->String);
194   for(int i = 1; i < Fct->NbrArguments; i++){
195     if((A+i)->Type != SCALAR){
196       Message::Error("Non-scalar argument in SetCumulativeVariable");
197       return;
198     }
199     char tmp2[256];
200     sprintf(tmp2, "_%g", (A+i)->Val[0]);
201     strcat(tmp, tmp2);
202   }
203 
204   struct Value V_saved ;
205   Cal_GetValueSaved(&V_saved, tmp);
206   Cal_AddValue(V, &V_saved, V);
207 
208   Cal_StoreInVariable(V, tmp);
209 }
210 
F_GetVariable(F_ARG)211 void F_GetVariable (F_ARG)
212 {
213   if(!Fct->String){
214     Message::Error("Missing runtime variable name: use GetVariable[...]{$name}");
215     return;
216   }
217 
218   char tmp[256];
219   strcpy(tmp, Fct->String);
220   for(int i = 0; i < Fct->NbrArguments; i++){
221     if((A+i)->Type != SCALAR){
222       Message::Error("Non-scalar argument in GetVariable");
223       return;
224     }
225     char tmp2[256];
226     sprintf(tmp2, "_%g", (A+i)->Val[0]);
227     strcat(tmp, tmp2);
228   }
229 
230   Cal_GetValueSaved(V, tmp);
231 }
232 
F_ValueFromTable(F_ARG)233 void F_ValueFromTable (F_ARG)
234 {
235   if(!Fct->String){
236     Message::Error("Missing ElementTable or NodeTable name: use ValueFromTable[...]{name}");
237     return;
238   }
239 
240   std::map<int, std::vector<double> > &table(GetDPNumbersMap[Fct->String]);
241   std::vector<double> &val(table[Current.NumEntity]);
242   if(val.size() == 1){
243     V->Val[0] = val[0];
244     V->Type = SCALAR ;
245     return;
246   }
247 
248   //if(GetDPNumbersMap.size()){
249   //  Message::Warning("No element or node table found with name %s, or "
250   //                   "no entity index %d in the table", Fct->String,
251   //                   Current.NumEntity);
252   //}
253 
254   for(int i = 0; i < Fct->NbrArguments; i++){
255     if(i != 0){
256       Message::Warning("Ignoring %d-th argument in ValueFromTable");
257       continue;
258     }
259     if((A+i)->Type != SCALAR){
260       Message::Error("Non-scalar default argument in ValueFromTable");
261       return;
262     }
263     Cal_CopyValue(A + i, V);
264     return;
265   }
266 
267   Message::Warning("Missing table data or default value in ValueFromTable");
268   V->Val[0] = 0. ;
269   V->Type = SCALAR ;
270 
271 }
272 
F_VirtualWork(F_ARG)273 void F_VirtualWork (F_ARG)
274 {
275   MATRIX3x3 Jac ;
276   double DetJac ;
277   double DetJac_dx[3], squF[6] ;
278   double s[3] = {0.,0.,0.};
279 
280   if(!Current.Element)
281     Message::Error("Uninitialized Element in 'F_VirtualWork'");
282 
283   Current.flagAssDiag = 1; /*+++prov*/
284 
285   int numNode = Current.NumEntity;
286 
287   int i = 0 ;
288   while (i < Current.Element->GeoElement->NbrNodes &&
289          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
290 
291   if (i < Current.Element->GeoElement->NbrNodes ) {
292     for(int j = 0; j < 3; j++) {
293       squF[j]   = A->Val[j] * A->Val[j] ;
294       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
295     }
296 
297     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
298     DetJac  = Current.Element->DetJac ;
299     Jac     = Current.Element->Jac ;
300 
301     DetJac_dx[0] =
302       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
303       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
304       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
305 
306     DetJac_dx[1] =
307       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
308       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
309       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
310 
311     DetJac_dx[2] =
312       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
313       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
314       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
315 
316     if(DetJac != 0){
317       s[0] = ( DetJac_dx[0] * ( - squF[0] + squF[1] + squF[2] )
318                -  2 * DetJac_dx[1] * squF[3]
319                -  2 * DetJac_dx[2] * squF[5])/DetJac ;
320 
321       s[1] = ( DetJac_dx[1] * ( squF[0] - squF[1] + squF[2] )
322                -  2 * DetJac_dx[0] * squF[3]
323                -  2 * DetJac_dx[2] * squF[4])/DetJac ;
324 
325       s[2] =  ( DetJac_dx[2] * ( squF[0] + squF[1] - squF[2] )
326                 -  2 * DetJac_dx[0] * squF[5]
327                 -  2 * DetJac_dx[1] * squF[4])/DetJac ;
328     }
329     else {
330       Message::Warning("Zero determinant in 'F_VirtualWork'") ;
331     }
332   }
333 
334   V->Type = VECTOR ;
335   V->Val[0] = s[0] ;
336   V->Val[1] = s[1] ;
337   V->Val[2] = s[2] ;
338 }
339 
F_NodeForceDensity(F_ARG)340 void F_NodeForceDensity (F_ARG)
341 {
342   MATRIX3x3 Jac ;
343   double DetJac ;
344   double Grad_n[3] ;
345   double s11 = 0., s12 = 0., s13 = 0. ;
346   double s21 = 0., s22 = 0., s23 = 0. ;
347   double s31 = 0., s32 = 0., s33 = 0. ;
348   double s[3] = {0.,0.,0.};
349 
350   if(!Current.Element)
351     Message::Error("Uninitialized Element in 'F_NodeForceDensity'");
352 
353   Current.flagAssDiag = 1; /*+++prov*/
354 
355   int numNode = Current.NumEntity;
356 
357   int i = 0 ;
358   while (i < Current.Element->GeoElement->NbrNodes &&
359          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
360 
361   if (i < Current.Element->GeoElement->NbrNodes ) {
362     if(A->Type == TENSOR_SYM){
363       s11 = A->Val[0] ;
364       s12 = A->Val[1] ;
365       s13 = A->Val[2] ;
366       s21 = s12;
367       s22 = A->Val[3] ;
368       s23 = A->Val[4] ;
369       s31 = s13;
370       s32 = s23;
371       s33 = A->Val[5] ;
372     }
373     else if(A->Type == TENSOR){
374       s11 = A->Val[0] ;
375       s12 = A->Val[1] ;
376       s13 = A->Val[2] ;
377       s21 = A->Val[3] ;
378       s22 = A->Val[4] ;
379       s23 = A->Val[5] ;
380       s31 = A->Val[6] ;
381       s32 = A->Val[7] ;
382       s33 = A->Val[8] ;
383     }
384     else{
385       Message::Error("Unknown tensor type in 'F_NodeForceDensity'") ;
386     }
387 
388     DetJac  = Current.Element->DetJac ;
389     Jac     = Current.Element->Jac ;
390 
391     Grad_n[0] =
392       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
393       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
394       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
395 
396     Grad_n[1] =
397       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
398       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
399       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
400 
401     Grad_n[2] =
402       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
403       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
404       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
405 
406     if(DetJac != 0){
407       s[0] = ( Grad_n[0] * s11 + Grad_n[1] * s12 + Grad_n[2] * s13 ) / DetJac ;
408       s[1] = ( Grad_n[0] * s21 + Grad_n[1] * s22 + Grad_n[2] * s23 ) / DetJac ;
409       s[2] = ( Grad_n[0] * s31 + Grad_n[1] * s32 + Grad_n[2] * s33 ) / DetJac ;
410     }
411     else {
412       Message::Warning("Zero determinant in 'F_NodeForceDensity'") ;
413     }
414   }
415 
416   V->Type = VECTOR ;
417   V->Val[0] = s[0] ;
418   V->Val[1] = s[1] ;
419   V->Val[2] = s[2] ;
420 }
421 
422 // Blex added 25/04/14 update 06/06/14
F_Felec(F_ARG)423 void F_Felec (F_ARG)
424 {
425   MATRIX3x3 Jac ;
426   double DetJac ;
427   double DetJac_dx[3], squF[6] ;
428   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
429   double s[3] = {0.,0.,0.};
430 
431   if(!Current.Element)
432     Message::Error("Uninitialized Element in 'F_Felec'");
433 
434   Current.flagAssDiag = 1; /*+++prov*/
435 
436   int numNode = Current.NumEntity;
437 
438   int i = 0 ;
439   while (i < Current.Element->GeoElement->NbrNodes &&
440          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
441 
442   if (i < Current.Element->GeoElement->NbrNodes ) {
443 
444     for(int j = 0; j < 3; j++) {
445       squF[j]   = A->Val[j] * A->Val[j] ;
446       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
447     }
448 
449     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
450     DetJac  = Current.Element->DetJac ;
451     Jac     = Current.Element->Jac ;
452 
453     DetJac_dx[0] =
454       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
455       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
456       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
457 
458     DetJac_dx[1] =
459       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
460       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
461       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
462 
463     DetJac_dx[2] =
464       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
465       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
466       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
467 
468     if(DetJac != 0){
469       dsdx[0] = DetJac_dx[0]/DetJac ;
470       dsdx[1] = DetJac_dx[1]/DetJac ;
471       dsdx[2] = DetJac_dx[2]/DetJac ;
472 
473       s[0] = - 0.5 * ( dsdx[0] * ( squF[0] - squF[1] - squF[2] )
474                +  2 * dsdx[1] * squF[3]
475                +  2 * dsdx[2] * squF[5]) ;
476 
477       s[1] = - 0.5 * ( dsdx[1] * ( - squF[0] + squF[1] - squF[2] )
478                +  2 * dsdx[0] * squF[3]
479                +  2 * dsdx[2] * squF[4]) ;
480 
481       s[2] = - 0.5 * ( dsdx[2] * ( - squF[0] - squF[1] + squF[2] )
482                +  2 * dsdx[0] * squF[5]
483                +  2 * dsdx[1] * squF[4]) ;
484     }
485     else {
486       Message::Warning("Zero determinant in 'F_Felec'") ;
487     }
488   }
489 
490   V->Type = VECTOR ;
491   V->Val[0] = s[0] ;
492   V->Val[1] = s[1] ;
493   V->Val[2] = s[2] ;
494 }
495 
496 
F_dFxdux(F_ARG)497 void F_dFxdux (F_ARG)
498 {
499   MATRIX3x3 Jac ;
500   double DetJac ;
501   double DetJac_dx[3], squF[6] ;
502   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
503   double s[3] = {0.,0.,0.};
504 
505   if(!Current.Element)
506     Message::Error("Uninitialized Element in 'F_dFxdux'");
507 
508 
509   int numNode = Current.NumEntity;
510 
511   int i = 0 ;
512   while (i < Current.Element->GeoElement->NbrNodes &&
513          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
514 
515   if (i < Current.Element->GeoElement->NbrNodes ) {
516 
517     for(int j = 0; j < 3; j++) {
518       squF[j]   = A->Val[j] * A->Val[j] ;
519       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
520     }
521 
522     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
523     DetJac  = Current.Element->DetJac ;
524     Jac     = Current.Element->Jac ;
525 
526     DetJac_dx[0] =
527       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
528       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
529       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
530 
531     DetJac_dx[1] =
532       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
533       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
534       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
535 
536     DetJac_dx[2] =
537       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
538       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
539       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
540 
541     if(DetJac != 0){
542       dsdx[0] = DetJac_dx[0]/DetJac ;
543       dsdx[1] = DetJac_dx[1]/DetJac ;
544       dsdx[2] = DetJac_dx[2]/DetJac ;
545 
546       s[0] = - squF[0] * dsdx[0] ;
547 
548       s[1] = - squF[0] * dsdx[1] ;
549 
550       s[2] = - squF[0] * dsdx[2] ;
551     }
552     else {
553       Message::Warning("Zero determinant in 'F_dFxdux'") ;
554     }
555   }
556 
557   V->Type = VECTOR ;
558   V->Val[0] = s[0] ;
559   V->Val[1] = s[1] ;
560   V->Val[2] = s[2] ;
561 }
562 
F_dFydux(F_ARG)563 void F_dFydux (F_ARG)
564 {
565   MATRIX3x3 Jac ;
566   double DetJac ;
567   double DetJac_dx[3], squF[6] ;
568   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
569   double s[3] = {0.,0.,0.};
570 
571   if(!Current.Element)
572     Message::Error("Uninitialized Element in 'F_dFydux'");
573 
574 
575   int numNode = Current.NumEntity;
576 
577   int i = 0 ;
578   while (i < Current.Element->GeoElement->NbrNodes &&
579          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
580 
581   if (i < Current.Element->GeoElement->NbrNodes ) {
582 
583     for(int j = 0; j < 3; j++) {
584       squF[j]   = A->Val[j] * A->Val[j] ;
585       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
586     }
587 
588     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
589     DetJac  = Current.Element->DetJac ;
590     Jac     = Current.Element->Jac ;
591 
592     DetJac_dx[0] =
593       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
594       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
595       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
596 
597     DetJac_dx[1] =
598       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
599       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
600       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
601 
602     DetJac_dx[2] =
603       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
604       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
605       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
606 
607     if(DetJac != 0){
608       dsdx[0] = DetJac_dx[0]/DetJac ;
609       dsdx[1] = DetJac_dx[1]/DetJac ;
610       dsdx[2] = DetJac_dx[2]/DetJac ;
611 
612       s[0] = - 0.5 * (2 * squF[3] * dsdx[0] + (- squF[0] - squF[1] + squF[2]) * dsdx[1] - 2 * squF[4] * dsdx[2])  ;
613 
614       s[1] = - 0.5 * ((squF[0] + squF[1] - squF[2]) * dsdx[0] + 2 * squF[3] * dsdx[1] + 2 * squF[5] * dsdx[2]) ;
615 
616       s[2] = - 0.5 * (2 * squF[4] * dsdx[0] - 2 * squF[5] * dsdx[1] + 2 * squF[3] * dsdx[2]) ;
617     }
618     else {
619       Message::Warning("Zero determinant in 'F_dFydux'") ;
620     }
621   }
622 
623   V->Type = VECTOR ;
624   V->Val[0] = s[0] ;
625   V->Val[1] = s[1] ;
626   V->Val[2] = s[2] ;
627 }
628 
629 
F_dFzdux(F_ARG)630 void F_dFzdux (F_ARG)
631 {
632   MATRIX3x3 Jac ;
633   double DetJac ;
634   double DetJac_dx[3], squF[6] ;
635   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
636   double s[3] = {0.,0.,0.};
637 
638   if(!Current.Element)
639     Message::Error("Uninitialized Element in 'F_dFzdux'");
640 
641 
642   int numNode = Current.NumEntity;
643 
644   int i = 0 ;
645   while (i < Current.Element->GeoElement->NbrNodes &&
646          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
647 
648   if (i < Current.Element->GeoElement->NbrNodes ) {
649 
650     for(int j = 0; j < 3; j++) {
651       squF[j]   = A->Val[j] * A->Val[j] ;
652       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
653     }
654 
655     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
656     DetJac  = Current.Element->DetJac ;
657     Jac     = Current.Element->Jac ;
658 
659     DetJac_dx[0] =
660       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
661       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
662       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
663 
664     DetJac_dx[1] =
665       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
666       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
667       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
668 
669     DetJac_dx[2] =
670       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
671       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
672       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
673 
674     if(DetJac != 0){
675       dsdx[0] = DetJac_dx[0]/DetJac ;
676       dsdx[1] = DetJac_dx[1]/DetJac ;
677       dsdx[2] = DetJac_dx[2]/DetJac ;
678 
679       s[0] = - 0.5 * (2 * squF[5] * dsdx[0] - 2 * squF[4] * dsdx[1] + (-squF[0] + squF[1] - squF[2]) * dsdx[2]) ;
680 
681       s[1] = - 0.5 * (2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] - 2 * squF[3] * dsdx[2]) ;
682 
683       s[2] = - 0.5 * ((squF[0] - squF[1] + squF[2]) * dsdx[0] + 2 * squF[3] * dsdx[1] + 2 * squF[5] * dsdx[2]) ;
684     }
685     else {
686       Message::Warning("Zero determinant in 'F_dFzdux'") ;
687     }
688   }
689 
690   V->Type = VECTOR ;
691   V->Val[0] = s[0] ;
692   V->Val[1] = s[1] ;
693   V->Val[2] = s[2] ;
694 }
695 
696 
F_dFxduy(F_ARG)697 void F_dFxduy (F_ARG)
698 {
699   MATRIX3x3 Jac ;
700   double DetJac ;
701   double DetJac_dx[3], squF[6] ;
702   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
703   double s[3] = {0.,0.,0.};
704 
705   if(!Current.Element)
706     Message::Error("Uninitialized Element in 'F_dFxduy'");
707 
708 
709   int numNode = Current.NumEntity;
710 
711   int i = 0 ;
712   while (i < Current.Element->GeoElement->NbrNodes &&
713          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
714 
715   if (i < Current.Element->GeoElement->NbrNodes ) {
716 
717     for(int j = 0; j < 3; j++) {
718       squF[j]   = A->Val[j] * A->Val[j] ;
719       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
720     }
721 
722     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
723     DetJac  = Current.Element->DetJac ;
724     Jac     = Current.Element->Jac ;
725 
726     DetJac_dx[0] =
727       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
728       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
729       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
730 
731     DetJac_dx[1] =
732       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
733       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
734       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
735 
736     DetJac_dx[2] =
737       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
738       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
739       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
740 
741     if(DetJac != 0){
742       dsdx[0] = DetJac_dx[0]/DetJac ;
743       dsdx[1] = DetJac_dx[1]/DetJac ;
744       dsdx[2] = DetJac_dx[2]/DetJac ;
745 
746       s[0] = - 0.5 * (2 * squF[3] * dsdx[0] + (squF[0] + squF[1] - squF[2]) * dsdx[1] + 2 * squF[4] * dsdx[2]) ;
747 
748       s[1] = - 0.5 * ((-squF[0] - squF[1] + squF[2]) * dsdx[0] + 2 * squF[3] * dsdx[1] - 2 * squF[5] * dsdx[2]) ;
749 
750       s[2] = - 0.5 * (- 2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] + 2 * squF[3] * dsdx[2]) ;
751     }
752     else {
753       Message::Warning("Zero determinant in 'F_dFxduy'") ;
754     }
755   }
756 
757   V->Type = VECTOR ;
758   V->Val[0] = s[0] ;
759   V->Val[1] = s[1] ;
760   V->Val[2] = s[2] ;
761 }
762 
763 
F_dFyduy(F_ARG)764 void F_dFyduy (F_ARG)
765 {
766   MATRIX3x3 Jac ;
767   double DetJac ;
768   double DetJac_dx[3], squF[6] ;
769   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
770   double s[3] = {0.,0.,0.};
771 
772   if(!Current.Element)
773     Message::Error("Uninitialized Element in 'F_dFyduy'");
774 
775 
776   int numNode = Current.NumEntity;
777 
778   int i = 0 ;
779   while (i < Current.Element->GeoElement->NbrNodes &&
780          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
781 
782   if (i < Current.Element->GeoElement->NbrNodes ) {
783 
784     for(int j = 0; j < 3; j++) {
785       squF[j]   = A->Val[j] * A->Val[j] ;
786       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
787     }
788 
789     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
790     DetJac  = Current.Element->DetJac ;
791     Jac     = Current.Element->Jac ;
792 
793     DetJac_dx[0] =
794       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
795       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
796       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
797 
798     DetJac_dx[1] =
799       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
800       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
801       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
802 
803     DetJac_dx[2] =
804       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
805       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
806       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
807 
808     if(DetJac != 0){
809       dsdx[0] = DetJac_dx[0]/DetJac ;
810       dsdx[1] = DetJac_dx[1]/DetJac ;
811       dsdx[2] = DetJac_dx[2]/DetJac ;
812 
813       s[0] = - squF[1] * dsdx[0] ;
814 
815       s[1] = - squF[1] * dsdx[1] ;
816 
817       s[2] = - squF[1] * dsdx[2] ;
818     }
819     else {
820       Message::Warning("Zero determinant in 'F_dFyduy'") ;
821     }
822   }
823 
824   V->Type = VECTOR ;
825   V->Val[0] = s[0] ;
826   V->Val[1] = s[1] ;
827   V->Val[2] = s[2] ;
828 }
829 
830 
F_dFzduy(F_ARG)831 void F_dFzduy (F_ARG)
832 {
833   MATRIX3x3 Jac ;
834   double DetJac ;
835   double DetJac_dx[3], squF[6] ;
836   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
837   double s[3] = {0.,0.,0.};
838 
839   if(!Current.Element)
840     Message::Error("Uninitialized Element in 'F_dFzduy'");
841 
842 
843   int numNode = Current.NumEntity;
844 
845   int i = 0 ;
846   while (i < Current.Element->GeoElement->NbrNodes &&
847          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
848 
849   if (i < Current.Element->GeoElement->NbrNodes ) {
850 
851     for(int j = 0; j < 3; j++) {
852       squF[j]   = A->Val[j] * A->Val[j] ;
853       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
854     }
855 
856     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
857     DetJac  = Current.Element->DetJac ;
858     Jac     = Current.Element->Jac ;
859 
860     DetJac_dx[0] =
861       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
862       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
863       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
864 
865     DetJac_dx[1] =
866       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
867       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
868       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
869 
870     DetJac_dx[2] =
871       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
872       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
873       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
874 
875     if(DetJac != 0){
876       dsdx[0] = DetJac_dx[0]/DetJac ;
877       dsdx[1] = DetJac_dx[1]/DetJac ;
878       dsdx[2] = DetJac_dx[2]/DetJac ;
879 
880       s[0] = - 0.5 * (2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] - 2 * squF[3] * dsdx[2])  ;
881 
882       s[1] = - 0.5 * (-2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] + (squF[0] - squF[1] - squF[2]) * dsdx[2]) ;
883 
884       s[2] = - 0.5 * (2 * squF[3] * dsdx[0] + (-squF[0] + squF[1] + squF[2]) * dsdx[1] + 2 * squF[4] * dsdx[2]) ;
885     }
886     else {
887       Message::Warning("Zero determinant in 'F_dFzduy'") ;
888     }
889   }
890 
891   V->Type = VECTOR ;
892   V->Val[0] = s[0] ;
893   V->Val[1] = s[1] ;
894   V->Val[2] = s[2] ;
895 }
896 
897 
F_dFxduz(F_ARG)898 void F_dFxduz (F_ARG)
899 {
900   MATRIX3x3 Jac ;
901   double DetJac ;
902   double DetJac_dx[3], squF[6] ;
903   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
904   double s[3] = {0.,0.,0.};
905 
906   if(!Current.Element)
907     Message::Error("Uninitialized Element in 'F_dFxduz'");
908 
909 
910   int numNode = Current.NumEntity;
911 
912   int i = 0 ;
913   while (i < Current.Element->GeoElement->NbrNodes &&
914          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
915 
916   if (i < Current.Element->GeoElement->NbrNodes ) {
917 
918     for(int j = 0; j < 3; j++) {
919       squF[j]   = A->Val[j] * A->Val[j] ;
920       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
921     }
922 
923     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
924     DetJac  = Current.Element->DetJac ;
925     Jac     = Current.Element->Jac ;
926 
927     DetJac_dx[0] =
928       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
929       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
930       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
931 
932     DetJac_dx[1] =
933       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
934       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
935       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
936 
937     DetJac_dx[2] =
938       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
939       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
940       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
941 
942     if(DetJac != 0){
943       dsdx[0] = DetJac_dx[0]/DetJac ;
944       dsdx[1] = DetJac_dx[1]/DetJac ;
945       dsdx[2] = DetJac_dx[2]/DetJac ;
946 
947       s[0] = - 0.5 * (2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] + (squF[0] - squF[1] + squF[2]) * dsdx[2]) ;
948 
949       s[1] = - 0.5 * (-2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] + 2 * squF[3] * dsdx[2]) ;
950 
951       s[2] = - 0.5 * ((-squF[0] + squF[1] - squF[2]) * dsdx[0] - 2 * squF[3] * dsdx[1] + 2 * squF[5] * dsdx[2]) ;
952     }
953     else {
954       Message::Warning("Zero determinant in 'F_dFxduz'") ;
955     }
956   }
957 
958   V->Type = VECTOR ;
959   V->Val[0] = s[0] ;
960   V->Val[1] = s[1] ;
961   V->Val[2] = s[2] ;
962 }
963 
964 
F_dFyduz(F_ARG)965 void F_dFyduz (F_ARG)
966 {
967   MATRIX3x3 Jac ;
968   double DetJac ;
969   double DetJac_dx[3], squF[6] ;
970   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
971   double s[3] = {0.,0.,0.};
972 
973   if(!Current.Element)
974     Message::Error("Uninitialized Element in 'F_dFyduz'");
975 
976 
977   int numNode = Current.NumEntity;
978 
979   int i = 0 ;
980   while (i < Current.Element->GeoElement->NbrNodes &&
981          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
982 
983   if (i < Current.Element->GeoElement->NbrNodes ) {
984 
985     for(int j = 0; j < 3; j++) {
986       squF[j]   = A->Val[j] * A->Val[j] ;
987       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
988     }
989 
990     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
991     DetJac  = Current.Element->DetJac ;
992     Jac     = Current.Element->Jac ;
993 
994     DetJac_dx[0] =
995       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
996       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
997       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
998 
999     DetJac_dx[1] =
1000       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1001       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1002       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1003 
1004     DetJac_dx[2] =
1005       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1006       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1007       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1008 
1009     if(DetJac != 0){
1010       dsdx[0] = DetJac_dx[0]/DetJac ;
1011       dsdx[1] = DetJac_dx[1]/DetJac ;
1012       dsdx[2] = DetJac_dx[2]/DetJac ;
1013 
1014       s[0] = - 0.5 * (2 * squF[4] * dsdx[0] - 2 * squF[5] * dsdx[1] + 2 * squF[3] * dsdx[2]) ;
1015 
1016       s[1] = - 0.5 * (2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] + (-squF[0] + squF[1] + squF[2]) * dsdx[2]) ;
1017 
1018       s[2] = - 0.5 * (-2 * squF[3] * dsdx[0] + (squF[0] - squF[1] - squF[2]) * dsdx[1] + 2 * squF[4] * dsdx[2]) ;
1019     }
1020     else {
1021       Message::Warning("Zero determinant in 'F_dFyduz'") ;
1022     }
1023   }
1024 
1025   V->Type = VECTOR ;
1026   V->Val[0] = s[0] ;
1027   V->Val[1] = s[1] ;
1028   V->Val[2] = s[2] ;
1029 }
1030 
1031 
F_dFzduz(F_ARG)1032 void F_dFzduz (F_ARG)
1033 {
1034   MATRIX3x3 Jac ;
1035   double DetJac ;
1036   double DetJac_dx[3], squF[6] ;
1037   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1038   double s[3] = {0.,0.,0.};
1039 
1040   if(!Current.Element)
1041     Message::Error("Uninitialized Element in 'F_dFzduz'");
1042 
1043 
1044   int numNode = Current.NumEntity;
1045 
1046   int i = 0 ;
1047   while (i < Current.Element->GeoElement->NbrNodes &&
1048          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1049 
1050   if (i < Current.Element->GeoElement->NbrNodes ) {
1051 
1052     for(int j = 0; j < 3; j++) {
1053       squF[j]   = A->Val[j] * A->Val[j] ;
1054       squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ;
1055     }
1056 
1057     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1058     DetJac  = Current.Element->DetJac ;
1059     Jac     = Current.Element->Jac ;
1060 
1061     DetJac_dx[0] =
1062       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1063       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1064       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1065 
1066     DetJac_dx[1] =
1067       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1068       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1069       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1070 
1071     DetJac_dx[2] =
1072       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1073       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1074       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1075 
1076     if(DetJac != 0){
1077       dsdx[0] = DetJac_dx[0]/DetJac ;
1078       dsdx[1] = DetJac_dx[1]/DetJac ;
1079       dsdx[2] = DetJac_dx[2]/DetJac ;
1080 
1081       s[0] = - squF[2] * dsdx[0]  ;
1082 
1083       s[1] = - squF[2] * dsdx[1] ;
1084 
1085       s[2] = - squF[2] * dsdx[2] ;
1086     }
1087     else {
1088       Message::Warning("Zero determinant in 'F_dFzduz'") ;
1089     }
1090   }
1091 
1092   V->Type = VECTOR ;
1093   V->Val[0] = s[0] ;
1094   V->Val[1] = s[1] ;
1095   V->Val[2] = s[2] ;
1096 }
1097 
1098 
1099 
F_dFxdv(F_ARG)1100 void F_dFxdv (F_ARG)
1101 {
1102   MATRIX3x3 Jac ;
1103   double DetJac ;
1104   double DetJac_dx[3], squF[3] ;
1105   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1106   double s[3] = {0.,0.,0.};
1107 
1108   if(!Current.Element)
1109     Message::Error("Uninitialized Element in 'F_dFxdv'");
1110 
1111 
1112   int numNode = Current.NumEntity;
1113 
1114   int i = 0 ;
1115   while (i < Current.Element->GeoElement->NbrNodes &&
1116          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1117 
1118   if (i < Current.Element->GeoElement->NbrNodes ) {
1119 
1120     for(int j = 0; j < 3; j++) {
1121       squF[j]   = A->Val[j] ;
1122     }
1123 
1124     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1125     DetJac  = Current.Element->DetJac ;
1126     Jac     = Current.Element->Jac ;
1127 
1128     DetJac_dx[0] =
1129       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1130       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1131       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1132 
1133     DetJac_dx[1] =
1134       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1135       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1136       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1137 
1138     DetJac_dx[2] =
1139       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1140       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1141       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1142 
1143     if(DetJac != 0){
1144       dsdx[0] = DetJac_dx[0]/DetJac ;
1145       dsdx[1] = DetJac_dx[1]/DetJac ;
1146       dsdx[2] = DetJac_dx[2]/DetJac ;
1147 
1148       s[0] = - squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2]  ;
1149 
1150       s[1] = squF[1] * dsdx[0] - squF[0] * dsdx[1] ;
1151 
1152       s[2] = squF[2] * dsdx[0] - squF[0] * dsdx[2] ;
1153     }
1154     else {
1155       Message::Warning("Zero determinant in 'F_dFxdv'") ;
1156     }
1157   }
1158 
1159   V->Type = VECTOR ;
1160   V->Val[0] = s[0] ;
1161   V->Val[1] = s[1] ;
1162   V->Val[2] = s[2] ;
1163 }
1164 
F_dFydv(F_ARG)1165 void F_dFydv (F_ARG)
1166 {
1167   MATRIX3x3 Jac ;
1168   double DetJac ;
1169   double DetJac_dx[3], squF[3] ;
1170   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1171   double s[3] = {0.,0.,0.};
1172 
1173   if(!Current.Element)
1174     Message::Error("Uninitialized Element in 'F_dFydv'");
1175 
1176 
1177   int numNode = Current.NumEntity;
1178 
1179   int i = 0 ;
1180   while (i < Current.Element->GeoElement->NbrNodes &&
1181          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1182 
1183   if (i < Current.Element->GeoElement->NbrNodes ) {
1184 
1185     for(int j = 0; j < 3; j++) {
1186       squF[j]   = A->Val[j] ;
1187     }
1188 
1189     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1190     DetJac  = Current.Element->DetJac ;
1191     Jac     = Current.Element->Jac ;
1192 
1193     DetJac_dx[0] =
1194       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1195       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1196       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1197 
1198     DetJac_dx[1] =
1199       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1200       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1201       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1202 
1203     DetJac_dx[2] =
1204       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1205       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1206       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1207 
1208     if(DetJac != 0){
1209       dsdx[0] = DetJac_dx[0]/DetJac ;
1210       dsdx[1] = DetJac_dx[1]/DetJac ;
1211       dsdx[2] = DetJac_dx[2]/DetJac ;
1212 
1213       s[0] = - squF[1] * dsdx[0] + squF[0] * dsdx[1]  ;
1214 
1215       s[1] = - squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2] ;
1216 
1217       s[2] = squF[2] * dsdx[1] - squF[1] * dsdx[2] ;
1218     }
1219     else {
1220       Message::Warning("Zero determinant in 'F_dFydv'") ;
1221     }
1222   }
1223 
1224   V->Type = VECTOR ;
1225   V->Val[0] = s[0] ;
1226   V->Val[1] = s[1] ;
1227   V->Val[2] = s[2] ;
1228 }
1229 
F_dFzdv(F_ARG)1230 void F_dFzdv (F_ARG)
1231 {
1232   MATRIX3x3 Jac ;
1233   double DetJac ;
1234   double DetJac_dx[3], squF[3] ;
1235   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1236   double s[3] = {0.,0.,0.};
1237 
1238   if(!Current.Element)
1239     Message::Error("Uninitialized Element in 'F_dFzdv'");
1240 
1241 
1242   int numNode = Current.NumEntity;
1243 
1244   int i = 0 ;
1245   while (i < Current.Element->GeoElement->NbrNodes &&
1246          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1247 
1248   if (i < Current.Element->GeoElement->NbrNodes ) {
1249 
1250     for(int j = 0; j < 3; j++) {
1251       squF[j]   = A->Val[j] ;
1252     }
1253 
1254     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1255     DetJac  = Current.Element->DetJac ;
1256     Jac     = Current.Element->Jac ;
1257 
1258     DetJac_dx[0] =
1259       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1260       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1261       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1262 
1263     DetJac_dx[1] =
1264       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1265       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1266       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1267 
1268     DetJac_dx[2] =
1269       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1270       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1271       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1272 
1273     if(DetJac != 0){
1274       dsdx[0] = DetJac_dx[0]/DetJac ;
1275       dsdx[1] = DetJac_dx[1]/DetJac ;
1276       dsdx[2] = DetJac_dx[2]/DetJac ;
1277 
1278       s[0] = - squF[2] * dsdx[0] + squF[0] * dsdx[2]  ;
1279 
1280       s[1] = - squF[2] * dsdx[1] + squF[1] * dsdx[2] ;
1281 
1282       s[2] = - squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2] ;
1283     }
1284     else {
1285       Message::Warning("Zero determinant in 'F_dFzdv'") ;
1286     }
1287   }
1288 
1289   V->Type = VECTOR ;
1290   V->Val[0] = s[0] ;
1291   V->Val[1] = s[1] ;
1292   V->Val[2] = s[2] ;
1293 }
1294 
F_dWedxdv(F_ARG)1295 void F_dWedxdv (F_ARG)
1296 {
1297   MATRIX3x3 Jac ;
1298   double DetJac ;
1299   double DetJac_dx[3], squF[3] ;
1300   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1301   double s[3] = {0.,0.,0.};
1302 
1303   if(!Current.Element)
1304     Message::Error("Uninitialized Element in 'F_dWedxdv'");
1305 
1306 
1307   int numNode = Current.NumEntity;
1308 
1309   int i = 0 ;
1310   while (i < Current.Element->GeoElement->NbrNodes &&
1311          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1312 
1313   if (i < Current.Element->GeoElement->NbrNodes ) {
1314 
1315     for(int j = 0; j < 3; j++) {
1316       squF[j]   = A->Val[j] ;
1317     }
1318 
1319     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1320     DetJac  = Current.Element->DetJac ;
1321     Jac     = Current.Element->Jac ;
1322 
1323     DetJac_dx[0] =
1324       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1325       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1326       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1327 
1328     DetJac_dx[1] =
1329       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1330       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1331       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1332 
1333     DetJac_dx[2] =
1334       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1335       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1336       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1337 
1338     if(DetJac != 0){
1339       dsdx[0] = DetJac_dx[0]/DetJac ;
1340       dsdx[1] = DetJac_dx[1]/DetJac ;
1341       dsdx[2] = DetJac_dx[2]/DetJac ;
1342 
1343       s[0] = -squF[0] * dsdx[0] + squF[1] * dsdx[1] + squF[2] * dsdx[2] ;
1344 
1345       s[1] = -squF[1] * dsdx[0] - squF[0] * dsdx[1] ;
1346 
1347       s[2] = -squF[2] * dsdx[0] - squF[0] * dsdx[2] ;
1348     }
1349     else {
1350       Message::Warning("Zero determinant in 'F_dWedxdv'") ;
1351     }
1352   }
1353 
1354   V->Type = VECTOR ;
1355   V->Val[0] = s[0] ;
1356   V->Val[1] = s[1] ;
1357   V->Val[2] = s[2] ;
1358 }
1359 
F_dWedydv(F_ARG)1360 void F_dWedydv (F_ARG)
1361 {
1362   MATRIX3x3 Jac ;
1363   double DetJac ;
1364   double DetJac_dx[3], squF[3] ;
1365   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1366   double s[3] = {0.,0.,0.};
1367 
1368   if(!Current.Element)
1369     Message::Error("Uninitialized Element in 'F_dWedydv'");
1370 
1371 
1372   int numNode = Current.NumEntity;
1373 
1374   int i = 0 ;
1375   while (i < Current.Element->GeoElement->NbrNodes &&
1376          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1377 
1378   if (i < Current.Element->GeoElement->NbrNodes ) {
1379 
1380     for(int j = 0; j < 3; j++) {
1381       squF[j]   = A->Val[j] ;
1382     }
1383 
1384     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1385     DetJac  = Current.Element->DetJac ;
1386     Jac     = Current.Element->Jac ;
1387 
1388     DetJac_dx[0] =
1389       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1390       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1391       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1392 
1393     DetJac_dx[1] =
1394       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1395       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1396       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1397 
1398     DetJac_dx[2] =
1399       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1400       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1401       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1402 
1403     if(DetJac != 0){
1404       dsdx[0] = DetJac_dx[0]/DetJac ;
1405       dsdx[1] = DetJac_dx[1]/DetJac ;
1406       dsdx[2] = DetJac_dx[2]/DetJac ;
1407 
1408       s[0] = -squF[1] * dsdx[0] - squF[0] * dsdx[1] ;
1409 
1410       s[1] = squF[0] * dsdx[0] - squF[1] * dsdx[1] + squF[2] * dsdx[2] ;
1411 
1412       s[2] = -squF[2] * dsdx[1] - squF[1] * dsdx[2] ;
1413     }
1414     else {
1415       Message::Warning("Zero determinant in 'F_dWedydv'") ;
1416     }
1417   }
1418 
1419   V->Type = VECTOR ;
1420   V->Val[0] = s[0] ;
1421   V->Val[1] = s[1] ;
1422   V->Val[2] = s[2] ;
1423 }
1424 
F_dWedzdv(F_ARG)1425 void F_dWedzdv (F_ARG)
1426 {
1427   MATRIX3x3 Jac ;
1428   double DetJac ;
1429   double DetJac_dx[3], squF[3] ;
1430   double dsdx[3] ; //Derivative of the base functions with respect to x, y and z
1431   double s[3] = {0.,0.,0.};
1432 
1433   if(!Current.Element)
1434     Message::Error("Uninitialized Element in 'F_dWedzdv'");
1435 
1436 
1437   int numNode = Current.NumEntity;
1438 
1439   int i = 0 ;
1440   while (i < Current.Element->GeoElement->NbrNodes &&
1441          Current.Element->GeoElement->NumNodes[i] != numNode)  i++;
1442 
1443   if (i < Current.Element->GeoElement->NbrNodes ) {
1444 
1445     for(int j = 0; j < 3; j++) {
1446       squF[j]   = A->Val[j] ;
1447     }
1448 
1449     //Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
1450     DetJac  = Current.Element->DetJac ;
1451     Jac     = Current.Element->Jac ;
1452 
1453     DetJac_dx[0] =
1454       Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 )
1455       - Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 )
1456       + Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 );
1457 
1458     DetJac_dx[1] =
1459       - Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 )
1460       + Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 )
1461       - Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 );
1462 
1463     DetJac_dx[2] =
1464       Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 )
1465       - Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 )
1466       + Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 );
1467 
1468     if(DetJac != 0){
1469       dsdx[0] = DetJac_dx[0]/DetJac ;
1470       dsdx[1] = DetJac_dx[1]/DetJac ;
1471       dsdx[2] = DetJac_dx[2]/DetJac ;
1472 
1473       s[0] = -squF[2] * dsdx[0] - squF[0] * dsdx[2] ;
1474 
1475       s[1] = -squF[2] * dsdx[1] - squF[1] * dsdx[2] ;
1476 
1477       s[2] = squF[0] * dsdx[0] + squF[1] * dsdx[1] - squF[2] * dsdx[2] ;
1478     }
1479     else {
1480       Message::Warning("Zero determinant in 'F_dWedzdv'") ;
1481     }
1482   }
1483 
1484   V->Type = VECTOR ;
1485   V->Val[0] = s[0] ;
1486   V->Val[1] = s[1] ;
1487   V->Val[2] = s[2] ;
1488 }
1489 // End Blex added
1490 
F_AssDiag(F_ARG)1491 void F_AssDiag(F_ARG)
1492 {
1493   int  k ;
1494 
1495   if (Fct->NbrParameters == 1)
1496     Current.flagAssDiag = Fct->Para[0];
1497   else
1498     Current.flagAssDiag = 2; /*+++prov*/
1499 
1500   V->Val[0] = 1.;
1501 
1502   if (Current.NbrHar != 1){
1503     V->Val[MAX_DIM] = 0. ;
1504     for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
1505       V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
1506   }
1507 
1508   V->Type = SCALAR ;
1509 }
1510 
F_AtIndex(F_ARG)1511 void F_AtIndex(F_ARG)
1512 {
1513   int index = 0;
1514   double ret = 0.;
1515   if(A->Type != SCALAR){
1516     Message::Error("Non scalar argument for function 'AtIndex");
1517   }
1518   else{
1519     index = (int)A->Val[0];
1520     if (index < 0 || index >= Fct->NbrParameters){
1521       Message::Error("Wrong index in function 'AtIndex': %d not in [0,%d[",
1522                      index, Fct->NbrParameters);
1523     }
1524     else{
1525       ret = Fct->Para[index];
1526     }
1527   }
1528 
1529   V->Type = SCALAR;
1530   V->Val[0] = ret;
1531   if (Current.NbrHar != 1){
1532     V->Val[MAX_DIM] = 0. ;
1533     for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
1534       V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
1535   }
1536 }
1537