8.2
general documentation
cs_math.h
Go to the documentation of this file.
1 #ifndef __CS_MATH_H__
2 #define __CS_MATH_H__
3 
4 /*============================================================================
5  * Mathematical base functions.
6  *============================================================================*/
7 
8 /*
9  This file is part of code_saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-2024 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 #if defined(__NVCC__) && defined(__CUDA_ARCH__)
40 #include <float.h>
41 #endif
42 
43 /*----------------------------------------------------------------------------
44  * Local headers
45  *----------------------------------------------------------------------------*/
46 
47 #if defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
48 #include "bft_error.h"
49 #endif
50 
51 /*----------------------------------------------------------------------------*/
52 
53 BEGIN_C_DECLS
54 
55 /*=============================================================================
56  * Local Macro definitions
57  *============================================================================*/
58 
59 /*============================================================================
60  * Type definition
61  *============================================================================*/
62 
63 /* Symmetric tensor component name */
64 
65 typedef enum {
66 
67  XX,
68  YY,
69  ZZ,
70  XY,
71  YZ,
72  XZ
73 
74 } cs_math_sym_tensor_component_t;
75 
76 /*============================================================================
77  * Global variables
78  *============================================================================*/
79 
80 /* Numerical constants */
81 
82 #if defined(__NVCC__) && defined(__CUDA_ARCH__)
83 
84 /* On GPU, global variables are usually not accessible. */
85 
86 #define cs_math_zero_threshold FLT_MIN
87 #define cs_math_epzero 1e-12
88 #define cs_math_infinite_r 1.e30
89 #define cs_math_big_r 1.e12
90 #define cs_math_pi 3.14159265358979323846
91 
92 #else
93 
94 /* General constants accessible on CPU */
95 
96 extern const cs_real_t cs_math_zero_threshold;
97 extern const cs_real_t cs_math_1ov3;
98 extern const cs_real_t cs_math_2ov3;
99 extern const cs_real_t cs_math_4ov3;
100 extern const cs_real_t cs_math_5ov3;
101 extern const cs_real_t cs_math_1ov6;
102 extern const cs_real_t cs_math_1ov12;
103 extern const cs_real_t cs_math_1ov24;
104 extern const cs_real_t cs_math_epzero;
105 extern const cs_real_t cs_math_infinite_r;
106 extern const cs_real_t cs_math_big_r;
107 extern const cs_real_t cs_math_pi;
108 
109 /* Identity matrix in dimension 3 */
110 static const cs_real_33_t cs_math_33_identity = {{1., 0., 0.,},
111  {0., 1., 0.},
112  {0., 0., 1.}};
113 static const cs_real_6_t cs_math_sym_33_identity = {1., 1., 1., 0. ,0., 0.};
114 
115 #endif
116 
117 /*=============================================================================
118  * Inline static functions
119  *============================================================================*/
120 
121 /*----------------------------------------------------------------------------*/
130 /*----------------------------------------------------------------------------*/
131 
132 static inline int
133 cs_math_binom(int n,
134  int k)
135 {
136  int ret = 1;
137  assert(n >= k);
138 
139  const int n_iter = (k > n-k) ? n-k : k;
140  for (int j = 1; j <= n_iter; j++, n--) {
141  if (n % j == 0)
142  ret *= n/j;
143  else if (ret % j == 0)
144  ret = ret/j*n;
145  else
146  ret = (ret*n)/j;
147  }
148 
149  return ret;
150 }
151 
152 /*----------------------------------------------------------------------------*/
160 /*----------------------------------------------------------------------------*/
161 
162 CS_F_HOST_DEVICE static inline cs_real_t
163 cs_math_fabs(cs_real_t x)
164 {
165  cs_real_t ret = (x < 0) ? -x : x;
166 
167  return ret;
168 }
169 
170 /*----------------------------------------------------------------------------*/
178 /*----------------------------------------------------------------------------*/
179 
180 CS_F_HOST_DEVICE static inline cs_real_t
181 cs_math_fmin(cs_real_t x,
182  cs_real_t y)
183 {
184  cs_real_t ret = (x < y) ? x : y;
185 
186  return ret;
187 }
188 
189 /*----------------------------------------------------------------------------*/
197 /*----------------------------------------------------------------------------*/
198 
199 CS_F_HOST_DEVICE static inline cs_real_t
200 cs_math_fmax(cs_real_t x,
201  cs_real_t y)
202 {
203  cs_real_t ret = (x < y) ? y : x;
204 
205  return ret;
206 }
207 
208 /*----------------------------------------------------------------------------*/
219 /*----------------------------------------------------------------------------*/
220 
221 CS_F_HOST_DEVICE static inline cs_real_t
222 cs_math_clamp(cs_real_t x,
223  cs_real_t xmin,
224  cs_real_t xmax)
225 {
226  cs_real_t ret = cs_math_fmin(xmax, cs_math_fmax(xmin, x));
227 
228  return ret;
229 }
230 
231 /*----------------------------------------------------------------------------*/
239 /*----------------------------------------------------------------------------*/
240 
241 CS_F_HOST_DEVICE static inline cs_real_t
242 cs_math_sq(cs_real_t x)
243 {
244  return x*x;
245 }
246 
247 /*----------------------------------------------------------------------------*/
255 /*----------------------------------------------------------------------------*/
256 
257 CS_F_HOST_DEVICE static inline cs_real_t
258 cs_math_pow2(cs_real_t x)
259 {
260  return x*x;
261 }
262 
263 /*----------------------------------------------------------------------------*/
271 /*----------------------------------------------------------------------------*/
272 
273 CS_F_HOST_DEVICE static inline cs_real_t
274 cs_math_pow3(cs_real_t x)
275 {
276  return x*x*x;
277 }
278 
279 /*----------------------------------------------------------------------------*/
287 /*----------------------------------------------------------------------------*/
288 
289 CS_F_HOST_DEVICE static inline cs_real_t
290 cs_math_pow4(cs_real_t x)
291 {
292  return (x*x)*(x*x);
293 }
294 
295 /*----------------------------------------------------------------------------*/
303 /*----------------------------------------------------------------------------*/
304 
305 CS_F_HOST_DEVICE static inline cs_real_t
306 cs_math_pow5(cs_real_t x)
307 {
308  return x*(x*x)*(x*x);
309 }
310 
311 /*----------------------------------------------------------------------------*/
321 /*----------------------------------------------------------------------------*/
322 
323 CS_F_HOST_DEVICE static inline cs_real_t
324 cs_math_3_distance(const cs_real_t xa[3],
325  const cs_real_t xb[3])
326 {
327  cs_real_t v[3];
328 
329  v[0] = xb[0] - xa[0];
330  v[1] = xb[1] - xa[1];
331  v[2] = xb[2] - xa[2];
332 
333  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
334 }
335 
336 /*----------------------------------------------------------------------------*/
346 /*----------------------------------------------------------------------------*/
347 
348 CS_F_HOST_DEVICE static inline cs_real_t
349 cs_math_3_distance_dot_product(const cs_real_t xa[3],
350  const cs_real_t xb[3],
351  const cs_real_t xc[3])
352 {
353  return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
354 }
355 
356 /*----------------------------------------------------------------------------*/
366 /*----------------------------------------------------------------------------*/
367 
368 CS_F_HOST_DEVICE static inline cs_real_t
369 cs_math_3_square_distance(const cs_real_t xa[3],
370  const cs_real_t xb[3])
371 {
372  cs_real_t v[3] = {xb[0] - xa[0],
373  xb[1] - xa[1],
374  xb[2] - xa[2]};
375 
376  return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
377 }
378 
379 /*----------------------------------------------------------------------------*/
388 /*----------------------------------------------------------------------------*/
389 
390 CS_F_HOST_DEVICE static inline cs_real_t
391 cs_math_3_dot_product(const cs_real_t u[3],
392  const cs_real_t v[3])
393 {
394  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
395 
396  return uv;
397 }
398 
399 /*----------------------------------------------------------------------------*/
409 /*----------------------------------------------------------------------------*/
410 
411 CS_F_HOST_DEVICE static inline cs_real_t
412 cs_math_3_33_3_dot_product(const cs_real_t n1[3],
413  const cs_real_t t[3][3],
414  const cs_real_t n2[3])
415 {
416  cs_real_t n_t_n
417  = ( n1[0]*t[0][0]*n2[0] + n1[1]*t[1][0]*n2[0] + n1[2]*t[2][0]*n2[0]
418  + n1[0]*t[0][1]*n2[1] + n1[1]*t[1][1]*n2[1] + n1[2]*t[2][1]*n2[1]
419  + n1[0]*t[0][2]*n2[2] + n1[1]*t[1][2]*n2[2] + n1[2]*t[2][2]*n2[2]);
420  return n_t_n;
421 }
422 
423 /*----------------------------------------------------------------------------*/
437 /*----------------------------------------------------------------------------*/
438 
439 CS_F_HOST_DEVICE static inline cs_real_t
440 cs_math_3_sym_33_3_dot_product(const cs_real_t n1[3],
441  const cs_real_t t[6],
442  const cs_real_t n2[3])
443 {
444  return ( n1[0] * (t[0]*n2[0] + t[3]*n2[1] + t[5]*n2[2])
445  + n1[1] * (t[3]*n2[0] + t[1]*n2[1] + t[4]*n2[2])
446  + n1[2] * (t[5]*n2[0] + t[4]*n2[1] + t[2]*n2[2]));
447 }
448 
449 /*----------------------------------------------------------------------------*/
457 /*----------------------------------------------------------------------------*/
458 
459 CS_F_HOST_DEVICE static inline cs_real_t
460 cs_math_3_norm(const cs_real_t v[3])
461 {
462  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
463 }
464 
465 /*----------------------------------------------------------------------------*/
473 /*----------------------------------------------------------------------------*/
474 
475 CS_F_HOST_DEVICE static inline cs_real_t
476 cs_math_3_square_norm(const cs_real_t v[3])
477 {
478  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
479 
480  return v2;
481 }
482 
483 /*----------------------------------------------------------------------------*/
492 /*----------------------------------------------------------------------------*/
493 
494 static inline void
495 cs_math_3_normalise(const cs_real_t vin[3],
496  cs_real_t vout[3])
497 {
498  cs_real_t norm = cs_math_3_norm(vin);
499 
500  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
501 
502  vout[0] = inv_norm * vin[0];
503  vout[1] = inv_norm * vin[1];
504  vout[2] = inv_norm * vin[2];
505 }
506 
507 /*----------------------------------------------------------------------------*/
516 /*----------------------------------------------------------------------------*/
517 
518 CS_F_HOST_DEVICE static inline void
519 cs_math_3_normalize(const cs_real_t vin[3],
520  cs_real_t vout[3])
521 {
522  cs_real_t norm = cs_math_3_norm(vin);
523 
524  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
525 
526  vout[0] = inv_norm * vin[0];
527  vout[1] = inv_norm * vin[1];
528  vout[2] = inv_norm * vin[2];
529 }
530 
531 /*----------------------------------------------------------------------------*/
542 /*----------------------------------------------------------------------------*/
543 
544 CS_F_HOST_DEVICE static inline void
545 cs_math_3_normalize_threshold(const cs_real_t vin[3],
546  const cs_real_t thres,
547  cs_real_t vout[3])
548 {
549  cs_real_t norm = cs_math_3_norm(vin);
550 
551  cs_real_t inv_norm = ((norm > thres) ? 1. / norm : 1. / thres);
552 
553  vout[0] = inv_norm * vin[0];
554  vout[1] = inv_norm * vin[1];
555  vout[2] = inv_norm * vin[2];
556 }
557 
558 /*----------------------------------------------------------------------------*/
567 /*----------------------------------------------------------------------------*/
568 
569 CS_F_HOST_DEVICE static inline void
570 cs_math_3_orthogonal_projection(const cs_real_t n[3],
571  const cs_real_t v[3],
572  cs_real_t vout[restrict 3])
573 {
574  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
575  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
576  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
577 }
578 
579 /*----------------------------------------------------------------------------*/
588 /*----------------------------------------------------------------------------*/
589 
590 CS_F_HOST_DEVICE static inline void
591 cs_math_3_normal_scaling(const cs_real_t n[3],
592  cs_real_t factor,
593  cs_real_t v[3])
594 {
595  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
596  for (int i = 0; i < 3; i++)
597  v[i] += v_dot_n * n[i];
598 }
599 
600 /*----------------------------------------------------------------------------*/
610 /*----------------------------------------------------------------------------*/
611 
612 CS_F_HOST_DEVICE static inline void
613 cs_math_33_normal_scaling_add(const cs_real_t n[3],
614  cs_real_t factor,
615  cs_real_t t[3][3])
616 {
617  cs_real_t n_t_n = (factor -1.) *
618  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
619  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
620  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
621  for (int i = 0; i < 3; i++) {
622  for (int j = 0; j < 3; j++)
623  t[i][j] += n_t_n * n[i] * n[j];
624  }
625 }
626 /*----------------------------------------------------------------------------*/
635 /*----------------------------------------------------------------------------*/
636 
637 CS_F_HOST_DEVICE static inline void
638 cs_math_33_3_product(const cs_real_t m[3][3],
639  const cs_real_t v[3],
640  cs_real_t mv[restrict 3])
641 {
642  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
643  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
644  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
645 }
646 
647 /*----------------------------------------------------------------------------*/
656 /*----------------------------------------------------------------------------*/
657 
658 CS_F_HOST_DEVICE static inline void
659 cs_math_33_3_product_add(const cs_real_t m[3][3],
660  const cs_real_t v[3],
661  cs_real_t mv[restrict 3])
662 {
663  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
664  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
665  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
666 }
667 
668 /*----------------------------------------------------------------------------*/
677 /*----------------------------------------------------------------------------*/
678 
679 CS_F_HOST_DEVICE static inline void
680 cs_math_33t_3_product(const cs_real_t m[3][3],
681  const cs_real_t v[3],
682  cs_real_t mv[restrict 3])
683 {
684  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
685  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
686  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
687 }
688 
689 /*----------------------------------------------------------------------------*/
699 /*----------------------------------------------------------------------------*/
700 
701 CS_F_HOST_DEVICE static inline void
702 cs_math_sym_33_3_product(const cs_real_t m[6],
703  const cs_real_t v[3],
704  cs_real_t mv[restrict 3])
705 {
706  mv[0] = m[0]*v[0] + m[3]*v[1] + m[5]*v[2];
707  mv[1] = m[3]*v[0] + m[1]*v[1] + m[4]*v[2];
708  mv[2] = m[5]*v[0] + m[4]*v[1] + m[2]*v[2];
709 }
710 
711 /*----------------------------------------------------------------------------*/
721 /*----------------------------------------------------------------------------*/
722 
723 CS_F_HOST_DEVICE static inline void
724 cs_math_sym_33_3_product_add(const cs_real_t m[6],
725  const cs_real_t v[3],
726  cs_real_t mv[restrict 3])
727 {
728  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
729  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
730  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
731 }
732 
733 /*----------------------------------------------------------------------------*/
743 /*----------------------------------------------------------------------------*/
744 
745 CS_F_HOST_DEVICE static inline cs_real_t
746 cs_math_sym_33_sym_33_product_trace(const cs_real_t m1[6],
747  const cs_real_t m2[6])
748 {
749  return m1[0]*m2[0] + 2.*m1[3]*m2[3] + 2.*m1[5]*m2[5]
750  + m1[1]*m2[1] + 2.*m1[4]*m2[4]
751  + m1[2]*m2[2];
752 }
753 
754 /*----------------------------------------------------------------------------*/
762 /*----------------------------------------------------------------------------*/
763 
764 CS_F_HOST_DEVICE static inline cs_real_t
765 cs_math_33_trace(const cs_real_t t[3][3])
766 {
767  return (t[0][0] + t[1][1] + t[2][2]);
768 }
769 
770 /*----------------------------------------------------------------------------*/
778 /*----------------------------------------------------------------------------*/
779 
780 CS_F_HOST_DEVICE static inline cs_real_t
781 cs_math_6_trace(const cs_real_t t[6])
782 {
783  return (t[0] + t[1] + t[2]);
784 }
785 
786 /*----------------------------------------------------------------------------*/
795 /*----------------------------------------------------------------------------*/
796 
797 CS_F_HOST_DEVICE static inline void
798 cs_math_66_6_product(const cs_real_t m[6][6],
799  const cs_real_t v[6],
800  cs_real_t mv[restrict 6])
801 {
802  for (int i = 0; i < 6; i++) {
803  for (int j = 0; j < 6; j++)
804  mv[i] = m[i][j] * v[j];
805  }
806 }
807 
808 /*----------------------------------------------------------------------------*/
817 /*----------------------------------------------------------------------------*/
818 
819 CS_F_HOST_DEVICE static inline void
820 cs_math_66_6_product_add(const cs_real_t m[6][6],
821  const cs_real_t v[6],
822  cs_real_t mv[restrict 6])
823 {
824  for (int i = 0; i < 6; i++) {
825  for (int j = 0; j < 6; j++)
826  mv[i] += m[i][j] * v[j];
827  }
828 }
829 
830 /*----------------------------------------------------------------------------*/
838 /*----------------------------------------------------------------------------*/
839 
840 CS_F_HOST_DEVICE static inline cs_real_t
841 cs_math_33_determinant(const cs_real_t m[3][3])
842 {
843  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
844  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
845  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
846 
847  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
848 }
849 
850 /*----------------------------------------------------------------------------*/
858 /*----------------------------------------------------------------------------*/
859 
860 CS_F_HOST_DEVICE static inline cs_real_t
861 cs_math_sym_33_determinant(const cs_real_6_t m)
862 {
863  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
864  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
865  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
866 
867  return m[0]*com0 + m[3]*com1 + m[5]*com2;
868 }
869 
870 /*----------------------------------------------------------------------------*/
878 /*----------------------------------------------------------------------------*/
879 
880 #if defined(__INTEL_COMPILER)
881 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
882 #endif
883 
884 CS_F_HOST_DEVICE static inline void
885 cs_math_3_cross_product(const cs_real_t u[3],
886  const cs_real_t v[3],
887  cs_real_t uv[restrict 3])
888 {
889  uv[0] = u[1]*v[2] - u[2]*v[1];
890  uv[1] = u[2]*v[0] - u[0]*v[2];
891  uv[2] = u[0]*v[1] - u[1]*v[0];
892 }
893 
894 /*----------------------------------------------------------------------------*/
904 /*----------------------------------------------------------------------------*/
905 
906 #if defined(__INTEL_COMPILER)
907 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
908 #endif
909 
910 CS_F_HOST_DEVICE static inline cs_real_t
911 cs_math_3_triple_product(const cs_real_t u[3],
912  const cs_real_t v[3],
913  const cs_real_t w[3])
914 {
915  return (u[1]*v[2] - u[2]*v[1]) * w[0]
916  + (u[2]*v[0] - u[0]*v[2]) * w[1]
917  + (u[0]*v[1] - u[1]*v[0]) * w[2];
918 }
919 
920 /*----------------------------------------------------------------------------*/
930 /*----------------------------------------------------------------------------*/
931 
932 CS_F_HOST_DEVICE static inline void
933 cs_math_3_orthonormal_basis(const cs_real_t vect[3],
934  cs_real_t axes[3][3])
935 {
936  assert(cs_math_3_norm(vect) > cs_math_zero_threshold);
937 
938  // Compute normal - third axis
939  cs_math_3_normalize(vect, axes[2]);
940 
941  // Compute base's first axis
942  // First test projection of Ox
943  cs_real_t Ox[3] = {1., 0., 0.};
944  cs_real_t w[3] = {0.};
945 
946  cs_math_3_orthogonal_projection(axes[2], Ox, w);
947 
948  // If Ox projection is null, project Oy
949  if (cs_math_3_norm(w) < cs_math_zero_threshold) {
950  cs_real_t Oy[3] = {0., 1., 0.};
951  cs_math_3_orthogonal_projection(axes[2], Oy, w);
952  }
953 
954  cs_math_3_normalize(w, axes[0]);
955 
956  // Compute base's second axis using cross product
957  cs_math_3_cross_product(axes[2], axes[0], axes[1]);
958 }
959 
960 /*----------------------------------------------------------------------------*/
967 /*----------------------------------------------------------------------------*/
968 
969 CS_F_HOST_DEVICE static inline void
970 cs_math_33_inv_cramer(const cs_real_t in[3][3],
971  cs_real_t out[3][3])
972 {
973  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
974  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
975  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
976 
977  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
978  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
979  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
980 
981  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
982  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
983  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
984 
985  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
986  const double invdet = 1./det;
987 
988  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
989  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
990  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
991 }
992 
993 /*----------------------------------------------------------------------------*/
999 /*----------------------------------------------------------------------------*/
1000 
1001 CS_F_HOST_DEVICE static inline void
1002 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
1003 {
1004  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
1005  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
1006  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
1007  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
1008  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
1009  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
1010  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
1011  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
1012  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
1013 
1014  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
1015 
1016  a[0][0] = a00 * det_inv;
1017  a[0][1] = a01 * det_inv;
1018  a[0][2] = a02 * det_inv;
1019  a[1][0] = a10 * det_inv;
1020  a[1][1] = a11 * det_inv;
1021  a[1][2] = a12 * det_inv;
1022  a[2][0] = a20 * det_inv;
1023  a[2][1] = a21 * det_inv;
1024  a[2][2] = a22 * det_inv;
1025 }
1026 
1027 /*----------------------------------------------------------------------------*/
1034 /*----------------------------------------------------------------------------*/
1035 
1036 CS_F_HOST_DEVICE static inline void
1037 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
1038 {
1039  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
1040  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
1041  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
1042  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
1043  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
1044  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
1045 
1046  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
1047 
1048  a[0][0] = a00 * det_inv;
1049  a[0][1] = a01 * det_inv;
1050  a[0][2] = a02 * det_inv;
1051  a[1][0] = a01 * det_inv;
1052  a[1][1] = a11 * det_inv;
1053  a[1][2] = a12 * det_inv;
1054  a[2][0] = a02 * det_inv;
1055  a[2][1] = a12 * det_inv;
1056  a[2][2] = a22 * det_inv;
1057 }
1058 
1059 /*----------------------------------------------------------------------------*/
1069 /*----------------------------------------------------------------------------*/
1070 
1071 CS_F_HOST_DEVICE static inline void
1072 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
1073  cs_real_t sout[restrict 6])
1074 {
1075  double detinv;
1076 
1077  sout[0] = s[1]*s[2] - s[4]*s[4];
1078  sout[1] = s[0]*s[2] - s[5]*s[5];
1079  sout[2] = s[0]*s[1] - s[3]*s[3];
1080  sout[3] = s[4]*s[5] - s[3]*s[2];
1081  sout[4] = s[3]*s[5] - s[0]*s[4];
1082  sout[5] = s[3]*s[4] - s[1]*s[5];
1083 
1084  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
1085 
1086  sout[0] *= detinv;
1087  sout[1] *= detinv;
1088  sout[2] *= detinv;
1089  sout[3] *= detinv;
1090  sout[4] *= detinv;
1091  sout[5] *= detinv;
1092 }
1093 
1094 /*----------------------------------------------------------------------------*/
1102 /*----------------------------------------------------------------------------*/
1103 
1104 CS_F_HOST_DEVICE static inline void
1105 cs_math_33_product(const cs_real_t m1[3][3],
1106  const cs_real_t m2[3][3],
1107  cs_real_t mout[3][3])
1108 {
1109  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1110  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1111  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1112 
1113  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1114  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1115  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1116 
1117  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1118  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1119  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1120 }
1121 
1122 /*----------------------------------------------------------------------------*/
1131 /*----------------------------------------------------------------------------*/
1132 
1133 CS_F_HOST_DEVICE static inline void
1134 cs_math_33_transform_r_to_a(const cs_real_t m[3][3],
1135  const cs_real_t q[3][3],
1136  cs_real_t mout[3][3])
1137 {
1138  /* _m = M.Q */
1139  cs_real_33_t _m;
1140  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[1][0] + m[0][2]*q[2][0];
1141  _m[0][1] = m[0][0]*q[0][1] + m[0][1]*q[1][1] + m[0][2]*q[2][1];
1142  _m[0][2] = m[0][0]*q[0][2] + m[0][1]*q[1][2] + m[0][2]*q[2][2];
1143 
1144  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[1][0] + m[1][2]*q[2][0];
1145  _m[1][1] = m[1][0]*q[0][1] + m[1][1]*q[1][1] + m[1][2]*q[2][1];
1146  _m[1][2] = m[1][0]*q[0][2] + m[1][1]*q[1][2] + m[1][2]*q[2][2];
1147 
1148  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[1][0] + m[2][2]*q[2][0];
1149  _m[2][1] = m[2][0]*q[0][1] + m[2][1]*q[1][1] + m[2][2]*q[2][1];
1150  _m[2][2] = m[2][0]*q[0][2] + m[2][1]*q[1][2] + m[2][2]*q[2][2];
1151 
1152  /* mout = Q^t _m */
1153  mout[0][0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1154  mout[0][1] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1155  mout[0][2] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1156 
1157  mout[1][0] = q[0][1]*_m[0][0] + q[1][1]*_m[1][0] + q[2][1]*_m[2][0];
1158  mout[1][1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1159  mout[1][2] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1160 
1161  mout[2][0] = q[0][2]*_m[0][0] + q[1][2]*_m[1][0] + q[2][2]*_m[2][0];
1162  mout[2][1] = q[0][2]*_m[0][1] + q[1][2]*_m[1][1] + q[2][2]*_m[2][1];
1163  mout[2][2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1164 }
1165 
1166 /*----------------------------------------------------------------------------*/
1175 /*----------------------------------------------------------------------------*/
1176 
1177 CS_F_HOST_DEVICE static inline void
1178 cs_math_sym_33_transform_r_to_a(const cs_real_t m[6],
1179  const cs_real_t q[3][3],
1180  cs_real_t mout[6])
1181 {
1182  /* _m = M.Q */
1183  cs_real_33_t _m;
1184  _m[0][0] = m[0]*q[0][0] + m[3]*q[1][0] + m[5]*q[2][0];
1185  _m[0][1] = m[0]*q[0][1] + m[3]*q[1][1] + m[5]*q[2][1];
1186  _m[0][2] = m[0]*q[0][2] + m[3]*q[1][2] + m[5]*q[2][2];
1187 
1188  _m[1][0] = m[3]*q[0][0] + m[1]*q[1][0] + m[4]*q[2][0];
1189  _m[1][1] = m[3]*q[0][1] + m[1]*q[1][1] + m[4]*q[2][1];
1190  _m[1][2] = m[3]*q[0][2] + m[1]*q[1][2] + m[4]*q[2][2];
1191 
1192  _m[2][0] = m[5]*q[0][0] + m[4]*q[1][0] + m[2]*q[2][0];
1193  _m[2][1] = m[5]*q[0][1] + m[4]*q[1][1] + m[2]*q[2][1];
1194  _m[2][2] = m[5]*q[0][2] + m[4]*q[1][2] + m[2]*q[2][2];
1195 
1196  /* mout = Q^t _m */
1197  mout[0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1198  mout[1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1199  mout[2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1200 
1201  mout[3] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1202  mout[4] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1203  mout[5] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1204 }
1205 
1206 /*----------------------------------------------------------------------------*/
1215 /*----------------------------------------------------------------------------*/
1216 
1217 CS_F_HOST_DEVICE static inline void
1218 cs_math_33_transform_a_to_r(const cs_real_t m[3][3],
1219  const cs_real_t q[3][3],
1220  cs_real_t mout[3][3])
1221 {
1222  /* _m = M.Q^t */
1223  cs_real_33_t _m;
1224  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[0][1] + m[0][2]*q[0][2];
1225  _m[0][1] = m[0][0]*q[1][0] + m[0][1]*q[1][1] + m[0][2]*q[1][2];
1226  _m[0][2] = m[0][0]*q[2][0] + m[0][1]*q[2][1] + m[0][2]*q[2][2];
1227 
1228  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[0][1] + m[1][2]*q[0][2];
1229  _m[1][1] = m[1][0]*q[1][0] + m[1][1]*q[1][1] + m[1][2]*q[1][2];
1230  _m[1][2] = m[1][0]*q[2][0] + m[1][1]*q[2][1] + m[1][2]*q[2][2];
1231 
1232  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[0][1] + m[2][2]*q[0][2];
1233  _m[2][1] = m[2][0]*q[1][0] + m[2][1]*q[1][1] + m[2][2]*q[1][2];
1234  _m[2][2] = m[2][0]*q[2][0] + m[2][1]*q[2][1] + m[2][2]*q[2][2];
1235 
1236  /* mout = Q _m */
1237  mout[0][0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1238  mout[0][1] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1239  mout[0][2] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1240 
1241  mout[1][0] = q[1][0]*_m[0][0] + q[1][1]*_m[1][0] + q[1][2]*_m[2][0];
1242  mout[1][1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1243  mout[1][2] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1244 
1245  mout[2][0] = q[2][0]*_m[0][0] + q[2][1]*_m[1][0] + q[2][2]*_m[2][0];
1246  mout[2][1] = q[2][0]*_m[0][1] + q[2][1]*_m[1][1] + q[2][2]*_m[2][1];
1247  mout[2][2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1248 }
1249 
1250 /*----------------------------------------------------------------------------*/
1259 /*----------------------------------------------------------------------------*/
1260 
1261 CS_F_HOST_DEVICE static inline void
1262 cs_math_sym_33_transform_a_to_r(const cs_real_t m[6],
1263  const cs_real_t q[3][3],
1264  cs_real_t mout[6])
1265 {
1266  /* _m = M.Q^t */
1267  cs_real_33_t _m;
1268  _m[0][0] = m[0]*q[0][0] + m[3]*q[0][1] + m[5]*q[0][2];
1269  _m[0][1] = m[0]*q[1][0] + m[3]*q[1][1] + m[5]*q[1][2];
1270  _m[0][2] = m[0]*q[2][0] + m[3]*q[2][1] + m[5]*q[2][2];
1271 
1272  _m[1][0] = m[3]*q[0][0] + m[1]*q[0][1] + m[4]*q[0][2];
1273  _m[1][1] = m[3]*q[1][0] + m[1]*q[1][1] + m[4]*q[1][2];
1274  _m[1][2] = m[3]*q[2][0] + m[1]*q[2][1] + m[4]*q[2][2];
1275 
1276  _m[2][0] = m[5]*q[0][0] + m[4]*q[0][1] + m[2]*q[0][2];
1277  _m[2][1] = m[5]*q[1][0] + m[4]*q[1][1] + m[2]*q[1][2];
1278  _m[2][2] = m[5]*q[2][0] + m[4]*q[2][1] + m[2]*q[2][2];
1279 
1280  /* mout = Q _m */
1281  mout[0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1282  mout[1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1283  mout[2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1284 
1285 
1286  mout[3] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1287  mout[4] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1288  mout[5] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1289 }
1290 
1291 /*----------------------------------------------------------------------------*/
1300 /*----------------------------------------------------------------------------*/
1301 
1302 CS_F_HOST_DEVICE static inline void
1303 cs_math_33_extract_sym_ant(const cs_real_t m[3][3],
1304  cs_real_t m_sym[3][3],
1305  cs_real_t m_ant[3][3])
1306 {
1307  /* sym = 0.5 (m + m_transpose) */
1308  m_sym[0][0] = 0.5 * (m[0][0] + m[0][0]);
1309  m_sym[0][1] = 0.5 * (m[0][1] + m[1][0]);
1310  m_sym[0][2] = 0.5 * (m[0][2] + m[2][0]);
1311  m_sym[1][0] = 0.5 * (m[1][0] + m[0][1]);
1312  m_sym[1][1] = 0.5 * (m[1][1] + m[1][1]);
1313  m_sym[1][2] = 0.5 * (m[1][2] + m[2][1]);
1314  m_sym[2][0] = 0.5 * (m[2][0] + m[0][2]);
1315  m_sym[2][1] = 0.5 * (m[2][1] + m[1][2]);
1316  m_sym[2][2] = 0.5 * (m[2][2] + m[2][2]);
1317 
1318  /* ant = 0.5 (m - m_transpose) */
1319  m_ant[0][0] = 0.5 * (m[0][0] - m[0][0]);
1320  m_ant[0][1] = 0.5 * (m[0][1] - m[1][0]);
1321  m_ant[0][2] = 0.5 * (m[0][2] - m[2][0]);
1322  m_ant[1][0] = 0.5 * (m[1][0] - m[0][1]);
1323  m_ant[1][1] = 0.5 * (m[1][1] - m[1][1]);
1324  m_ant[1][2] = 0.5 * (m[1][2] - m[2][1]);
1325  m_ant[2][0] = 0.5 * (m[2][0] - m[0][2]);
1326  m_ant[2][1] = 0.5 * (m[2][1] - m[1][2]);
1327  m_ant[2][2] = 0.5 * (m[2][2] - m[2][2]);
1328 }
1329 
1330 /*----------------------------------------------------------------------------*/
1339 /*----------------------------------------------------------------------------*/
1340 
1341 CS_F_HOST_DEVICE static inline cs_real_t
1342 cs_math_33_main_invariant_2(const cs_real_t m[3][3])
1343 {
1344  /* sym = 0.5 (m + m_transpose) */
1345  return cs_math_pow2(m[0][0])
1346  + cs_math_pow2(m[1][1])
1347  + cs_math_pow2(m[2][2])
1348  + 0.5 * ( cs_math_pow2(m[0][1] + m[1][0])
1349  + cs_math_pow2(m[0][2] + m[2][0])
1350  + cs_math_pow2(m[1][2] + m[2][1]));
1351 }
1352 
1353 /*----------------------------------------------------------------------------*/
1361 /*----------------------------------------------------------------------------*/
1362 
1363 CS_F_HOST_DEVICE static inline void
1364 cs_math_33_product_add(const cs_real_t m1[3][3],
1365  const cs_real_t m2[3][3],
1366  cs_real_t mout[restrict 3][3])
1367 {
1368  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1369  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1370  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1371 
1372  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1373  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1374  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1375 
1376  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1377  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1378  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1379 }
1380 
1381 /*----------------------------------------------------------------------------*/
1395 /*----------------------------------------------------------------------------*/
1396 
1397 CS_F_HOST_DEVICE static inline void
1398 cs_math_sym_33_product(const cs_real_t s1[6],
1399  const cs_real_t s2[6],
1400  cs_real_t sout[restrict 6])
1401 {
1402  /* S11 */
1403  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1404  /* S22 */
1405  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1406  /* S33 */
1407  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1408  /* S12 = S21 */
1409  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1410  /* S23 = S32 */
1411  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1412  /* S13 = S31 */
1413  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1414 }
1415 
1416 /*----------------------------------------------------------------------------*/
1424 /*----------------------------------------------------------------------------*/
1425 
1426 CS_F_HOST_DEVICE static inline void
1427 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1428  cs_real_t sout[restrict 6][6])
1429 {
1430  const int t2v[3][3] = {{0, 3, 5},
1431  {3, 1, 4},
1432  {5, 4, 2}};
1433 
1434  const int iv2t[6] = {0, 1, 2, 0, 1, 0};
1435  const int jv2t[6] = {0, 1, 2, 1, 2, 2};
1436 
1437  for (int i = 0; i < 6; i++) {
1438  for (int j = 0; j < 6; j++)
1439  sout[i][j] = 0;
1440  }
1441 
1442  /* Consider : W = s*R + R*s^t .
1443  * W_ij = Sum_{k<3} [s_ik*r_jk + s_jk*r_ik]
1444  * We look for A_(ij,pq) such as A*R = W
1445  *
1446  * so
1447  * A_(ij,jk) takes s_ik
1448  * and
1449  * A_(ij,ik) takes s_jk
1450  */
1451  for (int ij = 0; ij < 6; ij++) {
1452  int i = iv2t[ij];
1453  int j = jv2t[ij];
1454  for (int k = 0; k < 3; k++) {
1455  int ik = t2v[i][k];
1456  int jk = t2v[j][k];
1457 
1458  sout[ij][ik] += s[j][k];
1459  sout[ij][jk] += s[i][k];
1460  }
1461  }
1462 }
1463 
1464 /*----------------------------------------------------------------------------*/
1476 /*----------------------------------------------------------------------------*/
1477 
1478 CS_F_HOST_DEVICE static inline void
1479 cs_math_sym_33_double_product(const cs_real_t s1[6],
1480  const cs_real_t s2[6],
1481  const cs_real_t s3[6],
1482  cs_real_t sout[restrict 3][3])
1483 {
1484  cs_real_t _sout[3][3];
1485 
1486  /* S11 */
1487  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1488  /* S22 */
1489  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1490  /* S33 */
1491  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1492  /* S12 */
1493  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1494  /* S21 */
1495  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1496  /* S23 */
1497  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1498  /* S32 */
1499  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1500  /* S13 */
1501  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1502  /* S31 */
1503  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1504 
1505  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1506  /* S22 */
1507  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1508  /* S33 */
1509  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1510  /* S12 */
1511  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1512  /* S21 */
1513  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1514  /* S23 */
1515  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1516  /* S32 */
1517  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1518  /* S13 */
1519  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1520  /* S31 */
1521  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1522 }
1523 
1524 /*----------------------------------------------------------------------------*/
1531 /*----------------------------------------------------------------------------*/
1532 
1533 CS_F_HOST_DEVICE static inline void
1534 cs_nvec3(const cs_real_3_t v,
1535  cs_nvec3_t *qv)
1536 {
1537  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1538 
1539  qv->meas = magnitude;
1540  if (fabs(magnitude) > cs_math_zero_threshold) {
1541 
1542  const cs_real_t inv = 1/magnitude;
1543  qv->unitv[0] = inv * v[0];
1544  qv->unitv[1] = inv * v[1];
1545  qv->unitv[2] = inv * v[2];
1546 
1547  }
1548  else
1549  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1550 }
1551 
1552 /*=============================================================================
1553  * Public function prototypes
1554  *============================================================================*/
1555 
1556 /*----------------------------------------------------------------------------*/
1557 /*
1558  * \brief Compute the length (Euclidean norm) between two points xa and xb in
1559  * a Cartesian coordinate system of dimension 3
1560  *
1561  * \param[in] xa coordinate of the first extremity
1562  * \param[in] xb coordinate of the second extremity
1563  * \param[out] len pointer to the length of the vector va -> vb
1564  * \param[out] unitv unitary vector along xa -> xb
1565  */
1566 /*----------------------------------------------------------------------------*/
1567 
1568 CS_F_HOST_DEVICE void
1569 cs_math_3_length_unitv(const cs_real_t xa[3],
1570  const cs_real_t xb[3],
1571  cs_real_t *len,
1572  cs_real_3_t unitv);
1573 
1574 /*----------------------------------------------------------------------------*/
1586 /*----------------------------------------------------------------------------*/
1587 
1588 void
1589 cs_math_sym_33_eigen(const cs_real_t m[6],
1590  cs_real_t eig_vals[3]);
1591 
1592 /*----------------------------------------------------------------------------*/
1605 /*----------------------------------------------------------------------------*/
1606 
1607 void
1608 cs_math_33_eigen(const cs_real_t m[3][3],
1609  cs_real_t *eig_ratio,
1610  cs_real_t *eig_max);
1611 
1612 /*----------------------------------------------------------------------------*/
1623 /*----------------------------------------------------------------------------*/
1624 
1625 double
1626 cs_math_surftri(const cs_real_t xv[3],
1627  const cs_real_t xe[3],
1628  const cs_real_t xf[3]);
1629 
1630 /*----------------------------------------------------------------------------*/
1642 /*----------------------------------------------------------------------------*/
1643 
1644 double
1645 cs_math_voltet(const cs_real_t xv[3],
1646  const cs_real_t xe[3],
1647  const cs_real_t xf[3],
1648  const cs_real_t xc[3]);
1649 
1650 /*----------------------------------------------------------------------------*/
1663 /*----------------------------------------------------------------------------*/
1664 
1665 void
1666 cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1667  const cs_real_t tol_err,
1668  cs_real_t eig_val[restrict 3],
1669  cs_real_t eig_vec[restrict 3][3]);
1670 
1671 /*----------------------------------------------------------------------------*/
1681 /*----------------------------------------------------------------------------*/
1682 
1683 void
1684 cs_math_fact_lu(cs_lnum_t n_blocks,
1685  const int b_size,
1686  const cs_real_t *a,
1687  cs_real_t *a_lu);
1688 
1689 /*----------------------------------------------------------------------------*/
1699 /*----------------------------------------------------------------------------*/
1700 
1701 void
1702 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1703  const int n,
1704  cs_real_t x[],
1705  const cs_real_t b[]);
1706 
1707 /*----------------------------------------------------------------------------*/
1717 /*----------------------------------------------------------------------------*/
1718 
1719 void
1720 cs_math_sym_44_factor_ldlt(cs_real_t ldlt[10]);
1721 
1722 /*----------------------------------------------------------------------------*/
1735 /*----------------------------------------------------------------------------*/
1736 
1737 cs_real_t
1738 cs_math_sym_44_partial_solve_ldlt(const cs_real_t ldlt[10],
1739  const cs_real_t rhs[4]);
1740 
1741 /*----------------------------------------------------------------------------*/
1742 
1743 END_C_DECLS
1744 
1745 #endif /* __CS_MATH_H__ */
bft_error.h
cs_defs.h
restrict
#define restrict
Definition: cs_defs.h:141
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:528
CS_F_HOST_DEVICE
#define CS_F_HOST_DEVICE
Definition: cs_defs.h:546
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:347
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:349
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:529
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:356
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
t
@ t
Definition: cs_field_pointer.h:92
k
@ k
Definition: cs_field_pointer.h:70
cs_math_3_normalize
static CS_F_HOST_DEVICE 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:519
cs_math_1ov6
const cs_real_t cs_math_1ov6
cs_math_pow3
static CS_F_HOST_DEVICE cs_real_t cs_math_pow3(cs_real_t x)
Compute the cube of a real value.
Definition: cs_math.h:274
cs_math_3_sym_33_3_dot_product
static CS_F_HOST_DEVICE 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:440
cs_math_33t_3_product
static CS_F_HOST_DEVICE 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:680
cs_math_33_extract_sym_ant
static CS_F_HOST_DEVICE 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:1303
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.cpp:446
cs_math_sym_33_inv_cramer
static CS_F_HOST_DEVICE 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:1072
cs_math_33_inv_cramer_sym_in_place
static CS_F_HOST_DEVICE 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:1037
cs_math_33_3_product_add
static CS_F_HOST_DEVICE 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:659
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
cs_math_sym_33_double_product
static CS_F_HOST_DEVICE 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:1479
cs_math_pow2
static CS_F_HOST_DEVICE cs_real_t cs_math_pow2(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:258
cs_math_4ov3
const cs_real_t cs_math_4ov3
cs_math_3_cross_product
static CS_F_HOST_DEVICE 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:885
cs_math_3_square_distance
static CS_F_HOST_DEVICE 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:369
cs_math_sym_33_transform_r_to_a
static CS_F_HOST_DEVICE 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:1178
cs_math_33_transform_r_to_a
static CS_F_HOST_DEVICE 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:1134
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.cpp:234
cs_math_33_main_invariant_2
static CS_F_HOST_DEVICE 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:1342
cs_math_2ov3
const cs_real_t cs_math_2ov3
cs_math_pow4
static CS_F_HOST_DEVICE cs_real_t cs_math_pow4(cs_real_t x)
Compute the 4-th power of a real value.
Definition: cs_math.h:290
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.cpp:692
cs_math_1ov12
const cs_real_t cs_math_1ov12
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.cpp:634
cs_math_3_distance_dot_product
static CS_F_HOST_DEVICE 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:349
cs_math_3_square_norm
static CS_F_HOST_DEVICE 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:476
cs_math_33_normal_scaling_add
static CS_F_HOST_DEVICE 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:613
cs_math_3_normalize_threshold
static CS_F_HOST_DEVICE void cs_math_3_normalize_threshold(const cs_real_t vin[3], const cs_real_t thres, cs_real_t vout[3])
Normalise a vector of 3 real values and clip the norm using a threshold value.
Definition: cs_math.h:545
cs_math_66_6_product_add
static CS_F_HOST_DEVICE 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:820
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.cpp:323
cs_math_sym_33_product
static CS_F_HOST_DEVICE 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:1398
cs_math_1ov24
const cs_real_t cs_math_1ov24
cs_math_3_normal_scaling
static CS_F_HOST_DEVICE 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:591
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.cpp:796
cs_math_66_6_product
static CS_F_HOST_DEVICE 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:798
cs_math_33_inv_cramer_in_place
static CS_F_HOST_DEVICE 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:1002
cs_math_33_inv_cramer
static CS_F_HOST_DEVICE 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:970
cs_math_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition: cs_math.h:65
ZZ
@ ZZ
Definition: cs_math.h:69
XY
@ XY
Definition: cs_math.h:70
XZ
@ XZ
Definition: cs_math.h:72
YZ
@ YZ
Definition: cs_math.h:71
YY
@ YY
Definition: cs_math.h:68
XX
@ XX
Definition: cs_math.h:67
cs_math_3_length_unitv
CS_F_HOST_DEVICE 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.cpp:413
cs_math_33_3_product
static CS_F_HOST_DEVICE 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:638
cs_math_33_product_add
static CS_F_HOST_DEVICE 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:1364
cs_math_3_triple_product
static CS_F_HOST_DEVICE 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:911
cs_math_3_33_3_dot_product
static CS_F_HOST_DEVICE 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:412
cs_math_3_distance
static CS_F_HOST_DEVICE 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:324
cs_math_pow5
static CS_F_HOST_DEVICE cs_real_t cs_math_pow5(cs_real_t x)
Compute the 5-th power of a real value.
Definition: cs_math.h:306
cs_math_1ov3
const cs_real_t cs_math_1ov3
cs_math_sym_33_transform_a_to_r
static CS_F_HOST_DEVICE 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:1262
cs_math_sym_33_determinant
static CS_F_HOST_DEVICE 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:861
cs_math_3_dot_product
static CS_F_HOST_DEVICE 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:391
cs_math_5ov3
const cs_real_t cs_math_5ov3
cs_math_33_determinant
static CS_F_HOST_DEVICE 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:841
cs_math_sym_33_3_product_add
static CS_F_HOST_DEVICE 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:724
cs_math_epzero
const cs_real_t cs_math_epzero
cs_math_reduce_sym_prod_33_to_66
static CS_F_HOST_DEVICE 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:1427
cs_math_33_product
static CS_F_HOST_DEVICE 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:1105
cs_math_big_r
const cs_real_t cs_math_big_r
cs_math_33_transform_a_to_r
static CS_F_HOST_DEVICE 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:1218
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.cpp:476
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.cpp:510
cs_math_binom
static int cs_math_binom(int n, int k)
Computes the binomial coefficient of n and k.
Definition: cs_math.h:133
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.cpp:737
cs_nvec3
static CS_F_HOST_DEVICE 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:1534
cs_math_6_trace
static CS_F_HOST_DEVICE cs_real_t cs_math_6_trace(const cs_real_t t[6])
Compute the trace of a symmetric tensor.
Definition: cs_math.h:781
cs_math_sym_33_sym_33_product_trace
static CS_F_HOST_DEVICE 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:746
cs_math_pi
const cs_real_t cs_math_pi
cs_math_clamp
static CS_F_HOST_DEVICE 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:222
cs_math_33_identity
static const cs_real_33_t cs_math_33_identity
Definition: cs_math.h:110
cs_math_sym_33_identity
static const cs_real_6_t cs_math_sym_33_identity
Definition: cs_math.h:113
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:495
cs_math_sym_33_3_product
static CS_F_HOST_DEVICE 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:702
cs_math_3_norm
static CS_F_HOST_DEVICE 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:460
cs_math_fabs
static CS_F_HOST_DEVICE cs_real_t cs_math_fabs(cs_real_t x)
Compute the absolute value of a real value.
Definition: cs_math.h:163
cs_math_33_trace
static CS_F_HOST_DEVICE 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:765
cs_math_fmin
static CS_F_HOST_DEVICE 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:181
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
cs_math_sq
static CS_F_HOST_DEVICE cs_real_t cs_math_sq(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:242
cs_math_3_orthogonal_projection
static CS_F_HOST_DEVICE 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:570
cs_math_3_orthonormal_basis
static CS_F_HOST_DEVICE 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:933
cs_math_fmax
static CS_F_HOST_DEVICE 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:200
atimbr::u
double precision, dimension(:,:,:), allocatable u
Definition: atimbr.f90:112
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:113
numvar::ik
integer, save ik
turbulent kinetic energy
Definition: numvar.f90:75
cs_nvec3_t
Definition: cs_defs.h:386
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:388
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:389