7.3
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  return ( n1[0] * (t[0]*n2[0] + t[3]*n2[1] + t[5]*n2[2])
386  + n1[1] * (t[3]*n2[0] + t[1]*n2[1] + t[4]*n2[2])
387  + n1[2] * (t[5]*n2[0] + t[4]*n2[1] + t[2]*n2[2]));
388 }
389 
390 /*----------------------------------------------------------------------------*/
398 /*----------------------------------------------------------------------------*/
399 
400 static inline cs_real_t
401 cs_math_3_norm(const cs_real_t v[3])
402 {
403  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
404 }
405 
406 /*----------------------------------------------------------------------------*/
414 /*----------------------------------------------------------------------------*/
415 
416 static inline cs_real_t
417 cs_math_3_square_norm(const cs_real_t v[3])
418 {
419  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
420 
421  return v2;
422 }
423 
424 /*----------------------------------------------------------------------------*/
433 /*----------------------------------------------------------------------------*/
434 
435 static inline void
436 cs_math_3_normalise(const cs_real_t vin[3],
437  cs_real_t vout[3])
438 {
439  cs_real_t norm = cs_math_3_norm(vin);
440 
441  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
442 
443  vout[0] = inv_norm * vin[0];
444  vout[1] = inv_norm * vin[1];
445  vout[2] = inv_norm * vin[2];
446 }
447 
448 /*----------------------------------------------------------------------------*/
457 /*----------------------------------------------------------------------------*/
458 
459 static inline void
460 cs_math_3_normalize(const cs_real_t vin[3],
461  cs_real_t vout[3])
462 {
463  cs_real_t norm = cs_math_3_norm(vin);
464 
465  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
466 
467  vout[0] = inv_norm * vin[0];
468  vout[1] = inv_norm * vin[1];
469  vout[2] = inv_norm * vin[2];
470 }
471 
472 /*----------------------------------------------------------------------------*/
481 /*----------------------------------------------------------------------------*/
482 
483 static inline void
484 cs_math_3_orthogonal_projection(const cs_real_t n[3],
485  const cs_real_t v[3],
486  cs_real_t vout[restrict 3])
487 {
488  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
489  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
490  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
491 }
492 
493 /*----------------------------------------------------------------------------*/
502 /*----------------------------------------------------------------------------*/
503 
504 static inline void
505 cs_math_3_normal_scaling(const cs_real_t n[3],
506  cs_real_t factor,
507  cs_real_t v[3])
508 {
509  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
510  for (int i = 0; i < 3; i++)
511  v[i] += v_dot_n * n[i];
512 }
513 
514 /*----------------------------------------------------------------------------*/
524 /*----------------------------------------------------------------------------*/
525 
526 static inline void
527 cs_math_33_normal_scaling_add(const cs_real_t n[3],
528  cs_real_t factor,
529  cs_real_t t[3][3])
530 {
531  cs_real_t n_t_n = (factor -1.) *
532  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
533  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
534  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
535  for (int i = 0; i < 3; i++) {
536  for (int j = 0; j < 3; j++)
537  t[i][j] += n_t_n * n[i] * n[j];
538  }
539 }
540 /*----------------------------------------------------------------------------*/
549 /*----------------------------------------------------------------------------*/
550 
551 static inline void
552 cs_math_33_3_product(const cs_real_t m[3][3],
553  const cs_real_t v[3],
554  cs_real_t mv[restrict 3])
555 {
556  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
557  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
558  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
559 }
560 
561 /*----------------------------------------------------------------------------*/
570 /*----------------------------------------------------------------------------*/
571 
572 static inline void
573 cs_math_33_3_product_add(const cs_real_t m[3][3],
574  const cs_real_t v[3],
575  cs_real_t mv[restrict 3])
576 {
577  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
578  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
579  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
580 }
581 
582 /*----------------------------------------------------------------------------*/
591 /*----------------------------------------------------------------------------*/
592 
593 static inline void
594 cs_math_33t_3_product(const cs_real_t m[3][3],
595  const cs_real_t v[3],
596  cs_real_t mv[restrict 3])
597 {
598  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
599  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
600  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
601 }
602 
603 /*----------------------------------------------------------------------------*/
613 /*----------------------------------------------------------------------------*/
614 
615 static inline void
616 cs_math_sym_33_3_product(const cs_real_t m[6],
617  const cs_real_t v[3],
618  cs_real_t mv[restrict 3])
619 {
620  mv[0] = m[0]*v[0] + m[3]*v[1] + m[5]*v[2];
621  mv[1] = m[3]*v[0] + m[1]*v[1] + m[4]*v[2];
622  mv[2] = m[5]*v[0] + m[4]*v[1] + m[2]*v[2];
623 }
624 
625 /*----------------------------------------------------------------------------*/
635 /*----------------------------------------------------------------------------*/
636 
637 static inline void
638 cs_math_sym_33_3_product_add(const cs_real_t m[6],
639  const cs_real_t v[3],
640  cs_real_t mv[restrict 3])
641 {
642  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
643  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
644  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
645 }
646 
647 /*----------------------------------------------------------------------------*/
655 /*----------------------------------------------------------------------------*/
656 
657 static inline cs_real_t
658 cs_math_6_trace(const cs_real_t t[6])
659 {
660  return (t[0] + t[1] + t[2]);
661 }
662 
663 /*----------------------------------------------------------------------------*/
672 /*----------------------------------------------------------------------------*/
673 
674 static inline void
675 cs_math_66_6_product(const cs_real_t m[6][6],
676  const cs_real_t v[6],
677  cs_real_t mv[restrict 6])
678 {
679  for (int i = 0; i < 6; i++) {
680  for (int j = 0; j < 6; j++)
681  mv[i] = m[i][j] * v[j];
682  }
683 }
684 
685 /*----------------------------------------------------------------------------*/
694 /*----------------------------------------------------------------------------*/
695 
696 static inline void
697 cs_math_66_6_product_add(const cs_real_t m[6][6],
698  const cs_real_t v[6],
699  cs_real_t mv[restrict 6])
700 {
701  for (int i = 0; i < 6; i++) {
702  for (int j = 0; j < 6; j++)
703  mv[i] += m[i][j] * v[j];
704  }
705 }
706 
707 /*----------------------------------------------------------------------------*/
715 /*----------------------------------------------------------------------------*/
716 
717 static inline cs_real_t
718 cs_math_33_determinant(const cs_real_t m[3][3])
719 {
720  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
721  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
722  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
723 
724  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
725 }
726 
727 /*----------------------------------------------------------------------------*/
735 /*----------------------------------------------------------------------------*/
736 
737 static inline cs_real_t
738 cs_math_sym_33_determinant(const cs_real_6_t m)
739 {
740  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
741  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
742  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
743 
744  return m[0]*com0 + m[3]*com1 + m[5]*com2;
745 }
746 
747 /*----------------------------------------------------------------------------*/
755 /*----------------------------------------------------------------------------*/
756 
757 #if defined(__INTEL_COMPILER)
758 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
759 #endif
760 
761 static inline void
762 cs_math_3_cross_product(const cs_real_t u[3],
763  const cs_real_t v[3],
764  cs_real_t uv[restrict 3])
765 {
766  uv[0] = u[1]*v[2] - u[2]*v[1];
767  uv[1] = u[2]*v[0] - u[0]*v[2];
768  uv[2] = u[0]*v[1] - u[1]*v[0];
769 }
770 
771 /*----------------------------------------------------------------------------*/
781 /*----------------------------------------------------------------------------*/
782 
783 #if defined(__INTEL_COMPILER)
784 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
785 #endif
786 
787 static inline cs_real_t
788 cs_math_3_triple_product(const cs_real_t u[3],
789  const cs_real_t v[3],
790  const cs_real_t w[3])
791 {
792  return (u[1]*v[2] - u[2]*v[1]) * w[0]
793  + (u[2]*v[0] - u[0]*v[2]) * w[1]
794  + (u[0]*v[1] - u[1]*v[0]) * w[2];
795 }
796 
797 /*----------------------------------------------------------------------------*/
807 /*----------------------------------------------------------------------------*/
808 
809 static inline void
810 cs_math_3_orthonormal_basis(const cs_real_t vect[3],
811  cs_real_t axes[3][3])
812 {
813  assert(cs_math_3_norm(vect) > cs_math_zero_threshold);
814 
815  // Compute first axis
816  cs_math_3_normalize(vect, axes[0]);
817 
818  // Compute second axis
819  // First test projection of Ox
820  cs_real_t Ox[3] = {1., 0., 0.};
821  cs_real_t w[3] = {0.};
822 
823  cs_math_3_orthogonal_projection(axes[0], Ox, w);
824 
825  // If Ox projection is null, project Oy
826  if (cs_math_3_norm(w) < cs_math_zero_threshold) {
827  cs_real_t Oy[3] = {0., 1., 0.};
828  cs_math_3_orthogonal_projection(axes[0], Oy, w);
829  }
830 
831  cs_math_3_normalize(w, axes[1]);
832 
833  // Compute third axis using cross product
834  cs_math_3_cross_product(axes[0], axes[1], axes[2]);
835 }
836 
837 /*----------------------------------------------------------------------------*/
844 /*----------------------------------------------------------------------------*/
845 
846 static inline void
847 cs_math_33_inv_cramer(const cs_real_t in[3][3],
848  cs_real_t out[3][3])
849 {
850  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
851  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
852  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
853 
854  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
855  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
856  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
857 
858  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
859  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
860  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
861 
862  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
863  const double invdet = 1./det;
864 
865  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
866  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
867  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
868 }
869 
870 /*----------------------------------------------------------------------------*/
876 /*----------------------------------------------------------------------------*/
877 
878 static inline void
879 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
880 {
881  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
882  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
883  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
884  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
885  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
886  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
887  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
888  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
889  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
890 
891  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
892 
893  a[0][0] = a00 * det_inv;
894  a[0][1] = a01 * det_inv;
895  a[0][2] = a02 * det_inv;
896  a[1][0] = a10 * det_inv;
897  a[1][1] = a11 * det_inv;
898  a[1][2] = a12 * det_inv;
899  a[2][0] = a20 * det_inv;
900  a[2][1] = a21 * det_inv;
901  a[2][2] = a22 * det_inv;
902 }
903 
904 /*----------------------------------------------------------------------------*/
911 /*----------------------------------------------------------------------------*/
912 
913 static inline void
914 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
915 {
916  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
917  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
918  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
919  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
920  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
921  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
922 
923  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
924 
925  a[0][0] = a00 * det_inv;
926  a[0][1] = a01 * det_inv;
927  a[0][2] = a02 * det_inv;
928  a[1][0] = a01 * det_inv;
929  a[1][1] = a11 * det_inv;
930  a[1][2] = a12 * det_inv;
931  a[2][0] = a02 * det_inv;
932  a[2][1] = a12 * det_inv;
933  a[2][2] = a22 * det_inv;
934 }
935 
936 /*----------------------------------------------------------------------------*/
946 /*----------------------------------------------------------------------------*/
947 
948 static inline void
949 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
950  cs_real_t sout[restrict 6])
951 {
952  double detinv;
953 
954  sout[0] = s[1]*s[2] - s[4]*s[4];
955  sout[1] = s[0]*s[2] - s[5]*s[5];
956  sout[2] = s[0]*s[1] - s[3]*s[3];
957  sout[3] = s[4]*s[5] - s[3]*s[2];
958  sout[4] = s[3]*s[5] - s[0]*s[4];
959  sout[5] = s[3]*s[4] - s[1]*s[5];
960 
961  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
962 
963  sout[0] *= detinv;
964  sout[1] *= detinv;
965  sout[2] *= detinv;
966  sout[3] *= detinv;
967  sout[4] *= detinv;
968  sout[5] *= detinv;
969 }
970 
971 /*----------------------------------------------------------------------------*/
979 /*----------------------------------------------------------------------------*/
980 
981 static inline void
982 cs_math_33_product(const cs_real_t m1[3][3],
983  const cs_real_t m2[3][3],
984  cs_real_t mout[3][3])
985 {
986  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
987  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
988  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
989 
990  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
991  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
992  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
993 
994  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
995  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
996  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
997 }
998 
999 /*----------------------------------------------------------------------------*/
1008 /*----------------------------------------------------------------------------*/
1009 
1010 static inline void
1011 cs_math_33_transform_r_to_a(const cs_real_t m[3][3],
1012  const cs_real_t q[3][3],
1013  cs_real_t mout[3][3])
1014 {
1015  /* _m = M.Q */
1016  cs_real_33_t _m;
1017  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[1][0] + m[0][2]*q[2][0];
1018  _m[0][1] = m[0][0]*q[0][1] + m[0][1]*q[1][1] + m[0][2]*q[2][1];
1019  _m[0][2] = m[0][0]*q[0][2] + m[0][1]*q[1][2] + m[0][2]*q[2][2];
1020 
1021  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[1][0] + m[1][2]*q[2][0];
1022  _m[1][1] = m[1][0]*q[0][1] + m[1][1]*q[1][1] + m[1][2]*q[2][1];
1023  _m[1][2] = m[1][0]*q[0][2] + m[1][1]*q[1][2] + m[1][2]*q[2][2];
1024 
1025  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[1][0] + m[2][2]*q[2][0];
1026  _m[2][1] = m[2][0]*q[0][1] + m[2][1]*q[1][1] + m[2][2]*q[2][1];
1027  _m[2][2] = m[2][0]*q[0][2] + m[2][1]*q[1][2] + m[2][2]*q[2][2];
1028 
1029  /* mout = Q^t _m */
1030  mout[0][0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1031  mout[0][1] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1032  mout[0][2] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1033 
1034  mout[1][0] = q[0][1]*_m[0][0] + q[1][1]*_m[1][0] + q[2][1]*_m[2][0];
1035  mout[1][1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1036  mout[1][2] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1037 
1038  mout[2][0] = q[0][2]*_m[0][0] + q[1][2]*_m[1][0] + q[2][2]*_m[2][0];
1039  mout[2][1] = q[0][2]*_m[0][1] + q[1][2]*_m[1][1] + q[2][2]*_m[2][1];
1040  mout[2][2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1041 }
1042 
1043 /*----------------------------------------------------------------------------*/
1052 /*----------------------------------------------------------------------------*/
1053 
1054 static inline void
1055 cs_math_sym_33_transform_r_to_a(const cs_real_t m[6],
1056  const cs_real_t q[3][3],
1057  cs_real_t mout[6])
1058 {
1059  /* _m = M.Q */
1060  cs_real_33_t _m;
1061  _m[0][0] = m[0]*q[0][0] + m[3]*q[1][0] + m[5]*q[2][0];
1062  _m[0][1] = m[0]*q[0][1] + m[3]*q[1][1] + m[5]*q[2][1];
1063  _m[0][2] = m[0]*q[0][2] + m[3]*q[1][2] + m[5]*q[2][2];
1064 
1065  _m[1][0] = m[3]*q[0][0] + m[1]*q[1][0] + m[4]*q[2][0];
1066  _m[1][1] = m[3]*q[0][1] + m[1]*q[1][1] + m[4]*q[2][1];
1067  _m[1][2] = m[3]*q[0][2] + m[1]*q[1][2] + m[4]*q[2][2];
1068 
1069  _m[2][0] = m[5]*q[0][0] + m[4]*q[1][0] + m[2]*q[2][0];
1070  _m[2][1] = m[5]*q[0][1] + m[4]*q[1][1] + m[2]*q[2][1];
1071  _m[2][2] = m[5]*q[0][2] + m[4]*q[1][2] + m[2]*q[2][2];
1072 
1073  /* mout = Q^t _m */
1074  mout[0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1075  mout[1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1076  mout[2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1077 
1078  mout[3] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1079  mout[4] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1080  mout[5] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1081 
1082 }
1083 
1084 /*----------------------------------------------------------------------------*/
1093 /*----------------------------------------------------------------------------*/
1094 
1095 static inline void
1096 cs_math_33_transform_a_to_r(const cs_real_t m[3][3],
1097  const cs_real_t q[3][3],
1098  cs_real_t mout[3][3])
1099 {
1100  /* _m = M.Q^t */
1101  cs_real_33_t _m;
1102  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[0][1] + m[0][2]*q[0][2];
1103  _m[0][1] = m[0][0]*q[1][0] + m[0][1]*q[1][1] + m[0][2]*q[1][2];
1104  _m[0][2] = m[0][0]*q[2][0] + m[0][1]*q[2][1] + m[0][2]*q[2][2];
1105 
1106  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[0][1] + m[1][2]*q[0][2];
1107  _m[1][1] = m[1][0]*q[1][0] + m[1][1]*q[1][1] + m[1][2]*q[1][2];
1108  _m[1][2] = m[1][0]*q[2][0] + m[1][1]*q[2][1] + m[1][2]*q[2][2];
1109 
1110  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[0][1] + m[2][2]*q[0][2];
1111  _m[2][1] = m[2][0]*q[1][0] + m[2][1]*q[1][1] + m[2][2]*q[1][2];
1112  _m[2][2] = m[2][0]*q[2][0] + m[2][1]*q[2][1] + m[2][2]*q[2][2];
1113 
1114  /* mout = Q _m */
1115  mout[0][0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1116  mout[0][1] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1117  mout[0][2] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1118 
1119  mout[1][0] = q[1][0]*_m[0][0] + q[1][1]*_m[1][0] + q[1][2]*_m[2][0];
1120  mout[1][1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1121  mout[1][2] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1122 
1123  mout[2][0] = q[2][0]*_m[0][0] + q[2][1]*_m[1][0] + q[2][2]*_m[2][0];
1124  mout[2][1] = q[2][0]*_m[0][1] + q[2][1]*_m[1][1] + q[2][2]*_m[2][1];
1125  mout[2][2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1126 }
1127 
1128 /*----------------------------------------------------------------------------*/
1137 /*----------------------------------------------------------------------------*/
1138 
1139 static inline void
1140 cs_math_sym_33_transform_a_to_r(const cs_real_t m[6],
1141  const cs_real_t q[3][3],
1142  cs_real_t mout[6])
1143 {
1144  /* _m = M.Q^t */
1145  cs_real_33_t _m;
1146  _m[0][0] = m[0]*q[0][0] + m[3]*q[0][1] + m[5]*q[0][2];
1147  _m[0][1] = m[0]*q[1][0] + m[3]*q[1][1] + m[5]*q[1][2];
1148  _m[0][2] = m[0]*q[2][0] + m[3]*q[2][1] + m[5]*q[2][2];
1149 
1150  _m[1][0] = m[3]*q[0][0] + m[1]*q[0][1] + m[4]*q[0][2];
1151  _m[1][1] = m[3]*q[1][0] + m[1]*q[1][1] + m[4]*q[1][2];
1152  _m[1][2] = m[3]*q[2][0] + m[1]*q[2][1] + m[4]*q[2][2];
1153 
1154  _m[2][0] = m[5]*q[0][0] + m[4]*q[0][1] + m[2]*q[0][2];
1155  _m[2][1] = m[5]*q[1][0] + m[4]*q[1][1] + m[2]*q[1][2];
1156  _m[2][2] = m[5]*q[2][0] + m[4]*q[2][1] + m[2]*q[2][2];
1157 
1158  /* mout = Q _m */
1159  mout[0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1160  mout[1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1161  mout[2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1162 
1163 
1164  mout[3] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1165  mout[4] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1166  mout[5] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1167 }
1168 
1169 /*----------------------------------------------------------------------------*/
1178 /*----------------------------------------------------------------------------*/
1179 
1180 static inline void
1181 cs_math_33_extract_sym_ant(const cs_real_t m[3][3],
1182  cs_real_t m_sym[3][3],
1183  cs_real_t m_ant[3][3])
1184 {
1185  /* sym = 0.5 (m + m_transpose) */
1186  m_sym[0][0] = 0.5 * (m[0][0] + m[0][0]);
1187  m_sym[0][1] = 0.5 * (m[0][1] + m[1][0]);
1188  m_sym[0][2] = 0.5 * (m[0][2] + m[2][0]);
1189  m_sym[1][0] = 0.5 * (m[1][0] + m[0][1]);
1190  m_sym[1][1] = 0.5 * (m[1][1] + m[1][1]);
1191  m_sym[1][2] = 0.5 * (m[1][2] + m[2][1]);
1192  m_sym[2][0] = 0.5 * (m[2][0] + m[0][2]);
1193  m_sym[2][1] = 0.5 * (m[2][1] + m[1][2]);
1194  m_sym[2][2] = 0.5 * (m[2][2] + m[2][2]);
1195 
1196  /* ant = 0.5 (m - m_transpose) */
1197  m_ant[0][0] = 0.5 * (m[0][0] - m[0][0]);
1198  m_ant[0][1] = 0.5 * (m[0][1] - m[1][0]);
1199  m_ant[0][2] = 0.5 * (m[0][2] - m[2][0]);
1200  m_ant[1][0] = 0.5 * (m[1][0] - m[0][1]);
1201  m_ant[1][1] = 0.5 * (m[1][1] - m[1][1]);
1202  m_ant[1][2] = 0.5 * (m[1][2] - m[2][1]);
1203  m_ant[2][0] = 0.5 * (m[2][0] - m[0][2]);
1204  m_ant[2][1] = 0.5 * (m[2][1] - m[1][2]);
1205  m_ant[2][2] = 0.5 * (m[2][2] - m[2][2]);
1206 }
1207 
1208 /*----------------------------------------------------------------------------*/
1216 /*----------------------------------------------------------------------------*/
1217 
1218 static inline void
1219 cs_math_33_product_add(const cs_real_t m1[3][3],
1220  const cs_real_t m2[3][3],
1221  cs_real_t mout[restrict 3][3])
1222 {
1223  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1224  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1225  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1226 
1227  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1228  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1229  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1230 
1231  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1232  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1233  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1234 }
1235 
1236 /*----------------------------------------------------------------------------*/
1250 /*----------------------------------------------------------------------------*/
1251 
1252 static inline void
1253 cs_math_sym_33_product(const cs_real_t s1[6],
1254  const cs_real_t s2[6],
1255  cs_real_t sout[restrict 6])
1256 {
1257  /* S11 */
1258  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1259  /* S22 */
1260  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1261  /* S33 */
1262  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1263  /* S12 = S21 */
1264  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1265  /* S23 = S32 */
1266  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1267  /* S13 = S31 */
1268  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1269 }
1270 
1271 /*----------------------------------------------------------------------------*/
1279 /*----------------------------------------------------------------------------*/
1280 
1281 static inline void
1282 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1283  cs_real_t sout[restrict 6][6])
1284 {
1285  const int t2v[3][3] = {{0, 3, 5},
1286  {3, 1, 4},
1287  {5, 4, 2}};
1288 
1289  const int iv2t[6] = {0, 1, 2, 0, 1, 0};
1290  const int jv2t[6] = {0, 1, 2, 1, 2, 2};
1291 
1292  for (int i = 0; i < 6; i++) {
1293  for (int j = 0; j < 6; j++)
1294  sout[i][j] = 0;
1295  }
1296 
1297  /* Consider : W = R*s^t + s*R.
1298  * W_ij = Sum_{k<3} [s_jk*r_ik + s_ik*r_jk]
1299  * We look for A such as A*R = W
1300  */
1301  for (int i = 0; i < 6; i++) {
1302  int ii = iv2t[i];
1303  int jj = jv2t[i];
1304  for (int k = 0; k < 3; k++) {
1305  int ik = t2v[k][ii];
1306  int jk = t2v[k][jj];
1307 
1308  sout[ik][i] += s[k][jj];
1309  sout[jk][i] += s[k][ii];
1310  }
1311  }
1312 }
1313 
1314 /*----------------------------------------------------------------------------*/
1326 /*----------------------------------------------------------------------------*/
1327 
1328 static inline void
1329 cs_math_sym_33_double_product(const cs_real_t s1[6],
1330  const cs_real_t s2[6],
1331  const cs_real_t s3[6],
1332  cs_real_t sout[restrict 3][3])
1333 {
1334  cs_real_t _sout[3][3];
1335 
1336  /* S11 */
1337  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1338  /* S22 */
1339  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1340  /* S33 */
1341  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1342  /* S12 */
1343  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1344  /* S21 */
1345  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1346  /* S23 */
1347  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1348  /* S32 */
1349  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1350  /* S13 */
1351  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1352  /* S31 */
1353  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1354 
1355  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1356  /* S22 */
1357  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1358  /* S33 */
1359  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1360  /* S12 */
1361  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1362  /* S21 */
1363  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1364  /* S23 */
1365  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1366  /* S32 */
1367  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1368  /* S13 */
1369  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1370  /* S31 */
1371  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1372 }
1373 
1374 /*----------------------------------------------------------------------------*/
1381 /*----------------------------------------------------------------------------*/
1382 
1383 static inline void
1384 cs_nvec3(const cs_real_3_t v,
1385  cs_nvec3_t *qv)
1386 {
1387  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1388 
1389  qv->meas = magnitude;
1390  if (fabs(magnitude) > cs_math_zero_threshold) {
1391 
1392  const cs_real_t inv = 1/magnitude;
1393  qv->unitv[0] = inv * v[0];
1394  qv->unitv[1] = inv * v[1];
1395  qv->unitv[2] = inv * v[2];
1396 
1397  }
1398  else
1399  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1400 }
1401 
1402 /*=============================================================================
1403  * Public function prototypes
1404  *============================================================================*/
1405 
1406 /*----------------------------------------------------------------------------*/
1410 /*----------------------------------------------------------------------------*/
1411 
1412 void
1413 cs_math_set_machine_epsilon(void);
1414 
1415 /*----------------------------------------------------------------------------*/
1419 /*----------------------------------------------------------------------------*/
1420 
1421 double
1422 cs_math_get_machine_epsilon(void);
1423 
1424 /*----------------------------------------------------------------------------*/
1434 /*----------------------------------------------------------------------------*/
1435 
1436 void
1437 cs_math_3_length_unitv(const cs_real_t xa[3],
1438  const cs_real_t xb[3],
1439  cs_real_t *len,
1440  cs_real_3_t unitv);
1441 
1442 /*----------------------------------------------------------------------------*/
1454 /*----------------------------------------------------------------------------*/
1455 
1456 void
1457 cs_math_sym_33_eigen(const cs_real_t m[6],
1458  cs_real_t eig_vals[3]);
1459 
1460 /*----------------------------------------------------------------------------*/
1473 /*----------------------------------------------------------------------------*/
1474 
1475 void
1476 cs_math_33_eigen(const cs_real_t m[3][3],
1477  cs_real_t *eig_ratio,
1478  cs_real_t *eig_max);
1479 
1480 /*----------------------------------------------------------------------------*/
1491 /*----------------------------------------------------------------------------*/
1492 
1493 double
1494 cs_math_surftri(const cs_real_t xv[3],
1495  const cs_real_t xe[3],
1496  const cs_real_t xf[3]);
1497 
1498 /*----------------------------------------------------------------------------*/
1510 /*----------------------------------------------------------------------------*/
1511 
1512 double
1513 cs_math_voltet(const cs_real_t xv[3],
1514  const cs_real_t xe[3],
1515  const cs_real_t xf[3],
1516  const cs_real_t xc[3]);
1517 
1518 /*----------------------------------------------------------------------------*/
1531 /*----------------------------------------------------------------------------*/
1532 
1533 void
1534 cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1535  const cs_real_t tol_err,
1536  cs_real_t eig_val[restrict 3],
1537  cs_real_t eig_vec[restrict 3][3]);
1538 
1539 /*----------------------------------------------------------------------------*/
1549 /*----------------------------------------------------------------------------*/
1550 
1551 void
1552 cs_math_fact_lu(cs_lnum_t n_blocks,
1553  const int b_size,
1554  const cs_real_t *a,
1555  cs_real_t *a_lu);
1556 
1557 /*----------------------------------------------------------------------------*/
1567 /*----------------------------------------------------------------------------*/
1568 
1569 void
1570 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1571  const int n,
1572  cs_real_t x[],
1573  const cs_real_t b[]);
1574 
1575 /*----------------------------------------------------------------------------*/
1576 
1577 END_C_DECLS
1578 
1579 #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:949
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:658
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:638
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:594
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:675
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:879
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:436
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:527
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:788
cs_math_1ov6
const cs_real_t cs_math_1ov6
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:512
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:738
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:1329
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:370
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:417
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:460
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:1055
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:505
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:1253
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:1096
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:697
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:1219
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:914
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:1011
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:148
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:372
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:718
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:401
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:847
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:484
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:762
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:373
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:1384
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:513
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:982
cs_math_3_orthonormal_basis
static void cs_math_3_orthonormal_basis(const cs_real_t vect[3], cs_real_t axes[3][3])
Build an orthonormal basis based on a first vector "vect". axes[0] is vect normalized, while (axes[0], axes[1], axes[23]) is an orthonormal base.
Definition: cs_math.h:810
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:344
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:616
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:1181
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:573
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:552
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:1282
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:1140
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:148