8.3
general documentation
cs_math.h
Go to the documentation of this file.
1#ifndef __CS_MATH_H__
2#define __CS_MATH_H__
3
4/*============================================================================
5 * Mathematical base functions.
6 *============================================================================*/
7
8/*
9 This file is part of code_saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-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 || defined(SYCL_LANGUAGE_VERSION) \
41 || defined(HAVE_OPENMP_TARGET)
42#include <float.h>
43#endif
44
45/*----------------------------------------------------------------------------
46 * Local headers
47 *----------------------------------------------------------------------------*/
48
49#if defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
50#include "bft_error.h"
51#endif
52
53/*----------------------------------------------------------------------------*/
54
55BEGIN_C_DECLS
56
57/*=============================================================================
58 * Local Macro definitions
59 *============================================================================*/
60
61/*============================================================================
62 * Type definition
63 *============================================================================*/
64
65/* Symmetric tensor component name */
66
67typedef enum {
68
69 XX,
70 YY,
71 ZZ,
72 XY,
73 YZ,
74 XZ
75
76} cs_math_sym_tensor_component_t;
77
78/*============================================================================
79 * Global variables
80 *============================================================================*/
81
82/* Numerical constants */
83
84#if (defined(__NVCC__) && defined(__CUDA_ARCH__)) \
85 || defined(SYCL_LANGUAGE_VERSION) \
86 || defined(HAVE_OPENMP_TARGET)
87
88/* On GPU, global variables are usually not accessible. */
89
90#define cs_math_zero_threshold FLT_MIN
91#define cs_math_epzero 1e-12
92#define cs_math_infinite_r 1.e30
93#define cs_math_big_r 1.e12
94#define cs_math_pi 3.14159265358979323846
95
96#else
97
98/* General constants accessible on CPU */
99
100extern const cs_real_t cs_math_zero_threshold;
101extern const cs_real_t cs_math_epzero;
102extern const cs_real_t cs_math_infinite_r;
103extern const cs_real_t cs_math_big_r;
104extern const cs_real_t cs_math_pi;
105
106#endif
107
108#if !(defined(__NVCC__) && defined(__CUDA_ARCH__))
109
110extern const cs_real_t cs_math_1ov3;
111extern const cs_real_t cs_math_2ov3;
112extern const cs_real_t cs_math_4ov3;
113extern const cs_real_t cs_math_5ov3;
114extern const cs_real_t cs_math_1ov6;
115extern const cs_real_t cs_math_1ov12;
116extern const cs_real_t cs_math_1ov24;
117
118/* Identity matrix in dimension 3 */
119static const cs_real_33_t cs_math_33_identity = {{1., 0., 0.,},
120 {0., 1., 0.},
121 {0., 0., 1.}};
122static const cs_real_6_t cs_math_sym_33_identity = {1., 1., 1., 0. ,0., 0.};
123
124#endif
125
126/*----------------------------------------------------------------------------*/
127
128END_C_DECLS
129
130/*=============================================================================
131 * Templated inline functions
132 *============================================================================*/
133
134#if defined(__cplusplus)
135
136/*----------------------------------------------------------------------------*/
149/*----------------------------------------------------------------------------*/
150
151template <typename T, typename U>
152CS_F_HOST_DEVICE inline cs_real_t
153cs_math_3_distance_dot_product(const T xa[3],
154 const T xb[3],
155 const U xc[3])
156{
157 return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
158}
159
160/*----------------------------------------------------------------------------*/
172/*----------------------------------------------------------------------------*/
173
174template <typename T, typename U>
175CS_F_HOST_DEVICE inline cs_real_t
176cs_math_3_dot_product(const T u[3],
177 const U v[3])
178{
179 double uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
180
181 return uv;
182}
183
184/*----------------------------------------------------------------------------*/
194/*----------------------------------------------------------------------------*/
195
196template <typename T>
197CS_F_HOST_DEVICE inline cs_real_t
198cs_math_3_norm(const T v[3])
199{
200 return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
201}
202
203/*----------------------------------------------------------------------------*/
221/*----------------------------------------------------------------------------*/
222
223template <typename T, typename U, typename V>
224CS_F_HOST_DEVICE inline cs_real_t
225cs_math_3_sym_33_3_dot_product(const T n1[3],
226 const U t[6],
227 const V n2[3])
228{
229 return ( n1[0] * (t[0]*n2[0] + t[3]*n2[1] + t[5]*n2[2])
230 + n1[1] * (t[3]*n2[0] + t[1]*n2[1] + t[4]*n2[2])
231 + n1[2] * (t[5]*n2[0] + t[4]*n2[1] + t[2]*n2[2]));
232}
233
234/*----------------------------------------------------------------------------*/
244/*----------------------------------------------------------------------------*/
245
246template <typename T>
247CS_F_HOST_DEVICE inline cs_real_t
248cs_math_3_square_norm(const T v[3])
249{
250 cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
251
252 return v2;
253}
254
255/*----------------------------------------------------------------------------*/
267/*----------------------------------------------------------------------------*/
268
269template <typename T, typename U>
270CS_F_HOST_DEVICE static inline void
271cs_math_3_normalize(const T vin[3],
272 U vout[3])
273{
274 cs_real_t norm = cs_math_3_norm(vin);
275
276 cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
277
278 vout[0] = inv_norm * vin[0];
279 vout[1] = inv_norm * vin[1];
280 vout[2] = inv_norm * vin[2];
281}
282
283/*----------------------------------------------------------------------------*/
297/*----------------------------------------------------------------------------*/
298
299template <typename T, typename U, typename V>
300CS_F_HOST_DEVICE inline void
301cs_math_sym_33_3_product(const T m[6],
302 const U v[3],
303 V *restrict mv)
304{
305 mv[0] = m[0]*v[0] + m[3]*v[1] + m[5]*v[2];
306 mv[1] = m[3]*v[0] + m[1]*v[1] + m[4]*v[2];
307 mv[2] = m[5]*v[0] + m[4]*v[1] + m[2]*v[2];
308}
309
310/*----------------------------------------------------------------------------*/
323/*----------------------------------------------------------------------------*/
324
325template <typename T, typename U, typename V>
326CS_F_HOST_DEVICE inline void
327cs_math_3_normal_scaling(const T n[3],
328 U factor,
329 V v[3])
330{
331 cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
332 for (int i = 0; i < 3; i++)
333 v[i] += v_dot_n * n[i];
334}
335
336/*----------------------------------------------------------------------------*/
350/*----------------------------------------------------------------------------*/
351
352template <typename T, typename U, typename V>
353CS_F_HOST_DEVICE inline cs_real_t
354cs_math_3_33_3_dot_product(const T n1[3],
355 const U t[3][3],
356 const V n2[3])
357{
358 cs_real_t n_t_n
359 = ( n1[0]*t[0][0]*n2[0] + n1[1]*t[1][0]*n2[0] + n1[2]*t[2][0]*n2[0]
360 + n1[0]*t[0][1]*n2[1] + n1[1]*t[1][1]*n2[1] + n1[2]*t[2][1]*n2[1]
361 + n1[0]*t[0][2]*n2[2] + n1[1]*t[1][2]*n2[2] + n1[2]*t[2][2]*n2[2]);
362 return n_t_n;
363}
364
365/*----------------------------------------------------------------------------*/
378/*----------------------------------------------------------------------------*/
379
380template <typename T, typename U, typename V>
381CS_F_HOST_DEVICE inline void
382cs_math_3_orthogonal_projection(const T n[3],
383 const U v[3],
384 V *restrict vout)
385{
386 vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
387 vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
388 vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
389}
390
391/*----------------------------------------------------------------------------*/
403/*----------------------------------------------------------------------------*/
404
405template <typename T, typename U, typename V>
406CS_F_HOST_DEVICE static inline void
407cs_math_3_cross_product(const T u[3],
408 const U v[3],
409 V *restrict uv)
410{
411 uv[0] = u[1]*v[2] - u[2]*v[1];
412 uv[1] = u[2]*v[0] - u[0]*v[2];
413 uv[2] = u[0]*v[1] - u[1]*v[0];
414}
415
416#endif // defined(__cplusplus)
417
418BEGIN_C_DECLS
419
420/*=============================================================================
421 * Inline static functions
422 *============================================================================*/
423
424/*----------------------------------------------------------------------------*/
433/*----------------------------------------------------------------------------*/
434
435static inline int
436cs_math_binom(int n,
437 int k)
438{
439 int ret = 1;
440 assert(n >= k);
441
442 const int n_iter = (k > n-k) ? n-k : k;
443 for (int j = 1; j <= n_iter; j++, n--) {
444 if (n % j == 0)
445 ret *= n/j;
446 else if (ret % j == 0)
447 ret = ret/j*n;
448 else
449 ret = (ret*n)/j;
450 }
451
452 return ret;
453}
454
455/*----------------------------------------------------------------------------*/
463/*----------------------------------------------------------------------------*/
464
465CS_F_HOST_DEVICE static inline cs_real_t
466cs_math_fabs(cs_real_t x)
467{
468 cs_real_t ret = (x < 0) ? -x : x;
469
470 return ret;
471}
472
473/*----------------------------------------------------------------------------*/
481/*----------------------------------------------------------------------------*/
482
483CS_F_HOST_DEVICE static inline cs_real_t
484cs_math_fmin(cs_real_t x,
485 cs_real_t y)
486{
487 cs_real_t ret = (x < y) ? x : y;
488
489 return ret;
490}
491
492/*----------------------------------------------------------------------------*/
500/*----------------------------------------------------------------------------*/
501
502CS_F_HOST_DEVICE static inline cs_real_t
503cs_math_fmax(cs_real_t x,
504 cs_real_t y)
505{
506 cs_real_t ret = (x < y) ? y : x;
507
508 return ret;
509}
510
511/*----------------------------------------------------------------------------*/
522/*----------------------------------------------------------------------------*/
523
524CS_F_HOST_DEVICE static inline cs_real_t
525cs_math_clamp(cs_real_t x,
526 cs_real_t xmin,
527 cs_real_t xmax)
528{
529 cs_real_t ret = cs_math_fmin(xmax, cs_math_fmax(xmin, x));
530
531 return ret;
532}
533
534/*----------------------------------------------------------------------------*/
542/*----------------------------------------------------------------------------*/
543
544CS_F_HOST_DEVICE static inline cs_real_t
545cs_math_sq(cs_real_t x)
546{
547 return x*x;
548}
549
550/*----------------------------------------------------------------------------*/
558/*----------------------------------------------------------------------------*/
559
560CS_F_HOST_DEVICE static inline cs_real_t
561cs_math_pow2(cs_real_t x)
562{
563 return x*x;
564}
565
566/*----------------------------------------------------------------------------*/
574/*----------------------------------------------------------------------------*/
575
576CS_F_HOST_DEVICE static inline cs_real_t
577cs_math_pow3(cs_real_t x)
578{
579 return x*x*x;
580}
581
582/*----------------------------------------------------------------------------*/
590/*----------------------------------------------------------------------------*/
591
592CS_F_HOST_DEVICE static inline cs_real_t
593cs_math_pow4(cs_real_t x)
594{
595 cs_real_t x2 = x * x;
596 return x2 * x2;
597}
598
599/*----------------------------------------------------------------------------*/
607/*----------------------------------------------------------------------------*/
608
609CS_F_HOST_DEVICE static inline cs_real_t
610cs_math_pow5(cs_real_t x)
611{
612 cs_real_t x2 = x * x;
613 return x * x2 * x2;
614}
615
616/*----------------------------------------------------------------------------*/
626/*----------------------------------------------------------------------------*/
627
628CS_F_HOST_DEVICE static inline cs_real_t
629cs_math_3_distance(const cs_real_t xa[3],
630 const cs_real_t xb[3])
631{
632 cs_real_t v[3];
633
634 v[0] = xb[0] - xa[0];
635 v[1] = xb[1] - xa[1];
636 v[2] = xb[2] - xa[2];
637
638 return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
639}
640
641/*----------------------------------------------------------------------------*/
651/*----------------------------------------------------------------------------*/
652
653CS_F_HOST_DEVICE static inline cs_real_t
654cs_math_3_distance_dot_product(const cs_real_t xa[3],
655 const cs_real_t xb[3],
656 const cs_real_t xc[3])
657{
658 return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
659}
660
661/*----------------------------------------------------------------------------*/
671/*----------------------------------------------------------------------------*/
672
673CS_F_HOST_DEVICE static inline cs_real_t
674cs_math_3_square_distance(const cs_real_t xa[3],
675 const cs_real_t xb[3])
676{
677 cs_real_t v[3] = {xb[0] - xa[0],
678 xb[1] - xa[1],
679 xb[2] - xa[2]};
680
681 return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
682}
683
684/*----------------------------------------------------------------------------*/
693/*----------------------------------------------------------------------------*/
694
695CS_F_HOST_DEVICE static inline cs_real_t
696cs_math_3_dot_product(const cs_real_t u[3],
697 const cs_real_t v[3])
698{
699 cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
700
701 return uv;
702}
703
704/*----------------------------------------------------------------------------*/
714/*----------------------------------------------------------------------------*/
715
716CS_F_HOST_DEVICE static inline cs_real_t
717cs_math_3_33_3_dot_product(const cs_real_t n1[3],
718 const cs_real_t t[3][3],
719 const cs_real_t n2[3])
720{
721 cs_real_t n_t_n
722 = ( n1[0]*t[0][0]*n2[0] + n1[1]*t[1][0]*n2[0] + n1[2]*t[2][0]*n2[0]
723 + n1[0]*t[0][1]*n2[1] + n1[1]*t[1][1]*n2[1] + n1[2]*t[2][1]*n2[1]
724 + n1[0]*t[0][2]*n2[2] + n1[1]*t[1][2]*n2[2] + n1[2]*t[2][2]*n2[2]);
725 return n_t_n;
726}
727
728/*----------------------------------------------------------------------------*/
742/*----------------------------------------------------------------------------*/
743
744CS_F_HOST_DEVICE static inline cs_real_t
745cs_math_3_sym_33_3_dot_product(const cs_real_t n1[3],
746 const cs_real_t t[6],
747 const cs_real_t n2[3])
748{
749 return ( n1[0] * (t[0]*n2[0] + t[3]*n2[1] + t[5]*n2[2])
750 + n1[1] * (t[3]*n2[0] + t[1]*n2[1] + t[4]*n2[2])
751 + n1[2] * (t[5]*n2[0] + t[4]*n2[1] + t[2]*n2[2]));
752}
753
754/*----------------------------------------------------------------------------*/
762/*----------------------------------------------------------------------------*/
763
764CS_F_HOST_DEVICE static inline cs_real_t
765cs_math_3_norm(const cs_real_t v[3])
766{
767 return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
768}
769
770/*----------------------------------------------------------------------------*/
778/*----------------------------------------------------------------------------*/
779
780CS_F_HOST_DEVICE static inline cs_real_t
781cs_math_3_square_norm(const cs_real_t v[3])
782{
783 cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
784
785 return v2;
786}
787
788/*----------------------------------------------------------------------------*/
797/*----------------------------------------------------------------------------*/
798
799CS_F_HOST_DEVICE static inline void
800cs_math_3_normalize(const cs_real_t vin[3],
801 cs_real_t vout[3])
802{
803 cs_real_t norm = cs_math_3_norm(vin);
804
805 cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
806
807 vout[0] = inv_norm * vin[0];
808 vout[1] = inv_norm * vin[1];
809 vout[2] = inv_norm * vin[2];
810}
811
812/*----------------------------------------------------------------------------*/
823/*----------------------------------------------------------------------------*/
824
825CS_F_HOST_DEVICE static inline void
826cs_math_3_normalize_threshold(const cs_real_t vin[3],
827 const cs_real_t thres,
828 cs_real_t vout[3])
829{
830 cs_real_t norm = cs_math_3_norm(vin);
831
832 cs_real_t inv_norm = ((norm > thres) ? 1. / norm : 1. / thres);
833
834 vout[0] = inv_norm * vin[0];
835 vout[1] = inv_norm * vin[1];
836 vout[2] = inv_norm * vin[2];
837}
838
839/*----------------------------------------------------------------------------*/
848/*----------------------------------------------------------------------------*/
849
850CS_F_HOST_DEVICE static inline void
851cs_math_3_orthogonal_projection(const cs_real_t n[3],
852 const cs_real_t v[3],
853 cs_real_t *restrict vout)
854{
855 vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
856 vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
857 vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
858}
859
860/*----------------------------------------------------------------------------*/
869/*----------------------------------------------------------------------------*/
870
871CS_F_HOST_DEVICE static inline void
872cs_math_3_normal_scaling(const cs_real_t n[3],
873 cs_real_t factor,
874 cs_real_t v[3])
875{
876 cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
877 for (int i = 0; i < 3; i++)
878 v[i] += v_dot_n * n[i];
879}
880
881/*----------------------------------------------------------------------------*/
891/*----------------------------------------------------------------------------*/
892
893CS_F_HOST_DEVICE static inline void
894cs_math_33_normal_scaling_add(const cs_real_t n[3],
895 cs_real_t factor,
896 cs_real_t t[3][3])
897{
898 cs_real_t n_t_n = (factor -1.) *
899 ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
900 + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
901 + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
902 for (int i = 0; i < 3; i++) {
903 for (int j = 0; j < 3; j++)
904 t[i][j] += n_t_n * n[i] * n[j];
905 }
906}
907/*----------------------------------------------------------------------------*/
916/*----------------------------------------------------------------------------*/
917
918CS_F_HOST_DEVICE static inline void
919cs_math_33_3_product(const cs_real_t m[3][3],
920 const cs_real_t v[3],
921 cs_real_t *restrict mv)
922{
923 mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
924 mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
925 mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
926}
927
928/*----------------------------------------------------------------------------*/
937/*----------------------------------------------------------------------------*/
938
939CS_F_HOST_DEVICE static inline void
940cs_math_33_3_product_add(const cs_real_t m[3][3],
941 const cs_real_t v[3],
942 cs_real_t *restrict mv)
943{
944 mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
945 mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
946 mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
947}
948
949/*----------------------------------------------------------------------------*/
958/*----------------------------------------------------------------------------*/
959
960CS_F_HOST_DEVICE static inline void
961cs_math_33t_3_product(const cs_real_t m[3][3],
962 const cs_real_t v[3],
963 cs_real_t *restrict mv)
964{
965 mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
966 mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
967 mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
968}
969
970/*----------------------------------------------------------------------------*/
980/*----------------------------------------------------------------------------*/
981
982CS_F_HOST_DEVICE static inline void
983cs_math_sym_33_3_product(const cs_real_t m[6],
984 const cs_real_t v[3],
985 cs_real_t *restrict mv)
986{
987 mv[0] = m[0]*v[0] + m[3]*v[1] + m[5]*v[2];
988 mv[1] = m[3]*v[0] + m[1]*v[1] + m[4]*v[2];
989 mv[2] = m[5]*v[0] + m[4]*v[1] + m[2]*v[2];
990}
991
992/*----------------------------------------------------------------------------*/
1002/*----------------------------------------------------------------------------*/
1003
1004CS_F_HOST_DEVICE static inline void
1005cs_math_sym_33_3_product_add(const cs_real_t m[6],
1006 const cs_real_t v[3],
1007 cs_real_t *restrict mv)
1008{
1009 mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
1010 mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
1011 mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
1012}
1013
1014/*----------------------------------------------------------------------------*/
1024/*----------------------------------------------------------------------------*/
1025
1026CS_F_HOST_DEVICE static inline cs_real_t
1027cs_math_sym_33_sym_33_product_trace(const cs_real_t m1[6],
1028 const cs_real_t m2[6])
1029{
1030 return m1[0]*m2[0] + 2.*m1[3]*m2[3] + 2.*m1[5]*m2[5]
1031 + m1[1]*m2[1] + 2.*m1[4]*m2[4]
1032 + m1[2]*m2[2];
1033}
1034
1035/*----------------------------------------------------------------------------*/
1043/*----------------------------------------------------------------------------*/
1044
1045CS_F_HOST_DEVICE static inline cs_real_t
1046cs_math_33_trace(const cs_real_t t[3][3])
1047{
1048 return (t[0][0] + t[1][1] + t[2][2]);
1049}
1050
1051/*----------------------------------------------------------------------------*/
1059/*----------------------------------------------------------------------------*/
1060
1061CS_F_HOST_DEVICE static inline cs_real_t
1062cs_math_6_trace(const cs_real_t t[6])
1063{
1064 return (t[0] + t[1] + t[2]);
1065}
1066
1067/*----------------------------------------------------------------------------*/
1076/*----------------------------------------------------------------------------*/
1077
1078CS_F_HOST_DEVICE static inline void
1079cs_math_66_6_product(const cs_real_t m[6][6],
1080 const cs_real_t v[6],
1081 cs_real_t *restrict mv)
1082{
1083 for (int i = 0; i < 6; i++) {
1084 for (int j = 0; j < 6; j++)
1085 mv[i] = m[i][j] * v[j];
1086 }
1087}
1088
1089/*----------------------------------------------------------------------------*/
1098/*----------------------------------------------------------------------------*/
1099
1100CS_F_HOST_DEVICE static inline void
1101cs_math_66_6_product_add(const cs_real_t m[6][6],
1102 const cs_real_t v[6],
1103 cs_real_t *restrict mv)
1104{
1105 for (int i = 0; i < 6; i++) {
1106 for (int j = 0; j < 6; j++)
1107 mv[i] += m[i][j] * v[j];
1108 }
1109}
1110
1111/*----------------------------------------------------------------------------*/
1119/*----------------------------------------------------------------------------*/
1120
1121CS_F_HOST_DEVICE static inline cs_real_t
1122cs_math_33_determinant(const cs_real_t m[3][3])
1123{
1124 const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
1125 const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
1126 const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
1127
1128 return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
1129}
1130
1131/*----------------------------------------------------------------------------*/
1139/*----------------------------------------------------------------------------*/
1140
1141CS_F_HOST_DEVICE static inline cs_real_t
1142cs_math_sym_33_determinant(const cs_real_t m[6])
1143{
1144 const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
1145 const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
1146 const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
1147
1148 return m[0]*com0 + m[3]*com1 + m[5]*com2;
1149}
1150
1151/*----------------------------------------------------------------------------*/
1159/*----------------------------------------------------------------------------*/
1160
1161CS_F_HOST_DEVICE static inline void
1162cs_math_3_average(const cs_real_t u[3],
1163 const cs_real_t v[3],
1164 cs_real_t *restrict uv)
1165{
1166 uv[0] = (u[0] + v[0]) / 2.0;
1167 uv[1] = (u[1] + v[1]) / 2.0;
1168 uv[2] = (u[2] + v[2]) / 2.0;
1169}
1170
1171/*----------------------------------------------------------------------------*/
1179/*----------------------------------------------------------------------------*/
1180
1181#if defined(__INTEL_COMPILER)
1182#pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
1183#endif
1184
1185CS_F_HOST_DEVICE static inline void
1186cs_math_3_cross_product(const cs_real_t u[3],
1187 const cs_real_t v[3],
1188 cs_real_t *restrict uv)
1189{
1190 uv[0] = u[1]*v[2] - u[2]*v[1];
1191 uv[1] = u[2]*v[0] - u[0]*v[2];
1192 uv[2] = u[0]*v[1] - u[1]*v[0];
1193}
1194
1195/*----------------------------------------------------------------------------*/
1205/*----------------------------------------------------------------------------*/
1206
1207#if defined(__INTEL_COMPILER)
1208#pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
1209#endif
1210
1211CS_F_HOST_DEVICE static inline cs_real_t
1212cs_math_3_triple_product(const cs_real_t u[3],
1213 const cs_real_t v[3],
1214 const cs_real_t w[3])
1215{
1216 return (u[1]*v[2] - u[2]*v[1]) * w[0]
1217 + (u[2]*v[0] - u[0]*v[2]) * w[1]
1218 + (u[0]*v[1] - u[1]*v[0]) * w[2];
1219}
1220
1221/*----------------------------------------------------------------------------*/
1231/*----------------------------------------------------------------------------*/
1232
1233CS_F_HOST_DEVICE static inline void
1234cs_math_3_orthonormal_basis(const cs_real_t vect[3],
1235 cs_real_t axes[3][3])
1236{
1237 assert(cs_math_3_norm(vect) > cs_math_zero_threshold);
1238
1239 // Compute normal - third axis
1240 cs_math_3_normalize(vect, axes[2]);
1241
1242 // Compute base's first axis
1243 // First test projection of Ox
1244 cs_real_t Ox[3] = {1., 0., 0.};
1245 cs_real_t w[3] = {0.};
1246
1247 cs_math_3_orthogonal_projection(axes[2], Ox, w);
1248
1249 // If Ox projection is null, project Oy
1250 if (cs_math_3_norm(w) < cs_math_zero_threshold) {
1251 cs_real_t Oy[3] = {0., 1., 0.};
1252 cs_math_3_orthogonal_projection(axes[2], Oy, w);
1253 }
1254
1255 cs_math_3_normalize(w, axes[0]);
1256
1257 // Compute base's second axis using cross product
1258 cs_math_3_cross_product(axes[2], axes[0], axes[1]);
1259}
1260
1261/*----------------------------------------------------------------------------*/
1268/*----------------------------------------------------------------------------*/
1269
1270CS_F_HOST_DEVICE static inline void
1271cs_math_33_inv_cramer(const cs_real_t in[3][3],
1272 cs_real_t out[3][3])
1273{
1274 out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
1275 out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
1276 out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
1277
1278 out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
1279 out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
1280 out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
1281
1282 out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
1283 out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
1284 out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
1285
1286 const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
1287 const double invdet = 1./det;
1288
1289 out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
1290 out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
1291 out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
1292}
1293
1294/*----------------------------------------------------------------------------*/
1300/*----------------------------------------------------------------------------*/
1301
1302CS_F_HOST_DEVICE static inline void
1303cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
1304{
1305 cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
1306 cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
1307 cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
1308 cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
1309 cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
1310 cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
1311 cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
1312 cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
1313 cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
1314
1315 double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
1316
1317 a[0][0] = a00 * det_inv;
1318 a[0][1] = a01 * det_inv;
1319 a[0][2] = a02 * det_inv;
1320 a[1][0] = a10 * det_inv;
1321 a[1][1] = a11 * det_inv;
1322 a[1][2] = a12 * det_inv;
1323 a[2][0] = a20 * det_inv;
1324 a[2][1] = a21 * det_inv;
1325 a[2][2] = a22 * det_inv;
1326}
1327
1328/*----------------------------------------------------------------------------*/
1335/*----------------------------------------------------------------------------*/
1336
1337CS_F_HOST_DEVICE static inline void
1338cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
1339{
1340 cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
1341 cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
1342 cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
1343 cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
1344 cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
1345 cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
1346
1347 double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
1348
1349 a[0][0] = a00 * det_inv;
1350 a[0][1] = a01 * det_inv;
1351 a[0][2] = a02 * det_inv;
1352 a[1][0] = a01 * det_inv;
1353 a[1][1] = a11 * det_inv;
1354 a[1][2] = a12 * det_inv;
1355 a[2][0] = a02 * det_inv;
1356 a[2][1] = a12 * det_inv;
1357 a[2][2] = a22 * det_inv;
1358}
1359
1360/*----------------------------------------------------------------------------*/
1370/*----------------------------------------------------------------------------*/
1371
1372CS_F_HOST_DEVICE static inline void
1373cs_math_sym_33_inv_cramer(const cs_real_t s[6],
1374 cs_real_t *restrict sout)
1375{
1376 double detinv;
1377
1378 sout[0] = s[1]*s[2] - s[4]*s[4];
1379 sout[1] = s[0]*s[2] - s[5]*s[5];
1380 sout[2] = s[0]*s[1] - s[3]*s[3];
1381 sout[3] = s[4]*s[5] - s[3]*s[2];
1382 sout[4] = s[3]*s[5] - s[0]*s[4];
1383 sout[5] = s[3]*s[4] - s[1]*s[5];
1384
1385 detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
1386
1387 sout[0] *= detinv;
1388 sout[1] *= detinv;
1389 sout[2] *= detinv;
1390 sout[3] *= detinv;
1391 sout[4] *= detinv;
1392 sout[5] *= detinv;
1393}
1394
1395/*----------------------------------------------------------------------------*/
1403/*----------------------------------------------------------------------------*/
1404
1405CS_F_HOST_DEVICE static inline void
1406cs_math_33_product(const cs_real_t m1[3][3],
1407 const cs_real_t m2[3][3],
1408 cs_real_t mout[3][3])
1409{
1410 mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1411 mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1412 mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1413
1414 mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1415 mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1416 mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1417
1418 mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1419 mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1420 mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1421}
1422
1423/*----------------------------------------------------------------------------*/
1432/*----------------------------------------------------------------------------*/
1433
1434CS_F_HOST_DEVICE static inline void
1435cs_math_33_transform_r_to_a(const cs_real_t m[3][3],
1436 const cs_real_t q[3][3],
1437 cs_real_t mout[3][3])
1438{
1439 /* _m = M.Q */
1440 cs_real_33_t _m;
1441 _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[1][0] + m[0][2]*q[2][0];
1442 _m[0][1] = m[0][0]*q[0][1] + m[0][1]*q[1][1] + m[0][2]*q[2][1];
1443 _m[0][2] = m[0][0]*q[0][2] + m[0][1]*q[1][2] + m[0][2]*q[2][2];
1444
1445 _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[1][0] + m[1][2]*q[2][0];
1446 _m[1][1] = m[1][0]*q[0][1] + m[1][1]*q[1][1] + m[1][2]*q[2][1];
1447 _m[1][2] = m[1][0]*q[0][2] + m[1][1]*q[1][2] + m[1][2]*q[2][2];
1448
1449 _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[1][0] + m[2][2]*q[2][0];
1450 _m[2][1] = m[2][0]*q[0][1] + m[2][1]*q[1][1] + m[2][2]*q[2][1];
1451 _m[2][2] = m[2][0]*q[0][2] + m[2][1]*q[1][2] + m[2][2]*q[2][2];
1452
1453 /* mout = Q^t _m */
1454 mout[0][0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1455 mout[0][1] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1456 mout[0][2] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1457
1458 mout[1][0] = q[0][1]*_m[0][0] + q[1][1]*_m[1][0] + q[2][1]*_m[2][0];
1459 mout[1][1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1460 mout[1][2] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1461
1462 mout[2][0] = q[0][2]*_m[0][0] + q[1][2]*_m[1][0] + q[2][2]*_m[2][0];
1463 mout[2][1] = q[0][2]*_m[0][1] + q[1][2]*_m[1][1] + q[2][2]*_m[2][1];
1464 mout[2][2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1465}
1466
1467/*----------------------------------------------------------------------------*/
1476/*----------------------------------------------------------------------------*/
1477
1478CS_F_HOST_DEVICE static inline void
1479cs_math_sym_33_transform_r_to_a(const cs_real_t m[6],
1480 const cs_real_t q[3][3],
1481 cs_real_t mout[6])
1482{
1483 /* _m = M.Q */
1484 cs_real_33_t _m;
1485 _m[0][0] = m[0]*q[0][0] + m[3]*q[1][0] + m[5]*q[2][0];
1486 _m[0][1] = m[0]*q[0][1] + m[3]*q[1][1] + m[5]*q[2][1];
1487 _m[0][2] = m[0]*q[0][2] + m[3]*q[1][2] + m[5]*q[2][2];
1488
1489 _m[1][0] = m[3]*q[0][0] + m[1]*q[1][0] + m[4]*q[2][0];
1490 _m[1][1] = m[3]*q[0][1] + m[1]*q[1][1] + m[4]*q[2][1];
1491 _m[1][2] = m[3]*q[0][2] + m[1]*q[1][2] + m[4]*q[2][2];
1492
1493 _m[2][0] = m[5]*q[0][0] + m[4]*q[1][0] + m[2]*q[2][0];
1494 _m[2][1] = m[5]*q[0][1] + m[4]*q[1][1] + m[2]*q[2][1];
1495 _m[2][2] = m[5]*q[0][2] + m[4]*q[1][2] + m[2]*q[2][2];
1496
1497 /* mout = Q^t _m */
1498 mout[0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
1499 mout[1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
1500 mout[2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
1501
1502 mout[3] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
1503 mout[4] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
1504 mout[5] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
1505}
1506
1507/*----------------------------------------------------------------------------*/
1516/*----------------------------------------------------------------------------*/
1517
1518CS_F_HOST_DEVICE static inline void
1519cs_math_33_transform_a_to_r(const cs_real_t m[3][3],
1520 const cs_real_t q[3][3],
1521 cs_real_t mout[3][3])
1522{
1523 /* _m = M.Q^t */
1524 cs_real_33_t _m;
1525 _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[0][1] + m[0][2]*q[0][2];
1526 _m[0][1] = m[0][0]*q[1][0] + m[0][1]*q[1][1] + m[0][2]*q[1][2];
1527 _m[0][2] = m[0][0]*q[2][0] + m[0][1]*q[2][1] + m[0][2]*q[2][2];
1528
1529 _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[0][1] + m[1][2]*q[0][2];
1530 _m[1][1] = m[1][0]*q[1][0] + m[1][1]*q[1][1] + m[1][2]*q[1][2];
1531 _m[1][2] = m[1][0]*q[2][0] + m[1][1]*q[2][1] + m[1][2]*q[2][2];
1532
1533 _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[0][1] + m[2][2]*q[0][2];
1534 _m[2][1] = m[2][0]*q[1][0] + m[2][1]*q[1][1] + m[2][2]*q[1][2];
1535 _m[2][2] = m[2][0]*q[2][0] + m[2][1]*q[2][1] + m[2][2]*q[2][2];
1536
1537 /* mout = Q _m */
1538 mout[0][0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1539 mout[0][1] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1540 mout[0][2] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1541
1542 mout[1][0] = q[1][0]*_m[0][0] + q[1][1]*_m[1][0] + q[1][2]*_m[2][0];
1543 mout[1][1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1544 mout[1][2] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1545
1546 mout[2][0] = q[2][0]*_m[0][0] + q[2][1]*_m[1][0] + q[2][2]*_m[2][0];
1547 mout[2][1] = q[2][0]*_m[0][1] + q[2][1]*_m[1][1] + q[2][2]*_m[2][1];
1548 mout[2][2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1549}
1550
1551/*----------------------------------------------------------------------------*/
1560/*----------------------------------------------------------------------------*/
1561
1562CS_F_HOST_DEVICE static inline void
1563cs_math_sym_33_transform_a_to_r(const cs_real_t m[6],
1564 const cs_real_t q[3][3],
1565 cs_real_t mout[6])
1566{
1567 /* _m = M.Q^t */
1568 cs_real_33_t _m;
1569 _m[0][0] = m[0]*q[0][0] + m[3]*q[0][1] + m[5]*q[0][2];
1570 _m[0][1] = m[0]*q[1][0] + m[3]*q[1][1] + m[5]*q[1][2];
1571 _m[0][2] = m[0]*q[2][0] + m[3]*q[2][1] + m[5]*q[2][2];
1572
1573 _m[1][0] = m[3]*q[0][0] + m[1]*q[0][1] + m[4]*q[0][2];
1574 _m[1][1] = m[3]*q[1][0] + m[1]*q[1][1] + m[4]*q[1][2];
1575 _m[1][2] = m[3]*q[2][0] + m[1]*q[2][1] + m[4]*q[2][2];
1576
1577 _m[2][0] = m[5]*q[0][0] + m[4]*q[0][1] + m[2]*q[0][2];
1578 _m[2][1] = m[5]*q[1][0] + m[4]*q[1][1] + m[2]*q[1][2];
1579 _m[2][2] = m[5]*q[2][0] + m[4]*q[2][1] + m[2]*q[2][2];
1580
1581 /* mout = Q _m */
1582 mout[0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1583 mout[1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1584 mout[2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1585
1586
1587 mout[3] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1588 mout[4] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1589 mout[5] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1590}
1591
1592/*----------------------------------------------------------------------------*/
1601/*----------------------------------------------------------------------------*/
1602
1603CS_F_HOST_DEVICE static inline void
1604cs_math_33_extract_sym_ant(const cs_real_t m[3][3],
1605 cs_real_t m_sym[3][3],
1606 cs_real_t m_ant[3][3])
1607{
1608 /* sym = 0.5 (m + m_transpose) */
1609 m_sym[0][0] = 0.5 * (m[0][0] + m[0][0]);
1610 m_sym[0][1] = 0.5 * (m[0][1] + m[1][0]);
1611 m_sym[0][2] = 0.5 * (m[0][2] + m[2][0]);
1612 m_sym[1][0] = 0.5 * (m[1][0] + m[0][1]);
1613 m_sym[1][1] = 0.5 * (m[1][1] + m[1][1]);
1614 m_sym[1][2] = 0.5 * (m[1][2] + m[2][1]);
1615 m_sym[2][0] = 0.5 * (m[2][0] + m[0][2]);
1616 m_sym[2][1] = 0.5 * (m[2][1] + m[1][2]);
1617 m_sym[2][2] = 0.5 * (m[2][2] + m[2][2]);
1618
1619 /* ant = 0.5 (m - m_transpose) */
1620 m_ant[0][0] = 0.5 * (m[0][0] - m[0][0]);
1621 m_ant[0][1] = 0.5 * (m[0][1] - m[1][0]);
1622 m_ant[0][2] = 0.5 * (m[0][2] - m[2][0]);
1623 m_ant[1][0] = 0.5 * (m[1][0] - m[0][1]);
1624 m_ant[1][1] = 0.5 * (m[1][1] - m[1][1]);
1625 m_ant[1][2] = 0.5 * (m[1][2] - m[2][1]);
1626 m_ant[2][0] = 0.5 * (m[2][0] - m[0][2]);
1627 m_ant[2][1] = 0.5 * (m[2][1] - m[1][2]);
1628 m_ant[2][2] = 0.5 * (m[2][2] - m[2][2]);
1629}
1630
1631/*----------------------------------------------------------------------------*/
1640/*----------------------------------------------------------------------------*/
1641
1642CS_F_HOST_DEVICE static inline cs_real_t
1643cs_math_33_main_invariant_2(const cs_real_t m[3][3])
1644{
1645 /* sym = 0.5 (m + m_transpose) */
1646 return cs_math_pow2(m[0][0])
1647 + cs_math_pow2(m[1][1])
1648 + cs_math_pow2(m[2][2])
1649 + 0.5 * ( cs_math_pow2(m[0][1] + m[1][0])
1650 + cs_math_pow2(m[0][2] + m[2][0])
1651 + cs_math_pow2(m[1][2] + m[2][1]));
1652}
1653
1654/*----------------------------------------------------------------------------*/
1662/*----------------------------------------------------------------------------*/
1663
1664CS_F_HOST_DEVICE static inline void
1665cs_math_33_product_add(const cs_real_t m1[3][3],
1666 const cs_real_t m2[3][3],
1667 cs_real_t (*restrict mout)[3])
1668{
1669 mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1670 mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1671 mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1672
1673 mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1674 mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1675 mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1676
1677 mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1678 mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1679 mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1680}
1681
1682/*----------------------------------------------------------------------------*/
1696/*----------------------------------------------------------------------------*/
1697
1698CS_F_HOST_DEVICE static inline void
1699cs_math_sym_33_product(const cs_real_t s1[6],
1700 const cs_real_t s2[6],
1701 cs_real_t *restrict sout)
1702{
1703 /* S11 */
1704 sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1705 /* S22 */
1706 sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1707 /* S33 */
1708 sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1709 /* S12 = S21 */
1710 sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1711 /* S23 = S32 */
1712 sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1713 /* S13 = S31 */
1714 sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1715}
1716
1717/*----------------------------------------------------------------------------*/
1725/*----------------------------------------------------------------------------*/
1726
1727CS_F_HOST_DEVICE static inline void
1728cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1729 cs_real_t (*restrict sout)[6])
1730{
1731 const int t2v[3][3] = {{0, 3, 5},
1732 {3, 1, 4},
1733 {5, 4, 2}};
1734
1735 const int iv2t[6] = {0, 1, 2, 0, 1, 0};
1736 const int jv2t[6] = {0, 1, 2, 1, 2, 2};
1737
1738 for (int i = 0; i < 6; i++) {
1739 for (int j = 0; j < 6; j++)
1740 sout[i][j] = 0;
1741 }
1742
1743 /* Consider : W = s*R + R*s^t .
1744 * W_ij = Sum_{k<3} [s_ik*r_jk + s_jk*r_ik]
1745 * We look for A_(ij,pq) such as A*R = W
1746 *
1747 * so
1748 * A_(ij,jk) takes s_ik
1749 * and
1750 * A_(ij,ik) takes s_jk
1751 */
1752 for (int ij = 0; ij < 6; ij++) {
1753 int i = iv2t[ij];
1754 int j = jv2t[ij];
1755 for (int k = 0; k < 3; k++) {
1756 int ik = t2v[i][k];
1757 int jk = t2v[j][k];
1758
1759 sout[ij][ik] += s[j][k];
1760 sout[ij][jk] += s[i][k];
1761 }
1762 }
1763}
1764
1765/*----------------------------------------------------------------------------*/
1777/*----------------------------------------------------------------------------*/
1778
1779CS_F_HOST_DEVICE static inline void
1780cs_math_sym_33_double_product(const cs_real_t s1[6],
1781 const cs_real_t s2[6],
1782 const cs_real_t s3[6],
1783 cs_real_t (*restrict sout)[3])
1784{
1785 cs_real_t _sout[3][3];
1786
1787 /* S11 */
1788 _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1789 /* S22 */
1790 _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1791 /* S33 */
1792 _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1793 /* S12 */
1794 _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1795 /* S21 */
1796 _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1797 /* S23 */
1798 _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1799 /* S32 */
1800 _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1801 /* S13 */
1802 _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1803 /* S31 */
1804 _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1805
1806 sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1807 /* S22 */
1808 sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1809 /* S33 */
1810 sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1811 /* S12 */
1812 sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1813 /* S21 */
1814 sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1815 /* S23 */
1816 sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1817 /* S32 */
1818 sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1819 /* S13 */
1820 sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1821 /* S31 */
1822 sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1823}
1824
1825/*----------------------------------------------------------------------------*/
1832/*----------------------------------------------------------------------------*/
1833
1834CS_F_HOST_DEVICE static inline void
1835cs_nvec3(const cs_real_3_t v,
1836 cs_nvec3_t *qv)
1837{
1838 cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1839
1840 qv->meas = magnitude;
1841 if (fabs(magnitude) > cs_math_zero_threshold) {
1842
1843 const cs_real_t inv = 1/magnitude;
1844 qv->unitv[0] = inv * v[0];
1845 qv->unitv[1] = inv * v[1];
1846 qv->unitv[2] = inv * v[2];
1847
1848 }
1849 else
1850 qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1851}
1852
1853/*=============================================================================
1854 * Public function prototypes
1855 *============================================================================*/
1856
1857/*----------------------------------------------------------------------------*/
1858/*
1859 * \brief Compute the length (Euclidean norm) between two points xa and xb in
1860 * a Cartesian coordinate system of dimension 3
1861 *
1862 * \param[in] xa coordinate of the first extremity
1863 * \param[in] xb coordinate of the second extremity
1864 * \param[out] len pointer to the length of the vector va -> vb
1865 * \param[out] unitv unitary vector along xa -> xb
1866 */
1867/*----------------------------------------------------------------------------*/
1868
1869CS_F_HOST_DEVICE void
1870cs_math_3_length_unitv(const cs_real_t xa[3],
1871 const cs_real_t xb[3],
1872 cs_real_t *len,
1873 cs_real_3_t unitv);
1874
1875/*----------------------------------------------------------------------------*/
1887/*----------------------------------------------------------------------------*/
1888
1889void
1890cs_math_sym_33_eigen(const cs_real_t m[6],
1891 cs_real_t eig_vals[3]);
1892
1893/*----------------------------------------------------------------------------*/
1906/*----------------------------------------------------------------------------*/
1907
1908void
1909cs_math_33_eigen(const cs_real_t m[3][3],
1910 cs_real_t *eig_ratio,
1911 cs_real_t *eig_max);
1912
1913/*----------------------------------------------------------------------------*/
1924/*----------------------------------------------------------------------------*/
1925
1926double
1927cs_math_surftri(const cs_real_t xv[3],
1928 const cs_real_t xe[3],
1929 const cs_real_t xf[3]);
1930
1931/*----------------------------------------------------------------------------*/
1943/*----------------------------------------------------------------------------*/
1944
1945double
1946cs_math_voltet(const cs_real_t xv[3],
1947 const cs_real_t xe[3],
1948 const cs_real_t xf[3],
1949 const cs_real_t xc[3]);
1950
1951/*----------------------------------------------------------------------------*/
1964/*----------------------------------------------------------------------------*/
1965
1966void
1967cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1968 const cs_real_t tol_err,
1969 cs_real_t eig_val[3],
1970 cs_real_t eig_vec[3][3]);
1971
1972/*----------------------------------------------------------------------------*/
1982/*----------------------------------------------------------------------------*/
1983
1984void
1985cs_math_fact_lu(cs_lnum_t n_blocks,
1986 const int b_size,
1987 const cs_real_t *a,
1988 cs_real_t *a_lu);
1989
1990/*----------------------------------------------------------------------------*/
2000/*----------------------------------------------------------------------------*/
2001
2002void
2003cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
2004 const int n,
2005 cs_real_t x[],
2006 const cs_real_t b[]);
2007
2008/*----------------------------------------------------------------------------*/
2018/*----------------------------------------------------------------------------*/
2019
2020void
2021cs_math_sym_44_factor_ldlt(cs_real_t ldlt[10]);
2022
2023/*----------------------------------------------------------------------------*/
2036/*----------------------------------------------------------------------------*/
2037
2038cs_real_t
2039cs_math_sym_44_partial_solve_ldlt(const cs_real_t ldlt[10],
2040 const cs_real_t rhs[4]);
2041
2042/*----------------------------------------------------------------------------*/
2043
2044END_C_DECLS
2045
2046#endif /* __CS_MATH_H__ */
bft_error.h
cs_defs.h
restrict
#define restrict
Definition: cs_defs.h:145
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:542
CS_F_HOST_DEVICE
#define CS_F_HOST_DEVICE
Definition: cs_defs.h:561
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:359
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:361
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:543
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:368
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:335
t
@ t
Definition: cs_field_pointer.h:94
k
@ k
Definition: cs_field_pointer.h:72
x2
@ x2
Definition: cs_field_pointer.h:226
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:800
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:577
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:745
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 *restrict sout)
Compute the inverse of a symmetric matrix using Cramer's rule.
Definition: cs_math.h:1373
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[3], cs_real_t eig_vec[3][3])
Evaluate eigenvalues and eigenvectors of a real symmetric matrix m1[3,3]: m1*m2 = lambda*m2.
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:1604
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:442
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 *restrict mv)
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:1005
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:1338
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
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 *restrict mv)
Compute the product of a symmetric matrix of 3x3 real values by a vector of 3 real values....
Definition: cs_math.h:983
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:561
cs_math_4ov3
const cs_real_t cs_math_4ov3
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 *restrict mv)
Compute the product of the transpose of a matrix of 3x3 real values by a vector of 3 real values.
Definition: cs_math.h:961
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:674
cs_math_sym_33_determinant
static CS_F_HOST_DEVICE cs_real_t cs_math_sym_33_determinant(const cs_real_t m[6])
Compute the determinant of a 3x3 symmetric matrix.
Definition: cs_math.h:1142
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:1479
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 *restrict sout)
Compute the product of two symmetric matrices.
Definition: cs_math.h:1699
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:1435
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 *restrict mv)
Compute the product of a matrix of 6x6 real values by a vector of 6 real values.
Definition: cs_math.h:1079
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:230
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:1643
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:593
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:688
cs_math_1ov12
const cs_real_t cs_math_1ov12
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(*restrict sout)[3])
Compute the product of three symmetric matrices.
Definition: cs_math.h:1780
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 *restrict mv)
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:1101
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:630
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(*restrict sout)[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:1728
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:654
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:781
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:894
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:826
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 *restrict vout)
Orthogonal projection of a vector with respect to a normalised vector.
Definition: cs_math.h:851
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:319
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:872
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:792
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:1303
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:1271
cs_math_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition: cs_math.h:67
ZZ
@ ZZ
Definition: cs_math.h:71
XY
@ XY
Definition: cs_math.h:72
XZ
@ XZ
Definition: cs_math.h:74
YZ
@ YZ
Definition: cs_math.h:73
YY
@ YY
Definition: cs_math.h:70
XX
@ XX
Definition: cs_math.h:69
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:409
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:1212
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 *restrict mv)
Compute the product of a matrix of 3x3 real values by a vector of 3 real values add.
Definition: cs_math.h:940
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:717
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:629
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:610
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:1563
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:696
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:1122
cs_math_epzero
const cs_real_t cs_math_epzero
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:1406
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:1519
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:472
cs_math_binom
static int cs_math_binom(int n, int k)
Computes the binomial coefficient of n and k.
Definition: cs_math.h:436
cs_math_3_average
static CS_F_HOST_DEVICE void cs_math_3_average(const cs_real_t u[3], const cs_real_t v[3], cs_real_t *restrict uv)
Compute the average of two vector of dimension 3.
Definition: cs_math.h:1162
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(*restrict mout)[3])
Add the product of two 3x3 real matrices to a matrix.
Definition: cs_math.h:1665
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:733
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:1835
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:1062
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:1027
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:525
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 *restrict uv)
Compute the cross product of two vectors of 3 real values.
Definition: cs_math.h:1186
cs_math_33_identity
static const cs_real_33_t cs_math_33_identity
Definition: cs_math.h:119
cs_math_sym_33_identity
static const cs_real_6_t cs_math_sym_33_identity
Definition: cs_math.h:122
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:765
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:466
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:1046
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:484
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
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 *restrict mv)
Compute the product of a matrix of 3x3 real values by a vector of 3 real values.
Definition: cs_math.h:919
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:545
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:1234
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:503
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:400
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:402
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:403