7.2
general 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-2022 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 /* Symmetric tensor component name */
60 
61 typedef enum {
62 
63  XX,
64  YY,
65  ZZ,
66  XY,
67  YZ,
68  XZ
69 
70 } cs_math_sym_tensor_component_t;
71 
72 /*============================================================================
73  * Global variables
74  *============================================================================*/
75 
76 /* Numerical constants */
77 
78 extern const cs_real_t cs_math_zero_threshold;
79 extern const cs_real_t cs_math_1ov3;
80 extern const cs_real_t cs_math_2ov3;
81 extern const cs_real_t cs_math_4ov3;
82 extern const cs_real_t cs_math_5ov3;
83 extern const cs_real_t cs_math_1ov6;
84 extern const cs_real_t cs_math_1ov12;
85 extern const cs_real_t cs_math_1ov24;
86 extern const cs_real_t cs_math_epzero;
87 extern const cs_real_t cs_math_infinite_r;
88 extern const cs_real_t cs_math_big_r;
89 extern const cs_real_t cs_math_pi;
90 
91 /* Identity matrix in dimension 3 */
92 static const cs_real_33_t cs_math_33_identity = {{1., 0., 0.,},
93  {0., 1., 0.},
94  {0., 0., 1.}};
95 static const cs_real_6_t cs_math_sym_33_identity = {1., 1., 1., 0. ,0., 0.};
96 
97 /*=============================================================================
98  * Inline static function prototypes
99  *============================================================================*/
100 
101 /*----------------------------------------------------------------------------*/
110 /*----------------------------------------------------------------------------*/
111 
112 static inline int
113 cs_math_binom(int n,
114  int k)
115 {
116  int ret = 1;
117  assert(n >= k);
118 
119  const int n_iter = (k > n-k) ? n-k : k;
120  for (int j = 1; j <= n_iter; j++, n--) {
121  if (n % j == 0)
122  ret *= n/j;
123  else if (ret % j == 0)
124  ret = ret/j*n;
125  else
126  ret = (ret*n)/j;
127  }
128 
129  return ret;
130 }
131 
132 /*----------------------------------------------------------------------------*/
140 /*----------------------------------------------------------------------------*/
141 
142 static inline cs_real_t
143 cs_math_fabs(cs_real_t x)
144 {
145  cs_real_t ret = (x < 0) ? -x : x;
146 
147  return ret;
148 }
149 
150 /*----------------------------------------------------------------------------*/
158 /*----------------------------------------------------------------------------*/
159 
160 static inline cs_real_t
161 cs_math_fmin(cs_real_t x,
162  cs_real_t y)
163 {
164  cs_real_t ret = (x < y) ? x : y;
165 
166  return ret;
167 }
168 
169 /*----------------------------------------------------------------------------*/
177 /*----------------------------------------------------------------------------*/
178 
179 static inline cs_real_t
180 cs_math_fmax(cs_real_t x,
181  cs_real_t y)
182 {
183  cs_real_t ret = (x < y) ? y : x;
184 
185  return ret;
186 }
187 
188 /*----------------------------------------------------------------------------*/
196 /*----------------------------------------------------------------------------*/
197 
198 static inline cs_real_t
199 cs_math_sq(cs_real_t x)
200 {
201  return x*x;
202 }
203 
204 /*----------------------------------------------------------------------------*/
212 /*----------------------------------------------------------------------------*/
213 
214 static inline cs_real_t
215 cs_math_pow2(cs_real_t x)
216 {
217  return x*x;
218 }
219 
220 /*----------------------------------------------------------------------------*/
228 /*----------------------------------------------------------------------------*/
229 
230 static inline cs_real_t
231 cs_math_pow3(cs_real_t x)
232 {
233  return x*x*x;
234 }
235 
236 /*----------------------------------------------------------------------------*/
244 /*----------------------------------------------------------------------------*/
245 
246 static inline cs_real_t
247 cs_math_pow4(cs_real_t x)
248 {
249  return (x*x)*(x*x);
250 }
251 
252 /*----------------------------------------------------------------------------*/
262 /*----------------------------------------------------------------------------*/
263 
264 static inline cs_real_t
265 cs_math_3_distance(const cs_real_t xa[3],
266  const cs_real_t xb[3])
267 {
268  cs_real_t v[3];
269 
270  v[0] = xb[0] - xa[0];
271  v[1] = xb[1] - xa[1];
272  v[2] = xb[2] - xa[2];
273 
274  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
275 }
276 
277 /*----------------------------------------------------------------------------*/
287 /*----------------------------------------------------------------------------*/
288 
289 static inline cs_real_t
290 cs_math_3_distance_dot_product(const cs_real_t xa[3],
291  const cs_real_t xb[3],
292  const cs_real_t xc[3])
293 {
294  return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
295 }
296 
297 /*----------------------------------------------------------------------------*/
307 /*----------------------------------------------------------------------------*/
308 
309 static inline cs_real_t
310 cs_math_3_square_distance(const cs_real_t xa[3],
311  const cs_real_t xb[3])
312 {
313  cs_real_t v[3] = {xb[0] - xa[0],
314  xb[1] - xa[1],
315  xb[2] - xa[2]};
316 
317  return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
318 }
319 
320 /*----------------------------------------------------------------------------*/
329 /*----------------------------------------------------------------------------*/
330 
331 static inline cs_real_t
332 cs_math_3_dot_product(const cs_real_t u[3],
333  const cs_real_t v[3])
334 {
335  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
336 
337  return uv;
338 }
339 
340 /*----------------------------------------------------------------------------*/
350 /*----------------------------------------------------------------------------*/
351 
352 static inline cs_real_t
353 cs_math_3_33_3_dot_product(const cs_real_t n1[3],
354  const cs_real_t t[3][3],
355  const cs_real_t n2[3])
356 {
357  cs_real_t n_t_n
358  = ( n1[0]*t[0][0]*n2[0] + n1[1]*t[1][0]*n2[0] + n1[2]*t[2][0]*n2[0]
359  + n1[0]*t[0][1]*n2[1] + n1[1]*t[1][1]*n2[1] + n1[2]*t[2][1]*n2[1]
360  + n1[0]*t[0][2]*n2[2] + n1[1]*t[1][2]*n2[2] + n1[2]*t[2][2]*n2[2]);
361  return n_t_n;
362 }
363 
364 /*----------------------------------------------------------------------------*/
378 /*----------------------------------------------------------------------------*/
379 
380 static inline cs_real_t
381 cs_math_3_sym_33_3_dot_product(const cs_real_t n1[3],
382  const cs_real_t t[6],
383  const cs_real_t n2[3])
384 {
385  cs_real_t n_t_n
386  = ( n1[0]*t[0]*n2[0] + n1[1]*t[3]*n2[0] + n1[2]*t[5]*n2[0]
387  + n1[0]*t[3]*n2[1] + n1[1]*t[1]*n2[1] + n1[2]*t[4]*n2[1]
388  + n1[0]*t[5]*n2[2] + n1[1]*t[4]*n2[2] + n1[2]*t[2]*n2[2]);
389  return n_t_n;
390 }
391 
392 /*----------------------------------------------------------------------------*/
400 /*----------------------------------------------------------------------------*/
401 
402 static inline cs_real_t
403 cs_math_3_norm(const cs_real_t v[3])
404 {
405  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
406 }
407 
408 /*----------------------------------------------------------------------------*/
416 /*----------------------------------------------------------------------------*/
417 
418 static inline cs_real_t
419 cs_math_3_square_norm(const cs_real_t v[3])
420 {
421  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
422 
423  return v2;
424 }
425 
426 /*----------------------------------------------------------------------------*/
435 /*----------------------------------------------------------------------------*/
436 
437 static inline void
438 cs_math_3_normalise(const cs_real_t vin[3],
439  cs_real_t vout[3])
440 {
441  cs_real_t norm = cs_math_3_norm(vin);
442 
443  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
444 
445  vout[0] = inv_norm * vin[0];
446  vout[1] = inv_norm * vin[1];
447  vout[2] = inv_norm * vin[2];
448 }
449 
450 /*----------------------------------------------------------------------------*/
459 /*----------------------------------------------------------------------------*/
460 
461 static inline void
462 cs_math_3_normalize(const cs_real_t vin[3],
463  cs_real_t vout[3])
464 {
465  cs_real_t norm = cs_math_3_norm(vin);
466 
467  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
468 
469  vout[0] = inv_norm * vin[0];
470  vout[1] = inv_norm * vin[1];
471  vout[2] = inv_norm * vin[2];
472 }
473 
474 /*----------------------------------------------------------------------------*/
483 /*----------------------------------------------------------------------------*/
484 
485 static inline void
486 cs_math_3_orthogonal_projection(const cs_real_t n[3],
487  const cs_real_t v[3],
488  cs_real_t vout[restrict 3])
489 {
490  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
491  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
492  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
493 }
494 
495 /*----------------------------------------------------------------------------*/
504 /*----------------------------------------------------------------------------*/
505 
506 static inline void
507 cs_math_3_normal_scaling(const cs_real_t n[3],
508  cs_real_t factor,
509  cs_real_t v[3])
510 {
511  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
512  for (int i = 0; i < 3; i++)
513  v[i] += v_dot_n * n[i];
514 }
515 
516 /*----------------------------------------------------------------------------*/
526 /*----------------------------------------------------------------------------*/
527 
528 static inline void
529 cs_math_33_normal_scaling_add(const cs_real_t n[3],
530  cs_real_t factor,
531  cs_real_t t[3][3])
532 {
533  cs_real_t n_t_n = (factor -1.) *
534  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
535  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
536  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
537  for (int i = 0; i < 3; i++) {
538  for (int j = 0; j < 3; j++)
539  t[i][j] += n_t_n * n[i] * n[j];
540  }
541 }
542 /*----------------------------------------------------------------------------*/
551 /*----------------------------------------------------------------------------*/
552 
553 static inline void
554 cs_math_33_3_product(const cs_real_t m[3][3],
555  const cs_real_t v[3],
556  cs_real_t mv[restrict 3])
557 {
558  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
559  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
560  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
561 }
562 
563 /*----------------------------------------------------------------------------*/
572 /*----------------------------------------------------------------------------*/
573 
574 static inline void
575 cs_math_33_3_product_add(const cs_real_t m[3][3],
576  const cs_real_t v[3],
577  cs_real_t mv[restrict 3])
578 {
579  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
580  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
581  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
582 }
583 
584 /*----------------------------------------------------------------------------*/
593 /*----------------------------------------------------------------------------*/
594 
595 static inline void
596 cs_math_33t_3_product(const cs_real_t m[3][3],
597  const cs_real_t v[3],
598  cs_real_t mv[restrict 3])
599 {
600  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
601  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
602  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
603 }
604 
605 /*----------------------------------------------------------------------------*/
615 /*----------------------------------------------------------------------------*/
616 
617 static inline void
618 cs_math_sym_33_3_product(const cs_real_t m[6],
619  const cs_real_t v[3],
620  cs_real_t mv[restrict 3])
621 {
622  mv[0] = m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
623  mv[1] = m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
624  mv[2] = m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
625 }
626 
627 /*----------------------------------------------------------------------------*/
637 /*----------------------------------------------------------------------------*/
638 
639 static inline void
640 cs_math_sym_33_3_product_add(const cs_real_t m[6],
641  const cs_real_t v[3],
642  cs_real_t mv[restrict 3])
643 {
644  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
645  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
646  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
647 }
648 
649 /*----------------------------------------------------------------------------*/
657 /*----------------------------------------------------------------------------*/
658 
659 static inline cs_real_t
660 cs_math_6_trace(const cs_real_t t[6])
661 {
662  return (t[0] + t[1] + t[2]);
663 }
664 
665 /*----------------------------------------------------------------------------*/
674 /*----------------------------------------------------------------------------*/
675 
676 static inline void
677 cs_math_66_6_product(const cs_real_t m[6][6],
678  const cs_real_t v[6],
679  cs_real_t mv[restrict 6])
680 {
681  for (int i = 0; i < 6; i++) {
682  for (int j = 0; j < 6; j++)
683  mv[i] = m[i][j] * v[j];
684  }
685 }
686 
687 /*----------------------------------------------------------------------------*/
696 /*----------------------------------------------------------------------------*/
697 
698 static inline void
699 cs_math_66_6_product_add(const cs_real_t m[6][6],
700  const cs_real_t v[6],
701  cs_real_t mv[restrict 6])
702 {
703  for (int i = 0; i < 6; i++) {
704  for (int j = 0; j < 6; j++)
705  mv[i] += m[i][j] * v[j];
706  }
707 }
708 
709 /*----------------------------------------------------------------------------*/
717 /*----------------------------------------------------------------------------*/
718 
719 static inline cs_real_t
720 cs_math_33_determinant(const cs_real_t m[3][3])
721 {
722  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
723  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
724  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
725 
726  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
727 }
728 
729 /*----------------------------------------------------------------------------*/
737 /*----------------------------------------------------------------------------*/
738 
739 static inline cs_real_t
740 cs_math_sym_33_determinant(const cs_real_6_t m)
741 {
742  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
743  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
744  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
745 
746  return m[0]*com0 + m[3]*com1 + m[5]*com2;
747 }
748 
749 /*----------------------------------------------------------------------------*/
757 /*----------------------------------------------------------------------------*/
758 
759 #if defined(__INTEL_COMPILER)
760 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
761 #endif
762 
763 static inline void
764 cs_math_3_cross_product(const cs_real_t u[3],
765  const cs_real_t v[3],
766  cs_real_t uv[restrict 3])
767 {
768  uv[0] = u[1]*v[2] - u[2]*v[1];
769  uv[1] = u[2]*v[0] - u[0]*v[2];
770  uv[2] = u[0]*v[1] - u[1]*v[0];
771 }
772 
773 /*----------------------------------------------------------------------------*/
783 /*----------------------------------------------------------------------------*/
784 
785 #if defined(__INTEL_COMPILER)
786 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
787 #endif
788 
789 static inline cs_real_t
790 cs_math_3_triple_product(const cs_real_t u[3],
791  const cs_real_t v[3],
792  const cs_real_t w[3])
793 {
794  return (u[1]*v[2] - u[2]*v[1]) * w[0]
795  + (u[2]*v[0] - u[0]*v[2]) * w[1]
796  + (u[0]*v[1] - u[1]*v[0]) * w[2];
797 }
798 
799 /*----------------------------------------------------------------------------*/
806 /*----------------------------------------------------------------------------*/
807 
808 static inline void
809 cs_math_33_inv_cramer(const cs_real_t in[3][3],
810  cs_real_t out[3][3])
811 {
812  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
813  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
814  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
815 
816  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
817  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
818  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
819 
820  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
821  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
822  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
823 
824  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
825  const double invdet = 1./det;
826 
827  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
828  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
829  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
830 }
831 
832 /*----------------------------------------------------------------------------*/
838 /*----------------------------------------------------------------------------*/
839 
840 static inline void
841 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
842 {
843  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
844  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
845  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
846  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
847  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
848  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
849  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
850  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
851  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
852 
853  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
854 
855  a[0][0] = a00 * det_inv;
856  a[0][1] = a01 * det_inv;
857  a[0][2] = a02 * det_inv;
858  a[1][0] = a10 * det_inv;
859  a[1][1] = a11 * det_inv;
860  a[1][2] = a12 * det_inv;
861  a[2][0] = a20 * det_inv;
862  a[2][1] = a21 * det_inv;
863  a[2][2] = a22 * det_inv;
864 }
865 
866 /*----------------------------------------------------------------------------*/
873 /*----------------------------------------------------------------------------*/
874 
875 static inline void
876 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
877 {
878  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
879  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
880  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
881  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
882  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
883  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
884 
885  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
886 
887  a[0][0] = a00 * det_inv;
888  a[0][1] = a01 * det_inv;
889  a[0][2] = a02 * det_inv;
890  a[1][0] = a01 * det_inv;
891  a[1][1] = a11 * det_inv;
892  a[1][2] = a12 * det_inv;
893  a[2][0] = a02 * det_inv;
894  a[2][1] = a12 * det_inv;
895  a[2][2] = a22 * det_inv;
896 }
897 
898 /*----------------------------------------------------------------------------*/
908 /*----------------------------------------------------------------------------*/
909 
910 static inline void
911 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
912  cs_real_t sout[restrict 6])
913 {
914  double detinv;
915 
916  sout[0] = s[1]*s[2] - s[4]*s[4];
917  sout[1] = s[0]*s[2] - s[5]*s[5];
918  sout[2] = s[0]*s[1] - s[3]*s[3];
919  sout[3] = s[4]*s[5] - s[3]*s[2];
920  sout[4] = s[3]*s[5] - s[0]*s[4];
921  sout[5] = s[3]*s[4] - s[1]*s[5];
922 
923  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
924 
925  sout[0] *= detinv;
926  sout[1] *= detinv;
927  sout[2] *= detinv;
928  sout[3] *= detinv;
929  sout[4] *= detinv;
930  sout[5] *= detinv;
931 }
932 
933 /*----------------------------------------------------------------------------*/
941 /*----------------------------------------------------------------------------*/
942 
943 static inline void
944 cs_math_33_product(const cs_real_t m1[3][3],
945  const cs_real_t m2[3][3],
946  cs_real_t mout[3][3])
947 {
948  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
949  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
950  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
951 
952  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
953  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
954  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
955 
956  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
957  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
958  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
959 }
960 
961 /*----------------------------------------------------------------------------*/
970 /*----------------------------------------------------------------------------*/
971 
972 static inline void
973 cs_math_33_transform_r_to_a(const cs_real_t m[3][3],
974  const cs_real_t q[3][3],
975  cs_real_t mout[3][3])
976 {
977  /* _m = M.Q */
978  cs_real_33_t _m;
979  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[1][0] + m[0][2]*q[2][0];
980  _m[0][1] = m[0][0]*q[0][1] + m[0][1]*q[1][1] + m[0][2]*q[2][1];
981  _m[0][2] = m[0][0]*q[0][2] + m[0][1]*q[1][2] + m[0][2]*q[2][2];
982 
983  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[1][0] + m[1][2]*q[2][0];
984  _m[1][1] = m[1][0]*q[0][1] + m[1][1]*q[1][1] + m[1][2]*q[2][1];
985  _m[1][2] = m[1][0]*q[0][2] + m[1][1]*q[1][2] + m[1][2]*q[2][2];
986 
987  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[1][0] + m[2][2]*q[2][0];
988  _m[2][1] = m[2][0]*q[0][1] + m[2][1]*q[1][1] + m[2][2]*q[2][1];
989  _m[2][2] = m[2][0]*q[0][2] + m[2][1]*q[1][2] + m[2][2]*q[2][2];
990 
991  /* mout = Q^t _m */
992  mout[0][0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
993  mout[0][1] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
994  mout[0][2] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
995 
996  mout[1][0] = q[0][1]*_m[0][0] + q[1][1]*_m[1][0] + q[2][1]*_m[2][0];
997  mout[1][1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
998  mout[1][2] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
999 
1000  mout[2][0] = q[0][2]*_m[0][0] + q[1][2]*_m[1][0] + q[2][2]*_m[2][0];
1001  mout[2][1] = q[0][2]*_m[0][1] + q[1][2]*_m[1][1] + q[2][2]*_m[2][1];
1002  mout[2][2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1003 }
1004 
1005 /*----------------------------------------------------------------------------*/
1014 /*----------------------------------------------------------------------------*/
1015 
1016 static inline void
1017 cs_math_sym_33_transform_r_to_a(const cs_real_t m[6],
1018  const cs_real_t q[3][3],
1019  cs_real_t mout[6])
1020 {
1021  /* _m = M.Q */
1022  cs_real_33_t _m;
1023  _m[0][0] = m[0]*q[0][0] + m[3]*q[1][0] + m[5]*q[2][0];
1024  _m[0][1] = m[0]*q[0][1] + m[3]*q[1][1] + m[5]*q[2][1];
1025  _m[0][2] = m[0]*q[0][2] + m[3]*q[1][2] + m[5]*q[2][2];
1026 
1027  _m[1][0] = m[3]*q[0][0] + m[1]*q[1][0] + m[4]*q[2][0];
1028  _m[1][1] = m[3]*q[0][1] + m[1]*q[1][1] + m[4]*q[2][1];
1029  _m[1][2] = m[3]*q[0][2] + m[1]*q[1][2] + m[4]*q[2][2];
1030 
1031  _m[2][0] = m[5]*q[0][0] + m[4]*q[1][0] + m[2]*q[2][0];
1032  _m[2][1] = m[5]*q[0][1] + m[4]*q[1][1] + m[2]*q[2][1];
1033  _m[2][2] = m[5]*q[0][2] + m[4]*q[1][2] + m[2]*q[2][2];
1034 
1035  /* mout = Q^t _m */
1036  mout[0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1037  mout[1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1038  mout[2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1039 
1040  mout[3] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1041  mout[4] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1042  mout[5] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1043 
1044 }
1045 
1046 /*----------------------------------------------------------------------------*/
1055 /*----------------------------------------------------------------------------*/
1056 
1057 static inline void
1058 cs_math_33_transform_a_to_r(const cs_real_t m[3][3],
1059  const cs_real_t q[3][3],
1060  cs_real_t mout[3][3])
1061 {
1062  /* _m = M.Q^t */
1063  cs_real_33_t _m;
1064  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[0][1] + m[0][2]*q[0][2];
1065  _m[0][1] = m[0][0]*q[1][0] + m[0][1]*q[1][1] + m[0][2]*q[1][2];
1066  _m[0][2] = m[0][0]*q[2][0] + m[0][1]*q[2][1] + m[0][2]*q[2][2];
1067 
1068  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[0][1] + m[1][2]*q[0][2];
1069  _m[1][1] = m[1][0]*q[1][0] + m[1][1]*q[1][1] + m[1][2]*q[1][2];
1070  _m[1][2] = m[1][0]*q[2][0] + m[1][1]*q[2][1] + m[1][2]*q[2][2];
1071 
1072  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[0][1] + m[2][2]*q[0][2];
1073  _m[2][1] = m[2][0]*q[1][0] + m[2][1]*q[1][1] + m[2][2]*q[1][2];
1074  _m[2][2] = m[2][0]*q[2][0] + m[2][1]*q[2][1] + m[2][2]*q[2][2];
1075 
1076  /* mout = Q _m */
1077  mout[0][0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1078  mout[0][1] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1079  mout[0][2] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1080 
1081  mout[1][0] = q[1][0]*_m[0][0] + q[1][1]*_m[1][0] + q[1][2]*_m[2][0];
1082  mout[1][1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1083  mout[1][2] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1084 
1085  mout[2][0] = q[2][0]*_m[0][0] + q[2][1]*_m[1][0] + q[2][2]*_m[2][0];
1086  mout[2][1] = q[2][0]*_m[0][1] + q[2][1]*_m[1][1] + q[2][2]*_m[2][1];
1087  mout[2][2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1088 }
1089 
1090 /*----------------------------------------------------------------------------*/
1099 /*----------------------------------------------------------------------------*/
1100 
1101 static inline void
1102 cs_math_sym_33_transform_a_to_r(const cs_real_t m[6],
1103  const cs_real_t q[3][3],
1104  cs_real_t mout[6])
1105 {
1106  /* _m = M.Q^t */
1107  cs_real_33_t _m;
1108  _m[0][0] = m[0]*q[0][0] + m[3]*q[0][1] + m[5]*q[0][2];
1109  _m[0][1] = m[0]*q[1][0] + m[3]*q[1][1] + m[5]*q[1][2];
1110  _m[0][2] = m[0]*q[2][0] + m[3]*q[2][1] + m[5]*q[2][2];
1111 
1112  _m[1][0] = m[3]*q[0][0] + m[1]*q[0][1] + m[4]*q[0][2];
1113  _m[1][1] = m[3]*q[1][0] + m[1]*q[1][1] + m[4]*q[1][2];
1114  _m[1][2] = m[3]*q[2][0] + m[1]*q[2][1] + m[4]*q[2][2];
1115 
1116  _m[2][0] = m[5]*q[0][0] + m[4]*q[0][1] + m[2]*q[0][2];
1117  _m[2][1] = m[5]*q[1][0] + m[4]*q[1][1] + m[2]*q[1][2];
1118  _m[2][2] = m[5]*q[2][0] + m[4]*q[2][1] + m[2]*q[2][2];
1119 
1120  /* mout = Q _m */
1121  mout[0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1122  mout[1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1123  mout[2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1124 
1125 
1126  mout[3] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1127  mout[4] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1128  mout[5] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1129 }
1130 
1131 /*----------------------------------------------------------------------------*/
1140 /*----------------------------------------------------------------------------*/
1141 
1142 static inline void
1143 cs_math_33_extract_sym_ant(const cs_real_t m[3][3],
1144  cs_real_t m_sym[3][3],
1145  cs_real_t m_ant[3][3])
1146 {
1147  /* sym = 0.5 (m + m_transpose) */
1148  m_sym[0][0] = 0.5 * (m[0][0] + m[0][0]);
1149  m_sym[0][1] = 0.5 * (m[0][1] + m[1][0]);
1150  m_sym[0][2] = 0.5 * (m[0][2] + m[2][0]);
1151  m_sym[1][0] = 0.5 * (m[1][0] + m[0][1]);
1152  m_sym[1][1] = 0.5 * (m[1][1] + m[1][1]);
1153  m_sym[1][2] = 0.5 * (m[1][2] + m[2][1]);
1154  m_sym[2][0] = 0.5 * (m[2][0] + m[0][2]);
1155  m_sym[2][1] = 0.5 * (m[2][1] + m[1][2]);
1156  m_sym[2][2] = 0.5 * (m[2][2] + m[2][2]);
1157 
1158  /* ant = 0.5 (m - m_transpose) */
1159  m_ant[0][0] = 0.5 * (m[0][0] - m[0][0]);
1160  m_ant[0][1] = 0.5 * (m[0][1] - m[1][0]);
1161  m_ant[0][2] = 0.5 * (m[0][2] - m[2][0]);
1162  m_ant[1][0] = 0.5 * (m[1][0] - m[0][1]);
1163  m_ant[1][1] = 0.5 * (m[1][1] - m[1][1]);
1164  m_ant[1][2] = 0.5 * (m[1][2] - m[2][1]);
1165  m_ant[2][0] = 0.5 * (m[2][0] - m[0][2]);
1166  m_ant[2][1] = 0.5 * (m[2][1] - m[1][2]);
1167  m_ant[2][2] = 0.5 * (m[2][2] - m[2][2]);
1168 }
1169 
1170 /*----------------------------------------------------------------------------*/
1178 /*----------------------------------------------------------------------------*/
1179 
1180 static inline void
1181 cs_math_33_product_add(const cs_real_t m1[3][3],
1182  const cs_real_t m2[3][3],
1183  cs_real_t mout[restrict 3][3])
1184 {
1185  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1186  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1187  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1188 
1189  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1190  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1191  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1192 
1193  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1194  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1195  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1196 }
1197 
1198 /*----------------------------------------------------------------------------*/
1212 /*----------------------------------------------------------------------------*/
1213 
1214 static inline void
1215 cs_math_sym_33_product(const cs_real_t s1[6],
1216  const cs_real_t s2[6],
1217  cs_real_t sout[restrict 6])
1218 {
1219  /* S11 */
1220  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1221  /* S22 */
1222  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1223  /* S33 */
1224  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1225  /* S12 = S21 */
1226  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1227  /* S23 = S32 */
1228  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1229  /* S13 = S31 */
1230  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1231 }
1232 
1233 /*----------------------------------------------------------------------------*/
1241 /*----------------------------------------------------------------------------*/
1242 
1243 static inline void
1244 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1245  cs_real_t sout[restrict 6][6])
1246 {
1247  int tens2vect[3][3];
1248  int iindex[6], jindex[6];
1249 
1250  tens2vect[0][0] = 0; tens2vect[0][1] = 3; tens2vect[0][2] = 5;
1251  tens2vect[1][0] = 3; tens2vect[1][1] = 1; tens2vect[1][2] = 4;
1252  tens2vect[2][0] = 5; tens2vect[2][1] = 4; tens2vect[2][2] = 2;
1253 
1254  iindex[0] = 0; iindex[1] = 1; iindex[2] = 2;
1255  iindex[3] = 0; iindex[4] = 1; iindex[5] = 0;
1256 
1257  jindex[0] = 0; jindex[1] = 1; jindex[2] = 2;
1258  jindex[3] = 1; jindex[4] = 2; jindex[5] = 2;
1259 
1260  /* Consider : W = R*s^t + s*R.
1261  * W_ij = Sum_{k<3} [s_jk*r_ik + s_ik*r_jk]
1262  * We look for A such as A*R = W
1263  */
1264  for (int i = 0; i < 6; i++) {
1265  int ii = iindex[i];
1266  int jj = jindex[i];
1267  for (int k = 0; k < 3; k++) {
1268  int ik = tens2vect[k][ii];
1269  int jk = tens2vect[k][jj];
1270 
1271  sout[ik][i] += s[k][jj];
1272 
1273  sout[jk][i] += s[k][ii];
1274  }
1275  }
1276 }
1277 
1278 /*----------------------------------------------------------------------------*/
1290 /*----------------------------------------------------------------------------*/
1291 
1292 static inline void
1293 cs_math_sym_33_double_product(const cs_real_t s1[6],
1294  const cs_real_t s2[6],
1295  const cs_real_t s3[6],
1296  cs_real_t sout[restrict 3][3])
1297 {
1298  cs_real_t _sout[3][3];
1299 
1300  /* S11 */
1301  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1302  /* S22 */
1303  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1304  /* S33 */
1305  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1306  /* S12 */
1307  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1308  /* S21 */
1309  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1310  /* S23 */
1311  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1312  /* S32 */
1313  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1314  /* S13 */
1315  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1316  /* S31 */
1317  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1318 
1319  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1320  /* S22 */
1321  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1322  /* S33 */
1323  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1324  /* S12 */
1325  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1326  /* S21 */
1327  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1328  /* S23 */
1329  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1330  /* S32 */
1331  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1332  /* S13 */
1333  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1334  /* S31 */
1335  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1336 }
1337 
1338 /*----------------------------------------------------------------------------*/
1345 /*----------------------------------------------------------------------------*/
1346 
1347 static inline void
1348 cs_nvec3(const cs_real_3_t v,
1349  cs_nvec3_t *qv)
1350 {
1351  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1352 
1353  qv->meas = magnitude;
1354  if (fabs(magnitude) > cs_math_zero_threshold) {
1355 
1356  const cs_real_t inv = 1/magnitude;
1357  qv->unitv[0] = inv * v[0];
1358  qv->unitv[1] = inv * v[1];
1359  qv->unitv[2] = inv * v[2];
1360 
1361  }
1362  else
1363  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1364 }
1365 
1366 /*=============================================================================
1367  * Public function prototypes
1368  *============================================================================*/
1369 
1370 /*----------------------------------------------------------------------------*/
1374 /*----------------------------------------------------------------------------*/
1375 
1376 void
1377 cs_math_set_machine_epsilon(void);
1378 
1379 /*----------------------------------------------------------------------------*/
1383 /*----------------------------------------------------------------------------*/
1384 
1385 double
1386 cs_math_get_machine_epsilon(void);
1387 
1388 /*----------------------------------------------------------------------------*/
1398 /*----------------------------------------------------------------------------*/
1399 
1400 void
1401 cs_math_3_length_unitv(const cs_real_t xa[3],
1402  const cs_real_t xb[3],
1403  cs_real_t *len,
1404  cs_real_3_t unitv);
1405 
1406 /*----------------------------------------------------------------------------*/
1418 /*----------------------------------------------------------------------------*/
1419 
1420 void
1421 cs_math_sym_33_eigen(const cs_real_t m[6],
1422  cs_real_t eig_vals[3]);
1423 
1424 /*----------------------------------------------------------------------------*/
1437 /*----------------------------------------------------------------------------*/
1438 
1439 void
1440 cs_math_33_eigen(const cs_real_t m[3][3],
1441  cs_real_t *eig_ratio,
1442  cs_real_t *eig_max);
1443 
1444 /*----------------------------------------------------------------------------*/
1455 /*----------------------------------------------------------------------------*/
1456 
1457 double
1458 cs_math_surftri(const cs_real_t xv[3],
1459  const cs_real_t xe[3],
1460  const cs_real_t xf[3]);
1461 
1462 /*----------------------------------------------------------------------------*/
1474 /*----------------------------------------------------------------------------*/
1475 
1476 double
1477 cs_math_voltet(const cs_real_t xv[3],
1478  const cs_real_t xe[3],
1479  const cs_real_t xf[3],
1480  const cs_real_t xc[3]);
1481 
1482 /*----------------------------------------------------------------------------*/
1495 /*----------------------------------------------------------------------------*/
1496 
1497 void
1498 cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1499  const cs_real_t tol_err,
1500  cs_real_t eig_val[restrict 3],
1501  cs_real_t eig_vec[restrict 3][3]);
1502 
1503 /*----------------------------------------------------------------------------*/
1513 /*----------------------------------------------------------------------------*/
1514 
1515 void
1516 cs_math_fact_lu(cs_lnum_t n_blocks,
1517  const int b_size,
1518  const cs_real_t *a,
1519  cs_real_t *a_lu);
1520 
1521 /*----------------------------------------------------------------------------*/
1531 /*----------------------------------------------------------------------------*/
1532 
1533 void
1534 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1535  const int n,
1536  cs_real_t x[],
1537  const cs_real_t b[]);
1538 
1539 /*----------------------------------------------------------------------------*/
1540 
1541 END_C_DECLS
1542 
1543 #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:911
k
Definition: cs_field_pointer.h:70
numvar::ik
integer, save ik
Definition: numvar.f90:75
cs_math_fabs
static cs_real_t cs_math_fabs(cs_real_t x)
Compute the absolute value of a real value.
Definition: cs_math.h:143
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:290
restrict
#define restrict
Definition: cs_defs.h:142
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:332
cs_math_sym_33_identity
static const cs_real_6_t cs_math_sym_33_identity
Definition: cs_math.h:95
cs_math_6_trace
static cs_real_t cs_math_6_trace(const cs_real_t t[6])
Compute the trace of a symmetric tensor.
Definition: cs_math.h:660
cs_math_fmin
static cs_real_t cs_math_fmin(cs_real_t x, cs_real_t y)
Compute the min value of two real values.
Definition: cs_math.h:161
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:337
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:640
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_t mv[restrict 3])
Compute the product of the transpose of a matrix of 3x3 real values by a vector of 3 real values...
Definition: cs_math.h:596
cs_math_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition: cs_math.h:61
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:677
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:231
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:841
XZ
Definition: cs_math.h:68
cs_math_pi
const cs_real_t cs_math_pi
cs_math_3_normalise
static void cs_math_3_normalise(const cs_real_t vin[3], cs_real_t vout[3])
Normalize a vector of 3 real values.
Definition: cs_math.h:438
cs_defs.h
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_t t[3][3])
Add the dot product with a normal vector to the normal,normal component of a tensor: t += factor * n...
Definition: cs_math.h:529
cs_math_3_triple_product
static cs_real_t cs_math_3_triple_product(const cs_real_t u[3], const cs_real_t v[3], const cs_real_t w[3])
Compute the triple product.
Definition: cs_math.h:790
cs_math_1ov6
const cs_real_t cs_math_1ov6
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:510
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:715
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:740
cs_math_1ov24
const cs_real_t cs_math_1ov24
YY
Definition: cs_math.h:64
cs_math_1ov12
const cs_real_t cs_math_1ov12
cs_math_2ov3
const cs_real_t cs_math_2ov3
cs_math_4ov3
const cs_real_t cs_math_4ov3
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:265
ZZ
Definition: cs_math.h:65
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:351
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:1293
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:322
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.
Definition: cs_math.h:353
cs_nvec3_t
Definition: cs_defs.h:368
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:419
cs_math_3_normalize
static void cs_math_3_normalize(const cs_real_t vin[3], cs_real_t vout[3])
Normalise a vector of 3 real values.
Definition: cs_math.h:462
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:199
cs_math_binom
static int cs_math_binom(int n, int k)
Computes the binomial coefficient of n and k.
Definition: cs_math.h:113
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:471
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
cs_math_sym_33_transform_r_to_a
static void cs_math_sym_33_transform_r_to_a(const cs_real_t m[6], const cs_real_t q[3][3], cs_real_t mout[6])
Compute transformation from relative to absolute reference frame Q^t M Q.
Definition: cs_math.h:1017
cs_math_3_normal_scaling
static void cs_math_3_normal_scaling(const cs_real_t n[3], cs_real_t factor, cs_real_t v[3])
Add the dot product with a normal vector to the normal direction to a vector.
Definition: cs_math.h:507
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.
Definition: cs_math.h:1215
cs_math_1ov3
const cs_real_t cs_math_1ov3
cs_math_33_transform_a_to_r
static void cs_math_33_transform_a_to_r(const cs_real_t m[3][3], const cs_real_t q[3][3], cs_real_t mout[3][3])
Compute transformation from absolute to relative reference frame Q M Q^t.
Definition: cs_math.h:1058
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:699
cs_math_get_machine_epsilon
double cs_math_get_machine_epsilon(void)
Get the value related to the machine precision.
Definition: cs_math.c:245
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_t mout[restrict 3][3])
Add the product of two 3x3 real matrices to a matrix.
Definition: cs_math.h:1181
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:876
cs_math_33_transform_r_to_a
static void cs_math_33_transform_r_to_a(const cs_real_t m[3][3], const cs_real_t q[3][3], cs_real_t mout[3][3])
Compute transformation from relative to absolute reference frame Q^t M Q.
Definition: cs_math.h:973
cs_math_set_machine_epsilon
void cs_math_set_machine_epsilon(void)
Compute the value related to the machine precision.
Definition: cs_math.c:224
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:215
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:501
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:370
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:265
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:720
cs_math_5ov3
const cs_real_t cs_math_5ov3
atimbr::u
double precision, dimension(:,:,:), allocatable u
Definition: atimbr.f90:113
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:403
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:809
cs_math_33_eig_val_vec
void cs_math_33_eig_val_vec(const cs_real_t m_in[3][3], const cs_real_t tol_err, cs_real_t eig_val[restrict 3], cs_real_t eig_vec[restrict 3][3])
Evaluate eigenvalues and eigenvectors of a real symmetric matrix m1[3,3]: m1*m2 = lambda*m2...
Definition: cs_math.c:533
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:335
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:486
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:764
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
cs_math_33_identity
static const cs_real_33_t cs_math_33_identity
Definition: cs_math.h:92
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:371
XY
Definition: cs_math.h:66
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:316
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:247
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:310
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:1348
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:511
YZ
Definition: cs_math.h:67
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_t mout[3][3])
Compute the product of two 3x3 real valued matrices.
Definition: cs_math.h:944
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:342
cs_math_3_sym_33_3_dot_product
static cs_real_t cs_math_3_sym_33_3_dot_product(const cs_real_t n1[3], const cs_real_t t[6], const cs_real_t n2[3])
Compute the dot product of a symmetric tensor t with two vectors, n1 and n2.
Definition: cs_math.h:381
t
Definition: cs_field_pointer.h:92
XX
Definition: cs_math.h:63
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:618
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 (Euclidean norm) between two points xa and xb in a Cartesian coordinate system of ...
Definition: cs_math.c:438
cs_math_33_extract_sym_ant
static void cs_math_33_extract_sym_ant(const cs_real_t m[3][3], cs_real_t m_sym[3][3], cs_real_t m_ant[3][3])
Extract from the given matrix its symmetric and anti-symmetric part.
Definition: cs_math.h:1143
cs_math_fmax
static cs_real_t cs_math_fmax(cs_real_t x, cs_real_t y)
Compute the max value of two real values.
Definition: cs_math.h:180
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_t mv[restrict 3])
Compute the product of a matrix of 3x3 real values by a vector of 3 real values add.
Definition: cs_math.h:575
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
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_t mv[restrict 3])
Compute the product of a matrix of 3x3 real values by a vector of 3 real values.
Definition: cs_math.h:554
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:657
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:1244
cs_math_sym_33_transform_a_to_r
static void cs_math_sym_33_transform_a_to_r(const cs_real_t m[6], const cs_real_t q[3][3], cs_real_t mout[6])
Compute transformation from absolute to relative reference frame Q M Q^t.
Definition: cs_math.h:1102
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:146