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 // Contributor(s):
7 // Johan Gyselinck
8 //
9
10 #include <sstream>
11 #include <string.h>
12 #include <math.h>
13 #include "ProData.h"
14 #include "ProDefine.h"
15 #include "Cal_Value.h"
16 #include "Message.h"
17
18 extern struct CurrentData Current ;
19
20 #define SQU(a) ((a)*(a))
21
22 /* ------------------------------------------------------------------------ */
23 /* O p e r a t o r s o n V a l u e s */
24 /* ------------------------------------------------------------------------ */
25 /* Warning: the pointers V1 and R or V2 and R can be identical. You must */
26 /* use temporary variables in your computations: you can only */
27 /* affect to R at the very last time (when you're sure you will */
28 /* not use V1 or V2 any more). */
29 /* ------------------------------------------------------------------------ */
30
31 static double a0;
32 static double a1 [NBR_MAX_HARMONIC * MAX_DIM] ;
33 static double a2 [NBR_MAX_HARMONIC * MAX_DIM] ;
34 static double b1 [NBR_MAX_HARMONIC * MAX_DIM] ;
35 static double b2 [NBR_MAX_HARMONIC * MAX_DIM] ;
36 static double b3 [NBR_MAX_HARMONIC * MAX_DIM] ;
37 static double c1 [NBR_MAX_HARMONIC * MAX_DIM] ;
38 static double c2 [NBR_MAX_HARMONIC * MAX_DIM] ;
39 static double c3 [NBR_MAX_HARMONIC * MAX_DIM] ;
40 static double tmp[27][NBR_MAX_HARMONIC * MAX_DIM] ;
41
42 static int TENSOR_SYM_MAP[] = {0,1,2,1,3,4,2,4,5};
43 static int TENSOR_DIAG_MAP[] = {0,-1,-1,-1,1,-1,-1,-1,2};
44
Cal_ComplexProduct(double V1[],double V2[],double P[])45 void Cal_ComplexProduct(double V1[], double V2[], double P[])
46 {
47 P[0] = V1[0] * V2[0] - V1[MAX_DIM] * V2[MAX_DIM] ;
48 P[MAX_DIM] = V1[0] * V2[MAX_DIM] + V1[MAX_DIM] * V2[0] ;
49 }
50
Cal_ComplexDivision(double V1[],double V2[],double P[])51 void Cal_ComplexDivision(double V1[], double V2[], double P[])
52 {
53 double Mod2 ;
54
55 Mod2 = SQU(V2[0]) + SQU(V2[MAX_DIM]) ;
56 if(!Mod2){ Message::Error("Division by zero in 'Cal_ComplexDivision'"); return; }
57 P[0] = ( V1[0] * V2[0] + V1[MAX_DIM] * V2[MAX_DIM]) / Mod2 ;
58 P[MAX_DIM] = (- V1[0] * V2[MAX_DIM] + V1[MAX_DIM] * V2[0]) / Mod2 ;
59 }
60
Cal_ComplexInvert(double V1[],double P[])61 void Cal_ComplexInvert(double V1[], double P[])
62 {
63 double Mod2 ;
64
65 Mod2 = SQU(V1[0]) + SQU(V1[MAX_DIM]) ;
66 if(!Mod2){ Message::Error("Division by zero in 'Cal_ComplexInvert'"); return; }
67 P[0] = V1[0] / Mod2 ;
68 P[MAX_DIM] = - V1[MAX_DIM] / Mod2 ;
69 }
70
71 /* ------------------------------------------------------------------------
72 R <- V1
73 ------------------------------------------------------------------------ */
74
Cal_CopyValue(struct Value * V1,struct Value * R)75 void Cal_CopyValue(struct Value * V1, struct Value * R)
76 {
77 int k ;
78
79 if (V1->Type == SCALAR) {
80 R->Type = SCALAR ;
81 for (k = 0 ; k < Current.NbrHar ; k++)
82 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] ;
83 }
84 else if (V1->Type == VECTOR || V1->Type == TENSOR_DIAG){
85 R->Type = V1->Type ;
86 for (k = 0 ; k < Current.NbrHar ; k++) {
87 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] ;
88 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] ;
89 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] ;
90 }
91 }
92 else if (V1->Type == TENSOR_SYM){
93 R->Type = TENSOR_SYM ;
94 for (k = 0 ; k < Current.NbrHar ; k++) {
95 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] ;
96 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] ;
97 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] ;
98 R->Val[MAX_DIM*k+3] = V1->Val[MAX_DIM*k+3] ;
99 R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4] ;
100 R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+5] ;
101 }
102 }
103 else if (V1->Type == TENSOR){
104 R->Type = TENSOR ;
105 for (k = 0 ; k < Current.NbrHar ; k++) {
106 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] ;
107 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] ;
108 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] ;
109 R->Val[MAX_DIM*k+3] = V1->Val[MAX_DIM*k+3] ;
110 R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4] ;
111 R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+5] ;
112 R->Val[MAX_DIM*k+6] = V1->Val[MAX_DIM*k+6] ;
113 R->Val[MAX_DIM*k+7] = V1->Val[MAX_DIM*k+7] ;
114 R->Val[MAX_DIM*k+8] = V1->Val[MAX_DIM*k+8] ;
115 }
116 }
117 }
118
119 /* ------------------------------------------------------------------------
120 R <- 0
121 ------------------------------------------------------------------------ */
122
123 //static double VALUE_ZERO [NBR_MAX_HARMONIC * MAX_DIM] =
124 // {0.,0.,0., 0.,0.,0.,
125 // 0.,0.,0., 0.,0.,0.,
126 // 0.,0.,0., 0.,0.,0.} ;
127
Cal_ZeroValue(struct Value * R)128 void Cal_ZeroValue(struct Value * R)
129 {
130 //memcpy(R->Val, VALUE_ZERO, Current.NbrHar*MAX_DIM*sizeof(double));
131 for(int i = 0; i < Current.NbrHar*MAX_DIM; i++) R->Val[i] = 0.;
132 }
133
134 /* ------------------------------------------------------------------------
135 R <- V1 + V2
136 ------------------------------------------------------------------------ */
137
138 #define ADD(i) R->Val[i] = V1->Val[i] + V2->Val[i]
139 #define CADD(i) R->Val[MAX_DIM*k+i] = V1->Val[MAX_DIM*k+i] + V2->Val[MAX_DIM*k+i]
140
Cal_AddValue(struct Value * V1,struct Value * V2,struct Value * R)141 void Cal_AddValue(struct Value * V1, struct Value * V2, struct Value * R)
142 {
143 int i, k;
144 int i1,i2;
145 struct Value A;
146
147 if (V1->Type == SCALAR && V2->Type == SCALAR) {
148 if (Current.NbrHar == 1) {
149 ADD(0);
150 }
151 else {
152 for (k = 0 ; k < Current.NbrHar ; k++) {
153 CADD(0);
154 }
155 }
156 R->Type = SCALAR ;
157 }
158
159 else if ((V1->Type == VECTOR && V2->Type == VECTOR) ||
160 (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_DIAG)) {
161 if (Current.NbrHar == 1) {
162 ADD(0); ADD(1); ADD(2);
163 }
164 else {
165 for (k = 0 ; k < Current.NbrHar ; k++) {
166 CADD(0); CADD(1); CADD(2);
167 }
168 }
169 R->Type = V1->Type;
170 }
171
172 else if (V1->Type == TENSOR_SYM && V2->Type == TENSOR_SYM) {
173 if (Current.NbrHar == 1) {
174 ADD(0); ADD(1); ADD(2); ADD(3); ADD(4); ADD(5);
175 }
176 else {
177 for (k = 0 ; k < Current.NbrHar ; k++) {
178 CADD(0); CADD(1); CADD(2); CADD(3); CADD(4); CADD(5);
179 }
180 }
181 R->Type = TENSOR_SYM;
182 }
183
184 else if (V1->Type == TENSOR && V2->Type == TENSOR) {
185 if (Current.NbrHar == 1) {
186 ADD(0); ADD(1); ADD(2); ADD(3); ADD(4); ADD(5); ADD(6); ADD(7); ADD(8);
187 }
188 else {
189 for (k = 0 ; k < Current.NbrHar ; k++) {
190 CADD(0); CADD(1); CADD(2); CADD(3); CADD(4); CADD(5); CADD(6); CADD(7); CADD(8);
191 }
192 }
193 R->Type = TENSOR;
194 }
195
196 else if ((V1->Type == TENSOR && V2->Type == TENSOR_SYM) ||
197 (V1->Type == TENSOR && V2->Type == TENSOR_DIAG)||
198 (V1->Type == TENSOR_SYM && V2->Type == TENSOR_DIAG)){
199 A.Type = V1->Type;
200 for (k = 0 ; k < Current.NbrHar ; k++) {
201 for(i=0 ; i<9 ; i++){
202 i1 = (V1->Type==TENSOR)?i:TENSOR_SYM_MAP[i];
203 i2 = (V2->Type==TENSOR_SYM)?TENSOR_SYM_MAP[i]:TENSOR_DIAG_MAP[i];
204 A.Val[MAX_DIM*k+i1] =
205 V1->Val[MAX_DIM*k+i1] + ((i2<0)?0.0:V2->Val[MAX_DIM*k+i2]);
206 }
207 }
208 Cal_CopyValue(&A,R);
209 }
210
211 else if ((V1->Type == TENSOR_SYM && V2->Type == TENSOR) ||
212 (V1->Type == TENSOR_DIAG && V2->Type == TENSOR) ||
213 (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_SYM)){
214 A.Type = V2->Type;
215 for (k = 0 ; k < Current.NbrHar ; k++) {
216 for(i=0 ; i<9 ; i++){
217 i1 = (V1->Type==TENSOR_SYM)?TENSOR_SYM_MAP[i]:TENSOR_DIAG_MAP[i];
218 i2 = (V2->Type==TENSOR)?i:TENSOR_SYM_MAP[i];
219 A.Val[MAX_DIM*k+i2] =
220 ((i1<0)?0.0:V1->Val[MAX_DIM*k+i1]) + V2->Val[MAX_DIM*k+i2];
221 }
222 }
223 Cal_CopyValue(&A,R);
224 }
225
226 else {
227 Message::Error("Addition of different quantities: %s + %s",
228 Get_StringForDefine(Field_Type, V1->Type),
229 Get_StringForDefine(Field_Type, V2->Type));
230 }
231 }
232
233 #undef ADD
234 #undef CADD
235
236 /* ------------------------------------------------------------------------
237 R <- V1 * d , where d is a double
238 ------------------------------------------------------------------------ */
239
Cal_MultValue(struct Value * V1,double d,struct Value * R)240 void Cal_MultValue(struct Value * V1, double d, struct Value * R)
241 {
242 int k;
243
244 R->Type = V1->Type ;
245
246 switch(V1->Type){
247 case SCALAR :
248 if (Current.NbrHar == 1) {
249 R->Val[0] = V1->Val[0] * d;
250 }
251 else{
252 for (k = 0 ; k < Current.NbrHar ; k++) {
253 R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k] * d;
254 }
255 }
256 break;
257 case VECTOR :
258 case TENSOR_DIAG :
259 if (Current.NbrHar == 1) {
260 R->Val[0] = V1->Val[0] * d;
261 R->Val[1] = V1->Val[1] * d;
262 R->Val[2] = V1->Val[2] * d;
263 }
264 else{
265 for (k = 0 ; k < Current.NbrHar ; k++) {
266 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] * d;
267 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] * d;
268 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] * d;
269 }
270 }
271 break;
272 case TENSOR_SYM :
273 if (Current.NbrHar == 1) {
274 R->Val[0] = V1->Val[0] * d;
275 R->Val[1] = V1->Val[1] * d;
276 R->Val[2] = V1->Val[2] * d;
277 R->Val[3] = V1->Val[3] * d;
278 R->Val[4] = V1->Val[4] * d;
279 R->Val[5] = V1->Val[5] * d;
280 }
281 else{
282 for (k = 0 ; k < Current.NbrHar ; k++) {
283 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] * d;
284 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] * d;
285 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] * d;
286 R->Val[MAX_DIM*k+3] = V1->Val[MAX_DIM*k+3] * d;
287 R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4] * d;
288 R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+5] * d;
289 }
290 }
291 break;
292 case TENSOR :
293 if (Current.NbrHar == 1) {
294 R->Val[0] = V1->Val[0] * d;
295 R->Val[1] = V1->Val[1] * d;
296 R->Val[2] = V1->Val[2] * d;
297 R->Val[3] = V1->Val[3] * d;
298 R->Val[4] = V1->Val[4] * d;
299 R->Val[5] = V1->Val[5] * d;
300 R->Val[6] = V1->Val[6] * d;
301 R->Val[7] = V1->Val[7] * d;
302 R->Val[8] = V1->Val[8] * d;
303 }
304 else{
305 for (k = 0 ; k < Current.NbrHar ; k++) {
306 R->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] * d;
307 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] * d;
308 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] * d;
309 R->Val[MAX_DIM*k+3] = V1->Val[MAX_DIM*k+3] * d;
310 R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4] * d;
311 R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+5] * d;
312 R->Val[MAX_DIM*k+6] = V1->Val[MAX_DIM*k+6] * d;
313 R->Val[MAX_DIM*k+7] = V1->Val[MAX_DIM*k+7] * d;
314 R->Val[MAX_DIM*k+8] = V1->Val[MAX_DIM*k+8] * d;
315 }
316 }
317 break;
318 default :
319 Message::Error("Wrong argument type for 'Cal_MultValue'");
320 break;
321 }
322 }
323
324 /* ------------------------------------------------------------------------
325 R <- V1 + V2 * d , where d is a double
326 ------------------------------------------------------------------------ */
327
Cal_AddMultValue(struct Value * V1,struct Value * V2,double d,struct Value * R)328 void Cal_AddMultValue(struct Value * V1, struct Value * V2, double d, struct Value * R)
329 {
330 int k;
331 struct Value A ;
332
333 A.Type = V2->Type ;
334
335 switch(V2->Type){
336 case SCALAR :
337 if (Current.NbrHar == 1) {
338 A.Val[0] = V2->Val[0] * d;
339 }
340 else{
341 for (k = 0 ; k < Current.NbrHar ; k++) {
342 A.Val[MAX_DIM*k] = V2->Val[MAX_DIM*k] * d;
343 }
344 }
345 break;
346 case VECTOR :
347 case TENSOR_DIAG :
348 if (Current.NbrHar == 1) {
349 A.Val[0] = V2->Val[0] * d;
350 A.Val[1] = V2->Val[1] * d;
351 A.Val[2] = V2->Val[2] * d;
352 }
353 else{
354 for (k = 0 ; k < Current.NbrHar ; k++) {
355 A.Val[MAX_DIM*k ] = V2->Val[MAX_DIM*k ] * d;
356 A.Val[MAX_DIM*k+1] = V2->Val[MAX_DIM*k+1] * d;
357 A.Val[MAX_DIM*k+2] = V2->Val[MAX_DIM*k+2] * d;
358 }
359 }
360 break;
361 case TENSOR_SYM :
362 if (Current.NbrHar == 1) {
363 A.Val[0] = V2->Val[0] * d;
364 A.Val[1] = V2->Val[1] * d;
365 A.Val[2] = V2->Val[2] * d;
366 A.Val[3] = V2->Val[3] * d;
367 A.Val[4] = V2->Val[4] * d;
368 A.Val[5] = V2->Val[5] * d;
369 }
370 else{
371 for (k = 0 ; k < Current.NbrHar ; k++) {
372 A.Val[MAX_DIM*k ] = V2->Val[MAX_DIM*k ] * d;
373 A.Val[MAX_DIM*k+1] = V2->Val[MAX_DIM*k+1] * d;
374 A.Val[MAX_DIM*k+2] = V2->Val[MAX_DIM*k+2] * d;
375 A.Val[MAX_DIM*k+3] = V2->Val[MAX_DIM*k+3] * d;
376 A.Val[MAX_DIM*k+4] = V2->Val[MAX_DIM*k+4] * d;
377 A.Val[MAX_DIM*k+5] = V2->Val[MAX_DIM*k+5] * d;
378 }
379 }
380 break;
381 case TENSOR :
382 if (Current.NbrHar == 1) {
383 A.Val[0] = V2->Val[0] * d;
384 A.Val[1] = V2->Val[1] * d;
385 A.Val[2] = V2->Val[2] * d;
386 A.Val[3] = V2->Val[3] * d;
387 A.Val[4] = V2->Val[4] * d;
388 A.Val[5] = V2->Val[5] * d;
389 A.Val[6] = V2->Val[6] * d;
390 A.Val[7] = V2->Val[7] * d;
391 A.Val[8] = V2->Val[8] * d;
392 }
393 else{
394 for (k = 0 ; k < Current.NbrHar ; k++) {
395 A.Val[MAX_DIM*k ] = V2->Val[MAX_DIM*k ] * d;
396 A.Val[MAX_DIM*k+1] = V2->Val[MAX_DIM*k+1] * d;
397 A.Val[MAX_DIM*k+2] = V2->Val[MAX_DIM*k+2] * d;
398 A.Val[MAX_DIM*k+3] = V2->Val[MAX_DIM*k+3] * d;
399 A.Val[MAX_DIM*k+4] = V2->Val[MAX_DIM*k+4] * d;
400 A.Val[MAX_DIM*k+5] = V2->Val[MAX_DIM*k+5] * d;
401 A.Val[MAX_DIM*k+6] = V2->Val[MAX_DIM*k+6] * d;
402 A.Val[MAX_DIM*k+7] = V2->Val[MAX_DIM*k+7] * d;
403 A.Val[MAX_DIM*k+8] = V2->Val[MAX_DIM*k+8] * d;
404 }
405 }
406 break;
407 default :
408 Message::Error("Wrong argument type for 'Cal_AddMultValue'");
409 return;
410 }
411 Cal_AddValue(V1,&A,R);
412 }
413
414 /* ------------------------------------------------------------------------
415 V1 <- V1 * d1 + V2 * d2 , where d1, d2 are doubles
416 ------------------------------------------------------------------------ */
417
Cal_AddMultValue2(struct Value * V1,double d1,struct Value * V2,double d2)418 void Cal_AddMultValue2(struct Value * V1, double d1,
419 struct Value * V2, double d2)
420 {
421 int k;
422
423 switch(V1->Type){
424 case SCALAR :
425 if (Current.NbrHar == 1) {
426 V1->Val[0] = V1->Val[0] * d1 + V2->Val[0] * d2;
427 }
428 else{
429 for (k = 0 ; k < Current.NbrHar ; k++) {
430 V1->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k] * d1 + V2->Val[MAX_DIM*k] * d2;
431 }
432 }
433 break;
434 case VECTOR :
435 case TENSOR_DIAG :
436 if (Current.NbrHar == 1) {
437 V1->Val[0] = V1->Val[0] * d1 + V2->Val[0] * d2;
438 V1->Val[1] = V1->Val[1] * d1 + V2->Val[1] * d2;
439 V1->Val[2] = V1->Val[2] * d1 + V2->Val[2] * d2;
440 }
441 else{
442 for (k = 0 ; k < Current.NbrHar ; k++) {
443 V1->Val[MAX_DIM*k ] = V1->Val[MAX_DIM*k ] * d1 + V2->Val[MAX_DIM*k ] * d2;
444 V1->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] * d1 + V2->Val[MAX_DIM*k+1] * d2;
445 V1->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] * d1 + V2->Val[MAX_DIM*k+2] * d2;
446 }
447 }
448 break;
449 default :
450 Message::Error("Wrong argument type for 'Cal_AddMultValue2'");
451 break;
452 }
453 }
454
455 /* ------------------------------------------------------------------------
456 R <- V1 - V2
457 ------------------------------------------------------------------------ */
458
459 #define SUB(i) R->Val[i] = V1->Val[i] - V2->Val[i]
460 #define CSUB(i) R->Val[MAX_DIM*k+i] = V1->Val[MAX_DIM*k+i] - V2->Val[MAX_DIM*k+i]
461
Cal_SubstractValue(struct Value * V1,struct Value * V2,struct Value * R)462 void Cal_SubstractValue(struct Value * V1, struct Value * V2, struct Value * R)
463 {
464 int i, k;
465 int i1, i2;
466 struct Value A;
467
468 if (V1->Type == SCALAR && V2->Type == SCALAR) {
469 if (Current.NbrHar == 1) {
470 SUB(0);
471 }
472 else {
473 for (k = 0 ; k < Current.NbrHar ; k++) {
474 CSUB(0);
475 }
476 }
477 R->Type = SCALAR ;
478 }
479
480 else if ((V1->Type == VECTOR && V2->Type == VECTOR) ||
481 (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_DIAG)) {
482 if (Current.NbrHar == 1) {
483 SUB(0); SUB(1); SUB(2);
484 }
485 else {
486 for (k = 0 ; k < Current.NbrHar ; k++) {
487 CSUB(0); CSUB(1); CSUB(2);
488 }
489 }
490 R->Type = V1->Type ;
491 }
492
493 else if (V1->Type == TENSOR_SYM && V2->Type == TENSOR_SYM) {
494 if (Current.NbrHar == 1) {
495 SUB(0); SUB(1); SUB(2); SUB(3); SUB(4); SUB(5);
496 }
497 else {
498 for (k = 0 ; k < Current.NbrHar ; k++) {
499 CSUB(0); CSUB(1); CSUB(2); CSUB(3); CSUB(4); CSUB(5);
500 }
501 }
502 R->Type = TENSOR_SYM;
503 }
504
505 else if (V1->Type == TENSOR && V2->Type == TENSOR) {
506 if (Current.NbrHar == 1) {
507 SUB(0); SUB(1); SUB(2); SUB(3); SUB(4); SUB(5); SUB(6); SUB(7); SUB(8);
508 }
509 else {
510 for (k = 0 ; k < Current.NbrHar ; k++) {
511 CSUB(0); CSUB(1); CSUB(2); CSUB(3); CSUB(4); CSUB(5); CSUB(6); CSUB(7); CSUB(8);
512 }
513 }
514 R->Type = TENSOR;
515 }
516
517 else if ((V1->Type == TENSOR && V2->Type == TENSOR_SYM) ||
518 (V1->Type == TENSOR && V2->Type == TENSOR_DIAG)||
519 (V1->Type == TENSOR_SYM && V2->Type == TENSOR_DIAG)){
520 A.Type = V1->Type;
521 for (k = 0 ; k < Current.NbrHar ; k++) {
522 for(i=0 ; i<9 ; i++){
523 i1 = (V1->Type==TENSOR)?i:TENSOR_SYM_MAP[i];
524 i2 = (V2->Type==TENSOR_SYM)?TENSOR_SYM_MAP[i]:TENSOR_DIAG_MAP[i];
525 A.Val[MAX_DIM*k+i1] =
526 V1->Val[MAX_DIM*k+i1] - ((i2<0)?0.0:V2->Val[MAX_DIM*k+i2]);
527 }
528 }
529 Cal_CopyValue(&A,R);
530 }
531
532 else if ((V1->Type == TENSOR_SYM && V2->Type == TENSOR) ||
533 (V1->Type == TENSOR_DIAG && V2->Type == TENSOR) ||
534 (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_SYM)){
535 A.Type = V2->Type;
536 for (k = 0 ; k < Current.NbrHar ; k++) {
537 for(i=0 ; i<9 ; i++){
538 i1 = (V1->Type==TENSOR_SYM)?TENSOR_SYM_MAP[i]:TENSOR_DIAG_MAP[i];
539 i2 = (V2->Type==TENSOR)?i:TENSOR_SYM_MAP[i];
540 A.Val[MAX_DIM*k+i2] =
541 ((i1<0)?0.0:V1->Val[MAX_DIM*k+i1]) - V2->Val[MAX_DIM*k+i2];
542 }
543 }
544 Cal_CopyValue(&A,R);
545 }
546
547 else {
548 Message::Error("Substraction of different quantities: %s - %s",
549 Get_StringForDefine(Field_Type, V1->Type),
550 Get_StringForDefine(Field_Type, V2->Type));
551 }
552 }
553
554 #undef SUB
555 #undef CSUB
556
557 /* ------------------------------------------------------------------------
558 R <- V1 * V2
559 ------------------------------------------------------------------------ */
560
561 #define CMULT(a,b,c) \
562 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k+a]), &(V2->Val[MAX_DIM*k+b]), tmp[c])
563
564 #define CPUT(a) \
565 R->Val[MAX_DIM* k +a] = tmp[a][0] ; \
566 R->Val[MAX_DIM*(k+1)+a] = tmp[a][MAX_DIM]
567
568 #define CPUT3(a,b,c,d) \
569 R->Val[MAX_DIM* k +d] = tmp[a][0] +tmp[b][0] +tmp[c][0] ; \
570 R->Val[MAX_DIM*(k+1)+d] = tmp[a][MAX_DIM]+tmp[b][MAX_DIM]+tmp[c][MAX_DIM]
571
Cal_ProductValue(struct Value * V1,struct Value * V2,struct Value * R)572 void Cal_ProductValue(struct Value * V1, struct Value * V2, struct Value * R)
573 {
574 int k;
575
576 if (V1->Type == SCALAR && V2->Type == SCALAR) {
577 if (Current.NbrHar == 1) {
578 R->Val[0] = V1->Val[0]*V2->Val[0];
579 }
580 else {
581 for (k = 0 ; k < Current.NbrHar ; k += 2) {
582 CMULT(0,0,0);
583 CPUT(0);
584 }
585 }
586 R->Type = SCALAR ;
587 }
588
589 else if (V1->Type == SCALAR && (V2->Type == VECTOR || V2->Type == TENSOR_DIAG)) {
590 if (Current.NbrHar == 1) {
591 a0 = V1->Val[0] ;
592 R->Val[0] = a0 * V2->Val[0] ;
593 R->Val[1] = a0 * V2->Val[1] ;
594 R->Val[2] = a0 * V2->Val[2] ;
595 }
596 else {
597 for (k = 0 ; k < Current.NbrHar ; k += 2) {
598 CMULT(0,0,0); CMULT(0,1,1); CMULT(0,2,2);
599 CPUT(0); CPUT(1); CPUT(2);
600 }
601 }
602 R->Type = V2->Type ;
603 }
604
605 else if (V1->Type == SCALAR && V2->Type == TENSOR_SYM) {
606 if (Current.NbrHar == 1) {
607 a0 = V1->Val[0] ;
608 R->Val[0] = a0 * V2->Val[0] ;
609 R->Val[1] = a0 * V2->Val[1] ;
610 R->Val[2] = a0 * V2->Val[2] ;
611 R->Val[3] = a0 * V2->Val[3] ;
612 R->Val[4] = a0 * V2->Val[4] ;
613 R->Val[5] = a0 * V2->Val[5] ;
614 }
615 else {
616 for (k = 0 ; k < Current.NbrHar ; k += 2) {
617 CMULT(0,0,0); CMULT(0,1,1); CMULT(0,2,2); CMULT(0,3,3); CMULT(0,4,4); CMULT(0,5,5);
618 CPUT(0); CPUT(1); CPUT(2); CPUT(3); CPUT(4); CPUT(5);
619 }
620 }
621 R->Type = TENSOR_SYM ;
622 }
623
624 else if (V1->Type == SCALAR && V2->Type == TENSOR) {
625 if (Current.NbrHar == 1) {
626 a0 = V1->Val[0] ;
627 R->Val[0] = a0 * V2->Val[0] ;
628 R->Val[1] = a0 * V2->Val[1] ;
629 R->Val[2] = a0 * V2->Val[2] ;
630 R->Val[3] = a0 * V2->Val[3] ;
631 R->Val[4] = a0 * V2->Val[4] ;
632 R->Val[5] = a0 * V2->Val[5] ;
633 R->Val[6] = a0 * V2->Val[6] ;
634 R->Val[7] = a0 * V2->Val[7] ;
635 R->Val[8] = a0 * V2->Val[8] ;
636 }
637 else {
638 for (k = 0 ; k < Current.NbrHar ; k += 2) {
639 CMULT(0,0,0); CMULT(0,1,1); CMULT(0,2,2); CMULT(0,3,3); CMULT(0,4,4); CMULT(0,5,5);
640 CMULT(0,6,6); CMULT(0,7,7); CMULT(0,8,8);
641 CPUT(0); CPUT(1); CPUT(2); CPUT(3); CPUT(4); CPUT(5);
642 CPUT(6); CPUT(7); CPUT(8);
643 }
644 }
645 R->Type = TENSOR ;
646 }
647
648 else if ((V1->Type == VECTOR || V1->Type == TENSOR_DIAG) && V2->Type == SCALAR) {
649 if (Current.NbrHar == 1) {
650 a0 = V2->Val[0] ;
651 R->Val[0] = V1->Val[0] * a0 ;
652 R->Val[1] = V1->Val[1] * a0 ;
653 R->Val[2] = V1->Val[2] * a0 ;
654 }
655 else {
656 for (k = 0 ; k < Current.NbrHar ; k += 2) {
657 CMULT(0,0,0); CMULT(1,0,1); CMULT(2,0,2);
658 CPUT(0); CPUT(1); CPUT(2);
659 }
660 }
661 R->Type = V1->Type ;
662 }
663
664 else if (V1->Type == TENSOR_SYM && V2->Type == SCALAR) {
665 if (Current.NbrHar == 1) {
666 a0 = V2->Val[0] ;
667 R->Val[0] = V1->Val[0] * a0 ;
668 R->Val[1] = V1->Val[1] * a0 ;
669 R->Val[2] = V1->Val[2] * a0 ;
670 R->Val[3] = V1->Val[3] * a0 ;
671 R->Val[4] = V1->Val[4] * a0 ;
672 R->Val[5] = V1->Val[5] * a0 ;
673 }
674 else {
675 for (k = 0 ; k < Current.NbrHar ; k += 2) {
676 CMULT(0,0,0); CMULT(1,0,1); CMULT(2,0,2); CMULT(3,0,3); CMULT(4,0,4); CMULT(5,0,5);
677 CPUT(0); CPUT(1); CPUT(2); CPUT(3); CPUT(4); CPUT(5);
678 }
679 }
680 R->Type = TENSOR_SYM ;
681 }
682
683 else if (V1->Type == TENSOR && V2->Type == SCALAR) {
684 if (Current.NbrHar == 1) {
685 a0 = V2->Val[0] ;
686 R->Val[0] = V1->Val[0] * a0 ;
687 R->Val[1] = V1->Val[1] * a0 ;
688 R->Val[2] = V1->Val[2] * a0 ;
689 R->Val[3] = V1->Val[3] * a0 ;
690 R->Val[4] = V1->Val[4] * a0 ;
691 R->Val[5] = V1->Val[5] * a0 ;
692 R->Val[6] = V1->Val[6] * a0 ;
693 R->Val[7] = V1->Val[7] * a0 ;
694 R->Val[8] = V1->Val[8] * a0 ;
695 }
696 else {
697 for (k = 0 ; k < Current.NbrHar ; k += 2) {
698 CMULT(0,0,0); CMULT(1,0,1); CMULT(2,0,2); CMULT(3,0,3); CMULT(4,0,4); CMULT(5,0,5);
699 CMULT(6,0,6); CMULT(7,0,7); CMULT(8,0,8);
700 CPUT(0); CPUT(1); CPUT(2); CPUT(3); CPUT(4); CPUT(5);
701 CPUT(6); CPUT(7); CPUT(8);
702 }
703 }
704 R->Type = TENSOR ;
705 }
706
707 /* Scalar Product. See 'Cal_CrossProductValue' for the Vector Product */
708
709 else if (V1->Type == VECTOR && V2->Type == VECTOR) {
710 if (Current.NbrHar == 1) {
711 R->Val[0] = V1->Val[0] * V2->Val[0] +
712 V1->Val[1] * V2->Val[1] +
713 V1->Val[2] * V2->Val[2] ;
714 }
715 else {
716 for (k = 0 ; k < Current.NbrHar ; k += 2) {
717 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
718 CPUT3(0,1,2,0);
719 }
720 }
721 R->Type = SCALAR ;
722 }
723
724 else if ( (V1->Type == TENSOR_DIAG && V2->Type == VECTOR) ||
725 (V2->Type == TENSOR_DIAG && V1->Type == VECTOR) ) { /* WARNING WARNING! */
726 if (Current.NbrHar == 1) {
727 R->Val[0] = V1->Val[0]*V2->Val[0];
728 R->Val[1] = V1->Val[1]*V2->Val[1];
729 R->Val[2] = V1->Val[2]*V2->Val[2];
730 }
731 else {
732 for (k = 0 ; k < Current.NbrHar ; k += 2) {
733 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
734 CPUT(0); CPUT(1); CPUT(2);
735 }
736 }
737 R->Type = VECTOR ;
738 }
739
740 else if (V1->Type == TENSOR_SYM && V2->Type == VECTOR) {
741 if (Current.NbrHar == 1) {
742 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[1] + V1->Val[2]*V2->Val[2];
743 a1[1] = V1->Val[1]*V2->Val[0] + V1->Val[3]*V2->Val[1] + V1->Val[4]*V2->Val[2];
744 a1[2] = V1->Val[2]*V2->Val[0] + V1->Val[4]*V2->Val[1] + V1->Val[5]*V2->Val[2];
745 R->Val[0] = a1[0];
746 R->Val[1] = a1[1];
747 R->Val[2] = a1[2];
748 }
749 else {
750 for (k = 0 ; k < Current.NbrHar ; k += 2) {
751 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
752 CMULT(1,0,3); CMULT(3,1,4); CMULT(4,2,5);
753 CMULT(2,0,6); CMULT(4,1,7); CMULT(5,2,8);
754 CPUT3(0,1,2,0);
755 CPUT3(3,4,5,1);
756 CPUT3(6,7,8,2);
757 }
758 }
759 R->Type = VECTOR ;
760 }
761
762 else if (V1->Type == TENSOR && V2->Type == VECTOR) {
763 if (Current.NbrHar == 1) {
764 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[1] + V1->Val[2]*V2->Val[2];
765 a1[1] = V1->Val[3]*V2->Val[0] + V1->Val[4]*V2->Val[1] + V1->Val[5]*V2->Val[2];
766 a1[2] = V1->Val[6]*V2->Val[0] + V1->Val[7]*V2->Val[1] + V1->Val[8]*V2->Val[2];
767 R->Val[0] = a1[0];
768 R->Val[1] = a1[1];
769 R->Val[2] = a1[2];
770 }
771 else {
772 for (k = 0 ; k < Current.NbrHar ; k += 2) {
773 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
774 CMULT(3,0,3); CMULT(4,1,4); CMULT(5,2,5);
775 CMULT(6,0,6); CMULT(7,1,7); CMULT(8,2,8);
776 CPUT3(0,1,2,0);
777 CPUT3(3,4,5,1);
778 CPUT3(6,7,8,2);
779 }
780 }
781 R->Type = VECTOR ;
782 }
783
784 else if (V1->Type == VECTOR && V2->Type == TENSOR) {
785 if (Current.NbrHar == 1) {
786 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[3] + V1->Val[2]*V2->Val[6];
787 a1[1] = V1->Val[0]*V2->Val[1] + V1->Val[1]*V2->Val[4] + V1->Val[2]*V2->Val[7];
788 a1[2] = V1->Val[0]*V2->Val[2] + V1->Val[1]*V2->Val[5] + V1->Val[2]*V2->Val[8];
789 R->Val[0] = a1[0];
790 R->Val[1] = a1[1];
791 R->Val[2] = a1[2];
792 }
793 else {
794 for (k = 0 ; k < Current.NbrHar ; k += 2) {
795 CMULT(0,0,0); CMULT(0,1,1); CMULT(0,2,2);
796 CMULT(1,3,3); CMULT(1,4,4); CMULT(1,5,5);
797 CMULT(2,6,6); CMULT(2,7,7); CMULT(2,8,8);
798 CPUT3(0,3,6,0);
799 CPUT3(1,4,7,1);
800 CPUT3(2,5,8,2);
801 }
802 }
803 R->Type = VECTOR ;
804 }
805
806 else if (V1->Type == TENSOR && V2->Type == TENSOR_DIAG) {
807 if (Current.NbrHar == 1) {
808 a1[0] = V1->Val[0]*V2->Val[0];
809 a1[1] = V1->Val[1]*V2->Val[1];
810 a1[2] = V1->Val[2]*V2->Val[2];
811 a1[3] = V1->Val[3]*V2->Val[0];
812 a1[4] = V1->Val[4]*V2->Val[1];
813 a1[5] = V1->Val[5]*V2->Val[2];
814 a1[6] = V1->Val[6]*V2->Val[0];
815 a1[7] = V1->Val[7]*V2->Val[1];
816 a1[8] = V1->Val[8]*V2->Val[2];
817 R->Val[0] = a1[0]; R->Val[1] = a1[1]; R->Val[2] = a1[2];
818 R->Val[3] = a1[3]; R->Val[4] = a1[4]; R->Val[5] = a1[5];
819 R->Val[6] = a1[6]; R->Val[7] = a1[7]; R->Val[8] = a1[8];
820
821 }
822 else {
823 for (k = 0 ; k < Current.NbrHar ; k += 2) {
824 CMULT(0,0,0);
825 CMULT(1,1,1);
826 CMULT(2,2,2);
827 CMULT(3,0,3);
828 CMULT(4,1,4);
829 CMULT(5,2,5);
830 CMULT(6,0,6);
831 CMULT(7,1,7);
832 CMULT(8,2,8);
833 CPUT(0); CPUT(1); CPUT(2);
834 CPUT(3); CPUT(4); CPUT(5);
835 CPUT(6); CPUT(7); CPUT(8);
836 }
837 }
838 R->Type = TENSOR;
839 }
840
841 else if (V1->Type == TENSOR_DIAG && V2->Type == TENSOR) {
842 if (Current.NbrHar == 1) {
843 a1[0] = V1->Val[0]*V2->Val[0];
844 a1[1] = V1->Val[0]*V2->Val[1];
845 a1[2] = V1->Val[0]*V2->Val[2];
846 a1[3] = V1->Val[1]*V2->Val[3];
847 a1[4] = V1->Val[1]*V2->Val[4];
848 a1[5] = V1->Val[1]*V2->Val[5];
849 a1[6] = V1->Val[2]*V2->Val[6];
850 a1[7] = V1->Val[2]*V2->Val[7];
851 a1[8] = V1->Val[2]*V2->Val[8];
852 R->Val[0] = a1[0]; R->Val[1] = a1[1]; R->Val[2] = a1[2];
853 R->Val[3] = a1[3]; R->Val[4] = a1[4]; R->Val[5] = a1[5];
854 R->Val[6] = a1[6]; R->Val[7] = a1[7]; R->Val[8] = a1[8];
855
856 }
857 else {
858 for (k = 0 ; k < Current.NbrHar ; k += 2) {
859 CMULT(0,0,0);
860 CMULT(0,1,1);
861 CMULT(0,2,2);
862 CMULT(1,3,3);
863 CMULT(1,4,4);
864 CMULT(1,5,5);
865 CMULT(2,6,6);
866 CMULT(2,7,7);
867 CMULT(2,8,8);
868 CPUT(0); CPUT(1); CPUT(2);
869 CPUT(3); CPUT(4); CPUT(5);
870 CPUT(6); CPUT(7); CPUT(8);
871 }
872 }
873 R->Type = TENSOR;
874 }
875
876 else if (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_DIAG) {
877 if (Current.NbrHar == 1) {
878 R->Val[0] = V1->Val[0]*V2->Val[0];
879 R->Val[1] = V1->Val[1]*V2->Val[1];
880 R->Val[2] = V1->Val[2]*V2->Val[2];
881 }
882 else {
883 for (k = 0 ; k < Current.NbrHar ; k += 2) {
884 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
885 CPUT(0); CPUT(1); CPUT(2);
886 }
887 }
888 R->Type = TENSOR_DIAG;
889 }
890
891 else if (V1->Type == TENSOR_SYM && V2->Type == TENSOR_SYM) {
892 if (Current.NbrHar == 1) {
893 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[1] + V1->Val[2]*V2->Val[2];
894 a1[1] = V1->Val[0]*V2->Val[1] + V1->Val[1]*V2->Val[3] + V1->Val[2]*V2->Val[4];
895 a1[2] = V1->Val[0]*V2->Val[2] + V1->Val[1]*V2->Val[4] + V1->Val[2]*V2->Val[5];
896 a1[3] = V1->Val[1]*V2->Val[0] + V1->Val[3]*V2->Val[1] + V1->Val[4]*V2->Val[2];
897 a1[4] = V1->Val[1]*V2->Val[1] + V1->Val[3]*V2->Val[3] + V1->Val[4]*V2->Val[4];
898 a1[5] = V1->Val[1]*V2->Val[2] + V1->Val[3]*V2->Val[4] + V1->Val[4]*V2->Val[5];
899 a1[6] = V1->Val[2]*V2->Val[0] + V1->Val[4]*V2->Val[1] + V1->Val[5]*V2->Val[2];
900 a1[7] = V1->Val[2]*V2->Val[1] + V1->Val[4]*V2->Val[3] + V1->Val[5]*V2->Val[4];
901 a1[8] = V1->Val[2]*V2->Val[2] + V1->Val[4]*V2->Val[4] + V1->Val[5]*V2->Val[5];
902 R->Val[0] = a1[0]; R->Val[1] = a1[1]; R->Val[2] = a1[2];
903 R->Val[3] = a1[3]; R->Val[4] = a1[4]; R->Val[5] = a1[5];
904 R->Val[6] = a1[6]; R->Val[7] = a1[7]; R->Val[8] = a1[8];
905 }
906 else {
907 for (k = 0 ; k < Current.NbrHar ; k += 2) {
908 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
909 CMULT(0,1,3); CMULT(1,3,4); CMULT(2,4,5);
910 CMULT(0,2,6); CMULT(1,4,7); CMULT(2,5,8);
911 CMULT(1,0,9); CMULT(3,1,10); CMULT(4,2,11);
912 CMULT(1,1,12); CMULT(3,3,13); CMULT(4,4,14);
913 CMULT(1,2,15); CMULT(3,4,16); CMULT(4,5,17);
914 CMULT(2,0,18); CMULT(4,1,19); CMULT(5,2,20);
915 CMULT(2,1,21); CMULT(4,3,22); CMULT(5,4,23);
916 CMULT(2,2,24); CMULT(4,4,25); CMULT(5,5,26);
917 CPUT3(0,1,2,0); CPUT3(3,4,5,1); CPUT3(6,7,8,2);
918 CPUT3(9,10,11,3); CPUT3(12,13,14,4); CPUT3(15,16,17,5);
919 CPUT3(18,19,20,6); CPUT3(21,22,23,7); CPUT3(24,25,26,8);
920 }
921 }
922 R->Type = TENSOR;
923 }
924 else if (V1->Type == TENSOR && V2->Type == TENSOR) {
925 if (Current.NbrHar == 1) {
926 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[3] + V1->Val[2]*V2->Val[6];
927 a1[1] = V1->Val[0]*V2->Val[1] + V1->Val[1]*V2->Val[4] + V1->Val[2]*V2->Val[7];
928 a1[2] = V1->Val[0]*V2->Val[2] + V1->Val[1]*V2->Val[5] + V1->Val[2]*V2->Val[8];
929 a1[3] = V1->Val[3]*V2->Val[0] + V1->Val[4]*V2->Val[3] + V1->Val[5]*V2->Val[6];
930 a1[4] = V1->Val[3]*V2->Val[1] + V1->Val[4]*V2->Val[4] + V1->Val[5]*V2->Val[7];
931 a1[5] = V1->Val[3]*V2->Val[2] + V1->Val[4]*V2->Val[5] + V1->Val[5]*V2->Val[8];
932 a1[6] = V1->Val[6]*V2->Val[0] + V1->Val[7]*V2->Val[3] + V1->Val[8]*V2->Val[6];
933 a1[7] = V1->Val[6]*V2->Val[1] + V1->Val[7]*V2->Val[4] + V1->Val[8]*V2->Val[7];
934 a1[8] = V1->Val[6]*V2->Val[2] + V1->Val[7]*V2->Val[5] + V1->Val[8]*V2->Val[8];
935 R->Val[0] = a1[0]; R->Val[1] = a1[1]; R->Val[2] = a1[2];
936 R->Val[3] = a1[3]; R->Val[4] = a1[4]; R->Val[5] = a1[5];
937 R->Val[6] = a1[6]; R->Val[7] = a1[7]; R->Val[8] = a1[8];
938 }
939 else {
940 for (k = 0 ; k < Current.NbrHar ; k += 2) {
941 CMULT(0,0,0); CMULT(1,3,1); CMULT(2,6,2);
942 CMULT(0,1,3); CMULT(1,4,4); CMULT(2,7,5);
943 CMULT(0,2,6); CMULT(1,5,7); CMULT(2,8,8);
944 CMULT(3,0,9); CMULT(4,3,10); CMULT(5,6,11);
945 CMULT(3,1,12); CMULT(4,4,13); CMULT(5,7,14);
946 CMULT(3,2,15); CMULT(4,5,16); CMULT(5,8,17);
947 CMULT(6,0,18); CMULT(7,3,19); CMULT(8,6,20);
948 CMULT(6,1,21); CMULT(7,4,22); CMULT(8,7,23);
949 CMULT(6,2,24); CMULT(7,5,25); CMULT(8,8,26);
950 CPUT3(0,1,2,0); CPUT3(3,4,5,1); CPUT3(6,7,8,2);
951 CPUT3(9,10,11,3); CPUT3(12,13,14,4); CPUT3(15,16,17,5);
952 CPUT3(18,19,20,6); CPUT3(21,22,23,7); CPUT3(24,25,26,8);
953 }
954 }
955 R->Type = TENSOR;
956 }
957 else if (V1->Type == TENSOR_SYM && V2->Type == TENSOR) {
958 if (Current.NbrHar == 1) {
959 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[3] + V1->Val[2]*V2->Val[6];
960 a1[1] = V1->Val[0]*V2->Val[1] + V1->Val[1]*V2->Val[4] + V1->Val[2]*V2->Val[7];
961 a1[2] = V1->Val[0]*V2->Val[2] + V1->Val[1]*V2->Val[5] + V1->Val[2]*V2->Val[8];
962 a1[3] = V1->Val[1]*V2->Val[0] + V1->Val[3]*V2->Val[3] + V1->Val[4]*V2->Val[6];
963 a1[4] = V1->Val[1]*V2->Val[1] + V1->Val[3]*V2->Val[4] + V1->Val[4]*V2->Val[7];
964 a1[5] = V1->Val[1]*V2->Val[2] + V1->Val[3]*V2->Val[5] + V1->Val[4]*V2->Val[8];
965 a1[6] = V1->Val[2]*V2->Val[0] + V1->Val[4]*V2->Val[3] + V1->Val[5]*V2->Val[6];
966 a1[7] = V1->Val[2]*V2->Val[1] + V1->Val[4]*V2->Val[4] + V1->Val[5]*V2->Val[7];
967 a1[8] = V1->Val[2]*V2->Val[2] + V1->Val[4]*V2->Val[5] + V1->Val[5]*V2->Val[8];
968 R->Val[0] = a1[0]; R->Val[1] = a1[1]; R->Val[2] = a1[2];
969 R->Val[3] = a1[3]; R->Val[4] = a1[4]; R->Val[5] = a1[5];
970 R->Val[6] = a1[6]; R->Val[7] = a1[7]; R->Val[8] = a1[8];
971 }
972 else {
973 for (k = 0 ; k < Current.NbrHar ; k += 2) {
974 CMULT(0,0,0); CMULT(1,3,1); CMULT(2,6,2);
975 CMULT(0,1,3); CMULT(1,4,4); CMULT(2,7,5);
976 CMULT(0,2,6); CMULT(1,5,7); CMULT(2,8,8);
977 CMULT(1,0,9); CMULT(2,3,10); CMULT(4,6,11);
978 CMULT(1,1,12); CMULT(2,4,13); CMULT(4,7,14);
979 CMULT(1,2,15); CMULT(2,5,16); CMULT(4,8,17);
980 CMULT(3,0,18); CMULT(4,3,19); CMULT(5,6,20);
981 CMULT(3,1,21); CMULT(4,4,22); CMULT(5,7,23);
982 CMULT(3,2,24); CMULT(4,5,25); CMULT(5,8,26);
983 CPUT3(0,1,2,0); CPUT3(3,4,5,1); CPUT3(6,7,8,2);
984 CPUT3(9,10,11,3); CPUT3(12,13,14,4); CPUT3(15,16,17,5);
985 CPUT3(18,19,20,6); CPUT3(21,22,23,7); CPUT3(24,25,26,8);
986 }
987 }
988 R->Type = TENSOR;
989 }
990 else if (V1->Type == TENSOR && V2->Type == TENSOR_SYM) {
991 if (Current.NbrHar == 1) {
992 a1[0] = V1->Val[0]*V2->Val[0] + V1->Val[1]*V2->Val[1] + V1->Val[2]*V2->Val[2];
993 a1[1] = V1->Val[0]*V2->Val[1] + V1->Val[1]*V2->Val[3] + V1->Val[2]*V2->Val[4];
994 a1[2] = V1->Val[0]*V2->Val[2] + V1->Val[1]*V2->Val[4] + V1->Val[2]*V2->Val[5];
995 a1[3] = V1->Val[3]*V2->Val[0] + V1->Val[4]*V2->Val[1] + V1->Val[5]*V2->Val[2];
996 a1[4] = V1->Val[3]*V2->Val[1] + V1->Val[4]*V2->Val[3] + V1->Val[5]*V2->Val[4];
997 a1[5] = V1->Val[3]*V2->Val[2] + V1->Val[4]*V2->Val[4] + V1->Val[5]*V2->Val[5];
998 a1[6] = V1->Val[6]*V2->Val[0] + V1->Val[7]*V2->Val[1] + V1->Val[8]*V2->Val[2];
999 a1[7] = V1->Val[6]*V2->Val[1] + V1->Val[7]*V2->Val[3] + V1->Val[8]*V2->Val[4];
1000 a1[8] = V1->Val[6]*V2->Val[2] + V1->Val[7]*V2->Val[4] + V1->Val[8]*V2->Val[5];
1001 R->Val[0] = a1[0]; R->Val[1] = a1[1]; R->Val[2] = a1[2];
1002 R->Val[3] = a1[3]; R->Val[4] = a1[4]; R->Val[5] = a1[5];
1003 R->Val[6] = a1[6]; R->Val[7] = a1[7]; R->Val[8] = a1[8];
1004 }
1005 else {
1006 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1007 CMULT(0,0,0); CMULT(1,1,1); CMULT(2,2,2);
1008 CMULT(0,1,3); CMULT(1,3,4); CMULT(2,4,5);
1009 CMULT(0,2,6); CMULT(1,4,7); CMULT(2,5,8);
1010 CMULT(3,0,9); CMULT(4,1,10); CMULT(5,2,11);
1011 CMULT(3,1,12); CMULT(4,3,13); CMULT(5,4,14);
1012 CMULT(3,2,15); CMULT(4,4,16); CMULT(5,5,17);
1013 CMULT(6,0,18); CMULT(7,1,19); CMULT(8,2,20);
1014 CMULT(6,1,21); CMULT(7,3,22); CMULT(8,4,23);
1015 CMULT(6,2,24); CMULT(7,4,25); CMULT(8,5,26);
1016 CPUT3(0,1,2,0); CPUT3(3,4,5,1); CPUT3(6,7,8,2);
1017 CPUT3(9,10,11,3); CPUT3(12,13,14,4); CPUT3(15,16,17,5);
1018 CPUT3(18,19,20,6); CPUT3(21,22,23,7); CPUT3(24,25,26,8);
1019 }
1020 }
1021 R->Type = TENSOR;
1022 }
1023
1024 /* a faire: differents tenseurs entre eux */
1025
1026 else {
1027 Message::Error("Product of non adapted quantities: %s * %s",
1028 Get_StringForDefine(Field_Type, V1->Type),
1029 Get_StringForDefine(Field_Type, V2->Type));
1030 }
1031
1032 }
1033
1034 #undef CMULT
1035 #undef CPUT
1036 #undef CPUT3
1037
1038 /* ------------------------------------------------------------------------
1039 R <- V1 / V2
1040 ------------------------------------------------------------------------ */
1041
1042 #define CDIVI(a,b,c) \
1043 Cal_ComplexDivision(&(V1->Val[MAX_DIM*k+a]), &(V2->Val[MAX_DIM*k+b]), tmp[c])
1044
1045 #define CPUT(a) \
1046 R->Val[MAX_DIM* k +a] = tmp[a][0] ; \
1047 R->Val[MAX_DIM*(k+1)+a] = tmp[a][MAX_DIM]
1048
Cal_DivideValue(struct Value * V1,struct Value * V2,struct Value * R)1049 void Cal_DivideValue(struct Value * V1, struct Value * V2, struct Value * R)
1050 {
1051 int k ;
1052 struct Value V3 ;
1053
1054 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1055 if (Current.NbrHar == 1) {
1056 R->Val[0] = V1->Val[0]/V2->Val[0];
1057 }
1058 else {
1059 for (k = 0 ; k < Current.NbrHar ; k += 2) { /* meaning in multi-harmonics ??? */
1060 CDIVI(0,0,0);
1061 CPUT(0);
1062 }
1063 }
1064 R->Type = SCALAR ;
1065 }
1066
1067 else if ( (V1->Type == VECTOR || V1->Type == TENSOR_DIAG) && (V2->Type == SCALAR) ) {
1068 if (Current.NbrHar == 1) {
1069 a0 = V2->Val[0] ;
1070 R->Val[0] = V1->Val[0] / a0 ;
1071 R->Val[1] = V1->Val[1] / a0 ;
1072 R->Val[2] = V1->Val[2] / a0 ;
1073 }
1074 else {
1075 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1076 CDIVI(0,0,0); CDIVI(1,0,1); CDIVI(2,0,2);
1077 CPUT(0); CPUT(1); CPUT(2);
1078 }
1079 }
1080 R->Type = V1->Type ;
1081 }
1082
1083 else if ( (V1->Type == TENSOR_SYM) && (V2->Type == SCALAR) ) {
1084 if (Current.NbrHar == 1) {
1085 a0 = V2->Val[0] ;
1086 R->Val[0] = V1->Val[0] / a0 ;
1087 R->Val[1] = V1->Val[1] / a0 ;
1088 R->Val[2] = V1->Val[2] / a0 ;
1089 R->Val[3] = V1->Val[3] / a0 ;
1090 R->Val[4] = V1->Val[4] / a0 ;
1091 R->Val[5] = V1->Val[5] / a0 ;
1092 }
1093 else {
1094 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1095 CDIVI(0,0,0); CDIVI(1,0,1); CDIVI(2,0,2); CDIVI(3,0,3); CDIVI(4,0,4); CDIVI(5,0,5);
1096 CPUT(0); CPUT(1); CPUT(2); CPUT(3); CPUT(4); CPUT(5);
1097 }
1098 }
1099 R->Type = TENSOR_SYM ;
1100 }
1101
1102 else if ( (V1->Type == TENSOR) && (V2->Type == SCALAR) ) {
1103 if (Current.NbrHar == 1) {
1104 a0 = V2->Val[0] ;
1105 R->Val[0] = V1->Val[0] / a0 ;
1106 R->Val[1] = V1->Val[1] / a0 ;
1107 R->Val[2] = V1->Val[2] / a0 ;
1108 R->Val[3] = V1->Val[3] / a0 ;
1109 R->Val[4] = V1->Val[4] / a0 ;
1110 R->Val[5] = V1->Val[5] / a0 ;
1111 R->Val[6] = V1->Val[6] / a0 ;
1112 R->Val[7] = V1->Val[7] / a0 ;
1113 R->Val[8] = V1->Val[8] / a0 ;
1114 }
1115 else {
1116 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1117 CDIVI(0,0,0); CDIVI(1,0,1); CDIVI(2,0,2); CDIVI(3,0,3); CDIVI(4,0,4); CDIVI(5,0,5);
1118 CDIVI(6,0,6); CDIVI(7,0,7); CDIVI(8,0,8);
1119 CPUT(0); CPUT(1); CPUT(2); CPUT(3); CPUT(4); CPUT(5);
1120 CPUT(6); CPUT(7); CPUT(8);
1121 }
1122 }
1123 R->Type = TENSOR ;
1124 }
1125
1126 else if ( (V1->Type == SCALAR) &&
1127 (V2->Type == TENSOR || V2->Type == TENSOR_SYM || V2->Type == TENSOR_DIAG) ) {
1128 Cal_InvertValue(V2,&V3);
1129 Cal_ProductValue(V1,&V3,R);
1130 }
1131
1132 else {
1133 Message::Error("Division of non adapted quantities: %s / %s",
1134 Get_StringForDefine(Field_Type, V1->Type),
1135 Get_StringForDefine(Field_Type, V2->Type));
1136 }
1137
1138 }
1139
1140 #undef CDIVI
1141 #undef CPUT
1142
1143 /* ------------------------------------------------------------------------
1144 R <- V1 % V2
1145 ------------------------------------------------------------------------ */
1146
Cal_ModuloValue(struct Value * V1,struct Value * V2,struct Value * R)1147 void Cal_ModuloValue(struct Value * V1, struct Value * V2, struct Value * R)
1148 {
1149 int k ;
1150
1151 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1152 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1153 R->Val[MAX_DIM* k ] = (int)V1->Val[MAX_DIM*k] % (int)V2->Val[MAX_DIM*k] ;
1154 R->Val[MAX_DIM*(k+1)] = 0. ;
1155 }
1156 R->Type = SCALAR ;
1157 }
1158
1159 else if ( (V1->Type == VECTOR) && (V2->Type == SCALAR) ) {
1160 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1161 R->Val[MAX_DIM* k ] = (int)V1->Val[MAX_DIM*k ] % (int)V2->Val[MAX_DIM*k ] ;
1162 R->Val[MAX_DIM* k +1] = (int)V1->Val[MAX_DIM*k+1] % (int)V2->Val[MAX_DIM*k+1] ;
1163 R->Val[MAX_DIM* k +2] = (int)V1->Val[MAX_DIM*k+2] % (int)V2->Val[MAX_DIM*k+2] ;
1164 R->Val[MAX_DIM*(k+1)] = 0. ;
1165 R->Val[MAX_DIM*(k+1)+1] = 0. ;
1166 R->Val[MAX_DIM*(k+1)+2] = 0. ;
1167 }
1168 R->Type = VECTOR ;
1169 }
1170
1171 else {
1172 Message::Error("Modulo of non adapted quantities: %s %% %s",
1173 Get_StringForDefine(Field_Type, V1->Type),
1174 Get_StringForDefine(Field_Type, V2->Type));
1175 }
1176
1177 }
1178
1179 /* ------------------------------------------------------------------------
1180 R <- V1 X V2
1181 ------------------------------------------------------------------------ */
1182
Cal_CrossProductValue(struct Value * V1,struct Value * V2,struct Value * R)1183 void Cal_CrossProductValue (struct Value * V1, struct Value * V2, struct Value * R)
1184 {
1185 int k ;
1186
1187 if ( (V1->Type == VECTOR) && (V2->Type == VECTOR) ) {
1188 if (Current.NbrHar == 1) {
1189 a1[0] = V1->Val[1] * V2->Val[2] - V1->Val[2] * V2->Val[1] ;
1190 a1[1] = V1->Val[2] * V2->Val[0] - V1->Val[0] * V2->Val[2] ;
1191 a1[2] = V1->Val[0] * V2->Val[1] - V1->Val[1] * V2->Val[0] ;
1192 R->Val[0] = a1[0] ; R->Val[1] = a1[1] ; R->Val[2] = a1[2] ;
1193 }
1194 else {
1195 for (k = 0 ; k < Current.NbrHar ; k += 2) {
1196 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k+1]), &(V2->Val[MAX_DIM*k+2]), a1) ;
1197 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k+2]), &(V2->Val[MAX_DIM*k+1]), a2) ;
1198 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
1199
1200 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k+2]), &(V2->Val[MAX_DIM*k ]), a1) ;
1201 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k ]), &(V2->Val[MAX_DIM*k+2]), a2) ;
1202 b2[0] = a1[0] - a2[0] ; b2[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
1203
1204 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k ]), &(V2->Val[MAX_DIM*k+1]), a1) ;
1205 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k+1]), &(V2->Val[MAX_DIM*k ]), a2) ;
1206 b3[0] = a1[0] - a2[0] ; b3[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
1207
1208 R->Val[MAX_DIM*k ] = b1[0] ; R->Val[MAX_DIM*(k+1) ] = b1[MAX_DIM] ;
1209 R->Val[MAX_DIM*k+1] = b2[0] ; R->Val[MAX_DIM*(k+1)+1] = b2[MAX_DIM] ;
1210 R->Val[MAX_DIM*k+2] = b3[0] ; R->Val[MAX_DIM*(k+1)+2] = b3[MAX_DIM] ;
1211 }
1212 }
1213 R->Type = VECTOR ;
1214 }
1215
1216 else {
1217 Message::Error("Cross product of non vector quantities: %s /\\ %s",
1218 Get_StringForDefine(Field_Type, V1->Type),
1219 Get_StringForDefine(Field_Type, V2->Type));
1220 }
1221
1222 }
1223
1224 /* ------------------------------------------------------------------------
1225 R <- SQRT(V1)
1226 ------------------------------------------------------------------------ */
1227
Cal_SqrtValue(struct Value * V1,struct Value * R)1228 void Cal_SqrtValue(struct Value *V1, struct Value *R)
1229 {
1230 if( V1->Type == SCALAR ){
1231 struct Value P;
1232 P.Type = SCALAR;
1233 P.Val[0] = 0.5;
1234 Cal_PowerValue(V1, &P, R);
1235 }
1236 else {
1237 Message::Error("Square root of non scalar quantity: %s",
1238 Get_StringForDefine(Field_Type, V1->Type));
1239 }
1240 }
1241
1242 /* ------------------------------------------------------------------------
1243 R <- V1 ^ V2
1244 ------------------------------------------------------------------------ */
1245
Cal_PowerValue(struct Value * V1,struct Value * V2,struct Value * R)1246 void Cal_PowerValue(struct Value * V1, struct Value * V2, struct Value * R)
1247 {
1248 int k;
1249 double arg, abs ;
1250
1251 if ( V1->Type == SCALAR && V2->Type == SCALAR ){
1252
1253 if(V2->Val[0] == 1.){
1254 Cal_CopyValue(V1,R) ;
1255 }
1256 if(V2->Val[0] == 2.){
1257 if (Current.NbrHar == 1) {
1258 R->Val[0] = SQU(V1->Val[0]) ;
1259 }
1260 else{
1261 for (k = 0 ; k < Current.NbrHar ; k+=2) {
1262 Cal_ComplexProduct(&(V1->Val[MAX_DIM*k]), &(V1->Val[MAX_DIM*k]), a1) ;
1263 R->Val[MAX_DIM* k ] = a1[0];
1264 R->Val[MAX_DIM*(k+1)] = a1[MAX_DIM];
1265 }
1266 }
1267 }
1268 else{
1269 if (Current.NbrHar == 1) {
1270 R->Val[0] = pow(V1->Val[0],V2->Val[0]) ;
1271 }
1272 else{
1273 for (k = 0 ; k < Current.NbrHar ; k+=2) {
1274 abs = pow(sqrt(SQU(V1->Val[MAX_DIM*k])+SQU(V1->Val[MAX_DIM*(k+1)])),
1275 V2->Val[0]) ;
1276 arg = atan2(V1->Val[MAX_DIM*(k+1)], V1->Val[MAX_DIM*k]) ;
1277 R->Val[MAX_DIM* k ] = abs * cos(V2->Val[0] * arg) ;
1278 R->Val[MAX_DIM*(k+1)] = abs * sin(V2->Val[0] * arg) ;
1279 }
1280 }
1281 }
1282 R->Type = SCALAR ;
1283 }
1284
1285 else {
1286 Message::Error("Power of non scalar quantities: %s ^ %s",
1287 Get_StringForDefine(Field_Type, V1->Type),
1288 Get_StringForDefine(Field_Type, V2->Type));
1289 }
1290
1291 }
1292
1293 /* ------------------------------------------------------------------------
1294 R <- V1 < V2
1295 ------------------------------------------------------------------------ */
1296
Cal_LessValue(struct Value * V1,struct Value * V2,struct Value * R)1297 void Cal_LessValue(struct Value * V1, struct Value * V2, struct Value * R)
1298 {
1299 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1300 R->Val[0] = (V1->Val[0] < V2->Val[0]) ;
1301 R->Type = SCALAR ;
1302 }
1303 else {
1304 Message::Error("Comparison of non scalar quantities: %s < %s",
1305 Get_StringForDefine(Field_Type, V1->Type),
1306 Get_StringForDefine(Field_Type, V2->Type));
1307 }
1308 }
1309
1310 /* ------------------------------------------------------------------------
1311 R <- V1 <= V2
1312 ------------------------------------------------------------------------ */
1313
Cal_LessOrEqualValue(struct Value * V1,struct Value * V2,struct Value * R)1314 void Cal_LessOrEqualValue(struct Value * V1, struct Value * V2, struct Value * R)
1315 {
1316 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1317 R->Val[0] = (V1->Val[0] <= V2->Val[0]) ;
1318 R->Type = SCALAR ;
1319 }
1320 else {
1321 Message::Error("Comparison of non scalar quantities: %s <= %s",
1322 Get_StringForDefine(Field_Type, V1->Type),
1323 Get_StringForDefine(Field_Type, V2->Type));
1324 }
1325 }
1326
1327 /* ------------------------------------------------------------------------
1328 R <- V1 > V2
1329 ------------------------------------------------------------------------ */
1330
Cal_GreaterValue(struct Value * V1,struct Value * V2,struct Value * R)1331 void Cal_GreaterValue(struct Value * V1, struct Value * V2, struct Value * R)
1332 {
1333 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1334 R->Val[0] = (V1->Val[0] > V2->Val[0]) ;
1335 R->Type = SCALAR ;
1336 }
1337 else {
1338 Message::Error("Comparison of non scalar quantities: %s > %s",
1339 Get_StringForDefine(Field_Type, V1->Type),
1340 Get_StringForDefine(Field_Type, V2->Type));
1341 }
1342 }
1343
1344 /* ------------------------------------------------------------------------
1345 R <- V1 >= V2
1346 ------------------------------------------------------------------------ */
1347
Cal_GreaterOrEqualValue(struct Value * V1,struct Value * V2,struct Value * R)1348 void Cal_GreaterOrEqualValue(struct Value * V1, struct Value * V2, struct Value * R)
1349 {
1350 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1351 R->Val[0] = (V1->Val[0] >= V2->Val[0]) ;
1352 R->Type = SCALAR ;
1353 }
1354 else {
1355 Message::Error("Comparison of non scalar quantities: %s >= %s",
1356 Get_StringForDefine(Field_Type, V1->Type),
1357 Get_StringForDefine(Field_Type, V2->Type));
1358 }
1359 }
1360
1361 /* ------------------------------------------------------------------------
1362 R <- V1 == V2
1363 ------------------------------------------------------------------------ */
1364
Cal_EqualValue(struct Value * V1,struct Value * V2,struct Value * R)1365 void Cal_EqualValue(struct Value * V1, struct Value * V2, struct Value * R)
1366 {
1367 int k;
1368
1369 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1370 R->Val[0] = (V1->Val[0] == V2->Val[0]) ;
1371 for (k = 1 ; k < Current.NbrHar ; k++){
1372 if(!R->Val[0]) break;
1373 R->Val[0] = (V1->Val[MAX_DIM*k] == V2->Val[MAX_DIM*k]) ;
1374 }
1375 R->Type = SCALAR ;
1376 }
1377 else if ( ( (V1->Type == VECTOR) && (V2->Type == VECTOR) ) ||
1378 ( (V1->Type == TENSOR_DIAG) && (V2->Type == TENSOR_DIAG) ) ) {
1379 R->Val[0] = (V1->Val[0] == V2->Val[0] &&
1380 V1->Val[1] == V2->Val[1] &&
1381 V1->Val[2] == V2->Val[2]) ;
1382 for (k = 0 ; k < Current.NbrHar ; k++) {
1383 if(!R->Val[0]) break;
1384 R->Val[0] = (V1->Val[MAX_DIM*k ] == V2->Val[MAX_DIM*k ] &&
1385 V1->Val[MAX_DIM*k+1] == V2->Val[MAX_DIM*k+1] &&
1386 V1->Val[MAX_DIM*k+2] == V2->Val[MAX_DIM*k+2]) ;
1387 }
1388 R->Type = SCALAR ;
1389 }
1390 else if ( (V1->Type == TENSOR_SYM) && (V2->Type == TENSOR_SYM) ) {
1391 R->Val[0] = (V1->Val[0] == V2->Val[0] &&
1392 V1->Val[1] == V2->Val[1] &&
1393 V1->Val[2] == V2->Val[2] &&
1394 V1->Val[3] == V2->Val[3] &&
1395 V1->Val[4] == V2->Val[4] &&
1396 V1->Val[5] == V2->Val[5]) ;
1397 for (k = 0 ; k < Current.NbrHar ; k++) {
1398 if(!R->Val[0]) break;
1399 R->Val[0] = (V1->Val[MAX_DIM*k ] == V2->Val[MAX_DIM*k ] &&
1400 V1->Val[MAX_DIM*k+1] == V2->Val[MAX_DIM*k+1] &&
1401 V1->Val[MAX_DIM*k+2] == V2->Val[MAX_DIM*k+2] &&
1402 V1->Val[MAX_DIM*k+3] == V2->Val[MAX_DIM*k+3] &&
1403 V1->Val[MAX_DIM*k+4] == V2->Val[MAX_DIM*k+4] &&
1404 V1->Val[MAX_DIM*k+5] == V2->Val[MAX_DIM*k+5]) ;
1405 }
1406 R->Type = SCALAR ;
1407 }
1408 else if ( (V1->Type == TENSOR) && (V2->Type == TENSOR) ) {
1409 R->Val[0] = (V1->Val[0] == V2->Val[0] &&
1410 V1->Val[1] == V2->Val[1] &&
1411 V1->Val[2] == V2->Val[2] &&
1412 V1->Val[3] == V2->Val[3] &&
1413 V1->Val[4] == V2->Val[4] &&
1414 V1->Val[5] == V2->Val[5] &&
1415 V1->Val[6] == V2->Val[6] &&
1416 V1->Val[7] == V2->Val[7] &&
1417 V1->Val[8] == V2->Val[8]) ;
1418 for (k = 0 ; k < Current.NbrHar ; k++) {
1419 if(!R->Val[0]) break;
1420 R->Val[0] = (V1->Val[MAX_DIM*k ] == V2->Val[MAX_DIM*k ] &&
1421 V1->Val[MAX_DIM*k+1] == V2->Val[MAX_DIM*k+1] &&
1422 V1->Val[MAX_DIM*k+2] == V2->Val[MAX_DIM*k+2] &&
1423 V1->Val[MAX_DIM*k+3] == V2->Val[MAX_DIM*k+3] &&
1424 V1->Val[MAX_DIM*k+4] == V2->Val[MAX_DIM*k+4] &&
1425 V1->Val[MAX_DIM*k+5] == V2->Val[MAX_DIM*k+5] &&
1426 V1->Val[MAX_DIM*k+6] == V2->Val[MAX_DIM*k+6] &&
1427 V1->Val[MAX_DIM*k+7] == V2->Val[MAX_DIM*k+7] &&
1428 V1->Val[MAX_DIM*k+8] == V2->Val[MAX_DIM*k+8]) ;
1429 }
1430 R->Type = SCALAR ;
1431 }
1432 else {
1433 Message::Error("Comparison of different quantities: %s == %s",
1434 Get_StringForDefine(Field_Type, V1->Type),
1435 Get_StringForDefine(Field_Type, V2->Type));
1436 }
1437 }
1438
1439 /* ------------------------------------------------------------------------
1440 R <- V1 != V2
1441 ------------------------------------------------------------------------ */
1442
Cal_NotEqualValue(struct Value * V1,struct Value * V2,struct Value * R)1443 void Cal_NotEqualValue(struct Value * V1, struct Value * V2, struct Value * R)
1444 {
1445 int k;
1446
1447 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1448 R->Val[0] = (V1->Val[0] != V2->Val[0]) ;
1449 for (k = 1 ; k < Current.NbrHar ; k++){
1450 if(R->Val[0]) break;
1451 R->Val[0] = (V1->Val[MAX_DIM*k] != V2->Val[MAX_DIM*k]) ;
1452 }
1453 R->Type = SCALAR ;
1454 }
1455 else if ( ( (V1->Type == VECTOR) && (V2->Type == VECTOR) ) ||
1456 ( (V1->Type == TENSOR_DIAG) && (V2->Type == TENSOR_DIAG) ) ) {
1457 R->Val[0] = (V1->Val[0] != V2->Val[0] ||
1458 V1->Val[1] != V2->Val[1] ||
1459 V1->Val[2] != V2->Val[2]) ;
1460 for (k = 0 ; k < Current.NbrHar ; k++) {
1461 if(R->Val[0]) break;
1462 R->Val[0] = (V1->Val[MAX_DIM*k ] != V2->Val[MAX_DIM*k ] ||
1463 V1->Val[MAX_DIM*k+1] != V2->Val[MAX_DIM*k+1] ||
1464 V1->Val[MAX_DIM*k+2] != V2->Val[MAX_DIM*k+2]) ;
1465 }
1466 R->Type = SCALAR ;
1467 }
1468 else if ( (V1->Type == TENSOR_SYM) && (V2->Type == TENSOR_SYM) ) {
1469 R->Val[0] = (V1->Val[0] != V2->Val[0] ||
1470 V1->Val[1] != V2->Val[1] ||
1471 V1->Val[2] != V2->Val[2] ||
1472 V1->Val[3] != V2->Val[3] ||
1473 V1->Val[4] != V2->Val[4] ||
1474 V1->Val[5] != V2->Val[5]) ;
1475 for (k = 0 ; k < Current.NbrHar ; k++) {
1476 if(R->Val[0]) break;
1477 R->Val[0] = (V1->Val[MAX_DIM*k ] != V2->Val[MAX_DIM*k ] ||
1478 V1->Val[MAX_DIM*k+1] != V2->Val[MAX_DIM*k+1] ||
1479 V1->Val[MAX_DIM*k+2] != V2->Val[MAX_DIM*k+2] ||
1480 V1->Val[MAX_DIM*k+3] != V2->Val[MAX_DIM*k+3] ||
1481 V1->Val[MAX_DIM*k+4] != V2->Val[MAX_DIM*k+4] ||
1482 V1->Val[MAX_DIM*k+5] != V2->Val[MAX_DIM*k+5]) ;
1483 }
1484 R->Type = SCALAR ;
1485 }
1486 else if ( (V1->Type == TENSOR) && (V2->Type == TENSOR) ) {
1487 R->Val[0] = (V1->Val[0] != V2->Val[0] ||
1488 V1->Val[1] != V2->Val[1] ||
1489 V1->Val[2] != V2->Val[2] ||
1490 V1->Val[3] != V2->Val[3] ||
1491 V1->Val[4] != V2->Val[4] ||
1492 V1->Val[5] != V2->Val[5] ||
1493 V1->Val[6] != V2->Val[6] ||
1494 V1->Val[7] != V2->Val[7] ||
1495 V1->Val[8] != V2->Val[8]) ;
1496 for (k = 0 ; k < Current.NbrHar ; k++) {
1497 if(R->Val[0]) break;
1498 R->Val[0] = (V1->Val[MAX_DIM*k ] != V2->Val[MAX_DIM*k ] ||
1499 V1->Val[MAX_DIM*k+1] != V2->Val[MAX_DIM*k+1] ||
1500 V1->Val[MAX_DIM*k+2] != V2->Val[MAX_DIM*k+2] ||
1501 V1->Val[MAX_DIM*k+3] != V2->Val[MAX_DIM*k+3] ||
1502 V1->Val[MAX_DIM*k+4] != V2->Val[MAX_DIM*k+4] ||
1503 V1->Val[MAX_DIM*k+5] != V2->Val[MAX_DIM*k+5] ||
1504 V1->Val[MAX_DIM*k+6] != V2->Val[MAX_DIM*k+6] ||
1505 V1->Val[MAX_DIM*k+7] != V2->Val[MAX_DIM*k+7] ||
1506 V1->Val[MAX_DIM*k+8] != V2->Val[MAX_DIM*k+8]) ;
1507 }
1508 R->Type = SCALAR ;
1509 }
1510 else {
1511 Message::Error("Comparison of different quantities: %s != %s",
1512 Get_StringForDefine(Field_Type, V1->Type),
1513 Get_StringForDefine(Field_Type, V2->Type));
1514 }
1515 }
1516
1517 /* ------------------------------------------------------------------------
1518 R <- V1 ~= V2
1519 ------------------------------------------------------------------------ */
1520
Cal_ApproxEqualValue(struct Value * V1,struct Value * V2,struct Value * R)1521 void Cal_ApproxEqualValue(struct Value * V1, struct Value * V2, struct Value * R)
1522 {
1523 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1524 R->Val[0] = (fabs(V1->Val[0] - V2->Val[0]) < 1.e-10) ;
1525 R->Type = SCALAR ;
1526 }
1527 else {
1528 Message::Error("Comparison of non scalar quantities: %s ~= %s",
1529 Get_StringForDefine(Field_Type, V1->Type),
1530 Get_StringForDefine(Field_Type, V2->Type));
1531 }
1532 }
1533
1534 /* ------------------------------------------------------------------------
1535 R <- V1 && V2
1536 ------------------------------------------------------------------------ */
1537
Cal_AndValue(struct Value * V1,struct Value * V2,struct Value * R)1538 void Cal_AndValue(struct Value * V1, struct Value * V2, struct Value * R)
1539 {
1540 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1541 R->Val[0] = (V1->Val[0] && V2->Val[0]) ;
1542 R->Type = SCALAR ;
1543 }
1544 else {
1545 Message::Error("And of non scalar quantities: %s && %s",
1546 Get_StringForDefine(Field_Type, V1->Type),
1547 Get_StringForDefine(Field_Type, V2->Type));
1548 }
1549 }
1550
1551 /* ------------------------------------------------------------------------
1552 R <- V1 || V2
1553 ------------------------------------------------------------------------ */
1554
Cal_OrValue(struct Value * V1,struct Value * V2,struct Value * R)1555 void Cal_OrValue(struct Value * V1, struct Value * V2, struct Value * R)
1556 {
1557 if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
1558 R->Val[0] = (V1->Val[0] || V2->Val[0]) ;
1559 R->Type = SCALAR ;
1560 }
1561 else {
1562 Message::Error("Or of non scalar quantities: %s || %s",
1563 Get_StringForDefine(Field_Type, V1->Type),
1564 Get_StringForDefine(Field_Type, V2->Type));
1565 }
1566 }
1567
1568 /* ------------------------------------------------------------------------
1569 R <- -R
1570 ------------------------------------------------------------------------ */
1571
Cal_NegValue(struct Value * R)1572 void Cal_NegValue(struct Value * R)
1573 {
1574 int k ;
1575
1576 switch(R->Type) {
1577 case SCALAR :
1578 for (k = 0 ; k < Current.NbrHar ; k++){
1579 R->Val[MAX_DIM*k] = -R->Val[MAX_DIM*k] ;
1580 }
1581 break;
1582 case VECTOR :
1583 case TENSOR_DIAG :
1584 for (k = 0 ; k < Current.NbrHar ; k++){
1585 R->Val[MAX_DIM*k ] = -R->Val[MAX_DIM*k ] ;
1586 R->Val[MAX_DIM*k+1] = -R->Val[MAX_DIM*k+1] ;
1587 R->Val[MAX_DIM*k+2] = -R->Val[MAX_DIM*k+2] ;
1588 }
1589 break;
1590 case TENSOR_SYM :
1591 for (k = 0 ; k < Current.NbrHar ; k++){
1592 R->Val[MAX_DIM*k ] = -R->Val[MAX_DIM*k ] ;
1593 R->Val[MAX_DIM*k+1] = -R->Val[MAX_DIM*k+1] ;
1594 R->Val[MAX_DIM*k+2] = -R->Val[MAX_DIM*k+2] ;
1595 R->Val[MAX_DIM*k+3] = -R->Val[MAX_DIM*k+3] ;
1596 R->Val[MAX_DIM*k+4] = -R->Val[MAX_DIM*k+4] ;
1597 R->Val[MAX_DIM*k+5] = -R->Val[MAX_DIM*k+5] ;
1598 }
1599 break;
1600 case TENSOR :
1601 for (k = 0 ; k < Current.NbrHar ; k++){
1602 R->Val[MAX_DIM*k ] = -R->Val[MAX_DIM*k ] ;
1603 R->Val[MAX_DIM*k+1] = -R->Val[MAX_DIM*k+1] ;
1604 R->Val[MAX_DIM*k+2] = -R->Val[MAX_DIM*k+2] ;
1605 R->Val[MAX_DIM*k+3] = -R->Val[MAX_DIM*k+3] ;
1606 R->Val[MAX_DIM*k+4] = -R->Val[MAX_DIM*k+4] ;
1607 R->Val[MAX_DIM*k+5] = -R->Val[MAX_DIM*k+5] ;
1608 R->Val[MAX_DIM*k+6] = -R->Val[MAX_DIM*k+6] ;
1609 R->Val[MAX_DIM*k+7] = -R->Val[MAX_DIM*k+7] ;
1610 R->Val[MAX_DIM*k+8] = -R->Val[MAX_DIM*k+8] ;
1611 }
1612 break;
1613 default :
1614 Message::Error("Wrong argument type for Operator (-)");
1615 break;
1616 }
1617 }
1618
1619 /* ------------------------------------------------------------------------
1620 R <- !R
1621 ------------------------------------------------------------------------ */
1622
Cal_NotValue(struct Value * R)1623 void Cal_NotValue(struct Value * R)
1624 {
1625 if (R->Type == SCALAR){
1626 R->Val[0] = !R->Val[0] ;
1627 }
1628 else {
1629 Message::Error("Negation of non scalar quantity: ! %s",
1630 Get_StringForDefine(Field_Type, R->Type));
1631 }
1632 }
1633
1634 /* ------------------------------------------------------------------------
1635 R <- V1^T
1636 ------------------------------------------------------------------------ */
1637
Cal_TransposeValue(struct Value * V1,struct Value * R)1638 void Cal_TransposeValue(struct Value *V1, struct Value *R)
1639 {
1640 int k;
1641
1642 switch(V1->Type){
1643
1644 case TENSOR_DIAG :
1645 case TENSOR_SYM :
1646 Cal_CopyValue(V1,R);
1647 break;
1648
1649 case TENSOR :
1650 R->Type = TENSOR;
1651 if(Current.NbrHar==1){
1652 R->Val[0] = V1->Val[0];
1653 R->Val[4] = V1->Val[4];
1654 R->Val[8] = V1->Val[8];
1655 a1[0] = V1->Val[1];
1656 a1[1] = V1->Val[2];
1657 a1[2] = V1->Val[5];
1658 R->Val[1] = V1->Val[3];
1659 R->Val[2] = V1->Val[6];
1660 R->Val[5] = V1->Val[7];
1661 R->Val[3] = a1[0];
1662 R->Val[6] = a1[1];
1663 R->Val[7] = a1[2];
1664 }
1665 else{
1666 for(k=0 ; k<Current.NbrHar ; k+=1){
1667 R->Val[MAX_DIM*k+0] = V1->Val[MAX_DIM*k+0];
1668 R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4];
1669 R->Val[MAX_DIM*k+8] = V1->Val[MAX_DIM*k+8];
1670 a1[0] = V1->Val[MAX_DIM*k+1];
1671 a1[1] = V1->Val[MAX_DIM*k+2];
1672 a1[2] = V1->Val[MAX_DIM*k+5];
1673 R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+3];
1674 R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+6];
1675 R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+7];
1676 R->Val[MAX_DIM*k+3] = a1[0];
1677 R->Val[MAX_DIM*k+6] = a1[1];
1678 R->Val[MAX_DIM*k+7] = a1[2];
1679 }
1680 }
1681 break;
1682
1683 default:
1684 Message::Error("Wrong argument in 'Cal_TransposeValue'");
1685 break;
1686 }
1687 }
1688
1689 /* ------------------------------------------------------------------------
1690 R <- Trace(V1)
1691 ------------------------------------------------------------------------ */
1692
Cal_TraceValue(struct Value * V1,struct Value * R)1693 void Cal_TraceValue(struct Value *V1, struct Value *R)
1694 {
1695 int k;
1696
1697 switch(V1->Type){
1698
1699 case TENSOR_DIAG :
1700 if (Current.NbrHar == 1) {
1701 R->Val[0] = V1->Val[0]+V1->Val[1]+V1->Val[2];
1702 }
1703 else {
1704 for (k = 0 ; k < Current.NbrHar ; k++) {
1705 R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k ]+
1706 V1->Val[MAX_DIM*k+1]+
1707 V1->Val[MAX_DIM*k+2];
1708 }
1709 }
1710 R->Type = SCALAR ;
1711 break;
1712
1713 case TENSOR_SYM :
1714 if (Current.NbrHar == 1) {
1715 R->Val[0] = V1->Val[0]+V1->Val[3]+V1->Val[5];
1716 }
1717 else {
1718 for (k = 0 ; k < Current.NbrHar ; k++) {
1719 R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k ]+
1720 V1->Val[MAX_DIM*k+3]+
1721 V1->Val[MAX_DIM*k+5];
1722 }
1723 }
1724 R->Type = SCALAR ;
1725 break;
1726
1727 case TENSOR :
1728 if (Current.NbrHar == 1) {
1729 R->Val[0] = V1->Val[0]+V1->Val[4]+V1->Val[8];
1730 }
1731 else {
1732 for (k = 0 ; k < Current.NbrHar ; k++) {
1733 R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k ]+
1734 V1->Val[MAX_DIM*k+4]+
1735 V1->Val[MAX_DIM*k+8];
1736 }
1737 }
1738 R->Type = SCALAR ;
1739 break;
1740
1741 default:
1742 Message::Error("Wrong argument type in 'Cal_TraceValue'");
1743 break;
1744 }
1745 }
1746
1747 /* ------------------------------------------------------------------------
1748 R <- V1^T * V2 * V1 , V1 real
1749 ------------------------------------------------------------------------ */
1750
1751 #define A0 V1->Val[0]
1752 #define A1 V1->Val[1]
1753 #define A2 V1->Val[2]
1754 #define A3 V1->Val[3]
1755 #define A4 V1->Val[4]
1756 #define A5 V1->Val[5]
1757 #define A6 V1->Val[6]
1758 #define A7 V1->Val[7]
1759 #define A8 V1->Val[8]
1760
Cal_RotateValue(struct Value * V1,struct Value * V2,struct Value * R)1761 void Cal_RotateValue(struct Value *V1, struct Value *V2, struct Value *R)
1762 {
1763 int k;
1764 double t0,t1,t2,t3,t4,t5,t6,t7,t8;
1765 struct Value A;
1766
1767 switch(V2->Type){
1768
1769 case VECTOR :
1770 if(Current.NbrHar == 1){
1771 #define B0 V2->Val[0]
1772 #define B1 V2->Val[1]
1773 #define B2 V2->Val[2]
1774 A.Val[0]= A0*B0+A1*B1+A2*B2;
1775 A.Val[1]= A3*B0+A4*B1+A5*B2;
1776 A.Val[2]= A6*B0+A7*B1+A8*B2;
1777 A.Type = VECTOR;
1778 Cal_CopyValue(&A,R);
1779 #undef B0
1780 #undef B1
1781 #undef B2
1782 }
1783 else{ /* Attention: a modifier */
1784 #define B0 V2->Val[0]
1785 #define B1 V2->Val[1]
1786 #define B2 V2->Val[2]
1787 A.Val[0]= A0*B0+A1*B1+A2*B2;
1788 A.Val[1]= A3*B0+A4*B1+A5*B2;
1789 A.Val[2]= A6*B0+A7*B1+A8*B2;
1790 A.Type = VECTOR;
1791 Cal_CopyValue(&A,R);
1792 #undef B0
1793 #undef B1
1794 #undef B2
1795 }
1796 break ;
1797
1798 case TENSOR_DIAG :
1799 if(Current.NbrHar == 1){
1800 #define B0 V2->Val[0]
1801 #define B1 V2->Val[1]
1802 #define B2 V2->Val[2]
1803 A.Val[0]= A0*A0*B0+A3*A3*B1+A6*A6*B2;
1804 A.Val[1]= A1*A0*B0+A3*A4*B1+A6*A7*B2;
1805 A.Val[2]= A2*A0*B0+A3*A5*B1+A6*A8*B2;
1806 A.Val[3]= A1*A0*B0+A3*A4*B1+A6*A7*B2;
1807 A.Val[4]= A1*A1*B0+A4*A4*B1+A7*A7*B2;
1808 A.Val[5]= A2*A1*B0+A4*A5*B1+A7*A8*B2;
1809 A.Val[6]= A2*A0*B0+A3*A5*B1+A6*A8*B2;
1810 A.Val[7]= A2*A1*B0+A4*A5*B1+A7*A8*B2;
1811 A.Val[8]= A2*A2*B0+A5*A5*B1+A8*A8*B2;
1812 A.Type = TENSOR;
1813 Cal_CopyValue(&A,R);
1814 #undef B0
1815 #undef B1
1816 #undef B2
1817 }
1818 else{
1819 #define B0r V2->Val[MAX_DIM* k +0]
1820 #define B1r V2->Val[MAX_DIM* k +1]
1821 #define B2r V2->Val[MAX_DIM* k +2]
1822 #define B0i V2->Val[MAX_DIM*(k+1)+0]
1823 #define B1i V2->Val[MAX_DIM*(k+1)+1]
1824 #define B2i V2->Val[MAX_DIM*(k+1)+2]
1825 #define AFFECT(i) \
1826 A.Val[MAX_DIM* k +i] = t0*B0r+t1*B1r+t2*B2r; \
1827 A.Val[MAX_DIM*(k+1)+i] = t0*B0i+t1*B1i+t2*B2i
1828 for(k=0 ; k<Current.NbrHar ; k+=2){
1829 t0=A0*A0; t1=A3*A3; t2=A6*A6; AFFECT(0);
1830 t0=A1*A0; t1=A3*A4; t2=A6*A7; AFFECT(1);
1831 t0=A2*A0; t1=A3*A5; t2=A6*A8; AFFECT(2);
1832 t0=A1*A0; t1=A3*A4; t2=A6*A7; AFFECT(3);
1833 t0=A1*A1; t1=A4*A4; t2=A7*A7; AFFECT(4);
1834 t0=A2*A1; t1=A4*A5; t2=A7*A8; AFFECT(5);
1835 t0=A2*A0; t1=A3*A5; t2=A6*A8; AFFECT(6);
1836 t0=A2*A1; t1=A4*A5; t2=A7*A8; AFFECT(7);
1837 t0=A2*A2; t1=A5*A5; t2=A8*A8; AFFECT(8);
1838 }
1839 A.Type = TENSOR;
1840 Cal_CopyValue(&A,R);
1841 #undef B0r
1842 #undef B1r
1843 #undef B2r
1844 #undef B0i
1845 #undef B1i
1846 #undef B2i
1847 #undef AFFECT
1848 }
1849 break;
1850
1851 #define COMPUTE_A \
1852 A.Val[0]= A0*A0*B0+A0*A3*B3+A0*A6*B6+A3*A0*B1+A3*A3*B4+A3*A6*B7+A6*A0*B2+A6*A3*B5+A6*A6*B8; \
1853 A.Val[1]= A1*A0*B0+A1*A3*B3+A1*A6*B6+A4*A0*B1+A4*A3*B4+A4*A6*B7+A7*A0*B2+A7*A3*B5+A7*A6*B8; \
1854 A.Val[2]= A2*A0*B0+A2*A3*B3+A2*A6*B6+A5*A0*B1+A5*A3*B4+A5*A6*B7+A8*A0*B2+A8*A3*B5+A8*A6*B8; \
1855 A.Val[3]= A1*A0*B0+A0*A4*B3+A0*A7*B6+A3*A1*B1+A4*A3*B4+A3*A7*B7+A6*A1*B2+A6*A4*B5+A7*A6*B8; \
1856 A.Val[4]= A1*A1*B0+A1*A4*B3+A1*A7*B6+A4*A1*B1+A4*A4*B4+A4*A7*B7+A7*A1*B2+A7*A4*B5+A7*A7*B8; \
1857 A.Val[5]= A2*A1*B0+A2*A4*B3+A2*A7*B6+A5*A1*B1+A5*A4*B4+A5*A7*B7+A8*A1*B2+A8*A4*B5+A8*A7*B8; \
1858 A.Val[6]= A2*A0*B0+A0*A5*B3+A0*A8*B6+A3*A2*B1+A5*A3*B4+A3*A8*B7+A6*A2*B2+A6*A5*B5+A8*A6*B8; \
1859 A.Val[7]= A2*A1*B0+A1*A5*B3+A1*A8*B6+A4*A2*B1+A5*A4*B4+A4*A8*B7+A7*A2*B2+A7*A5*B5+A8*A7*B8; \
1860 A.Val[8]= A2*A2*B0+A2*A5*B3+A2*A8*B6+A5*A2*B1+A5*A5*B4+A5*A8*B7+A8*A2*B2+A8*A5*B5+A8*A8*B8
1861
1862 #define COMPLEX_COMPUTE_A \
1863 t0=A0*A0; t1=A0*A3; t2=A0*A6; t3=A3*A0; t4=A3*A3; t5=A3*A6; t6=A6*A0; t7=A6*A3; t8=A6*A6; \
1864 AFFECT(0); \
1865 t0=A1*A0; t1=A1*A3; t2=A1*A6; t3=A4*A0; t4=A4*A3; t5=A4*A6; t6=A7*A0; t7=A7*A3; t8=A7*A6; \
1866 AFFECT(1); \
1867 t0=A2*A0; t1=A2*A3; t2=A2*A6; t3=A5*A0; t4=A5*A3; t5=A5*A6; t6=A8*A0; t7=A8*A3; t8=A8*A6; \
1868 AFFECT(2); \
1869 t0=A1*A0; t1=A0*A4; t2=A0*A7; t3=A3*A1; t4=A4*A3; t5=A3*A7; t6=A6*A1; t7=A6*A4; t8=A7*A6; \
1870 AFFECT(3); \
1871 t0=A1*A1; t1=A1*A4; t2=A1*A7; t3=A4*A1; t4=A4*A4; t5=A4*A7; t6=A7*A1; t7=A7*A4; t8=A7*A7; \
1872 AFFECT(4); \
1873 t0=A2*A1; t1=A2*A4; t2=A2*A7; t3=A5*A1; t4=A5*A4; t5=A5*A7; t6=A8*A1; t7=A8*A4; t8=A8*A7; \
1874 AFFECT(5); \
1875 t0=A2*A0; t1=A0*A5; t2=A0*A8; t3=A3*A2; t4=A5*A3; t5=A3*A8; t6=A6*A2; t7=A6*A5; t8=A8*A6; \
1876 AFFECT(6); \
1877 t0=A2*A1; t1=A1*A5; t2=A1*A8; t3=A4*A2; t4=A5*A4; t5=A4*A8; t6=A7*A2; t7=A7*A5; t8=A8*A7; \
1878 AFFECT(7); \
1879 t0=A2*A2; t1=A2*A5; t2=A2*A8; t3=A5*A2; t4=A5*A5; t5=A5*A8; t6=A8*A2; t7=A8*A5; t8=A8*A8; \
1880 AFFECT(8)
1881
1882 case TENSOR_SYM :
1883 if(Current.NbrHar == 1){
1884 #define B0 V2->Val[0]
1885 #define B1 V2->Val[1]
1886 #define B2 V2->Val[2]
1887 #define B3 V2->Val[1]
1888 #define B4 V2->Val[3]
1889 #define B5 V2->Val[4]
1890 #define B6 V2->Val[2]
1891 #define B7 V2->Val[4]
1892 #define B8 V2->Val[5]
1893 COMPUTE_A;
1894 A.Type = TENSOR;
1895 Cal_CopyValue(&A,R);
1896 #undef B0
1897 #undef B1
1898 #undef B2
1899 #undef B3
1900 #undef B4
1901 #undef B5
1902 #undef B6
1903 #undef B7
1904 #undef B8
1905 }
1906 else{
1907 #define B0r V2->Val[MAX_DIM* k +0]
1908 #define B1r V2->Val[MAX_DIM* k +1]
1909 #define B2r V2->Val[MAX_DIM* k +2]
1910 #define B3r V2->Val[MAX_DIM* k +1]
1911 #define B4r V2->Val[MAX_DIM* k +3]
1912 #define B5r V2->Val[MAX_DIM* k +4]
1913 #define B6r V2->Val[MAX_DIM* k +2]
1914 #define B7r V2->Val[MAX_DIM* k +4]
1915 #define B8r V2->Val[MAX_DIM* k +5]
1916 #define B0i V2->Val[MAX_DIM*(k+1)+0]
1917 #define B1i V2->Val[MAX_DIM*(k+1)+1]
1918 #define B2i V2->Val[MAX_DIM*(k+1)+2]
1919 #define B3i V2->Val[MAX_DIM*(k+1)+1]
1920 #define B4i V2->Val[MAX_DIM*(k+1)+3]
1921 #define B5i V2->Val[MAX_DIM*(k+1)+4]
1922 #define B6i V2->Val[MAX_DIM*(k+1)+2]
1923 #define B7i V2->Val[MAX_DIM*(k+1)+4]
1924 #define B8i V2->Val[MAX_DIM*(k+1)+5]
1925 #define AFFECT(i) \
1926 A.Val[MAX_DIM* k +i] = t0*B0r+t1*B3r+t2*B6r+t3*B1r+t4*B4r+t5*B7r+t6*B2r+t7*B5r+t8*B8r; \
1927 A.Val[MAX_DIM*(k+1)+i] = t0*B0i+t1*B3i+t2*B6i+t3*B1i+t4*B4i+t5*B7i+t6*B2i+t7*B5i+t8*B8i
1928 for(k=0 ; k<Current.NbrHar ; k+=2){
1929 COMPLEX_COMPUTE_A;
1930 }
1931 A.Type = TENSOR;
1932 Cal_CopyValue(&A,R);
1933 #undef B0r
1934 #undef B1r
1935 #undef B2r
1936 #undef B3r
1937 #undef B4r
1938 #undef B5r
1939 #undef B6r
1940 #undef B7r
1941 #undef B8r
1942 #undef B0i
1943 #undef B1i
1944 #undef B2i
1945 #undef B3i
1946 #undef B4i
1947 #undef B5i
1948 #undef B6i
1949 #undef B7i
1950 #undef B8i
1951 #undef AFFECT
1952 }
1953 break;
1954
1955 case TENSOR :
1956 if(Current.NbrHar == 1){
1957 #define B0 V2->Val[0]
1958 #define B1 V2->Val[1]
1959 #define B2 V2->Val[2]
1960 #define B3 V2->Val[3]
1961 #define B4 V2->Val[4]
1962 #define B5 V2->Val[5]
1963 #define B6 V2->Val[6]
1964 #define B7 V2->Val[7]
1965 #define B8 V2->Val[8]
1966 COMPUTE_A;
1967 A.Type = TENSOR;
1968 Cal_CopyValue(&A,R);
1969 #undef B0
1970 #undef B1
1971 #undef B2
1972 #undef B3
1973 #undef B4
1974 #undef B5
1975 #undef B6
1976 #undef B7
1977 #undef B8
1978 }
1979 else{
1980 #define B0r V2->Val[MAX_DIM* k +0]
1981 #define B1r V2->Val[MAX_DIM* k +1]
1982 #define B2r V2->Val[MAX_DIM* k +2]
1983 #define B3r V2->Val[MAX_DIM* k +3]
1984 #define B4r V2->Val[MAX_DIM* k +4]
1985 #define B5r V2->Val[MAX_DIM* k +5]
1986 #define B6r V2->Val[MAX_DIM* k +6]
1987 #define B7r V2->Val[MAX_DIM* k +7]
1988 #define B8r V2->Val[MAX_DIM* k +8]
1989 #define B0i V2->Val[MAX_DIM*(k+1)+0]
1990 #define B1i V2->Val[MAX_DIM*(k+1)+1]
1991 #define B2i V2->Val[MAX_DIM*(k+1)+2]
1992 #define B3i V2->Val[MAX_DIM*(k+1)+3]
1993 #define B4i V2->Val[MAX_DIM*(k+1)+4]
1994 #define B5i V2->Val[MAX_DIM*(k+1)+5]
1995 #define B6i V2->Val[MAX_DIM*(k+1)+6]
1996 #define B7i V2->Val[MAX_DIM*(k+1)+7]
1997 #define B8i V2->Val[MAX_DIM*(k+1)+8]
1998 #define AFFECT(i) \
1999 A.Val[MAX_DIM* k +i] = t0*B0r+t1*B3r+t2*B6r+t3*B1r+t4*B4r+t5*B7r+t6*B2r+t7*B5r+t8*B8r; \
2000 A.Val[MAX_DIM*(k+1)+i] = t0*B0i+t1*B3i+t2*B6i+t3*B1i+t4*B4i+t5*B7i+t6*B2i+t7*B5i+t8*B8i
2001 for(k=0 ; k<Current.NbrHar ; k+=2){
2002 COMPLEX_COMPUTE_A;
2003 }
2004 A.Type = TENSOR;
2005 Cal_CopyValue(&A,R);
2006 #undef B0r
2007 #undef B1r
2008 #undef B2r
2009 #undef B3r
2010 #undef B4r
2011 #undef B5r
2012 #undef B6r
2013 #undef B7r
2014 #undef B8r
2015 #undef B0i
2016 #undef B1i
2017 #undef B2i
2018 #undef B3i
2019 #undef B4i
2020 #undef B5i
2021 #undef B6i
2022 #undef B7i
2023 #undef B8i
2024 #undef AFFECT
2025 }
2026 break;
2027
2028 #undef COMPUTE_A
2029 #undef COMPLEX_COMPUTE_A
2030
2031 default :
2032 Message::Error("Wrong argument type in 'Cal_RotateValue'");
2033 break;
2034 }
2035 }
2036
2037 #undef A0
2038 #undef A1
2039 #undef A2
2040 #undef A3
2041 #undef A4
2042 #undef A5
2043 #undef A6
2044 #undef A7
2045 #undef A8
2046
2047 /* ------------------------------------------------------------------------
2048 R <- Det(V1)
2049 ------------------------------------------------------------------------ */
2050
Cal_DetValue(struct Value * V1,struct Value * R)2051 void Cal_DetValue(struct Value *V1, struct Value *R)
2052 {
2053 int k;
2054 int V1Type;
2055
2056 V1Type = V1->Type;
2057 R->Type = SCALAR;
2058
2059 switch(V1Type){
2060
2061 case TENSOR_DIAG :
2062 if(Current.NbrHar==1){
2063 R->Val[0] = V1->Val[0]*V1->Val[1]*V1->Val[2];
2064 }
2065 else{
2066 for(k=0 ; k<Current.NbrHar ; k+=2){
2067 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+0], &V1->Val[MAX_DIM*k+1], a1);
2068 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+2], a1, a2);
2069 R->Val[MAX_DIM* k ] = a2[0];
2070 R->Val[MAX_DIM*(k+1)] = a2[MAX_DIM];
2071 }
2072 }
2073 break;
2074
2075 case TENSOR_SYM :
2076 if(Current.NbrHar==1){
2077 R->Val[0] = V1->Val[0]*(V1->Val[3]*V1->Val[5]-V1->Val[4]*V1->Val[4])
2078 - V1->Val[1]*(V1->Val[1]*V1->Val[5]-V1->Val[2]*V1->Val[4])
2079 + V1->Val[2]*(V1->Val[1]*V1->Val[4]-V1->Val[2]*V1->Val[3]);
2080 }
2081 else{
2082 for(k=0 ; k<Current.NbrHar ; k+=2){
2083 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+3], &V1->Val[MAX_DIM*k+5], a1);
2084 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+4], &V1->Val[MAX_DIM*k+4], a2);
2085 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
2086 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+0], b1, c1);
2087
2088 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+1], &V1->Val[MAX_DIM*k+5], a1);
2089 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+2], &V1->Val[MAX_DIM*k+4], a2);
2090 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
2091 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+1], b1, c2);
2092
2093 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+1], &V1->Val[MAX_DIM*k+4], a1);
2094 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+2], &V1->Val[MAX_DIM*k+3], a2);
2095 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
2096 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+2], b1, c3);
2097
2098 R->Val[MAX_DIM* k ] = c1[0] -c2[0] +c3[0];
2099 R->Val[MAX_DIM*(k+1)] = c1[MAX_DIM]-c2[MAX_DIM]+c3[MAX_DIM];
2100 }
2101 }
2102 break;
2103
2104 case TENSOR :
2105 if(Current.NbrHar==1){
2106 R->Val[0] = V1->Val[0]*(V1->Val[4]*V1->Val[8]-V1->Val[7]*V1->Val[5])
2107 - V1->Val[1]*(V1->Val[3]*V1->Val[8]-V1->Val[6]*V1->Val[5])
2108 + V1->Val[2]*(V1->Val[3]*V1->Val[7]-V1->Val[6]*V1->Val[4]);
2109 }
2110 else{
2111 for(k=0 ; k<Current.NbrHar ; k+=2){
2112 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+4], &V1->Val[MAX_DIM*k+8], a1);
2113 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+7], &V1->Val[MAX_DIM*k+5], a2);
2114 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
2115 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+0], b1, c1);
2116
2117 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+3], &V1->Val[MAX_DIM*k+8], a1);
2118 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+6], &V1->Val[MAX_DIM*k+5], a2);
2119 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
2120 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+1], b1, c2);
2121
2122 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+3], &V1->Val[MAX_DIM*k+7], a1);
2123 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+6], &V1->Val[MAX_DIM*k+4], a2);
2124 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM] ;
2125 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+2], b1, c3);
2126
2127 R->Val[MAX_DIM* k ] = c1[0] -c2[0] +c3[0];
2128 R->Val[MAX_DIM*(k+1)] = c1[MAX_DIM]-c2[MAX_DIM]+c3[MAX_DIM];
2129 }
2130 }
2131 break;
2132
2133 default:
2134 Message::Error("Wrong argument type in 'Cal_DetValue'");
2135 break;
2136 }
2137
2138 }
2139
2140 /* ------------------------------------------------------------------------
2141 R <- 1/V1
2142 ------------------------------------------------------------------------ */
2143
Cal_InvertValue(struct Value * V1,struct Value * R)2144 void Cal_InvertValue(struct Value *V1, struct Value *R)
2145 {
2146 int k;
2147 struct Value Det,A;
2148
2149 switch(V1->Type){
2150
2151 case SCALAR :
2152 R->Type = SCALAR;
2153
2154 if(Current.NbrHar==1){
2155 if(!V1->Val[0]){
2156 Message::Error("Division by zero in 'Cal_InvertValue'");
2157 return;
2158 }
2159 R->Val[0] = 1./V1->Val[0];
2160 }
2161 else{
2162 for(k=0 ; k<Current.NbrHar ; k+=2){
2163 Cal_ComplexInvert(&V1->Val[MAX_DIM*k], &R->Val[MAX_DIM*k]);
2164 }
2165 }
2166 break;
2167
2168 case TENSOR_DIAG :
2169 R->Type = TENSOR_DIAG;
2170
2171 if(Current.NbrHar==1){
2172 if(V1->Val[0] && V1->Val[1] && V1->Val[2]){
2173 R->Val[0] = 1./V1->Val[0];
2174 R->Val[1] = 1./V1->Val[1];
2175 R->Val[2] = 1./V1->Val[2];
2176 }
2177 else{
2178 Message::Error("Null determinant in 'Cal_InvertValue'");
2179 return;
2180 }
2181 }
2182 else{
2183 for(k=0 ; k<Current.NbrHar ; k+=2){
2184 Cal_ComplexInvert(&V1->Val[MAX_DIM*k ], &R->Val[MAX_DIM*k ]);
2185 Cal_ComplexInvert(&V1->Val[MAX_DIM*k+1], &R->Val[MAX_DIM*k+1]);
2186 Cal_ComplexInvert(&V1->Val[MAX_DIM*k+2], &R->Val[MAX_DIM*k+2]);
2187 }
2188 }
2189 break;
2190
2191 case TENSOR_SYM :
2192 Cal_DetValue(V1,&Det);
2193
2194 if(!Det.Val[0]){
2195 Message::Error("Null determinant in 'Cal_InvertValue'");
2196 return;
2197 }
2198
2199 if(Current.NbrHar==1){
2200 A.Val[0] = (V1->Val[3]*V1->Val[5]-V1->Val[4]*V1->Val[4])/Det.Val[0];
2201 A.Val[1] = -(V1->Val[1]*V1->Val[5]-V1->Val[4]*V1->Val[2])/Det.Val[0];
2202 A.Val[2] = (V1->Val[1]*V1->Val[4]-V1->Val[3]*V1->Val[2])/Det.Val[0];
2203 A.Val[3] = -(V1->Val[1]*V1->Val[5]-V1->Val[2]*V1->Val[4])/Det.Val[0];
2204 A.Val[4] = (V1->Val[0]*V1->Val[5]-V1->Val[2]*V1->Val[2])/Det.Val[0];
2205 A.Val[5] = -(V1->Val[0]*V1->Val[4]-V1->Val[2]*V1->Val[1])/Det.Val[0];
2206 A.Val[6] = (V1->Val[1]*V1->Val[4]-V1->Val[2]*V1->Val[3])/Det.Val[0];
2207 A.Val[7] = -(V1->Val[0]*V1->Val[4]-V1->Val[1]*V1->Val[2])/Det.Val[0];
2208 A.Val[8] = (V1->Val[0]*V1->Val[3]-V1->Val[1]*V1->Val[1])/Det.Val[0];
2209 A.Type = TENSOR;
2210 Cal_CopyValue(&A,R);
2211 }
2212 else{
2213 #define PRODSUBDIV(a,b,c,d,e) \
2214 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+a], &V1->Val[MAX_DIM*k+b], a1); \
2215 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+c], &V1->Val[MAX_DIM*k+d], a2); \
2216 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM]; \
2217 Cal_ComplexDivision(b1, &Det.Val[MAX_DIM*k], &A.Val[e])
2218
2219 #define ASSIGN1(i) \
2220 R->Val[MAX_DIM* k +i] = A.Val[MAX_DIM* k +i]; \
2221 R->Val[MAX_DIM*(k+1)+i] = A.Val[MAX_DIM*(k+1)+i]
2222
2223 #define ASSIGN2(i) \
2224 R->Val[MAX_DIM* k +i] = -A.Val[MAX_DIM* k +i]; \
2225 R->Val[MAX_DIM*(k+1)+i] = -A.Val[MAX_DIM*(k+1)+i]
2226 for(k=0 ; k<Current.NbrHar ; k+=2){
2227 PRODSUBDIV(3,5,4,4,0); PRODSUBDIV(1,5,4,2,1); PRODSUBDIV(1,4,3,2,2);
2228 PRODSUBDIV(1,5,2,4,3); PRODSUBDIV(0,5,2,2,4); PRODSUBDIV(0,4,2,1,5);
2229 PRODSUBDIV(1,4,2,3,6); PRODSUBDIV(0,4,1,2,7); PRODSUBDIV(0,3,1,1,8);
2230 ASSIGN1(0); ASSIGN2(1); ASSIGN1(2);
2231 ASSIGN2(3); ASSIGN1(4); ASSIGN2(5);
2232 ASSIGN1(6); ASSIGN2(7); ASSIGN1(8);
2233 }
2234 R->Type = TENSOR;
2235 #undef PRODSUBDIV
2236 #undef ASSIGN1
2237 #undef ASSIGN2
2238 }
2239 break;
2240
2241 case TENSOR :
2242 Cal_DetValue(V1,&Det);
2243
2244 if(!Det.Val[0]){
2245 Message::Error("Null determinant in 'Cal_InvertValue'");
2246 return;
2247 }
2248
2249 if(Current.NbrHar==1){
2250 A.Val[0] = (V1->Val[4]*V1->Val[8]-V1->Val[5]*V1->Val[7])/Det.Val[0];
2251 A.Val[1] = -(V1->Val[1]*V1->Val[8]-V1->Val[7]*V1->Val[2])/Det.Val[0];
2252 A.Val[2] = (V1->Val[1]*V1->Val[5]-V1->Val[4]*V1->Val[2])/Det.Val[0];
2253 A.Val[3] = -(V1->Val[3]*V1->Val[8]-V1->Val[6]*V1->Val[5])/Det.Val[0];
2254 A.Val[4] = (V1->Val[0]*V1->Val[8]-V1->Val[2]*V1->Val[6])/Det.Val[0];
2255 A.Val[5] = -(V1->Val[0]*V1->Val[5]-V1->Val[2]*V1->Val[3])/Det.Val[0];
2256 A.Val[6] = (V1->Val[3]*V1->Val[7]-V1->Val[6]*V1->Val[4])/Det.Val[0];
2257 A.Val[7] = -(V1->Val[0]*V1->Val[7]-V1->Val[1]*V1->Val[6])/Det.Val[0];
2258 A.Val[8] = (V1->Val[0]*V1->Val[4]-V1->Val[1]*V1->Val[3])/Det.Val[0];
2259 A.Type = TENSOR;
2260 Cal_CopyValue(&A,R);
2261 }
2262 else{
2263 #define PRODSUBDIV(a,b,c,d,e) \
2264 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+a], &V1->Val[MAX_DIM*k+b], a1); \
2265 Cal_ComplexProduct(&V1->Val[MAX_DIM*k+c], &V1->Val[MAX_DIM*k+d], a2); \
2266 b1[0] = a1[0] - a2[0] ; b1[MAX_DIM] = a1[MAX_DIM] - a2[MAX_DIM]; \
2267 Cal_ComplexDivision(b1, &Det.Val[MAX_DIM*k], &A.Val[e])
2268
2269 #define ASSIGN1(i) \
2270 R->Val[MAX_DIM* k +i] = A.Val[MAX_DIM* k +i]; \
2271 R->Val[MAX_DIM*(k+1)+i] = A.Val[MAX_DIM*(k+1)+i]
2272
2273 #define ASSIGN2(i) \
2274 R->Val[MAX_DIM* k +i] = -A.Val[MAX_DIM* k +i]; \
2275 R->Val[MAX_DIM*(k+1)+i] = -A.Val[MAX_DIM*(k+1)+i]
2276 for(k=0 ; k<Current.NbrHar ; k+=2){
2277 PRODSUBDIV(4,8,5,7,0); PRODSUBDIV(1,8,7,2,1); PRODSUBDIV(1,5,4,2,2);
2278 PRODSUBDIV(3,8,6,5,3); PRODSUBDIV(0,8,2,6,4); PRODSUBDIV(0,5,2,3,5);
2279 PRODSUBDIV(3,7,6,4,6); PRODSUBDIV(0,7,1,6,7); PRODSUBDIV(0,4,1,3,8);
2280 ASSIGN1(0); ASSIGN2(1); ASSIGN1(2);
2281 ASSIGN2(3); ASSIGN1(4); ASSIGN2(5);
2282 ASSIGN1(6); ASSIGN2(7); ASSIGN1(8);
2283 }
2284 R->Type = TENSOR;
2285 #undef PRODSUBDIV
2286 #undef ASSIGN1
2287 #undef ASSIGN2
2288 }
2289 break;
2290
2291 default :
2292 Message::Error("Wrong type of argument in 'Cal_InvertValue'");
2293 break;
2294 }
2295
2296 }
2297
2298 /* ------------------------------------------------------- */
2299 /* --> P r i n t _ V a l u e */
2300 /* ------------------------------------------------------- */
2301
Print_Value_ToString(struct Value * A)2302 std::string Print_Value_ToString(struct Value *A)
2303 {
2304 int i, j, k, index = 0;
2305 std::ostringstream sstream;
2306
2307 sstream.precision(16);
2308
2309 switch(A->Type){
2310
2311 case SCALAR :
2312 if(Current.NbrHar>1) sstream << "(";
2313 for (k = 0 ; k < Current.NbrHar ; k++) {
2314 if(k) sstream << ",";
2315 sstream << A->Val[MAX_DIM*k];
2316 }
2317 if(Current.NbrHar>1) sstream << ")";
2318 break;
2319
2320 case VECTOR :
2321 sstream << "[";
2322 for (i = 0 ; i < 3 ; i++) {
2323 if(i) sstream << " ";
2324 if(Current.NbrHar>1) sstream << "(";
2325 for (k = 0 ; k < Current.NbrHar ; k++) {
2326 if(k) sstream << ",";
2327 sstream << A->Val[MAX_DIM*k+i];
2328 }
2329 if(Current.NbrHar>1) sstream << ")";
2330 }
2331 sstream << "]";
2332 break;
2333
2334 case TENSOR_DIAG :
2335 case TENSOR_SYM :
2336 case TENSOR :
2337 sstream << "[[";
2338 for (i = 0 ; i < 3 ; i++) {
2339 if(i) sstream << "][";
2340 for (j = 0 ; j < 3 ; j++) {
2341 if(j) sstream << " ";
2342 if(Current.NbrHar>1) sstream << "(";
2343 switch(A->Type){
2344 case TENSOR_DIAG : index = TENSOR_DIAG_MAP[3*i+j]; break;
2345 case TENSOR_SYM : index = TENSOR_SYM_MAP[3*i+j]; break;
2346 case TENSOR : index = 3*i+j; break;
2347 }
2348 for (k = 0 ; k < Current.NbrHar ; k++) {
2349 if(k) sstream << ",";
2350 if(index<0)
2351 sstream << "0";
2352 else
2353 sstream << A->Val[MAX_DIM*k+index];
2354 }
2355 if(Current.NbrHar>1) sstream << ")";
2356 }
2357 }
2358 sstream << "]]";
2359 break;
2360
2361 default :
2362 Message::Error("Unknown type of argument in function 'Printf'");
2363 break;
2364 }
2365
2366 std::string ret(sstream.str());
2367 return ret;
2368 }
2369
Print_Value(struct Value * A,FILE * fp)2370 void Print_Value(struct Value *A, FILE *fp)
2371 {
2372 std::string s = Print_Value_ToString(A);
2373 if(fp && fp != stdout && fp != stderr)
2374 fprintf(fp, "%s\n", s.c_str());
2375 else
2376 Message::Direct("%s", s.c_str());
2377 }
2378
2379 /* ------------------------------------------------------------------------
2380 Complete harmonics
2381 ------------------------------------------------------------------------ */
2382
Cal_SetHarmonicValue(struct Value * R)2383 void Cal_SetHarmonicValue(struct Value *R)
2384 {
2385 int k ;
2386
2387 switch(R->Type){
2388
2389 case SCALAR :
2390 R->Val[MAX_DIM] = 0. ;
2391 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
2392 R->Val[MAX_DIM*k ] = R->Val[0] ;
2393 R->Val[MAX_DIM*(k+1)] = 0. ;
2394 }
2395 break;
2396
2397 case VECTOR :
2398 R->Val[MAX_DIM] = R->Val[MAX_DIM+1] = R->Val[MAX_DIM+2] = 0. ;
2399 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
2400 R->Val[MAX_DIM*k ] = R->Val[0] ;
2401 R->Val[MAX_DIM*k+1] = R->Val[1] ;
2402 R->Val[MAX_DIM*k+2] = R->Val[2] ;
2403 R->Val[MAX_DIM*(k+1) ] = 0. ;
2404 R->Val[MAX_DIM*(k+1)+1] = 0. ;
2405 R->Val[MAX_DIM*(k+1)+2] = 0. ;
2406 }
2407 break;
2408
2409 default :
2410 Message::Error("Unknown type of argument in function 'Cal_SetHarmonicValue'");
2411 }
2412
2413 }
2414
2415 /* ------------------------------------------------------------------------
2416 Set superfluous harmonics to zero (in case of CASTing)
2417 ------------------------------------------------------------------------ */
2418
Cal_SetZeroHarmonicValue(struct Value * R,int Save_NbrHar)2419 void Cal_SetZeroHarmonicValue(struct Value *R, int Save_NbrHar)
2420 {
2421 int k ;
2422
2423 switch(R->Type) {
2424 case SCALAR :
2425 for (k = Current.NbrHar ; k < Save_NbrHar ; k++) {
2426 R->Val[k*MAX_DIM ] = 0. ;
2427 }
2428 break ;
2429 case VECTOR :
2430 case TENSOR_DIAG :
2431 for (k = Current.NbrHar ; k < Save_NbrHar ; k++) {
2432 R->Val[k*MAX_DIM ] = 0. ;
2433 R->Val[k*MAX_DIM+1] = 0. ;
2434 R->Val[k*MAX_DIM+2] = 0. ;
2435 }
2436 break ;
2437 case TENSOR_SYM :
2438 for (k = Current.NbrHar ; k < Save_NbrHar ; k++) {
2439 R->Val[k*MAX_DIM ] = 0. ;
2440 R->Val[k*MAX_DIM+1] = 0. ;
2441 R->Val[k*MAX_DIM+2] = 0. ;
2442 R->Val[k*MAX_DIM+3] = 0. ;
2443 R->Val[k*MAX_DIM+4] = 0. ;
2444 R->Val[k*MAX_DIM+5] = 0. ;
2445 }
2446 break ;
2447 case TENSOR :
2448 for (k = Current.NbrHar ; k < Save_NbrHar ; k++) {
2449 R->Val[k*MAX_DIM ] = 0. ;
2450 R->Val[k*MAX_DIM+1] = 0. ;
2451 R->Val[k*MAX_DIM+2] = 0. ;
2452 R->Val[k*MAX_DIM+3] = 0. ;
2453 R->Val[k*MAX_DIM+4] = 0. ;
2454 R->Val[k*MAX_DIM+5] = 0. ;
2455 R->Val[k*MAX_DIM+6] = 0. ;
2456 R->Val[k*MAX_DIM+7] = 0. ;
2457 R->Val[k*MAX_DIM+8] = 0. ;
2458 }
2459 break ;
2460 default :
2461 Message::Error("Unknown type of argument in function 'Cal_SetZeroHarmonicValue'");
2462 }
2463
2464 }
2465
2466 /* ------------------------------------------------------- */
2467 /* --> S h o w _ V a l u e */
2468 /* ------------------------------------------------------- */
2469
2470 #define W(i,j) A->Val[MAX_DIM*(i)+j]
2471
NonZeroHar(int NbrComp,double Val[])2472 int NonZeroHar(int NbrComp, double Val[])
2473 {
2474 int iComp, nz, nh;
2475
2476 nh=Current.NbrHar-1;
2477 while( nh >= 0 ){
2478 nz=0;
2479 for (iComp=0 ; iComp<NbrComp ; iComp++)
2480 if(Val[MAX_DIM*nh+iComp])nz++;
2481 if(nz)break;
2482 nh--;
2483 }
2484 return nh+1;
2485 }
2486
Show_Value(struct Value * A)2487 void Show_Value(struct Value *A)
2488 {
2489 int k, nzh;
2490
2491 switch(A->Type){
2492
2493 case SCALAR :
2494 if((nzh=NonZeroHar(1,A->Val)) == 0){
2495 fprintf(stderr, "zero scalar \n") ;
2496 }
2497 else if(nzh == 1){
2498 fprintf(stderr, "real scalar %e \n", W(0,0) ) ;
2499 }
2500 else if (nzh == 2){
2501 fprintf(stderr, "complex scalar %e +j %e \n", W(0,0), W(1,0) ) ;
2502 }
2503 else {
2504 fprintf(stderr, "multi-freq scalar ");
2505 for (k = 0 ; k < Current.NbrHar ; k += 2)
2506 fprintf(stderr, " Freq %d : %e + j %e",k/2+1, W(k,0), W(k+1,0) ) ;
2507 fprintf(stderr, "\n");
2508 }
2509 break;
2510
2511 case VECTOR :
2512 if((nzh=NonZeroHar(3,A->Val)) == 0){
2513 fprintf(stderr, "zero vector \n") ;
2514 }
2515 else if (nzh == 1){
2516 fprintf(stderr, "real vector x %e y %e z %e \n", W(0,0), W(0,1), W(0,2));
2517 }
2518 else if (nzh == 2){
2519 fprintf(stderr, "complex vector x %e +j %e y %e +j %e z %e +j %e \n",
2520 W(0,0), W(1,0), W(0,1), W(1,1), W(0,2), W(1,2) );
2521 }
2522 else{
2523 fprintf(stderr, "multi-freq vector ");
2524 for (k = 0 ; k < Current.NbrHar ; k += 2)
2525 fprintf(stderr, " Freq %d : x %e +j %e y %e +j %e z %e +j %e", k/2+1,
2526 W(k,0), W(k+1,0), W(k,1), W(k+1,1), W(k,2), W(k+1,2) );
2527 fprintf(stderr, "\n");
2528 }
2529 break;
2530
2531 case TENSOR :
2532 if((nzh=NonZeroHar(9,A->Val)) == 0){
2533 fprintf(stderr, "zero tensor \n") ;
2534 }
2535 else if (nzh == 1){
2536 fprintf(stderr, "real tensor "
2537 " xx %e xy %e xz %e "
2538 " yx %e yy %e yz %e "
2539 " zx %e zy %e zz %e \n",
2540 W(0,0), W(0,1), W(0,2), W(0,3), W(0,4), W(0,5), W(0,6), W(0,7), W(0,8));
2541 }
2542 else if (nzh == 2){
2543 fprintf(stderr, "complex tensor "
2544 " xx %e +j %e xy %e +j %e xz %e +j %e "
2545 " yx %e +j %e yy %e +j %e yz %e +j %e "
2546 " zx %e +j %e zy %e +j %e zz %e +j %e\n",
2547 W(0,0), W(1,0), W(0,1), W(1,1), W(0,2), W(1,2),
2548 W(0,3), W(1,3), W(0,4), W(1,4), W(0,5), W(1,5),
2549 W(0,6), W(1,6), W(0,7), W(1,7), W(0,8), W(1,8));
2550 }
2551 else {
2552 fprintf(stderr, "multi-freq tensor ");
2553 for (k = 0 ; k < Current.NbrHar ; k += 2)
2554 fprintf(stderr, " Freq %d : "
2555 " xx %e +j %e xy %e +j %e xz %e +j %e "
2556 " yx %e +j %e yy %e +j %e yz %e +j %e "
2557 " zx %e +j %e zy %e +j %e zz %e +j %e", k/2+1,
2558 W(k,0), W(k+1,0), W(k,1), W(k+1,1), W(k,2), W(k+1,2),
2559 W(k,3), W(k+1,3), W(k,4), W(k+1,4), W(k,5), W(k+1,5),
2560 W(k,6), W(k+1,6), W(k,7), W(k+1,7), W(k,8), W(k+1,8));
2561 fprintf(stderr, "\n");
2562 }
2563 break;
2564
2565 case TENSOR_SYM :
2566 if((nzh=NonZeroHar(6,A->Val)) == 0){
2567 fprintf(stderr, "zero sym tensor \n") ;
2568 }
2569 else if (nzh == 1){
2570 fprintf(stderr, "real sym tensor "
2571 " xx %e xy %e xz %e "
2572 " yy %e yz %e zz %e \n",
2573 W(0,0), W(0,1), W(0,2), W(0,3), W(0,4), W(0,5));
2574 }
2575 else if (nzh == 2){
2576 fprintf(stderr, "complex sym tensor "
2577 " xx %e +j %e xy %e +j %e xz %e +j %e "
2578 " yy %e +j %e yz %e +j %e zz %e +j %e\n",
2579 W(0,0), W(1,0), W(0,1), W(1,1), W(0,2), W(1,2),
2580 W(0,3), W(1,3), W(0,4), W(1,4), W(0,5), W(1,5));
2581 }
2582 else {
2583 fprintf(stderr, "multi-freq sym tensor ");
2584 for (k = 0 ; k < Current.NbrHar ; k += 2)
2585 fprintf(stderr, " Freq %d : "
2586 " xx %e +j %e xy %e +j %e xz %e +j %e "
2587 " yy %e +j %e yz %e +j %e zz %e +j %e", k/2+1,
2588 W(k,0), W(k+1,0), W(k,1), W(k+1,1), W(k,2), W(k+1,2),
2589 W(k,3), W(k+1,3), W(k,4), W(k+1,4), W(k,5), W(k+1,5));
2590 fprintf(stderr, "\n");
2591 }
2592 break;
2593
2594 case TENSOR_DIAG :
2595 if((nzh=NonZeroHar(3,A->Val)) == 0){
2596 fprintf(stderr, "zero sym tensor \n") ;
2597 }
2598 else if (nzh == 1){
2599 fprintf(stderr, "real diag tensor xx %e yy %e zz %e \n",
2600 W(0,0), W(0,1), W(0,2));
2601 }
2602 else if (nzh == 2){
2603 fprintf(stderr, "complex diag tensor xx %e +j %e yy %e +j %e zz %e +j %e",
2604 W(0,0), W(1,0), W(0,1), W(1,1), W(0,2), W(1,2));
2605 }
2606 else {
2607 fprintf(stderr, "multi-freq diag tensor ");
2608 for (k = 0 ; k < Current.NbrHar ; k += 2)
2609 fprintf(stderr, " Freq %d : xx %e +j %e yy %e +j %e zz %e +j %e", k/2+1,
2610 W(k,0), W(k+1,0), W(k,1), W(k+1,1), W(k,2), W(k+1,2));
2611 fprintf(stderr, "\n");
2612 }
2613 break;
2614
2615 default :
2616 Message::Error("Unknown value type in Show_Value");
2617 }
2618 }
2619
Export_Value(struct Value * A,std::vector<double> & out,List_T * harmonics,bool append)2620 void Export_Value(struct Value *A, std::vector<double> &out, List_T *harmonics,
2621 bool append)
2622 {
2623 std::vector<int> har;
2624 if(harmonics){
2625 for(int i = 0; i < List_Nbr(harmonics); i++)
2626 har.push_back((int)*(double*)List_Pointer(harmonics, i));
2627 }
2628 else{
2629 for(int i = 0; i < Current.NbrHar; i++)
2630 har.push_back(i);
2631 }
2632
2633 if(!append) out.clear();
2634
2635 switch(A->Type){
2636
2637 case SCALAR :
2638 for(unsigned int k = 0; k < har.size(); k++)
2639 out.push_back(W(har[k], 0));
2640 break;
2641
2642 case VECTOR :
2643 case TENSOR_DIAG :
2644 for(unsigned int k = 0; k < har.size(); k++)
2645 for(int i = 0; i < 3; i++)
2646 out.push_back(W(har[k], i));
2647 break;
2648
2649 case TENSOR :
2650 for(unsigned int k = 0; k < har.size(); k++)
2651 for(int i = 0; i < 9; i++)
2652 out.push_back(W(har[k], i));
2653 break;
2654
2655 case TENSOR_SYM :
2656 for(unsigned int k = 0; k < har.size(); k++)
2657 for(int i = 0; i < 6; i++)
2658 out.push_back(W(har[k], i));
2659 break;
2660
2661 default :
2662 Message::Error("Unknown value type in Export_Value");
2663 }
2664 }
2665
2666 #undef W
2667