8.1
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-2023 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 /*----------------------------------------------------------------------------*/
199 /*----------------------------------------------------------------------------*/
200 
201 static inline cs_real_t
202 cs_math_clamp(cs_real_t x,
203  cs_real_t xmin,
204  cs_real_t xmax)
205 {
206  cs_real_t ret = cs_math_fmin(xmax, cs_math_fmax(xmin, x));
207 
208  return ret;
209 }
210 
211 /*----------------------------------------------------------------------------*/
219 /*----------------------------------------------------------------------------*/
220 
221 static inline cs_real_t
222 cs_math_sq(cs_real_t x)
223 {
224  return x*x;
225 }
226 
227 /*----------------------------------------------------------------------------*/
235 /*----------------------------------------------------------------------------*/
236 
237 static inline cs_real_t
238 cs_math_pow2(cs_real_t x)
239 {
240  return x*x;
241 }
242 
243 /*----------------------------------------------------------------------------*/
251 /*----------------------------------------------------------------------------*/
252 
253 static inline cs_real_t
254 cs_math_pow3(cs_real_t x)
255 {
256  return x*x*x;
257 }
258 
259 /*----------------------------------------------------------------------------*/
267 /*----------------------------------------------------------------------------*/
268 
269 static inline cs_real_t
270 cs_math_pow4(cs_real_t x)
271 {
272  return (x*x)*(x*x);
273 }
274 
275 /*----------------------------------------------------------------------------*/
283 /*----------------------------------------------------------------------------*/
284 
285 static inline cs_real_t
286 cs_math_pow5(cs_real_t x)
287 {
288  return x*(x*x)*(x*x);
289 }
290 
291 /*----------------------------------------------------------------------------*/
301 /*----------------------------------------------------------------------------*/
302 
303 static inline cs_real_t
304 cs_math_3_distance(const cs_real_t xa[3],
305  const cs_real_t xb[3])
306 {
307  cs_real_t v[3];
308 
309  v[0] = xb[0] - xa[0];
310  v[1] = xb[1] - xa[1];
311  v[2] = xb[2] - xa[2];
312 
313  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
314 }
315 
316 /*----------------------------------------------------------------------------*/
326 /*----------------------------------------------------------------------------*/
327 
328 static inline cs_real_t
329 cs_math_3_distance_dot_product(const cs_real_t xa[3],
330  const cs_real_t xb[3],
331  const cs_real_t xc[3])
332 {
333  return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
334 }
335 
336 /*----------------------------------------------------------------------------*/
346 /*----------------------------------------------------------------------------*/
347 
348 static inline cs_real_t
349 cs_math_3_square_distance(const cs_real_t xa[3],
350  const cs_real_t xb[3])
351 {
352  cs_real_t v[3] = {xb[0] - xa[0],
353  xb[1] - xa[1],
354  xb[2] - xa[2]};
355 
356  return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
357 }
358 
359 /*----------------------------------------------------------------------------*/
368 /*----------------------------------------------------------------------------*/
369 
370 static inline cs_real_t
371 cs_math_3_dot_product(const cs_real_t u[3],
372  const cs_real_t v[3])
373 {
374  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
375 
376  return uv;
377 }
378 
379 /*----------------------------------------------------------------------------*/
389 /*----------------------------------------------------------------------------*/
390 
391 static inline cs_real_t
392 cs_math_3_33_3_dot_product(const cs_real_t n1[3],
393  const cs_real_t t[3][3],
394  const cs_real_t n2[3])
395 {
396  cs_real_t n_t_n
397  = ( n1[0]*t[0][0]*n2[0] + n1[1]*t[1][0]*n2[0] + n1[2]*t[2][0]*n2[0]
398  + n1[0]*t[0][1]*n2[1] + n1[1]*t[1][1]*n2[1] + n1[2]*t[2][1]*n2[1]
399  + n1[0]*t[0][2]*n2[2] + n1[1]*t[1][2]*n2[2] + n1[2]*t[2][2]*n2[2]);
400  return n_t_n;
401 }
402 
403 /*----------------------------------------------------------------------------*/
417 /*----------------------------------------------------------------------------*/
418 
419 static inline cs_real_t
420 cs_math_3_sym_33_3_dot_product(const cs_real_t n1[3],
421  const cs_real_t t[6],
422  const cs_real_t n2[3])
423 {
424  return ( n1[0] * (t[0]*n2[0] + t[3]*n2[1] + t[5]*n2[2])
425  + n1[1] * (t[3]*n2[0] + t[1]*n2[1] + t[4]*n2[2])
426  + n1[2] * (t[5]*n2[0] + t[4]*n2[1] + t[2]*n2[2]));
427 }
428 
429 /*----------------------------------------------------------------------------*/
437 /*----------------------------------------------------------------------------*/
438 
439 static inline cs_real_t
440 cs_math_3_norm(const cs_real_t v[3])
441 {
442  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
443 }
444 
445 /*----------------------------------------------------------------------------*/
453 /*----------------------------------------------------------------------------*/
454 
455 static inline cs_real_t
456 cs_math_3_square_norm(const cs_real_t v[3])
457 {
458  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
459 
460  return v2;
461 }
462 
463 /*----------------------------------------------------------------------------*/
472 /*----------------------------------------------------------------------------*/
473 
474 static inline void
475 cs_math_3_normalise(const cs_real_t vin[3],
476  cs_real_t vout[3])
477 {
478  cs_real_t norm = cs_math_3_norm(vin);
479 
480  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
481 
482  vout[0] = inv_norm * vin[0];
483  vout[1] = inv_norm * vin[1];
484  vout[2] = inv_norm * vin[2];
485 }
486 
487 /*----------------------------------------------------------------------------*/
496 /*----------------------------------------------------------------------------*/
497 
498 static inline void
499 cs_math_3_normalize(const cs_real_t vin[3],
500  cs_real_t vout[3])
501 {
502  cs_real_t norm = cs_math_3_norm(vin);
503 
504  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
505 
506  vout[0] = inv_norm * vin[0];
507  vout[1] = inv_norm * vin[1];
508  vout[2] = inv_norm * vin[2];
509 }
510 
511 /*----------------------------------------------------------------------------*/
520 /*----------------------------------------------------------------------------*/
521 
522 static inline void
523 cs_math_3_orthogonal_projection(const cs_real_t n[3],
524  const cs_real_t v[3],
525  cs_real_t vout[restrict 3])
526 {
527  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
528  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
529  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
530 }
531 
532 /*----------------------------------------------------------------------------*/
541 /*----------------------------------------------------------------------------*/
542 
543 static inline void
544 cs_math_3_normal_scaling(const cs_real_t n[3],
545  cs_real_t factor,
546  cs_real_t v[3])
547 {
548  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
549  for (int i = 0; i < 3; i++)
550  v[i] += v_dot_n * n[i];
551 }
552 
553 /*----------------------------------------------------------------------------*/
563 /*----------------------------------------------------------------------------*/
564 
565 static inline void
566 cs_math_33_normal_scaling_add(const cs_real_t n[3],
567  cs_real_t factor,
568  cs_real_t t[3][3])
569 {
570  cs_real_t n_t_n = (factor -1.) *
571  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
572  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
573  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
574  for (int i = 0; i < 3; i++) {
575  for (int j = 0; j < 3; j++)
576  t[i][j] += n_t_n * n[i] * n[j];
577  }
578 }
579 /*----------------------------------------------------------------------------*/
588 /*----------------------------------------------------------------------------*/
589 
590 static inline void
591 cs_math_33_3_product(const cs_real_t m[3][3],
592  const cs_real_t v[3],
593  cs_real_t mv[restrict 3])
594 {
595  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
596  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
597  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
598 }
599 
600 /*----------------------------------------------------------------------------*/
609 /*----------------------------------------------------------------------------*/
610 
611 static inline void
612 cs_math_33_3_product_add(const cs_real_t m[3][3],
613  const cs_real_t v[3],
614  cs_real_t mv[restrict 3])
615 {
616  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
617  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
618  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
619 }
620 
621 /*----------------------------------------------------------------------------*/
630 /*----------------------------------------------------------------------------*/
631 
632 static inline void
633 cs_math_33t_3_product(const cs_real_t m[3][3],
634  const cs_real_t v[3],
635  cs_real_t mv[restrict 3])
636 {
637  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
638  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
639  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
640 }
641 
642 /*----------------------------------------------------------------------------*/
652 /*----------------------------------------------------------------------------*/
653 
654 static inline void
655 cs_math_sym_33_3_product(const cs_real_t m[6],
656  const cs_real_t v[3],
657  cs_real_t mv[restrict 3])
658 {
659  mv[0] = m[0]*v[0] + m[3]*v[1] + m[5]*v[2];
660  mv[1] = m[3]*v[0] + m[1]*v[1] + m[4]*v[2];
661  mv[2] = m[5]*v[0] + m[4]*v[1] + m[2]*v[2];
662 }
663 
664 /*----------------------------------------------------------------------------*/
674 /*----------------------------------------------------------------------------*/
675 
676 static inline void
677 cs_math_sym_33_3_product_add(const cs_real_t m[6],
678  const cs_real_t v[3],
679  cs_real_t mv[restrict 3])
680 {
681  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
682  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
683  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
684 }
685 
686 /*----------------------------------------------------------------------------*/
696 /*----------------------------------------------------------------------------*/
697 
698 static inline cs_real_t
699 cs_math_sym_33_sym_33_product_trace(const cs_real_t m1[6],
700  const cs_real_t m2[6])
701 {
702  return m1[0]*m2[0] + 2.*m1[3]*m2[3] + 2.*m1[5]*m2[5]
703  + m1[1]*m2[1] + 2.*m1[4]*m2[4]
704  + m1[2]*m2[2];
705 }
706 
707 /*----------------------------------------------------------------------------*/
715 /*----------------------------------------------------------------------------*/
716 
717 static inline cs_real_t
718 cs_math_33_trace(const cs_real_t t[3][3])
719 {
720  return (t[0][0] + t[1][1] + t[2][2]);
721 }
722 
723 /*----------------------------------------------------------------------------*/
731 /*----------------------------------------------------------------------------*/
732 
733 static inline cs_real_t
734 cs_math_6_trace(const cs_real_t t[6])
735 {
736  return (t[0] + t[1] + t[2]);
737 }
738 
739 /*----------------------------------------------------------------------------*/
748 /*----------------------------------------------------------------------------*/
749 
750 static inline void
751 cs_math_66_6_product(const cs_real_t m[6][6],
752  const cs_real_t v[6],
753  cs_real_t mv[restrict 6])
754 {
755  for (int i = 0; i < 6; i++) {
756  for (int j = 0; j < 6; j++)
757  mv[i] = m[i][j] * v[j];
758  }
759 }
760 
761 /*----------------------------------------------------------------------------*/
770 /*----------------------------------------------------------------------------*/
771 
772 static inline void
773 cs_math_66_6_product_add(const cs_real_t m[6][6],
774  const cs_real_t v[6],
775  cs_real_t mv[restrict 6])
776 {
777  for (int i = 0; i < 6; i++) {
778  for (int j = 0; j < 6; j++)
779  mv[i] += m[i][j] * v[j];
780  }
781 }
782 
783 /*----------------------------------------------------------------------------*/
791 /*----------------------------------------------------------------------------*/
792 
793 static inline cs_real_t
794 cs_math_33_determinant(const cs_real_t m[3][3])
795 {
796  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
797  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
798  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
799 
800  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
801 }
802 
803 /*----------------------------------------------------------------------------*/
811 /*----------------------------------------------------------------------------*/
812 
813 static inline cs_real_t
814 cs_math_sym_33_determinant(const cs_real_6_t m)
815 {
816  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
817  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
818  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
819 
820  return m[0]*com0 + m[3]*com1 + m[5]*com2;
821 }
822 
823 /*----------------------------------------------------------------------------*/
831 /*----------------------------------------------------------------------------*/
832 
833 #if defined(__INTEL_COMPILER)
834 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
835 #endif
836 
837 static inline void
838 cs_math_3_cross_product(const cs_real_t u[3],
839  const cs_real_t v[3],
840  cs_real_t uv[restrict 3])
841 {
842  uv[0] = u[1]*v[2] - u[2]*v[1];
843  uv[1] = u[2]*v[0] - u[0]*v[2];
844  uv[2] = u[0]*v[1] - u[1]*v[0];
845 }
846 
847 /*----------------------------------------------------------------------------*/
857 /*----------------------------------------------------------------------------*/
858 
859 #if defined(__INTEL_COMPILER)
860 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
861 #endif
862 
863 static inline cs_real_t
864 cs_math_3_triple_product(const cs_real_t u[3],
865  const cs_real_t v[3],
866  const cs_real_t w[3])
867 {
868  return (u[1]*v[2] - u[2]*v[1]) * w[0]
869  + (u[2]*v[0] - u[0]*v[2]) * w[1]
870  + (u[0]*v[1] - u[1]*v[0]) * w[2];
871 }
872 
873 /*----------------------------------------------------------------------------*/
883 /*----------------------------------------------------------------------------*/
884 
885 static inline void
886 cs_math_3_orthonormal_basis(const cs_real_t vect[3],
887  cs_real_t axes[3][3])
888 {
889  assert(cs_math_3_norm(vect) > cs_math_zero_threshold);
890 
891  // Compute normal - third axis
892  cs_math_3_normalize(vect, axes[2]);
893 
894  // Compute base's first axis
895  // First test projection of Ox
896  cs_real_t Ox[3] = {1., 0., 0.};
897  cs_real_t w[3] = {0.};
898 
899  cs_math_3_orthogonal_projection(axes[2], Ox, w);
900 
901  // If Ox projection is null, project Oy
902  if (cs_math_3_norm(w) < cs_math_zero_threshold) {
903  cs_real_t Oy[3] = {0., 1., 0.};
904  cs_math_3_orthogonal_projection(axes[2], Oy, w);
905  }
906 
907  cs_math_3_normalize(w, axes[0]);
908 
909  // Compute base's second axis using cross product
910  cs_math_3_cross_product(axes[2], axes[0], axes[1]);
911 }
912 
913 /*----------------------------------------------------------------------------*/
920 /*----------------------------------------------------------------------------*/
921 
922 static inline void
923 cs_math_33_inv_cramer(const cs_real_t in[3][3],
924  cs_real_t out[3][3])
925 {
926  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
927  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
928  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
929 
930  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
931  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
932  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
933 
934  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
935  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
936  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
937 
938  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
939  const double invdet = 1./det;
940 
941  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
942  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
943  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
944 }
945 
946 /*----------------------------------------------------------------------------*/
952 /*----------------------------------------------------------------------------*/
953 
954 static inline void
955 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
956 {
957  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
958  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
959  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
960  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
961  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
962  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
963  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
964  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
965  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
966 
967  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
968 
969  a[0][0] = a00 * det_inv;
970  a[0][1] = a01 * det_inv;
971  a[0][2] = a02 * det_inv;
972  a[1][0] = a10 * det_inv;
973  a[1][1] = a11 * det_inv;
974  a[1][2] = a12 * det_inv;
975  a[2][0] = a20 * det_inv;
976  a[2][1] = a21 * det_inv;
977  a[2][2] = a22 * det_inv;
978 }
979 
980 /*----------------------------------------------------------------------------*/
987 /*----------------------------------------------------------------------------*/
988 
989 static inline void
990 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
991 {
992  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
993  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
994  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
995  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
996  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
997  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
998 
999  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
1000 
1001  a[0][0] = a00 * det_inv;
1002  a[0][1] = a01 * det_inv;
1003  a[0][2] = a02 * det_inv;
1004  a[1][0] = a01 * det_inv;
1005  a[1][1] = a11 * det_inv;
1006  a[1][2] = a12 * det_inv;
1007  a[2][0] = a02 * det_inv;
1008  a[2][1] = a12 * det_inv;
1009  a[2][2] = a22 * det_inv;
1010 }
1011 
1012 /*----------------------------------------------------------------------------*/
1022 /*----------------------------------------------------------------------------*/
1023 
1024 static inline void
1025 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
1026  cs_real_t sout[restrict 6])
1027 {
1028  double detinv;
1029 
1030  sout[0] = s[1]*s[2] - s[4]*s[4];
1031  sout[1] = s[0]*s[2] - s[5]*s[5];
1032  sout[2] = s[0]*s[1] - s[3]*s[3];
1033  sout[3] = s[4]*s[5] - s[3]*s[2];
1034  sout[4] = s[3]*s[5] - s[0]*s[4];
1035  sout[5] = s[3]*s[4] - s[1]*s[5];
1036 
1037  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
1038 
1039  sout[0] *= detinv;
1040  sout[1] *= detinv;
1041  sout[2] *= detinv;
1042  sout[3] *= detinv;
1043  sout[4] *= detinv;
1044  sout[5] *= detinv;
1045 }
1046 
1047 /*----------------------------------------------------------------------------*/
1055 /*----------------------------------------------------------------------------*/
1056 
1057 static inline void
1058 cs_math_33_product(const cs_real_t m1[3][3],
1059  const cs_real_t m2[3][3],
1060  cs_real_t mout[3][3])
1061 {
1062  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1063  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1064  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1065 
1066  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1067  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1068  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1069 
1070  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1071  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1072  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1073 }
1074 
1075 /*----------------------------------------------------------------------------*/
1084 /*----------------------------------------------------------------------------*/
1085 
1086 static inline void
1087 cs_math_33_transform_r_to_a(const cs_real_t m[3][3],
1088  const cs_real_t q[3][3],
1089  cs_real_t mout[3][3])
1090 {
1091  /* _m = M.Q */
1092  cs_real_33_t _m;
1093  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[1][0] + m[0][2]*q[2][0];
1094  _m[0][1] = m[0][0]*q[0][1] + m[0][1]*q[1][1] + m[0][2]*q[2][1];
1095  _m[0][2] = m[0][0]*q[0][2] + m[0][1]*q[1][2] + m[0][2]*q[2][2];
1096 
1097  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[1][0] + m[1][2]*q[2][0];
1098  _m[1][1] = m[1][0]*q[0][1] + m[1][1]*q[1][1] + m[1][2]*q[2][1];
1099  _m[1][2] = m[1][0]*q[0][2] + m[1][1]*q[1][2] + m[1][2]*q[2][2];
1100 
1101  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[1][0] + m[2][2]*q[2][0];
1102  _m[2][1] = m[2][0]*q[0][1] + m[2][1]*q[1][1] + m[2][2]*q[2][1];
1103  _m[2][2] = m[2][0]*q[0][2] + m[2][1]*q[1][2] + m[2][2]*q[2][2];
1104 
1105  /* mout = Q^t _m */
1106  mout[0][0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1107  mout[0][1] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1108  mout[0][2] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1109 
1110  mout[1][0] = q[0][1]*_m[0][0] + q[1][1]*_m[1][0] + q[2][1]*_m[2][0];
1111  mout[1][1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1112  mout[1][2] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1113 
1114  mout[2][0] = q[0][2]*_m[0][0] + q[1][2]*_m[1][0] + q[2][2]*_m[2][0];
1115  mout[2][1] = q[0][2]*_m[0][1] + q[1][2]*_m[1][1] + q[2][2]*_m[2][1];
1116  mout[2][2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1117 }
1118 
1119 /*----------------------------------------------------------------------------*/
1128 /*----------------------------------------------------------------------------*/
1129 
1130 static inline void
1131 cs_math_sym_33_transform_r_to_a(const cs_real_t m[6],
1132  const cs_real_t q[3][3],
1133  cs_real_t mout[6])
1134 {
1135  /* _m = M.Q */
1136  cs_real_33_t _m;
1137  _m[0][0] = m[0]*q[0][0] + m[3]*q[1][0] + m[5]*q[2][0];
1138  _m[0][1] = m[0]*q[0][1] + m[3]*q[1][1] + m[5]*q[2][1];
1139  _m[0][2] = m[0]*q[0][2] + m[3]*q[1][2] + m[5]*q[2][2];
1140 
1141  _m[1][0] = m[3]*q[0][0] + m[1]*q[1][0] + m[4]*q[2][0];
1142  _m[1][1] = m[3]*q[0][1] + m[1]*q[1][1] + m[4]*q[2][1];
1143  _m[1][2] = m[3]*q[0][2] + m[1]*q[1][2] + m[4]*q[2][2];
1144 
1145  _m[2][0] = m[5]*q[0][0] + m[4]*q[1][0] + m[2]*q[2][0];
1146  _m[2][1] = m[5]*q[0][1] + m[4]*q[1][1] + m[2]*q[2][1];
1147  _m[2][2] = m[5]*q[0][2] + m[4]*q[1][2] + m[2]*q[2][2];
1148 
1149  /* mout = Q^t _m */
1150  mout[0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1151  mout[1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1152  mout[2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1153 
1154  mout[3] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1155  mout[4] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1156  mout[5] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1157 
1158 }
1159 
1160 /*----------------------------------------------------------------------------*/
1169 /*----------------------------------------------------------------------------*/
1170 
1171 static inline void
1172 cs_math_33_transform_a_to_r(const cs_real_t m[3][3],
1173  const cs_real_t q[3][3],
1174  cs_real_t mout[3][3])
1175 {
1176  /* _m = M.Q^t */
1177  cs_real_33_t _m;
1178  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[0][1] + m[0][2]*q[0][2];
1179  _m[0][1] = m[0][0]*q[1][0] + m[0][1]*q[1][1] + m[0][2]*q[1][2];
1180  _m[0][2] = m[0][0]*q[2][0] + m[0][1]*q[2][1] + m[0][2]*q[2][2];
1181 
1182  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[0][1] + m[1][2]*q[0][2];
1183  _m[1][1] = m[1][0]*q[1][0] + m[1][1]*q[1][1] + m[1][2]*q[1][2];
1184  _m[1][2] = m[1][0]*q[2][0] + m[1][1]*q[2][1] + m[1][2]*q[2][2];
1185 
1186  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[0][1] + m[2][2]*q[0][2];
1187  _m[2][1] = m[2][0]*q[1][0] + m[2][1]*q[1][1] + m[2][2]*q[1][2];
1188  _m[2][2] = m[2][0]*q[2][0] + m[2][1]*q[2][1] + m[2][2]*q[2][2];
1189 
1190  /* mout = Q _m */
1191  mout[0][0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1192  mout[0][1] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1193  mout[0][2] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1194 
1195  mout[1][0] = q[1][0]*_m[0][0] + q[1][1]*_m[1][0] + q[1][2]*_m[2][0];
1196  mout[1][1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1197  mout[1][2] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1198 
1199  mout[2][0] = q[2][0]*_m[0][0] + q[2][1]*_m[1][0] + q[2][2]*_m[2][0];
1200  mout[2][1] = q[2][0]*_m[0][1] + q[2][1]*_m[1][1] + q[2][2]*_m[2][1];
1201  mout[2][2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1202 }
1203 
1204 /*----------------------------------------------------------------------------*/
1213 /*----------------------------------------------------------------------------*/
1214 
1215 static inline void
1216 cs_math_sym_33_transform_a_to_r(const cs_real_t m[6],
1217  const cs_real_t q[3][3],
1218  cs_real_t mout[6])
1219 {
1220  /* _m = M.Q^t */
1221  cs_real_33_t _m;
1222  _m[0][0] = m[0]*q[0][0] + m[3]*q[0][1] + m[5]*q[0][2];
1223  _m[0][1] = m[0]*q[1][0] + m[3]*q[1][1] + m[5]*q[1][2];
1224  _m[0][2] = m[0]*q[2][0] + m[3]*q[2][1] + m[5]*q[2][2];
1225 
1226  _m[1][0] = m[3]*q[0][0] + m[1]*q[0][1] + m[4]*q[0][2];
1227  _m[1][1] = m[3]*q[1][0] + m[1]*q[1][1] + m[4]*q[1][2];
1228  _m[1][2] = m[3]*q[2][0] + m[1]*q[2][1] + m[4]*q[2][2];
1229 
1230  _m[2][0] = m[5]*q[0][0] + m[4]*q[0][1] + m[2]*q[0][2];
1231  _m[2][1] = m[5]*q[1][0] + m[4]*q[1][1] + m[2]*q[1][2];
1232  _m[2][2] = m[5]*q[2][0] + m[4]*q[2][1] + m[2]*q[2][2];
1233 
1234  /* mout = Q _m */
1235  mout[0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1236  mout[1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1237  mout[2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1238 
1239 
1240  mout[3] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1241  mout[4] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1242  mout[5] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1243 }
1244 
1245 /*----------------------------------------------------------------------------*/
1254 /*----------------------------------------------------------------------------*/
1255 
1256 static inline void
1257 cs_math_33_extract_sym_ant(const cs_real_t m[3][3],
1258  cs_real_t m_sym[3][3],
1259  cs_real_t m_ant[3][3])
1260 {
1261  /* sym = 0.5 (m + m_transpose) */
1262  m_sym[0][0] = 0.5 * (m[0][0] + m[0][0]);
1263  m_sym[0][1] = 0.5 * (m[0][1] + m[1][0]);
1264  m_sym[0][2] = 0.5 * (m[0][2] + m[2][0]);
1265  m_sym[1][0] = 0.5 * (m[1][0] + m[0][1]);
1266  m_sym[1][1] = 0.5 * (m[1][1] + m[1][1]);
1267  m_sym[1][2] = 0.5 * (m[1][2] + m[2][1]);
1268  m_sym[2][0] = 0.5 * (m[2][0] + m[0][2]);
1269  m_sym[2][1] = 0.5 * (m[2][1] + m[1][2]);
1270  m_sym[2][2] = 0.5 * (m[2][2] + m[2][2]);
1271 
1272  /* ant = 0.5 (m - m_transpose) */
1273  m_ant[0][0] = 0.5 * (m[0][0] - m[0][0]);
1274  m_ant[0][1] = 0.5 * (m[0][1] - m[1][0]);
1275  m_ant[0][2] = 0.5 * (m[0][2] - m[2][0]);
1276  m_ant[1][0] = 0.5 * (m[1][0] - m[0][1]);
1277  m_ant[1][1] = 0.5 * (m[1][1] - m[1][1]);
1278  m_ant[1][2] = 0.5 * (m[1][2] - m[2][1]);
1279  m_ant[2][0] = 0.5 * (m[2][0] - m[0][2]);
1280  m_ant[2][1] = 0.5 * (m[2][1] - m[1][2]);
1281  m_ant[2][2] = 0.5 * (m[2][2] - m[2][2]);
1282 }
1283 
1284 /*----------------------------------------------------------------------------*/
1293 /*----------------------------------------------------------------------------*/
1294 
1295 static inline cs_real_t
1296 cs_math_33_main_invariant_2(const cs_real_t m[3][3])
1297 {
1298  /* sym = 0.5 (m + m_transpose) */
1299  return cs_math_pow2(m[0][0])
1300  + cs_math_pow2(m[1][1])
1301  + cs_math_pow2(m[2][2])
1302  + 0.5 * ( cs_math_pow2(m[0][1] + m[1][0])
1303  + cs_math_pow2(m[0][2] + m[2][0])
1304  + cs_math_pow2(m[1][2] + m[2][1]));
1305 }
1306 
1307 /*----------------------------------------------------------------------------*/
1315 /*----------------------------------------------------------------------------*/
1316 
1317 static inline void
1318 cs_math_33_product_add(const cs_real_t m1[3][3],
1319  const cs_real_t m2[3][3],
1320  cs_real_t mout[restrict 3][3])
1321 {
1322  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1323  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1324  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1325 
1326  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1327  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1328  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1329 
1330  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1331  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1332  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1333 }
1334 
1335 /*----------------------------------------------------------------------------*/
1349 /*----------------------------------------------------------------------------*/
1350 
1351 static inline void
1352 cs_math_sym_33_product(const cs_real_t s1[6],
1353  const cs_real_t s2[6],
1354  cs_real_t sout[restrict 6])
1355 {
1356  /* S11 */
1357  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1358  /* S22 */
1359  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1360  /* S33 */
1361  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1362  /* S12 = S21 */
1363  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1364  /* S23 = S32 */
1365  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1366  /* S13 = S31 */
1367  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1368 }
1369 
1370 /*----------------------------------------------------------------------------*/
1378 /*----------------------------------------------------------------------------*/
1379 
1380 static inline void
1381 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1382  cs_real_t sout[restrict 6][6])
1383 {
1384  const int t2v[3][3] = {{0, 3, 5},
1385  {3, 1, 4},
1386  {5, 4, 2}};
1387 
1388  const int iv2t[6] = {0, 1, 2, 0, 1, 0};
1389  const int jv2t[6] = {0, 1, 2, 1, 2, 2};
1390 
1391  for (int i = 0; i < 6; i++) {
1392  for (int j = 0; j < 6; j++)
1393  sout[i][j] = 0;
1394  }
1395 
1396  /* Consider : W = s*R + R*s^t .
1397  * W_ij = Sum_{k<3} [s_ik*r_jk + s_jk*r_ik]
1398  * We look for A_(ij,pq) such as A*R = W
1399  *
1400  * so
1401  * A_(ij,jk) takes s_ik
1402  * and
1403  * A_(ij,ik) takes s_jk
1404  */
1405  for (int ij = 0; ij < 6; ij++) {
1406  int i = iv2t[ij];
1407  int j = jv2t[ij];
1408  for (int k = 0; k < 3; k++) {
1409  int ik = t2v[i][k];
1410  int jk = t2v[j][k];
1411 
1412  sout[ij][ik] += s[j][k];
1413  sout[ij][jk] += s[i][k];
1414  }
1415  }
1416 }
1417 
1418 /*----------------------------------------------------------------------------*/
1430 /*----------------------------------------------------------------------------*/
1431 
1432 static inline void
1433 cs_math_sym_33_double_product(const cs_real_t s1[6],
1434  const cs_real_t s2[6],
1435  const cs_real_t s3[6],
1436  cs_real_t sout[restrict 3][3])
1437 {
1438  cs_real_t _sout[3][3];
1439 
1440  /* S11 */
1441  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1442  /* S22 */
1443  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1444  /* S33 */
1445  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1446  /* S12 */
1447  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1448  /* S21 */
1449  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1450  /* S23 */
1451  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1452  /* S32 */
1453  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1454  /* S13 */
1455  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1456  /* S31 */
1457  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1458 
1459  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1460  /* S22 */
1461  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1462  /* S33 */
1463  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1464  /* S12 */
1465  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1466  /* S21 */
1467  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1468  /* S23 */
1469  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1470  /* S32 */
1471  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1472  /* S13 */
1473  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1474  /* S31 */
1475  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1476 }
1477 
1478 /*----------------------------------------------------------------------------*/
1485 /*----------------------------------------------------------------------------*/
1486 
1487 static inline void
1488 cs_nvec3(const cs_real_3_t v,
1489  cs_nvec3_t *qv)
1490 {
1491  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1492 
1493  qv->meas = magnitude;
1494  if (fabs(magnitude) > cs_math_zero_threshold) {
1495 
1496  const cs_real_t inv = 1/magnitude;
1497  qv->unitv[0] = inv * v[0];
1498  qv->unitv[1] = inv * v[1];
1499  qv->unitv[2] = inv * v[2];
1500 
1501  }
1502  else
1503  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1504 }
1505 
1506 /*=============================================================================
1507  * Public function prototypes
1508  *============================================================================*/
1509 
1510 /*----------------------------------------------------------------------------*/
1511 /*
1512  * \brief Compute the length (Euclidean norm) between two points xa and xb in
1513  * a Cartesian coordinate system of dimension 3
1514  *
1515  * \param[in] xa coordinate of the first extremity
1516  * \param[in] xb coordinate of the second extremity
1517  * \param[out] len pointer to the length of the vector va -> vb
1518  * \param[out] unitv unitary vector along xa -> xb
1519  */
1520 /*----------------------------------------------------------------------------*/
1521 
1522 void
1523 cs_math_3_length_unitv(const cs_real_t xa[3],
1524  const cs_real_t xb[3],
1525  cs_real_t *len,
1526  cs_real_3_t unitv);
1527 
1528 /*----------------------------------------------------------------------------*/
1540 /*----------------------------------------------------------------------------*/
1541 
1542 void
1543 cs_math_sym_33_eigen(const cs_real_t m[6],
1544  cs_real_t eig_vals[3]);
1545 
1546 /*----------------------------------------------------------------------------*/
1559 /*----------------------------------------------------------------------------*/
1560 
1561 void
1562 cs_math_33_eigen(const cs_real_t m[3][3],
1563  cs_real_t *eig_ratio,
1564  cs_real_t *eig_max);
1565 
1566 /*----------------------------------------------------------------------------*/
1577 /*----------------------------------------------------------------------------*/
1578 
1579 double
1580 cs_math_surftri(const cs_real_t xv[3],
1581  const cs_real_t xe[3],
1582  const cs_real_t xf[3]);
1583 
1584 /*----------------------------------------------------------------------------*/
1596 /*----------------------------------------------------------------------------*/
1597 
1598 double
1599 cs_math_voltet(const cs_real_t xv[3],
1600  const cs_real_t xe[3],
1601  const cs_real_t xf[3],
1602  const cs_real_t xc[3]);
1603 
1604 /*----------------------------------------------------------------------------*/
1617 /*----------------------------------------------------------------------------*/
1618 
1619 void
1620 cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1621  const cs_real_t tol_err,
1622  cs_real_t eig_val[restrict 3],
1623  cs_real_t eig_vec[restrict 3][3]);
1624 
1625 /*----------------------------------------------------------------------------*/
1635 /*----------------------------------------------------------------------------*/
1636 
1637 void
1638 cs_math_fact_lu(cs_lnum_t n_blocks,
1639  const int b_size,
1640  const cs_real_t *a,
1641  cs_real_t *a_lu);
1642 
1643 /*----------------------------------------------------------------------------*/
1653 /*----------------------------------------------------------------------------*/
1654 
1655 void
1656 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1657  const int n,
1658  cs_real_t x[],
1659  const cs_real_t b[]);
1660 
1661 /*----------------------------------------------------------------------------*/
1671 /*----------------------------------------------------------------------------*/
1672 
1673 void
1674 cs_math_sym_44_factor_ldlt(cs_real_t ldlt[10]);
1675 
1676 /*----------------------------------------------------------------------------*/
1689 /*----------------------------------------------------------------------------*/
1690 
1691 cs_real_t
1692 cs_math_sym_44_partial_solve_ldlt(const cs_real_t ldlt[10],
1693  const cs_real_t rhs[4]);
1694 
1695 /*----------------------------------------------------------------------------*/
1696 
1697 END_C_DECLS
1698 
1699 #endif /* __CS_MATH_H__ */
bft_error.h
cs_defs.h
restrict
#define restrict
Definition: cs_defs.h:139
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:514
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:319
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:334
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:336
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:515
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:343
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:313
t
@ t
Definition: cs_field_pointer.h:92
k
@ k
Definition: cs_field_pointer.h:70
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_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:814
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:349
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:403
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_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:544
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:1433
cs_math_1ov6
const cs_real_t cs_math_1ov6
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:566
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_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:655
cs_math_clamp
static cs_real_t cs_math_clamp(cs_real_t x, cs_real_t xmin, cs_real_t xmax)
Clamp function for a given scalar value.
Definition: cs_math.h:202
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:304
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:677
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:734
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,...
Definition: cs_math.h:886
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:436
cs_math_33_trace
static cs_real_t cs_math_33_trace(const cs_real_t t[3][3])
Compute the trace of a 3x3 tensor.
Definition: cs_math.h:718
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
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:440
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:456
cs_math_4ov3
const cs_real_t cs_math_4ov3
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:1172
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:238
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:1318
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's rule.
Definition: cs_math.h:955
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:230
cs_math_2ov3
const cs_real_t cs_math_2ov3
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:499
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:1216
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:680
cs_math_1ov12
const cs_real_t cs_math_1ov12
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:254
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:622
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:1352
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:1257
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:633
cs_math_pow5
static cs_real_t cs_math_pow5(cs_real_t x)
Compute the 5-th power of a real value.
Definition: cs_math.h:286
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's rule.
Definition: cs_math.h:990
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:316
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:923
cs_math_1ov24
const cs_real_t cs_math_1ov24
cs_math_sym_44_partial_solve_ldlt
cs_real_t cs_math_sym_44_partial_solve_ldlt(const cs_real_t ldlt[10], const cs_real_t rhs[4])
LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http...
Definition: cs_math.c:784
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:371
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:864
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:838
cs_math_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition: cs_math.h:61
ZZ
@ ZZ
Definition: cs_math.h:65
XY
@ XY
Definition: cs_math.h:66
XZ
@ XZ
Definition: cs_math.h:68
YZ
@ YZ
Definition: cs_math.h:67
YY
@ YY
Definition: cs_math.h:64
XX
@ XX
Definition: cs_math.h:63
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:523
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's rule.
Definition: cs_math.h:1025
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:591
cs_math_1ov3
const cs_real_t cs_math_1ov3
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:329
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:1131
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:612
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:1087
cs_math_5ov3
const cs_real_t cs_math_5ov3
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:222
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:773
cs_math_epzero
const cs_real_t cs_math_epzero
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:420
cs_math_big_r
const cs_real_t cs_math_big_r
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:466
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:498
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_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:1488
cs_math_sym_44_factor_ldlt
void cs_math_sym_44_factor_ldlt(cs_real_t ldlt[10])
LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http...
Definition: cs_math.c:725
cs_math_33_main_invariant_2
static cs_real_t cs_math_33_main_invariant_2(const cs_real_t m[3][3])
Compute the second main invariant of the symmetric part of a 3x3 tensor.
Definition: cs_math.h:1296
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:1381
cs_math_pi
const cs_real_t cs_math_pi
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:751
cs_math_33_identity
static const cs_real_33_t cs_math_33_identity
Definition: cs_math.h:92
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:270
cs_math_sym_33_identity
static const cs_real_6_t cs_math_sym_33_identity
Definition: cs_math.h:95
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:475
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:794
cs_math_sym_33_sym_33_product_trace
static cs_real_t cs_math_sym_33_sym_33_product_trace(const cs_real_t m1[6], const cs_real_t m2[6])
Compute the product of two symmetric matrices of 3x3 real values and take the trace....
Definition: cs_math.h:699
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:1058
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:392
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
atimbr::u
double precision, dimension(:,:,:), allocatable u
Definition: atimbr.f90:113
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
numvar::ik
integer, save ik
turbulent kinetic energy
Definition: numvar.f90:75
cs_nvec3_t
Definition: cs_defs.h:372
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:374
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:375