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 // Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
7
8 #include <math.h>
9 #include "GetDPConfig.h"
10 #include "ProData.h"
11 #include "F.h"
12 #include "Cal_Value.h"
13 #include "Message.h"
14
15 //
16 // Differential geometry operations on 3-dimensional manifold
17 //
18 // The user has the responsibility that the input values are in fact
19 // coefficients of k-vector fields or differential k-forms on
20 // 3-dimensional manifold in appropiate context
21 //
22
23 // Hodge operator of a k-form:
24 // Arguments:
25 // k-form coefficient vector
26 // metric tensor coefficient matrix
27 //
28 // Parameters:
29 // degree k of the form
30 //
31 // Output:
32 // (3-k)-form coefficient vector
33 //
F_Hodge(F_ARG)34 void F_Hodge(F_ARG)
35 {
36 int k;
37 struct Value detS;
38 struct Value *S;
39
40 k = Fct->Para[0];
41 S = A+1;
42
43 if( (A->Type != SCALAR && A->Type != VECTOR) ||
44 (S->Type != TENSOR_DIAG &&
45 S->Type != TENSOR_SYM &&
46 S->Type != TENSOR))
47 Message::Error("Wrong type of arguments for function 'Hodge'");
48
49 Cal_DetValue(S, &detS);
50 detS.Val[0] = sqrt(fabs(detS.Val[0]));
51
52 switch(k) {
53 case 0:
54 Cal_ProductValue(&detS, A, V);
55 break;
56 case 1:
57 Cal_InvertValue(S, S);
58 Cal_ProductValue(S, A, V);
59 Cal_ProductValue(&detS, V, V);
60 break;
61 case 2:
62 Cal_InvertValue(&detS, &detS);
63 Cal_ProductValue(S, A, V);
64 Cal_ProductValue(&detS, V, V);
65 break;
66 case 3:
67 Cal_InvertValue(&detS, &detS);
68 Cal_ProductValue(&detS, A, V);
69 break;
70 default:
71 Message::Error("Invalid parameter for function 'Hodge'");
72 break;
73 }
74 }
75
76 // Inner product of k-forms:
77 // Arguments:
78 // k-form coefficient vector
79 // k-form coefficient vector
80 // metric tensor coefficient matrix
81 //
82 // Parameters:
83 // degree k of the forms
84 //
85 // Output:
86 // scalar
87 //
F_InnerProduct(F_ARG)88 void F_InnerProduct(F_ARG)
89 {
90 int k;
91 struct Value detS;
92 struct Value *S;
93 struct Value *V1;
94 struct Value *V2;
95
96 k = Fct->Para[0];
97
98 V1 = A;
99 V2 = A+1;
100 S = A+2;
101
102 if( (V1->Type != SCALAR && V1->Type != VECTOR) ||
103 (V2->Type != SCALAR && V2->Type != VECTOR) ||
104 (V2->Type != V1->Type) ||
105 (S->Type != TENSOR_DIAG &&
106 S->Type != TENSOR_SYM &&
107 S->Type != TENSOR))
108 Message::Error("Wrong type of arguments for function 'InnerProduct'");
109
110 switch(k) {
111 case 0:
112 Cal_CopyValue(V2, V);
113 break;
114 case 1:
115 Cal_InvertValue(S, S);
116 Cal_ProductValue(S, V2, V);
117 break;
118 case 2:
119 Cal_InvertValue(&detS, &detS);
120 Cal_ProductValue(S, V2, V);
121 Cal_ProductValue(&detS, V, V);
122 break;
123 case 3:
124 Cal_InvertValue(&detS, &detS);
125 Cal_ProductValue(&detS, V2, V);
126 break;
127 default:
128 Message::Error("Invalid parameter for function 'InnerProduct'");
129 break;
130 }
131 Cal_ProductValue(V1, V, V);
132 }
133
134 // Sharp operator of a k-form:
135 // Arguments:
136 // k-form coefficient vector
137 // metric tensor coefficient matrix
138 //
139 // Parameters:
140 // degree k of the form
141 //
142 // Output:
143 // k-vector coefficient vector
144 //
F_Sharp(F_ARG)145 void F_Sharp(F_ARG)
146 {
147 if( (A->Type != SCALAR && A->Type != VECTOR) ||
148 ((A+1)->Type != TENSOR_DIAG &&
149 (A+1)->Type != TENSOR_SYM &&
150 (A+1)->Type != TENSOR))
151 Message::Error("Wrong type of arguments for function 'Sharp'");
152
153 if( Fct->Para[0] > 3 || Fct->Para[0] < 0 )
154 Message::Error("Invalid parameter for function 'Sharp'");
155
156 Cal_InvertValue(A+1, A+1);
157 F_Flat(Fct, A, V);
158 }
159
160 // Flat operator of a k-vector:
161 // Arguments:
162 // k-vector coefficient vector
163 // metric tensor coefficient matrix
164 //
165 // Parameters:
166 // degree k of the vector
167 //
168 // Output:
169 // k-form coefficient vector
170 //
F_Flat(F_ARG)171 void F_Flat(F_ARG)
172 {
173 int k;
174 struct Value detS;
175 struct Value *S;
176
177 k = Fct->Para[0];
178 S = A+1;
179
180 if( (A->Type != SCALAR && A->Type != VECTOR) ||
181 (S->Type != TENSOR_DIAG &&
182 S->Type != TENSOR_SYM &&
183 S->Type != TENSOR))
184 Message::Error("Wrong type of arguments for function 'Flat'");
185
186 Cal_DetValue(S, &detS);
187 detS.Val[0] = sqrt(fabs(detS.Val[0]));
188
189 switch(k) {
190 case 0:
191 Cal_CopyValue(A, V);
192 break;
193 case 1:
194 Cal_ProductValue(S, A, V);
195 break;
196 case 2:
197 Cal_InvertValue(S, S);
198 Cal_ProductValue(S, A, V);
199 Cal_ProductValue(&detS, V, V);
200 break;
201 case 3:
202 Cal_ProductValue(&detS, A, V);
203 break;
204 default:
205 Message::Error("Invalid parameter for function 'Flat'");
206 break;
207 }
208 }
209
210 // Wedge product of k-forms or k-vectors:
211 // Arguments:
212 // k1-form or k1-vector coefficient vector
213 // k2-form or k2-vector coefficient vector
214 //
215 // Parameters:
216 // degree k1 of the first argument
217 // degree k2 of the second argument
218 //
219 // Output:
220 // (k1+k2)-form or (k1+k2)-vector coefficient vector
221 //
F_WedgeProduct(F_ARG)222 void F_WedgeProduct(F_ARG)
223 {
224 int k1,k2;
225 struct Value *V1;
226 struct Value *V2;
227
228 k1 = Fct->Para[0];
229 k2 = Fct->Para[0];
230
231 V1 = A;
232 V2 = A+1;
233
234 if( (V1->Type != SCALAR && V1->Type != VECTOR) ||
235 (V2->Type != SCALAR && V2->Type != VECTOR) )
236 Message::Error("Wrong type of arguments for function 'WedgeProduct'");
237 if(k1 < 0 || k1 > 3 || k2 < 0 || k2 > 3)
238 Message::Error("Invalid parameter for function 'WedgeProduct'");
239
240 if( k1 == 0 || k2 == 0 || (k1 == 1 && k2 == 2) || (k1 == 2 && k2 == 1) )
241 Cal_ProductValue(V1, V2, V);
242 else if( k1 == 1 && k2 == 1 )
243 Cal_CrossProductValue(V1, V2, V);
244 else if( k1 + k2 > 3 )
245 Cal_ZeroValue(V);
246 else
247 Message::Error("Missing implementation in 'WedgeProduct'");
248 }
249
F_TensorProduct(F_ARG)250 void F_TensorProduct(F_ARG)
251 {
252 Message::Error("'TensorProduct' not implemented");
253 }
254
255 // Interior product of a 1-vector and a k-form:
256 // Arguments:
257 // 1-vector coefficient vector
258 // k-form coefficient vector
259 //
260 // Parameters:
261 // degree k of the k-form
262 //
263 // Output:
264 // (k-1)-form coefficient vector
265 //
F_InteriorProduct(F_ARG)266 void F_InteriorProduct(F_ARG)
267 {
268 int k;
269 struct Value *V1;
270 struct Value *V2;
271
272 k = Fct->Para[0];
273
274 V1 = A;
275 V2 = A+1;
276
277 if( V1->Type != VECTOR ||
278 (V2->Type != SCALAR && V2->Type != VECTOR) )
279 Message::Error("Wrong type of arguments for function 'InteriorProduct'");
280
281 switch(k) {
282 case 1:
283 case 3:
284 Cal_ProductValue(V1, V2, V);
285 break;
286 case 2:
287 Cal_CrossProductValue(V1, V2, V);
288 break;
289 default:
290 Message::Error("Invalid parameter for function 'InteriorProduct'");
291 break;
292 }
293 }
294
295 // Pullback of a k-form:
296 // Arguments:
297 // k-form coefficient vector
298 // Jacobian matrix of the transition map between charts of a manifold
299 //
300 // Parameters:
301 // degree k of the form
302 //
303 // Output:
304 // k-form coefficient vector
305 //
F_PullBack(F_ARG)306 void F_PullBack(F_ARG)
307 {
308 if( (A->Type != SCALAR && A->Type != VECTOR) ||
309 ((A+1)->Type != TENSOR_DIAG &&
310 (A+1)->Type != TENSOR_SYM &&
311 (A+1)->Type != TENSOR))
312 Message::Error("Wrong type of arguments for function 'PullBack'");
313
314 if( Fct->Para[0] < 0 || Fct->Para[0] > 3 )
315 Message::Error("Invalid parameter for function 'PullBack'");
316
317 Cal_TransposeValue(A+1, A+1);
318 F_PushForward(Fct, A, V);
319 }
320
321 // Pullback of a metric tensor:
322 // Arguments:
323 // metric tensor coefficient matrix
324 // Jacobian matrix of the transition map between charts of a manifold
325 //
326 // Parameters:
327 // none
328 //
329 // Output:
330 // metric tensor coefficient matrix
331 //
F_PullBackMetric(F_ARG)332 void F_PullBackMetric(F_ARG)
333 {
334 struct Value *S;
335 struct Value *J;
336
337 J = A+1;
338 S = A;
339
340 if( (S->Type != TENSOR_DIAG &&
341 S->Type != TENSOR_SYM &&
342 S->Type != TENSOR) ||
343 (J->Type != TENSOR_DIAG &&
344 J->Type != TENSOR_SYM &&
345 J->Type != TENSOR))
346 Message::Error("Wrong type of arguments for function 'PullBackMetric'");
347
348 Cal_ProductValue(S, J, V);
349 Cal_TransposeValue(J, J);
350 Cal_ProductValue(J, V, V);
351 }
352
353 // Inverse pullback of a k-form:
354 // Arguments:
355 // k-form coefficient vector
356 // Jacobian matrix of the transition map between charts of a manifold
357 //
358 // Parameters:
359 // degree k of the form
360 //
361 // Output:
362 // k-form coefficient vector
363 //
F_InvPullBack(F_ARG)364 void F_InvPullBack(F_ARG)
365 {
366 if( (A->Type != SCALAR && A->Type != VECTOR) ||
367 ((A+1)->Type != TENSOR_DIAG &&
368 (A+1)->Type != TENSOR_SYM &&
369 (A+1)->Type != TENSOR))
370 Message::Error("Wrong type of arguments for function 'InvPullBack'");
371
372 if( Fct->Para[0] < 0 || Fct->Para[0] > 3 )
373 Message::Error("Invalid parameter for function 'InvPullBack'");
374
375 Cal_InvertValue(A+1, A+1);
376 Cal_TransposeValue(A+1, A+1);
377 F_PushForward(Fct, A, V);
378 }
379
380 // Pushforward of a k-vector:
381 // Arguments:
382 // k-vector coefficient vector
383 // Jacobian matrix of the transition map between charts of a manifold
384 //
385 // Parameters:
386 // degree k of the k-vector
387 //
388 // Output:
389 // k-vector coefficient vector
390 //
F_PushForward(F_ARG)391 void F_PushForward(F_ARG)
392 {
393 int k;
394 struct Value *J;
395 struct Value detJ;
396
397 k = Fct->Para[0];
398 J = A+1;
399
400 if( (A->Type != SCALAR && A->Type != VECTOR) ||
401 (J->Type != TENSOR_DIAG &&
402 J->Type != TENSOR_SYM &&
403 J->Type != TENSOR))
404 Message::Error("Wrong type of arguments for function 'PushForward'");
405
406 switch(k) {
407 case 0:
408 Cal_CopyValue(A, V);
409 break;
410 case 1:
411 Cal_ProductValue(J, A, V);
412 break;
413 case 2:
414 Cal_InvertValue(J, J);
415 Cal_TransposeValue(J, J);
416 Cal_DetValue(J, &detJ);
417 Cal_ProductValue(J, A, V);
418 Cal_ProductValue(&detJ, V, V);
419 break;
420 case 3:
421 Cal_DetValue(J, &detJ);
422 Cal_ProductValue(&detJ, A, V);
423 break;
424 default:
425 Message::Error("Invalid parameter for function 'PushForward'");
426 break;
427 }
428 }
429
430 // Inverse pushforward of a k-vector:
431 // Arguments:
432 // k-vector coefficient vector
433 // Jacobian matrix of the transition map between charts of a manifold
434 //
435 // Parameters:
436 // degree k of the k-vector
437 //
438 // Output:
439 // k-vector coefficient vector
440 //
F_InvPushForward(F_ARG)441 void F_InvPushForward(F_ARG)
442 {
443 if( (A->Type != SCALAR && A->Type != VECTOR) ||
444 ((A+1)->Type != TENSOR_DIAG &&
445 (A+1)->Type != TENSOR_SYM &&
446 (A+1)->Type != TENSOR))
447 Message::Error("Wrong type of arguments for function 'InvPushForward'");
448
449 if( Fct->Para[0] < 0 || Fct->Para[0] > 3 )
450 Message::Error("Invalid parameter for function 'InvPushForward'");
451
452 Cal_InvertValue(A+1, A+1);
453 F_PushForward(Fct, A, V);
454 }
455