7.0
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-2021 EDF S.A.
12 
13  This program is free software; you can redistribute it and/or modify it under
14  the terms of the GNU General Public License as published by the Free Software
15  Foundation; either version 2 of the License, or (at your option) any later
16  version.
17 
18  This program is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21  details.
22 
23  You should have received a copy of the GNU General Public License along with
24  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25  Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27 
28 /*----------------------------------------------------------------------------*/
29 
30 #include "cs_defs.h"
31 
32 /*----------------------------------------------------------------------------
33  * Standard C library headers
34  *----------------------------------------------------------------------------*/
35 
36 #include <assert.h>
37 #include <math.h>
38 
39 /*----------------------------------------------------------------------------
40  * Local headers
41  *----------------------------------------------------------------------------*/
42 
43 #if defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
44 #include "bft_error.h"
45 #endif
46 
47 /*----------------------------------------------------------------------------*/
48 
49 BEGIN_C_DECLS
50 
51 /*=============================================================================
52  * Local Macro definitions
53  *============================================================================*/
54 
55 /*============================================================================
56  * Type definition
57  *============================================================================*/
58 
59 /* Symetric tensor component name */
60 
61 typedef enum {
62 
63  XX,
64  YY,
65  ZZ,
66  XY,
67  YZ,
68  XZ
69 
70 } cs_math_sym_tensor_component_t;
71 
72 /*============================================================================
73  * Global variables
74  *============================================================================*/
75 
76 /* Numerical constants */
77 
78 extern const cs_real_t cs_math_zero_threshold;
79 extern const cs_real_t cs_math_1ov3;
80 extern const cs_real_t cs_math_2ov3;
81 extern const cs_real_t cs_math_4ov3;
82 extern const cs_real_t cs_math_5ov3;
83 extern const cs_real_t cs_math_1ov6;
84 extern const cs_real_t cs_math_1ov12;
85 extern const cs_real_t cs_math_1ov24;
86 extern const cs_real_t cs_math_epzero;
87 extern const cs_real_t cs_math_infinite_r;
88 extern const cs_real_t cs_math_big_r;
89 extern const cs_real_t cs_math_pi;
90 
91 /*=============================================================================
92  * Inline static function prototypes
93  *============================================================================*/
94 
95 /*----------------------------------------------------------------------------*/
104 /*----------------------------------------------------------------------------*/
105 
106 static inline int
107 cs_math_binom(int n,
108  int k)
109 {
110  int ret = 1;
111  assert(n >= k);
112 
113  const int n_iter = (k > n-k) ? n-k : k;
114  for (int j = 1; j <= n_iter; j++, n--) {
115  if (n % j == 0)
116  ret *= n/j;
117  else if (ret % j == 0)
118  ret = ret/j*n;
119  else
120  ret = (ret*n)/j;
121  }
122 
123  return ret;
124 }
125 
126 /*----------------------------------------------------------------------------*/
134 /*----------------------------------------------------------------------------*/
135 
136 static inline cs_real_t
137 cs_math_fabs(cs_real_t x)
138 {
139  cs_real_t ret = (x < 0) ? -x : x;
140 
141  return ret;
142 }
143 
144 /*----------------------------------------------------------------------------*/
152 /*----------------------------------------------------------------------------*/
153 
154 static inline cs_real_t
155 cs_math_fmin(cs_real_t x,
156  cs_real_t y)
157 {
158  cs_real_t ret = (x < y) ? x : y;
159 
160  return ret;
161 }
162 
163 /*----------------------------------------------------------------------------*/
171 /*----------------------------------------------------------------------------*/
172 
173 static inline cs_real_t
174 cs_math_fmax(cs_real_t x,
175  cs_real_t y)
176 {
177  cs_real_t ret = (x < y) ? y : x;
178 
179  return ret;
180 }
181 
182 /*----------------------------------------------------------------------------*/
190 /*----------------------------------------------------------------------------*/
191 
192 static inline cs_real_t
193 cs_math_sq(cs_real_t x)
194 {
195  return x*x;
196 }
197 
198 /*----------------------------------------------------------------------------*/
206 /*----------------------------------------------------------------------------*/
207 
208 static inline cs_real_t
209 cs_math_pow2(cs_real_t x)
210 {
211  return x*x;
212 }
213 
214 /*----------------------------------------------------------------------------*/
222 /*----------------------------------------------------------------------------*/
223 
224 static inline cs_real_t
225 cs_math_pow3(cs_real_t x)
226 {
227  return x*x*x;
228 }
229 
230 /*----------------------------------------------------------------------------*/
238 /*----------------------------------------------------------------------------*/
239 
240 static inline cs_real_t
241 cs_math_pow4(cs_real_t x)
242 {
243  return (x*x)*(x*x);
244 }
245 
246 /*----------------------------------------------------------------------------*/
256 /*----------------------------------------------------------------------------*/
257 
258 static inline cs_real_t
259 cs_math_3_distance(const cs_real_t xa[3],
260  const cs_real_t xb[3])
261 {
262  cs_real_t v[3];
263 
264  v[0] = xb[0] - xa[0];
265  v[1] = xb[1] - xa[1];
266  v[2] = xb[2] - xa[2];
267 
268  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
269 }
270 
271 /*----------------------------------------------------------------------------*/
281 /*----------------------------------------------------------------------------*/
282 
283 static inline cs_real_t
284 cs_math_3_distance_dot_product(const cs_real_t xa[3],
285  const cs_real_t xb[3],
286  const cs_real_t xc[3])
287 {
288  return ((xb[0] - xa[0])*xc[0]+(xb[1] - xa[1])*xc[1]+(xb[2] - xa[2])*xc[2]);
289 }
290 
291 /*----------------------------------------------------------------------------*/
301 /*----------------------------------------------------------------------------*/
302 
303 static inline cs_real_t
304 cs_math_3_square_distance(const cs_real_t xa[3],
305  const cs_real_t xb[3])
306 {
307  cs_real_t v[3] = {xb[0] - xa[0],
308  xb[1] - xa[1],
309  xb[2] - xa[2]};
310 
311  return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
312 }
313 
314 /*----------------------------------------------------------------------------*/
323 /*----------------------------------------------------------------------------*/
324 
325 static inline cs_real_t
326 cs_math_3_dot_product(const cs_real_t u[3],
327  const cs_real_t v[3])
328 {
329  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
330 
331  return uv;
332 }
333 
334 /*----------------------------------------------------------------------------*/
344 /*----------------------------------------------------------------------------*/
345 
346 static inline cs_real_t
347 cs_math_3_33_3_dot_product(const cs_real_t n1[3],
348  const cs_real_t t[3][3],
349  const cs_real_t n2[3])
350 {
351  cs_real_t n_t_n
352  = ( n1[0]*t[0][0]*n2[0] + n1[1]*t[1][0]*n2[0] + n1[2]*t[2][0]*n2[0]
353  + n1[0]*t[0][1]*n2[1] + n1[1]*t[1][1]*n2[1] + n1[2]*t[2][1]*n2[1]
354  + n1[0]*t[0][2]*n2[2] + n1[1]*t[1][2]*n2[2] + n1[2]*t[2][2]*n2[2]);
355  return n_t_n;
356 }
357 
358 /*----------------------------------------------------------------------------*/
372 /*----------------------------------------------------------------------------*/
373 
374 static inline cs_real_t
375 cs_math_3_sym_33_3_dot_product(const cs_real_t n1[3],
376  const cs_real_t t[6],
377  const cs_real_t n2[3])
378 {
379  cs_real_t n_t_n
380  = ( n1[0]*t[0]*n2[0] + n1[1]*t[3]*n2[0] + n1[2]*t[5]*n2[0]
381  + n1[0]*t[3]*n2[1] + n1[1]*t[1]*n2[1] + n1[2]*t[4]*n2[1]
382  + n1[0]*t[5]*n2[2] + n1[1]*t[4]*n2[2] + n1[2]*t[2]*n2[2]);
383  return n_t_n;
384 }
385 
386 /*----------------------------------------------------------------------------*/
394 /*----------------------------------------------------------------------------*/
395 
396 static inline cs_real_t
397 cs_math_3_norm(const cs_real_t v[3])
398 {
399  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
400 }
401 
402 /*----------------------------------------------------------------------------*/
410 /*----------------------------------------------------------------------------*/
411 
412 static inline cs_real_t
413 cs_math_3_square_norm(const cs_real_t v[3])
414 {
415  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
416 
417  return v2;
418 }
419 
420 /*----------------------------------------------------------------------------*/
429 /*----------------------------------------------------------------------------*/
430 
431 static inline void
432 cs_math_3_normalise(const cs_real_t vin[3],
433  cs_real_t vout[3])
434 {
435  cs_real_t norm = cs_math_3_norm(vin);
436 
437  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
438 
439  vout[0] = inv_norm * vin[0];
440  vout[1] = inv_norm * vin[1];
441  vout[2] = inv_norm * vin[2];
442 }
443 
444 /*----------------------------------------------------------------------------*/
453 /*----------------------------------------------------------------------------*/
454 
455 static inline void
456 cs_math_3_normalize(const cs_real_t vin[3],
457  cs_real_t vout[3])
458 {
459  cs_real_t norm = cs_math_3_norm(vin);
460 
461  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
462 
463  vout[0] = inv_norm * vin[0];
464  vout[1] = inv_norm * vin[1];
465  vout[2] = inv_norm * vin[2];
466 }
467 
468 /*----------------------------------------------------------------------------*/
477 /*----------------------------------------------------------------------------*/
478 
479 static inline void
480 cs_math_3_orthogonal_projection(const cs_real_t n[3],
481  const cs_real_t v[3],
482  cs_real_t vout[restrict 3])
483 {
484  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
485  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
486  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
487 }
488 
489 /*----------------------------------------------------------------------------*/
498 /*----------------------------------------------------------------------------*/
499 
500 static inline void
501 cs_math_3_normal_scaling(const cs_real_t n[3],
502  cs_real_t factor,
503  cs_real_t v[3])
504 {
505  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
506  for (int i = 0; i < 3; i++)
507  v[i] += v_dot_n * n[i];
508 }
509 
510 /*----------------------------------------------------------------------------*/
520 /*----------------------------------------------------------------------------*/
521 
522 static inline void
523 cs_math_33_normal_scaling_add(const cs_real_t n[3],
524  cs_real_t factor,
525  cs_real_t t[3][3])
526 {
527  cs_real_t n_t_n = (factor -1.) *
528  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
529  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
530  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
531  for (int i = 0; i < 3; i++) {
532  for (int j = 0; j < 3; j++)
533  t[i][j] += n_t_n * n[i] * n[j];
534  }
535 }
536 /*----------------------------------------------------------------------------*/
545 /*----------------------------------------------------------------------------*/
546 
547 static inline void
548 cs_math_33_3_product(const cs_real_t m[3][3],
549  const cs_real_t v[3],
550  cs_real_t mv[restrict 3])
551 {
552  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
553  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
554  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
555 }
556 
557 /*----------------------------------------------------------------------------*/
566 /*----------------------------------------------------------------------------*/
567 
568 static inline void
569 cs_math_33_3_product_add(const cs_real_t m[3][3],
570  const cs_real_t v[3],
571  cs_real_t mv[restrict 3])
572 {
573  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
574  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
575  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
576 }
577 
578 /*----------------------------------------------------------------------------*/
587 /*----------------------------------------------------------------------------*/
588 
589 static inline void
590 cs_math_33t_3_product(const cs_real_t m[3][3],
591  const cs_real_t v[3],
592  cs_real_t mv[restrict 3])
593 {
594  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
595  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
596  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
597 }
598 
599 /*----------------------------------------------------------------------------*/
609 /*----------------------------------------------------------------------------*/
610 
611 static inline void
612 cs_math_sym_33_3_product(const cs_real_t m[6],
613  const cs_real_t v[3],
614  cs_real_t mv[restrict 3])
615 {
616  mv[0] = m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
617  mv[1] = m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
618  mv[2] = m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
619 }
620 
621 /*----------------------------------------------------------------------------*/
631 /*----------------------------------------------------------------------------*/
632 
633 static inline void
634 cs_math_sym_33_3_product_add(const cs_real_t m[6],
635  const cs_real_t v[3],
636  cs_real_t mv[restrict 3])
637 {
638  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
639  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
640  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
641 }
642 
643 /*----------------------------------------------------------------------------*/
651 /*----------------------------------------------------------------------------*/
652 
653 static inline cs_real_t
654 cs_math_6_trace(const cs_real_t t[6])
655 {
656  return (t[0] + t[1] + t[2]);
657 }
658 
659 /*----------------------------------------------------------------------------*/
668 /*----------------------------------------------------------------------------*/
669 
670 static inline void
671 cs_math_66_6_product(const cs_real_t m[6][6],
672  const cs_real_t v[6],
673  cs_real_t mv[restrict 6])
674 {
675  for (int i = 0; i < 6; i++) {
676  for (int j = 0; j < 6; j++)
677  mv[i] = m[i][j] * v[j];
678  }
679 }
680 
681 /*----------------------------------------------------------------------------*/
690 /*----------------------------------------------------------------------------*/
691 
692 static inline void
693 cs_math_66_6_product_add(const cs_real_t m[6][6],
694  const cs_real_t v[6],
695  cs_real_t mv[restrict 6])
696 {
697  for (int i = 0; i < 6; i++) {
698  for (int j = 0; j < 6; j++)
699  mv[i] += m[i][j] * v[j];
700  }
701 }
702 
703 /*----------------------------------------------------------------------------*/
711 /*----------------------------------------------------------------------------*/
712 
713 static inline cs_real_t
714 cs_math_33_determinant(const cs_real_t m[3][3])
715 {
716  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
717  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
718  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
719 
720  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
721 }
722 
723 /*----------------------------------------------------------------------------*/
731 /*----------------------------------------------------------------------------*/
732 
733 static inline cs_real_t
734 cs_math_sym_33_determinant(const cs_real_6_t m)
735 {
736  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
737  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
738  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
739 
740  return m[0]*com0 + m[3]*com1 + m[5]*com2;
741 }
742 
743 /*----------------------------------------------------------------------------*/
751 /*----------------------------------------------------------------------------*/
752 
753 #if defined(__INTEL_COMPILER)
754 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
755 #endif
756 
757 static inline void
758 cs_math_3_cross_product(const cs_real_t u[3],
759  const cs_real_t v[3],
760  cs_real_t uv[restrict 3])
761 {
762  uv[0] = u[1]*v[2] - u[2]*v[1];
763  uv[1] = u[2]*v[0] - u[0]*v[2];
764  uv[2] = u[0]*v[1] - u[1]*v[0];
765 }
766 
767 /*----------------------------------------------------------------------------*/
777 /*----------------------------------------------------------------------------*/
778 
779 #if defined(__INTEL_COMPILER)
780 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
781 #endif
782 
783 static inline cs_real_t
784 cs_math_3_triple_product(const cs_real_t u[3],
785  const cs_real_t v[3],
786  const cs_real_t w[3])
787 {
788  return (u[1]*v[2] - u[2]*v[1]) * w[0]
789  + (u[2]*v[0] - u[0]*v[2]) * w[1]
790  + (u[0]*v[1] - u[1]*v[0]) * w[2];
791 }
792 
793 /*----------------------------------------------------------------------------*/
800 /*----------------------------------------------------------------------------*/
801 
802 static inline void
803 cs_math_33_inv_cramer(const cs_real_t in[3][3],
804  cs_real_t out[3][3])
805 {
806  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
807  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
808  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
809 
810  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
811  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
812  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
813 
814  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
815  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
816  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
817 
818  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
819  const double invdet = 1/det;
820 
821  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
822  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
823  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
824 }
825 
826 /*----------------------------------------------------------------------------*/
832 /*----------------------------------------------------------------------------*/
833 
834 static inline void
835 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
836 {
837  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
838  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
839  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
840  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
841  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
842  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
843  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
844  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
845  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
846 
847  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
848 
849  a[0][0] = a00 * det_inv;
850  a[0][1] = a01 * det_inv;
851  a[0][2] = a02 * det_inv;
852  a[1][0] = a10 * det_inv;
853  a[1][1] = a11 * det_inv;
854  a[1][2] = a12 * det_inv;
855  a[2][0] = a20 * det_inv;
856  a[2][1] = a21 * det_inv;
857  a[2][2] = a22 * det_inv;
858 }
859 
860 /*----------------------------------------------------------------------------*/
867 /*----------------------------------------------------------------------------*/
868 
869 static inline void
870 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
871 {
872  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
873  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
874  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
875  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
876  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
877  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
878 
879  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
880 
881  a[0][0] = a00 * det_inv;
882  a[0][1] = a01 * det_inv;
883  a[0][2] = a02 * det_inv;
884  a[1][0] = a01 * det_inv;
885  a[1][1] = a11 * det_inv;
886  a[1][2] = a12 * det_inv;
887  a[2][0] = a02 * det_inv;
888  a[2][1] = a12 * det_inv;
889  a[2][2] = a22 * det_inv;
890 }
891 
892 /*----------------------------------------------------------------------------*/
902 /*----------------------------------------------------------------------------*/
903 
904 static inline void
905 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
906  cs_real_t sout[restrict 6])
907 {
908  double detinv;
909 
910  sout[0] = s[1]*s[2] - s[4]*s[4];
911  sout[1] = s[0]*s[2] - s[5]*s[5];
912  sout[2] = s[0]*s[1] - s[3]*s[3];
913  sout[3] = s[4]*s[5] - s[3]*s[2];
914  sout[4] = s[3]*s[5] - s[0]*s[4];
915  sout[5] = s[3]*s[4] - s[1]*s[5];
916 
917  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
918 
919  sout[0] *= detinv;
920  sout[1] *= detinv;
921  sout[2] *= detinv;
922  sout[3] *= detinv;
923  sout[4] *= detinv;
924  sout[5] *= detinv;
925 }
926 
927 /*----------------------------------------------------------------------------*/
935 /*----------------------------------------------------------------------------*/
936 
937 static inline void
938 cs_math_33_product(const cs_real_t m1[3][3],
939  const cs_real_t m2[3][3],
940  cs_real_t mout[3][3])
941 {
942  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
943  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
944  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
945 
946  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
947  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
948  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
949 
950  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
951  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
952  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
953 }
954 
955 /*----------------------------------------------------------------------------*/
964 /*----------------------------------------------------------------------------*/
965 
966 static inline void
967 cs_math_33_transform_r_to_a(const cs_real_t m[3][3],
968  const cs_real_t q[3][3],
969  cs_real_t mout[3][3])
970 {
971  /* _m = M.Q */
972  cs_real_33_t _m;
973  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[1][0] + m[0][2]*q[2][0];
974  _m[0][1] = m[0][0]*q[0][1] + m[0][1]*q[1][1] + m[0][2]*q[2][1];
975  _m[0][2] = m[0][0]*q[0][2] + m[0][1]*q[1][2] + m[0][2]*q[2][2];
976 
977  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[1][0] + m[1][2]*q[2][0];
978  _m[1][1] = m[1][0]*q[0][1] + m[1][1]*q[1][1] + m[1][2]*q[2][1];
979  _m[1][2] = m[1][0]*q[0][2] + m[1][1]*q[1][2] + m[1][2]*q[2][2];
980 
981  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[1][0] + m[2][2]*q[2][0];
982  _m[2][1] = m[2][0]*q[0][1] + m[2][1]*q[1][1] + m[2][2]*q[2][1];
983  _m[2][2] = m[2][0]*q[0][2] + m[2][1]*q[1][2] + m[2][2]*q[2][2];
984 
985  /* mout = Q^t _m */
986  mout[0][0] = q[0][0]*_m[0][0] + q[1][0]*_m[1][0] + q[2][0]*_m[2][0];
987  mout[0][1] = q[0][0]*_m[0][1] + q[1][0]*_m[1][1] + q[2][0]*_m[2][1];
988  mout[0][2] = q[0][0]*_m[0][2] + q[1][0]*_m[1][2] + q[2][0]*_m[2][2];
989 
990  mout[1][0] = q[0][1]*_m[0][0] + q[1][1]*_m[1][0] + q[2][1]*_m[2][0];
991  mout[1][1] = q[0][1]*_m[0][1] + q[1][1]*_m[1][1] + q[2][1]*_m[2][1];
992  mout[1][2] = q[0][1]*_m[0][2] + q[1][1]*_m[1][2] + q[2][1]*_m[2][2];
993 
994  mout[2][0] = q[0][2]*_m[0][0] + q[1][2]*_m[1][0] + q[2][2]*_m[2][0];
995  mout[2][1] = q[0][2]*_m[0][1] + q[1][2]*_m[1][1] + q[2][2]*_m[2][1];
996  mout[2][2] = q[0][2]*_m[0][2] + q[1][2]*_m[1][2] + q[2][2]*_m[2][2];
997 }
998 
999 /*----------------------------------------------------------------------------*/
1008 /*----------------------------------------------------------------------------*/
1009 
1010 static inline void
1011 cs_math_33_transform_a_to_r(const cs_real_t m[3][3],
1012  const cs_real_t q[3][3],
1013  cs_real_t mout[3][3])
1014 {
1015  /* _m = M.Q^t */
1016  cs_real_33_t _m;
1017  _m[0][0] = m[0][0]*q[0][0] + m[0][1]*q[0][1] + m[0][2]*q[0][2];
1018  _m[0][1] = m[0][0]*q[1][0] + m[0][1]*q[1][1] + m[0][2]*q[1][2];
1019  _m[0][2] = m[0][0]*q[2][0] + m[0][1]*q[2][1] + m[0][2]*q[2][2];
1020 
1021  _m[1][0] = m[1][0]*q[0][0] + m[1][1]*q[0][1] + m[1][2]*q[0][2];
1022  _m[1][1] = m[1][0]*q[1][0] + m[1][1]*q[1][1] + m[1][2]*q[1][2];
1023  _m[1][2] = m[1][0]*q[2][0] + m[1][1]*q[2][1] + m[1][2]*q[2][2];
1024 
1025  _m[2][0] = m[2][0]*q[0][0] + m[2][1]*q[0][1] + m[2][2]*q[0][2];
1026  _m[2][1] = m[2][0]*q[1][0] + m[2][1]*q[1][1] + m[2][2]*q[1][2];
1027  _m[2][2] = m[2][0]*q[2][0] + m[2][1]*q[2][1] + m[2][2]*q[2][2];
1028 
1029  /* mout = Q _m */
1030  mout[0][0] = q[0][0]*_m[0][0] + q[0][1]*_m[1][0] + q[0][2]*_m[2][0];
1031  mout[0][1] = q[0][0]*_m[0][1] + q[0][1]*_m[1][1] + q[0][2]*_m[2][1];
1032  mout[0][2] = q[0][0]*_m[0][2] + q[0][1]*_m[1][2] + q[0][2]*_m[2][2];
1033 
1034  mout[1][0] = q[1][0]*_m[0][0] + q[1][1]*_m[1][0] + q[1][2]*_m[2][0];
1035  mout[1][1] = q[1][0]*_m[0][1] + q[1][1]*_m[1][1] + q[1][2]*_m[2][1];
1036  mout[1][2] = q[1][0]*_m[0][2] + q[1][1]*_m[1][2] + q[1][2]*_m[2][2];
1037 
1038  mout[2][0] = q[2][0]*_m[0][0] + q[2][1]*_m[1][0] + q[2][2]*_m[2][0];
1039  mout[2][1] = q[2][0]*_m[0][1] + q[2][1]*_m[1][1] + q[2][2]*_m[2][1];
1040  mout[2][2] = q[2][0]*_m[0][2] + q[2][1]*_m[1][2] + q[2][2]*_m[2][2];
1041 }
1042 
1043 /*----------------------------------------------------------------------------*/
1052 /*----------------------------------------------------------------------------*/
1053 
1054 static inline void
1055 cs_math_33_extract_sym_ant(const cs_real_t m[3][3],
1056  cs_real_t m_sym[3][3],
1057  cs_real_t m_ant[3][3])
1058 {
1059  /* sym = 0.5 (m + m_transpose) */
1060  m_sym[0][0] = 0.5 * (m[0][0] + m[0][0]);
1061  m_sym[0][1] = 0.5 * (m[0][1] + m[1][0]);
1062  m_sym[0][2] = 0.5 * (m[0][2] + m[2][0]);
1063  m_sym[1][0] = 0.5 * (m[1][0] + m[0][1]);
1064  m_sym[1][1] = 0.5 * (m[1][1] + m[1][1]);
1065  m_sym[1][2] = 0.5 * (m[1][2] + m[2][1]);
1066  m_sym[2][0] = 0.5 * (m[2][0] + m[0][2]);
1067  m_sym[2][1] = 0.5 * (m[2][1] + m[1][2]);
1068  m_sym[2][2] = 0.5 * (m[2][2] + m[2][2]);
1069 
1070  /* ant = 0.5 (m - m_transpose) */
1071  m_ant[0][0] = 0.5 * (m[0][0] - m[0][0]);
1072  m_ant[0][1] = 0.5 * (m[0][1] - m[1][0]);
1073  m_ant[0][2] = 0.5 * (m[0][2] - m[2][0]);
1074  m_ant[1][0] = 0.5 * (m[1][0] - m[0][1]);
1075  m_ant[1][1] = 0.5 * (m[1][1] - m[1][1]);
1076  m_ant[1][2] = 0.5 * (m[1][2] - m[2][1]);
1077  m_ant[2][0] = 0.5 * (m[2][0] - m[0][2]);
1078  m_ant[2][1] = 0.5 * (m[2][1] - m[1][2]);
1079  m_ant[2][2] = 0.5 * (m[2][2] - m[2][2]);
1080 }
1081 
1082 /*----------------------------------------------------------------------------*/
1090 /*----------------------------------------------------------------------------*/
1091 
1092 static inline void
1093 cs_math_33_product_add(const cs_real_t m1[3][3],
1094  const cs_real_t m2[3][3],
1095  cs_real_t mout[restrict 3][3])
1096 {
1097  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
1098  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
1099  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
1100 
1101  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
1102  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
1103  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
1104 
1105  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
1106  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
1107  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
1108 }
1109 
1110 
1111 /*----------------------------------------------------------------------------*/
1125 /*----------------------------------------------------------------------------*/
1126 
1127 static inline void
1128 cs_math_sym_33_product(const cs_real_t s1[6],
1129  const cs_real_t s2[6],
1130  cs_real_t sout[restrict 6])
1131 {
1132  /* S11 */
1133  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1134  /* S22 */
1135  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1136  /* S33 */
1137  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1138  /* S12 = S21 */
1139  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1140  /* S23 = S32 */
1141  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1142  /* S13 = S31 */
1143  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1144 }
1145 
1146 /*----------------------------------------------------------------------------*/
1154 /*----------------------------------------------------------------------------*/
1155 
1156 static inline void
1157 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1158  cs_real_t sout[restrict 6][6])
1159 {
1160  int tens2vect[3][3];
1161  int iindex[6], jindex[6];
1162 
1163  tens2vect[0][0] = 0; tens2vect[0][1] = 3; tens2vect[0][2] = 5;
1164  tens2vect[1][0] = 3; tens2vect[1][1] = 1; tens2vect[1][2] = 4;
1165  tens2vect[2][0] = 5; tens2vect[2][1] = 4; tens2vect[2][2] = 2;
1166 
1167  iindex[0] = 0; iindex[1] = 1; iindex[2] = 2;
1168  iindex[3] = 0; iindex[4] = 1; iindex[5] = 0;
1169 
1170  jindex[0] = 0; jindex[1] = 1; jindex[2] = 2;
1171  jindex[3] = 1; jindex[4] = 2; jindex[5] = 2;
1172 
1173  /* Consider : W = R*s^t + s*R.
1174  * W_ij = Sum_{k<3} [s_jk*r_ik + s_ik*r_jk]
1175  * We look for A such as A*R = W
1176  */
1177  for (int i = 0; i < 6; i++) {
1178  int ii = iindex[i];
1179  int jj = jindex[i];
1180  for (int k = 0; k < 3; k++) {
1181  int ik = tens2vect[k][ii];
1182  int jk = tens2vect[k][jj];
1183 
1184  sout[ik][i] += s[k][jj];
1185 
1186  sout[jk][i] += s[k][ii];
1187  }
1188  }
1189 }
1190 
1191 /*----------------------------------------------------------------------------*/
1203 /*----------------------------------------------------------------------------*/
1204 
1205 static inline void
1206 cs_math_sym_33_double_product(const cs_real_t s1[6],
1207  const cs_real_t s2[6],
1208  const cs_real_t s3[6],
1209  cs_real_t sout[restrict 3][3])
1210 {
1211  cs_real_t _sout[3][3];
1212 
1213  /* S11 */
1214  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1215  /* S22 */
1216  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1217  /* S33 */
1218  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1219  /* S12 */
1220  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1221  /* S21 */
1222  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1223  /* S23 */
1224  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1225  /* S32 */
1226  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1227  /* S13 */
1228  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1229  /* S31 */
1230  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1231 
1232  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1233  /* S22 */
1234  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1235  /* S33 */
1236  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1237  /* S12 */
1238  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1239  /* S21 */
1240  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1241  /* S23 */
1242  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1243  /* S32 */
1244  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1245  /* S13 */
1246  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1247  /* S31 */
1248  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1249 }
1250 
1251 /*----------------------------------------------------------------------------*/
1258 /*----------------------------------------------------------------------------*/
1259 
1260 static inline void
1261 cs_nvec3(const cs_real_3_t v,
1262  cs_nvec3_t *qv)
1263 {
1264  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1265 
1266  qv->meas = magnitude;
1267  if (fabs(magnitude) > cs_math_zero_threshold) {
1268 
1269  const cs_real_t inv = 1/magnitude;
1270  qv->unitv[0] = inv * v[0];
1271  qv->unitv[1] = inv * v[1];
1272  qv->unitv[2] = inv * v[2];
1273 
1274  }
1275  else
1276  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1277 
1278 }
1279 
1280 /*=============================================================================
1281  * Public function prototypes
1282  *============================================================================*/
1283 
1284 /*----------------------------------------------------------------------------*/
1288 /*----------------------------------------------------------------------------*/
1289 
1290 void
1291 cs_math_set_machine_epsilon(void);
1292 
1293 /*----------------------------------------------------------------------------*/
1297 /*----------------------------------------------------------------------------*/
1298 
1299 double
1300 cs_math_get_machine_epsilon(void);
1301 
1302 /*----------------------------------------------------------------------------*/
1312 /*----------------------------------------------------------------------------*/
1313 
1314 void
1315 cs_math_3_length_unitv(const cs_real_t xa[3],
1316  const cs_real_t xb[3],
1317  cs_real_t *len,
1318  cs_real_3_t unitv);
1319 
1320 /*----------------------------------------------------------------------------*/
1332 /*----------------------------------------------------------------------------*/
1333 
1334 void
1335 cs_math_sym_33_eigen(const cs_real_t m[6],
1336  cs_real_t eig_vals[3]);
1337 
1338 /*----------------------------------------------------------------------------*/
1351 /*----------------------------------------------------------------------------*/
1352 
1353 void
1354 cs_math_33_eigen(const cs_real_t m[3][3],
1355  cs_real_t *eig_ratio,
1356  cs_real_t *eig_max);
1357 
1358 /*----------------------------------------------------------------------------*/
1369 /*----------------------------------------------------------------------------*/
1370 
1371 double
1372 cs_math_surftri(const cs_real_t xv[3],
1373  const cs_real_t xe[3],
1374  const cs_real_t xf[3]);
1375 
1376 /*----------------------------------------------------------------------------*/
1388 /*----------------------------------------------------------------------------*/
1389 
1390 double
1391 cs_math_voltet(const cs_real_t xv[3],
1392  const cs_real_t xe[3],
1393  const cs_real_t xf[3],
1394  const cs_real_t xc[3]);
1395 
1396 /*----------------------------------------------------------------------------*/
1409 /*----------------------------------------------------------------------------*/
1410 
1411 void
1412 cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1413  const cs_real_t tol_err,
1414  cs_real_t eig_val[restrict 3],
1415  cs_real_t eig_vec[restrict 3][3]);
1416 
1417 /*----------------------------------------------------------------------------*/
1427 /*----------------------------------------------------------------------------*/
1428 
1429 void
1430 cs_math_fact_lu(cs_lnum_t n_blocks,
1431  const int b_size,
1432  const cs_real_t *a,
1433  cs_real_t *a_lu);
1434 
1435 /*----------------------------------------------------------------------------*/
1445 /*----------------------------------------------------------------------------*/
1446 
1447 void
1448 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1449  const int n,
1450  cs_real_t x[],
1451  const cs_real_t b[]);
1452 
1453 /*----------------------------------------------------------------------------*/
1454 
1455 END_C_DECLS
1456 
1457 #endif /* __CS_MATH_H__ */
k
Definition: cs_field_pointer.h:70
numvar::ik
integer, save ik
Definition: numvar.f90:75
restrict
#define restrict
Definition: cs_defs.h:127
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:322
cs_math_big_r
const cs_real_t cs_math_big_r
cs_math_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition: cs_math.h:61
XZ
Definition: cs_math.h:68
cs_math_pi
const cs_real_t cs_math_pi
cs_defs.h
cs_math_1ov6
const cs_real_t cs_math_1ov6
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:495
cs_math_fw_and_bw_lu
void cs_math_fw_and_bw_lu(const cs_real_t a_lu[], const int n, cs_real_t x[], const cs_real_t b[])
Block Jacobi utilities. Compute forward and backward to solve an LU P*P system.
Definition: cs_math.c:700
cs_math_epzero
const cs_real_t cs_math_epzero
cs_math_1ov24
const cs_real_t cs_math_1ov24
YY
Definition: cs_math.h:64
cs_math_1ov12
const cs_real_t cs_math_1ov12
cs_math_2ov3
const cs_real_t cs_math_2ov3
cs_math_4ov3
const cs_real_t cs_math_4ov3
ZZ
Definition: cs_math.h:65
cs_math_33_eigen
void cs_math_33_eigen(const cs_real_t m[3][3], cs_real_t *eig_ratio, cs_real_t *eig_max)
Compute max/min eigenvalues ratio and max. eigenvalue of a 3x3 symmetric matrix with non-symmetric st...
Definition: cs_math.c:336
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:307
cs_nvec3_t
Definition: cs_defs.h:353
cs_math_surftri
double cs_math_surftri(const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3])
Compute the area of the convex_hull generated by 3 points. This corresponds to the computation of the...
Definition: cs_math.c:456
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
cs_math_1ov3
const cs_real_t cs_math_1ov3
cs_math_get_machine_epsilon
double cs_math_get_machine_epsilon(void)
Get the value related to the machine precision.
Definition: cs_math.c:230
cs_math_set_machine_epsilon
void cs_math_set_machine_epsilon(void)
Compute the value related to the machine precision.
Definition: cs_math.c:209
bft_error.h
cs_fuel_incl::a
double precision, save a
Definition: cs_fuel_incl.f90:146
cs_math_voltet
double cs_math_voltet(const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3], const cs_real_t xc[3])
Compute the volume of the convex_hull generated by 4 points. This is equivalent to the computation of...
Definition: cs_math.c:486
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:355
cs_math_sym_33_eigen
void cs_math_sym_33_eigen(const cs_real_t m[6], cs_real_t eig_vals[3])
Compute all eigenvalues of a 3x3 symmetric matrix with symmetric storage.
Definition: cs_math.c:250
cs_math_5ov3
const cs_real_t cs_math_5ov3
atimbr::u
double precision, dimension(:,:,:), allocatable u
Definition: atimbr.f90:113
cs_math_33_eig_val_vec
void cs_math_33_eig_val_vec(const cs_real_t m_in[3][3], const cs_real_t tol_err, cs_real_t eig_val[restrict 3], cs_real_t eig_vec[restrict 3][3])
Evaluate eigenvalues and eigenvectors of a real symmetric matrix m1[3,3]: m1*m2 = lambda*m2...
Definition: cs_math.c:518
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:320
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:356
XY
Definition: cs_math.h:66
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:301
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:496
YZ
Definition: cs_math.h:67
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:327
t
Definition: cs_field_pointer.h:98
XX
Definition: cs_math.h:63
cs_math_3_length_unitv
void cs_math_3_length_unitv(const cs_real_t xa[3], const cs_real_t xb[3], cs_real_t *len, cs_real_3_t unitv)
Compute the length (euclidien norm) between two points xa and xb in a cartesian coordinate system of ...
Definition: cs_math.c:423
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
cs_math_fact_lu
void cs_math_fact_lu(cs_lnum_t n_blocks, const int b_size, const cs_real_t *a, cs_real_t *a_lu)
Compute LU factorization of an array of dense matrices of identical size.
Definition: cs_math.c:642
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:146