programmer's documentation
cs_math.h
Go to the documentation of this file.
1 #ifndef __CS_MATH_H__
2 #define __CS_MATH_H__
3 
4 /*============================================================================
5  * Mathematical base functions.
6  *============================================================================*/
7 
8 /*
9  This file is part of Code_Saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-2018 EDF S.A.
12 
13  This program is free software; you can redistribute it and/or modify it under
14  the terms of the GNU General Public License as published by the Free Software
15  Foundation; either version 2 of the License, or (at your option) any later
16  version.
17 
18  This program is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21  details.
22 
23  You should have received a copy of the GNU General Public License along with
24  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25  Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27 
28 /*----------------------------------------------------------------------------*/
29 
30 #include "cs_defs.h"
31 
32 /*----------------------------------------------------------------------------
33  * Standard C library headers
34  *----------------------------------------------------------------------------*/
35 
36 #include <assert.h>
37 #include <math.h>
38 
39 /*----------------------------------------------------------------------------
40  * Local headers
41  *----------------------------------------------------------------------------*/
42 
43 #if defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
44 #include "bft_error.h"
45 #endif
46 
47 /*----------------------------------------------------------------------------*/
48 
49 BEGIN_C_DECLS
50 
51 /*=============================================================================
52  * Local Macro definitions
53  *============================================================================*/
54 
55 /*============================================================================
56  * Type definition
57  *============================================================================*/
58 
59 /*============================================================================
60  * Global variables
61  *============================================================================*/
62 
63 /* Numerical constants */
64 
65 extern const cs_real_t cs_math_zero_threshold;
66 extern const cs_real_t cs_math_onethird;
67 extern const cs_real_t cs_math_onesix;
68 extern const cs_real_t cs_math_onetwelve;
69 extern const cs_real_t cs_math_one24;
70 extern const cs_real_t cs_math_epzero;
71 extern const cs_real_t cs_math_infinite_r;
72 extern const cs_real_t cs_math_big_r;
73 extern const cs_real_t cs_math_pi;
74 
75 /*=============================================================================
76  * Inline static function prototypes
77  *============================================================================*/
78 
79 /*----------------------------------------------------------------------------*/
88 /*----------------------------------------------------------------------------*/
89 
90 static inline int
91 cs_math_binom(short int n,
92  short int k)
93 {
94  int ret = 1;
95  assert(n >= k);
96 
97  const short int n_iter = (k > n-k) ? n-k : k;
98  for (short int j = 1; j <= n_iter; j++, n--) {
99  if (n % j == 0)
100  ret *= n/j;
101  else if (ret % j == 0)
102  ret = ret/j*n;
103  else
104  ret = (ret*n)/j;
105  }
106 
107  return ret;
108 }
109 
110 /*----------------------------------------------------------------------------*/
118 /*----------------------------------------------------------------------------*/
119 
120 static inline cs_real_t
121 cs_math_sq(cs_real_t x)
122 {
123  return x*x;
124 }
125 
126 /*----------------------------------------------------------------------------*/
134 /*----------------------------------------------------------------------------*/
135 
136 static inline cs_real_t
137 cs_math_pow2(cs_real_t x)
138 {
139  return x*x;
140 }
141 
142 /*----------------------------------------------------------------------------*/
150 /*----------------------------------------------------------------------------*/
151 
152 static inline cs_real_t
153 cs_math_pow3(cs_real_t x)
154 {
155  return x*x*x;
156 }
157 
158 /*----------------------------------------------------------------------------*/
166 /*----------------------------------------------------------------------------*/
167 
168 static inline cs_real_t
169 cs_math_pow4(cs_real_t x)
170 {
171  return (x*x)*(x*x);
172 }
173 
174 /*----------------------------------------------------------------------------*/
184 /*----------------------------------------------------------------------------*/
185 
186 static inline cs_real_t
187 cs_math_3_distance(const cs_real_t xa[3],
188  const cs_real_t xb[3])
189 {
190  cs_real_t v[3];
191 
192  v[0] = xb[0] - xa[0];
193  v[1] = xb[1] - xa[1];
194  v[2] = xb[2] - xa[2];
195 
196  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
197 }
198 
199 /*----------------------------------------------------------------------------*/
209 /*----------------------------------------------------------------------------*/
210 
211 static inline cs_real_t
212 cs_math_3_distance_dot_product(const cs_real_t xa[3],
213  const cs_real_t xb[3],
214  const cs_real_t xc[3])
215 {
216  return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
217 }
218 
219 /*----------------------------------------------------------------------------*/
229 /*----------------------------------------------------------------------------*/
230 
231 static inline cs_real_t
232 cs_math_3_square_distance(const cs_real_t xa[3],
233  const cs_real_t xb[3])
234 {
235  cs_real_t v[3] = {xb[0] - xa[0],
236  xb[1] - xa[1],
237  xb[2] - xa[2]};
238 
239  return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
240 }
241 
242 /*----------------------------------------------------------------------------*/
251 /*----------------------------------------------------------------------------*/
252 
253 static inline cs_real_t
254 cs_math_3_dot_product(const cs_real_t u[3],
255  const cs_real_t v[3])
256 {
257  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
258 
259  return uv;
260 }
261 
262 /*----------------------------------------------------------------------------*/
274 /*----------------------------------------------------------------------------*/
275 
276 static inline cs_real_t
277 cs_math_3_33_3_dot_product(const cs_real_t n1[3],
278  const cs_real_t t[3][3],
279  const cs_real_t n2[3])
280 {
281  cs_real_t n_t_n =
282  ( n1[0] * t[0][0] * n2[0] + n1[1] * t[1][0] * n2[0] + n1[2] * t[2][0] * n2[0]
283  + n1[0] * t[0][1] * n2[1] + n1[1] * t[1][1] * n2[1] + n1[2] * t[2][1] * n2[1]
284  + n1[0] * t[0][2] * n2[2] + n1[1] * t[1][2] * n2[2] + n1[2] * t[2][2] * n2[2]);
285  return n_t_n;
286 }
287 
288 
289 /*----------------------------------------------------------------------------*/
297 /*----------------------------------------------------------------------------*/
298 
299 static inline cs_real_t
300 cs_math_3_norm(const cs_real_t v[3])
301 {
302  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
303 }
304 
305 /*----------------------------------------------------------------------------*/
313 /*----------------------------------------------------------------------------*/
314 
315 static inline cs_real_t
316 cs_math_3_square_norm(const cs_real_t v[3])
317 {
318  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
319 
320  return v2;
321 }
322 
323 /*----------------------------------------------------------------------------*/
330 /*----------------------------------------------------------------------------*/
331 
332 static inline void
333 cs_math_3_normalise(const cs_real_t vin[3],
334  cs_real_t vout[restrict 3])
335 {
336  cs_real_t norm = cs_math_3_norm(vin);
337 
338  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
339 
340  vout[0] = inv_norm * vin[0];
341  vout[1] = inv_norm * vin[1];
342  vout[2] = inv_norm * vin[2];
343 }
344 
345 /*----------------------------------------------------------------------------*/
354 /*----------------------------------------------------------------------------*/
355 
356 static inline void
357 cs_math_3_orthogonal_projection(const cs_real_t n[3],
358  const cs_real_t v[3],
359  cs_real_t vout[restrict 3])
360 {
361  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
362  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
363  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
364 }
365 
366 /*----------------------------------------------------------------------------*/
375 /*----------------------------------------------------------------------------*/
376 
377 static inline void
378 cs_math_3_normal_scaling(const cs_real_t n[3],
379  cs_real_t factor,
380  cs_real_3_t v)
381 {
382  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
383  for (int i = 0; i < 3; i++)
384  v[i] += v_dot_n * n[i];
385 }
386 
387 /*----------------------------------------------------------------------------*/
397 /*----------------------------------------------------------------------------*/
398 
399 static inline void
400 cs_math_33_normal_scaling_add(const cs_real_t n[3],
401  cs_real_t factor,
402  cs_real_33_t t)
403 {
404  cs_real_t n_t_n = (factor -1.) *
405  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
406  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
407  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
408  for (int i = 0; i < 3; i++) {
409  for (int j = 0; j < 3; j++)
410  t[i][j] += n_t_n * n[i] * n[j];
411  }
412 }
413 /*----------------------------------------------------------------------------*/
422 /*----------------------------------------------------------------------------*/
423 
424 static inline void
425 cs_math_33_3_product(const cs_real_t m[3][3],
426  const cs_real_t v[3],
427  cs_real_3_t mv)
428 {
429  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
430  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
431  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
432 }
433 
434 /*----------------------------------------------------------------------------*/
443 /*----------------------------------------------------------------------------*/
444 
445 static inline void
446 cs_math_33_3_product_add(const cs_real_t m[3][3],
447  const cs_real_t v[3],
448  cs_real_3_t mv)
449 {
450  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
451  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
452  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
453 }
454 
455 /*----------------------------------------------------------------------------*/
464 /*----------------------------------------------------------------------------*/
465 
466 static inline void
467 cs_math_33t_3_product(const cs_real_t m[3][3],
468  const cs_real_t v[3],
469  cs_real_3_t mv)
470 {
471  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
472  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
473  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
474 }
475 
476 /*----------------------------------------------------------------------------*/
486 /*----------------------------------------------------------------------------*/
487 
488 static inline void
489 cs_math_sym_33_3_product(const cs_real_t m[6],
490  const cs_real_t v[3],
491  cs_real_t mv[restrict 3])
492 {
493  mv[0] = m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
494  mv[1] = m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
495  mv[2] = m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
496 }
497 
498 /*----------------------------------------------------------------------------*/
508 /*----------------------------------------------------------------------------*/
509 
510 static inline void
511 cs_math_sym_33_3_product_add(const cs_real_t m[6],
512  const cs_real_t v[3],
513  cs_real_t mv[restrict 3])
514 {
515  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
516  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
517  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
518 }
519 
520 /*----------------------------------------------------------------------------*/
529 /*----------------------------------------------------------------------------*/
530 
531 static inline void
532 cs_math_66_6_product(const cs_real_t m[6][6],
533  const cs_real_t v[6],
534  cs_real_t mv[restrict 6])
535 {
536  for (int i = 0; i < 6; i++) {
537  for (int j = 0; j < 6; j++)
538  mv[i] = m[i][j] * v[j];
539  }
540 }
541 
542 /*----------------------------------------------------------------------------*/
551 /*----------------------------------------------------------------------------*/
552 
553 static inline void
554 cs_math_66_6_product_add(const cs_real_t m[6][6],
555  const cs_real_t v[6],
556  cs_real_t mv[restrict 6])
557 {
558  for (int i = 0; i < 6; i++) {
559  for (int j = 0; j < 6; j++)
560  mv[i] += m[i][j] * v[j];
561  }
562 }
563 
564 /*----------------------------------------------------------------------------*/
572 /*----------------------------------------------------------------------------*/
573 
574 static inline cs_real_t
575 cs_math_33_determinant(const cs_real_t m[3][3])
576 {
577  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
578  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
579  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
580 
581  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
582 }
583 
584 /*----------------------------------------------------------------------------*/
592 /*----------------------------------------------------------------------------*/
593 
594 static inline cs_real_t
595 cs_math_sym_33_determinant(const cs_real_6_t m)
596 {
597  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
598  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
599  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
600 
601  return m[0]*com0 + m[3]*com1 + m[5]*com2;
602 }
603 
604 /*----------------------------------------------------------------------------*/
612 /*----------------------------------------------------------------------------*/
613 
614 #if defined(__INTEL_COMPILER)
615 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
616 #endif
617 
618 static inline void
619 cs_math_3_cross_product(const cs_real_t u[3],
620  const cs_real_t v[3],
621  cs_real_t uv[restrict 3])
622 {
623  uv[0] = u[1]*v[2] - u[2]*v[1];
624  uv[1] = u[2]*v[0] - u[0]*v[2];
625  uv[2] = u[0]*v[1] - u[1]*v[0];
626 }
627 
628 /*----------------------------------------------------------------------------*/
635 /*----------------------------------------------------------------------------*/
636 
637 static inline void
638 cs_math_33_inv_cramer(const cs_real_t in[3][3],
639  cs_real_t out[3][3])
640 {
641  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
642  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
643  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
644 
645  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
646  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
647  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
648 
649  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
650  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
651  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
652 
653  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
654  const double invdet = 1/det;
655 
656  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
657  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
658  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
659 }
660 
661 /*----------------------------------------------------------------------------*/
667 /*----------------------------------------------------------------------------*/
668 
669 static inline void
670 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
671 {
672  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
673  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
674  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
675  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
676  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
677  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
678  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
679  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
680  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
681 
682  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
683 
684  a[0][0] = a00 * det_inv;
685  a[0][1] = a01 * det_inv;
686  a[0][2] = a02 * det_inv;
687  a[1][0] = a10 * det_inv;
688  a[1][1] = a11 * det_inv;
689  a[1][2] = a12 * det_inv;
690  a[2][0] = a20 * det_inv;
691  a[2][1] = a21 * det_inv;
692  a[2][2] = a22 * det_inv;
693 }
694 
695 /*----------------------------------------------------------------------------*/
702 /*----------------------------------------------------------------------------*/
703 
704 static inline void
705 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
706 {
707  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
708  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
709  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
710  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
711  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
712  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
713 
714  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
715 
716  a[0][0] = a00 * det_inv;
717  a[0][1] = a01 * det_inv;
718  a[0][2] = a02 * det_inv;
719  a[1][0] = a01 * det_inv;
720  a[1][1] = a11 * det_inv;
721  a[1][2] = a12 * det_inv;
722  a[2][0] = a02 * det_inv;
723  a[2][1] = a12 * det_inv;
724  a[2][2] = a22 * det_inv;
725 }
726 
727 /*----------------------------------------------------------------------------*/
737 /*----------------------------------------------------------------------------*/
738 
739 static inline void
740 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
741  cs_real_t sout[restrict 6])
742 {
743  double detinv;
744 
745  sout[0] = s[1]*s[2] - s[4]*s[4];
746  sout[1] = s[0]*s[2] - s[5]*s[5];
747  sout[2] = s[0]*s[1] - s[3]*s[3];
748  sout[3] = s[4]*s[5] - s[3]*s[2];
749  sout[4] = s[3]*s[5] - s[0]*s[4];
750  sout[5] = s[3]*s[4] - s[1]*s[5];
751 
752  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
753 
754  sout[0] *= detinv;
755  sout[1] *= detinv;
756  sout[2] *= detinv;
757  sout[3] *= detinv;
758  sout[4] *= detinv;
759  sout[5] *= detinv;
760 }
761 
762 /*----------------------------------------------------------------------------*/
771 /*----------------------------------------------------------------------------*/
772 
773 static inline void
774 cs_math_33_product(const cs_real_t m1[3][3],
775  const cs_real_t m2[3][3],
776  cs_real_33_t mout)
777 {
778  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
779  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
780  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
781 
782  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
783  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
784  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
785 
786  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
787  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
788  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
789 }
790 
791 /*----------------------------------------------------------------------------*/
800 /*----------------------------------------------------------------------------*/
801 
802 static inline void
803 cs_math_33_product_add(const cs_real_t m1[3][3],
804  const cs_real_t m2[3][3],
805  cs_real_33_t mout)
806 {
807  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
808  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
809  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
810 
811  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
812  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
813  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
814 
815  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
816  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
817  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
818 }
819 
820 
821 /*----------------------------------------------------------------------------*/
834 /*----------------------------------------------------------------------------*/
835 
836 static inline void
837 cs_math_sym_33_product(const cs_real_t s1[6],
838  const cs_real_t s2[6],
839  cs_real_t sout[restrict 6])
840 {
841  /* S11 */
842  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
843  /* S22 */
844  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
845  /* S33 */
846  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
847  /* S12 = S21 */
848  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
849  /* S23 = S32 */
850  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
851  /* S13 = S31 */
852  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
853 }
854 
855 /*----------------------------------------------------------------------------*/
863 /*----------------------------------------------------------------------------*/
864 
865 static inline void
866 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
867  cs_real_t sout[restrict 6][6])
868 {
869  int tens2vect[3][3];
870  int iindex[6], jindex[6];
871 
872  tens2vect[0][0] = 0; tens2vect[0][1] = 3; tens2vect[0][2] = 5;
873  tens2vect[1][0] = 3; tens2vect[1][1] = 1; tens2vect[1][2] = 4;
874  tens2vect[2][0] = 5; tens2vect[2][1] = 4; tens2vect[2][2] = 2;
875 
876  iindex[0] = 0; iindex[1] = 1; iindex[2] = 2;
877  iindex[3] = 0; iindex[4] = 1; iindex[5] = 0;
878 
879  jindex[0] = 0; jindex[1] = 1; jindex[2] = 2;
880  jindex[3] = 1; jindex[4] = 2; jindex[5] = 2;
881 
882  /* Consider : W = R*s^t + s*R.
883  * W_ij = Sum_{k<3} [s_jk*r_ik + s_ik*r_jk]
884  * We look for A such as A*R = W
885  */
886  for (int i = 0; i < 6; i++) {
887  int ii = iindex[i];
888  int jj = jindex[i];
889  for (int k = 0; k < 3; k++) {
890  int ik = tens2vect[k][ii];
891  int jk = tens2vect[k][jj];
892 
893  sout[ik][i] += s[k][jj];
894 
895  sout[jk][i] += s[k][ii];
896  }
897  }
898 }
899 
900 /*----------------------------------------------------------------------------*/
912 /*----------------------------------------------------------------------------*/
913 
914 static inline void
915 cs_math_sym_33_double_product(const cs_real_t s1[6],
916  const cs_real_t s2[6],
917  const cs_real_t s3[6],
918  cs_real_t sout[restrict 3][3])
919 {
920  cs_real_33_t _sout;
921 
922  /* S11 */
923  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
924  /* S22 */
925  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
926  /* S33 */
927  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
928  /* S12 */
929  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
930  /* S21 */
931  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
932  /* S23 */
933  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
934  /* S32 */
935  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
936  /* S13 */
937  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
938  /* S31 */
939  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
940 
941  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
942  /* S22 */
943  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
944  /* S33 */
945  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
946  /* S12 */
947  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
948  /* S21 */
949  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
950  /* S23 */
951  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
952  /* S32 */
953  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
954  /* S13 */
955  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
956  /* S31 */
957  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
958 }
959 
960 /*----------------------------------------------------------------------------*/
967 /*----------------------------------------------------------------------------*/
968 
969 static inline void
970 cs_nvec3(const cs_real_3_t v,
971  cs_nvec3_t *qv)
972 {
973  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
974 
975  qv->meas = magnitude;
976  if (fabs(magnitude) > cs_math_zero_threshold) {
977 
978  const cs_real_t inv = 1/magnitude;
979  qv->unitv[0] = inv * v[0];
980  qv->unitv[1] = inv * v[1];
981  qv->unitv[2] = inv * v[2];
982 
983  }
984  else
985  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
986 
987 }
988 
989 /*=============================================================================
990  * Public function prototypes
991  *============================================================================*/
992 
993 /*----------------------------------------------------------------------------*/
997 /*----------------------------------------------------------------------------*/
998 
999 void
1000 cs_math_set_machine_epsilon(void);
1001 
1002 /*----------------------------------------------------------------------------*/
1006 /*----------------------------------------------------------------------------*/
1007 
1008 double
1009 cs_math_get_machine_epsilon(void);
1010 
1011 /*----------------------------------------------------------------------------*/
1021 /*----------------------------------------------------------------------------*/
1022 
1023 void
1024 cs_math_3_length_unitv(const cs_real_t xa[3],
1025  const cs_real_t xb[3],
1026  cs_real_t *len,
1027  cs_real_3_t unitv);
1028 
1029 /*----------------------------------------------------------------------------*/
1041 /*----------------------------------------------------------------------------*/
1042 
1043 void
1044 cs_math_sym_33_eigen(const cs_real_t m[6],
1045  cs_real_t eig_vals[3]);
1046 
1047 /*----------------------------------------------------------------------------*/
1060 /*----------------------------------------------------------------------------*/
1061 
1062 void
1063 cs_math_33_eigen(const cs_real_t m[3][3],
1064  cs_real_t *eig_ratio,
1065  cs_real_t *eig_max);
1066 
1067 /*----------------------------------------------------------------------------*/
1078 /*----------------------------------------------------------------------------*/
1079 
1080 double
1081 cs_math_surftri(const cs_real_t xv[3],
1082  const cs_real_t xe[3],
1083  const cs_real_t xf[3]);
1084 
1085 /*----------------------------------------------------------------------------*/
1097 /*----------------------------------------------------------------------------*/
1098 
1099 double
1100 cs_math_voltet(const cs_real_t xv[3],
1101  const cs_real_t xe[3],
1102  const cs_real_t xf[3],
1103  const cs_real_t xc[3]);
1104 
1105 /*----------------------------------------------------------------------------*/
1115 /*----------------------------------------------------------------------------*/
1116 
1117 void
1118 cs_math_fact_lu(cs_lnum_t n_blocks,
1119  const int b_size,
1120  const cs_real_t *a,
1121  cs_real_t *a_lu);
1122 
1123 /*----------------------------------------------------------------------------*/
1133 /*----------------------------------------------------------------------------*/
1134 
1135 void
1136 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1137  const int n,
1138  cs_real_t x[],
1139  const cs_real_t b[]);
1140 
1141 /*----------------------------------------------------------------------------*/
1142 
1143 END_C_DECLS
1144 
1145 #endif /* __CS_MATH_H__ */
cs_math_sym_33_inv_cramer
static void cs_math_sym_33_inv_cramer(const cs_real_t s[6], cs_real_t sout[restrict 6])
Compute the inverse of a symmetric matrix using Cramer&#39;s rule.
Definition: cs_math.h:740
k
Definition: cs_field_pointer.h:70
numvar::ik
integer, save ik
Definition: numvar.f90:75
cs_math_3_distance_dot_product
static cs_real_t cs_math_3_distance_dot_product(const cs_real_t xa[3], const cs_real_t xb[3], const cs_real_t xc[3])
Compute .
Definition: cs_math.h:212
restrict
#define restrict
Definition: cs_defs.h:122
cs_math_3_dot_product
static cs_real_t cs_math_3_dot_product(const cs_real_t u[3], const cs_real_t v[3])
Compute the dot product of two vectors of 3 real values.
Definition: cs_math.h:254
cs_math_onesix
const cs_real_t cs_math_onesix
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:312
len
size_t len
Definition: mei_scanner.c:560
cs_math_3_normalise
static void cs_math_3_normalise(const cs_real_t vin[3], cs_real_t vout[restrict 3])
Normalise a vector of 3 real values.
Definition: cs_math.h:333
cs_math_big_r
const cs_real_t cs_math_big_r
cs_math_sym_33_3_product_add
static void cs_math_sym_33_3_product_add(const cs_real_t m[6], const cs_real_t v[3], cs_real_t mv[restrict 3])
Compute the product of a symmetric matrix of 3x3 real values by a vector of 3 real values and add it ...
Definition: cs_math.h:511
cs_math_66_6_product
static void cs_math_66_6_product(const cs_real_t m[6][6], const cs_real_t v[6], cs_real_t mv[restrict 6])
Compute the product of a matrix of 6x6 real values by a vector of 6 real values.
Definition: cs_math.h:532
cs_math_pow3
static cs_real_t cs_math_pow3(cs_real_t x)
Compute the cube of a real value.
Definition: cs_math.h:153
cs_math_33_inv_cramer_in_place
static void cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
Inverse a 3x3 matrix in place, using Cramer&#39;s rule.
Definition: cs_math.h:670
cs_math_33_3_product
static void cs_math_33_3_product(const cs_real_t m[3][3], const cs_real_t v[3], cs_real_3_t mv)
Compute the product of a matrix of 3x3 real values by a vector of 3 real values.
Definition: cs_math.h:425
cs_math_pi
const cs_real_t cs_math_pi
cs_defs.h
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:462
cs_math_fw_and_bw_lu
void cs_math_fw_and_bw_lu(const cs_real_t a_lu[], const int n, cs_real_t x[], const cs_real_t b[])
Block Jacobi utilities. Compute forward and backward to solve an LU P*P system.
Definition: cs_math.c:537
cs_math_epzero
const cs_real_t cs_math_epzero
cs_math_sym_33_determinant
static cs_real_t cs_math_sym_33_determinant(const cs_real_6_t m)
Compute the determinant of a 3x3 symmetric matrix.
Definition: cs_math.h:595
cs_math_3_distance
static cs_real_t cs_math_3_distance(const cs_real_t xa[3], const cs_real_t xb[3])
Compute the (euclidean) distance between two points xa and xb in a cartesian coordinate system of dim...
Definition: cs_math.h:187
cs_math_33_normal_scaling_add
static void cs_math_33_normal_scaling_add(const cs_real_t n[3], cs_real_t factor, cs_real_33_t t)
Add the dot product with a normal vector to the normal,normal component of a tensor: t += factor * n...
Definition: cs_math.h:400
cs_math_33_eigen
void cs_math_33_eigen(const cs_real_t m[3][3], cs_real_t *eig_ratio, cs_real_t *eig_max)
Compute max/min eigenvalues ratio and max. eigenvalue of a 3x3 symmetric matrix with non-symmetric st...
Definition: cs_math.c:300
cs_math_sym_33_double_product
static void cs_math_sym_33_double_product(const cs_real_t s1[6], const cs_real_t s2[6], const cs_real_t s3[6], cs_real_t sout[restrict 3][3])
Compute the product of three symmetric matrices.
Definition: cs_math.h:915
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:297
cs_math_3_33_3_dot_product
static cs_real_t cs_math_3_33_3_dot_product(const cs_real_t n1[3], const cs_real_t t[3][3], const cs_real_t n2[3])
Compute the dot product of a tensor t with two vectors n1, and n2 n1 t n2.
Definition: cs_math.h:277
cs_nvec3_t
Definition: cs_defs.h:338
u
Definition: cs_field_pointer.h:68
cs_math_one24
const cs_real_t cs_math_one24
cs_math_3_square_norm
static cs_real_t cs_math_3_square_norm(const cs_real_t v[3])
Compute the square norm of a vector of 3 real values.
Definition: cs_math.h:316
cs_math_sq
static cs_real_t cs_math_sq(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:121
cs_math_surftri
double cs_math_surftri(const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3])
Compute the area of the convex_hull generated by 3 points. This corresponds to the computation of the...
Definition: cs_math.c:420
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
cs_math_sym_33_product
static void cs_math_sym_33_product(const cs_real_t s1[6], const cs_real_t s2[6], cs_real_t sout[restrict 6])
Compute the product of two symmetric matrices. Warning: this is valid if and only if s1 and s2 commut...
Definition: cs_math.h:837
cs_math_66_6_product_add
static void cs_math_66_6_product_add(const cs_real_t m[6][6], const cs_real_t v[6], cs_real_t mv[restrict 6])
Compute the product of a matrix of 6x6 real values by a vector of 6 real values and add it to the vec...
Definition: cs_math.h:554
cs_math_get_machine_epsilon
double cs_math_get_machine_epsilon(void)
Get the value related to the machine precision.
Definition: cs_math.c:195
cs_math_onetwelve
const cs_real_t cs_math_onetwelve
cs_math_33_inv_cramer_sym_in_place
static void cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
Inverse a 3x3 symmetric matrix (with non-symmetric storage) in place, using Cramer&#39;s rule...
Definition: cs_math.h:705
cs_math_set_machine_epsilon
void cs_math_set_machine_epsilon(void)
Compute the value related to the machine precision.
Definition: cs_math.c:174
bft_error.h
cs_math_pow2
static cs_real_t cs_math_pow2(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:137
cs_fuel_incl::a
double precision, save a
Definition: cs_fuel_incl.f90:146
cs_math_voltet
double cs_math_voltet(const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3], const cs_real_t xc[3])
Compute the volume of the convex_hull generated by 4 points. This is equivalent to the computation of...
Definition: cs_math.c:450
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:340
cs_math_sym_33_eigen
void cs_math_sym_33_eigen(const cs_real_t m[6], cs_real_t eig_vals[3])
Compute all eigenvalues of a 3x3 symmetric matrix with symmetric storage.
Definition: cs_math.c:215
cs_math_33_determinant
static cs_real_t cs_math_33_determinant(const cs_real_t m[3][3])
Compute the determinant of a 3x3 matrix.
Definition: cs_math.h:575
cs_math_33_product
static void cs_math_33_product(const cs_real_t m1[3][3], const cs_real_t m2[3][3], cs_real_33_t mout)
Compute the product of a matrix of 3x3 real values by a matrix of 3x3 real values.
Definition: cs_math.h:774
cs_math_3_norm
static cs_real_t cs_math_3_norm(const cs_real_t v[3])
Compute the euclidean norm of a vector of dimension 3.
Definition: cs_math.h:300
cs_math_33_inv_cramer
static void cs_math_33_inv_cramer(const cs_real_t in[3][3], cs_real_t out[3][3])
Inverse a 3x3 matrix.
Definition: cs_math.h:638
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:310
cs_math_3_orthogonal_projection
static void cs_math_3_orthogonal_projection(const cs_real_t n[3], const cs_real_t v[3], cs_real_t vout[restrict 3])
Orthogonal projection of a vector with respect to a normalised vector.
Definition: cs_math.h:357
cs_math_3_cross_product
static void cs_math_3_cross_product(const cs_real_t u[3], const cs_real_t v[3], cs_real_t uv[restrict 3])
Compute the cross product of two vectors of 3 real values.
Definition: cs_math.h:619
cs_math_33_3_product_add
static void cs_math_33_3_product_add(const cs_real_t m[3][3], const cs_real_t v[3], cs_real_3_t mv)
Compute the product of a matrix of 3x3 real values by a vector of 3 real values add.
Definition: cs_math.h:446
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:341
cs_math_onethird
const cs_real_t cs_math_onethird
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:293
cs_math_pow4
static cs_real_t cs_math_pow4(cs_real_t x)
Compute the 4-th power of a real value.
Definition: cs_math.h:169
cs_math_binom
static int cs_math_binom(short int n, short int k)
Computes the binomial coefficient of n and k.
Definition: cs_math.h:91
cs_math_3_normal_scaling
static void cs_math_3_normal_scaling(const cs_real_t n[3], cs_real_t factor, cs_real_3_t v)
Add the dot product with a normal vector to the normal direction to a vector.
Definition: cs_math.h:378
cs_math_3_square_distance
static cs_real_t cs_math_3_square_distance(const cs_real_t xa[3], const cs_real_t xb[3])
Compute the squared distance between two points xa and xb in a cartesian coordinate system of dimensi...
Definition: cs_math.h:232
cs_nvec3
static void cs_nvec3(const cs_real_3_t v, cs_nvec3_t *qv)
Define a cs_nvec3_t structure from a cs_real_3_t.
Definition: cs_math.h:970
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:463
cs_math_33t_3_product
static void cs_math_33t_3_product(const cs_real_t m[3][3], const cs_real_t v[3], cs_real_3_t mv)
Compute the product of the transpose of a matrix of 3x3 real values by a vector of 3 real values...
Definition: cs_math.h:467
cs_math_33_product_add
static void cs_math_33_product_add(const cs_real_t m1[3][3], const cs_real_t m2[3][3], cs_real_33_t mout)
Compute the product of a matrix of 3x3 real values by a matrix of 3x3 real values and add...
Definition: cs_math.h:803
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:316
t
Definition: cs_field_pointer.h:96
cs_math_sym_33_3_product
static void cs_math_sym_33_3_product(const cs_real_t m[6], const cs_real_t v[3], cs_real_t mv[restrict 3])
Compute the product of a symmetric matrix of 3x3 real values by a vector of 3 real values...
Definition: cs_math.h:489
cs_math_3_length_unitv
void cs_math_3_length_unitv(const cs_real_t xa[3], const cs_real_t xb[3], cs_real_t *len, cs_real_3_t unitv)
Compute the length (euclidien norm) between two points xa and xb in a cartesian coordinate system of ...
Definition: cs_math.c:387
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
cs_math_fact_lu
void cs_math_fact_lu(cs_lnum_t n_blocks, const int b_size, const cs_real_t *a, cs_real_t *a_lu)
Compute LU factorization of an array of dense matrices of identical size.
Definition: cs_math.c:479
cs_math_reduce_sym_prod_33_to_66
static void cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3], cs_real_t sout[restrict 6][6])
Compute a 6x6 matrix A, equivalent to a 3x3 matrix s, such as: A*R_6 = R*s^t + s*R.
Definition: cs_math.h:866
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:146
Generated on Sat Feb 2 2019 15:43:25 by   doxygen 1.8.13