7.1
general documentation
cs_math.h
Go to the documentation of this file.
1 #ifndef __CS_MATH_H__
2 #define __CS_MATH_H__
3 
4 /*============================================================================
5  * Mathematical base functions.
6  *============================================================================*/
7 
8 /*
9  This file is part of Code_Saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-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 /* Symmetric tensor component name */
60 
61 typedef enum {
62 
63  XX,
64  YY,
65  ZZ,
66  XY,
67  YZ,
68  XZ
69 
70 } cs_math_sym_tensor_component_t;
71 
72 /*============================================================================
73  * Global variables
74  *============================================================================*/
75 
76 /* Numerical constants */
77 
78 extern const cs_real_t cs_math_zero_threshold;
79 extern const cs_real_t cs_math_1ov3;
80 extern const cs_real_t cs_math_2ov3;
81 extern const cs_real_t cs_math_4ov3;
82 extern const cs_real_t cs_math_5ov3;
83 extern const cs_real_t cs_math_1ov6;
84 extern const cs_real_t cs_math_1ov12;
85 extern const cs_real_t cs_math_1ov24;
86 extern const cs_real_t cs_math_epzero;
87 extern const cs_real_t cs_math_infinite_r;
88 extern const cs_real_t cs_math_big_r;
89 extern const cs_real_t cs_math_pi;
90 
91 /*=============================================================================
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 /*----------------------------------------------------------------------------*/
1124 /*----------------------------------------------------------------------------*/
1125 
1126 static inline void
1127 cs_math_sym_33_product(const cs_real_t s1[6],
1128  const cs_real_t s2[6],
1129  cs_real_t sout[restrict 6])
1130 {
1131  /* S11 */
1132  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1133  /* S22 */
1134  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1135  /* S33 */
1136  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1137  /* S12 = S21 */
1138  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1139  /* S23 = S32 */
1140  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1141  /* S13 = S31 */
1142  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1143 }
1144 
1145 /*----------------------------------------------------------------------------*/
1153 /*----------------------------------------------------------------------------*/
1154 
1155 static inline void
1156 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
1157  cs_real_t sout[restrict 6][6])
1158 {
1159  int tens2vect[3][3];
1160  int iindex[6], jindex[6];
1161 
1162  tens2vect[0][0] = 0; tens2vect[0][1] = 3; tens2vect[0][2] = 5;
1163  tens2vect[1][0] = 3; tens2vect[1][1] = 1; tens2vect[1][2] = 4;
1164  tens2vect[2][0] = 5; tens2vect[2][1] = 4; tens2vect[2][2] = 2;
1165 
1166  iindex[0] = 0; iindex[1] = 1; iindex[2] = 2;
1167  iindex[3] = 0; iindex[4] = 1; iindex[5] = 0;
1168 
1169  jindex[0] = 0; jindex[1] = 1; jindex[2] = 2;
1170  jindex[3] = 1; jindex[4] = 2; jindex[5] = 2;
1171 
1172  /* Consider : W = R*s^t + s*R.
1173  * W_ij = Sum_{k<3} [s_jk*r_ik + s_ik*r_jk]
1174  * We look for A such as A*R = W
1175  */
1176  for (int i = 0; i < 6; i++) {
1177  int ii = iindex[i];
1178  int jj = jindex[i];
1179  for (int k = 0; k < 3; k++) {
1180  int ik = tens2vect[k][ii];
1181  int jk = tens2vect[k][jj];
1182 
1183  sout[ik][i] += s[k][jj];
1184 
1185  sout[jk][i] += s[k][ii];
1186  }
1187  }
1188 }
1189 
1190 /*----------------------------------------------------------------------------*/
1202 /*----------------------------------------------------------------------------*/
1203 
1204 static inline void
1205 cs_math_sym_33_double_product(const cs_real_t s1[6],
1206  const cs_real_t s2[6],
1207  const cs_real_t s3[6],
1208  cs_real_t sout[restrict 3][3])
1209 {
1210  cs_real_t _sout[3][3];
1211 
1212  /* S11 */
1213  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1214  /* S22 */
1215  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1216  /* S33 */
1217  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1218  /* S12 */
1219  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1220  /* S21 */
1221  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1222  /* S23 */
1223  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1224  /* S32 */
1225  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1226  /* S13 */
1227  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1228  /* S31 */
1229  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1230 
1231  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1232  /* S22 */
1233  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1234  /* S33 */
1235  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1236  /* S12 */
1237  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1238  /* S21 */
1239  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1240  /* S23 */
1241  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1242  /* S32 */
1243  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1244  /* S13 */
1245  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1246  /* S31 */
1247  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1248 }
1249 
1250 /*----------------------------------------------------------------------------*/
1257 /*----------------------------------------------------------------------------*/
1258 
1259 static inline void
1260 cs_nvec3(const cs_real_3_t v,
1261  cs_nvec3_t *qv)
1262 {
1263  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1264 
1265  qv->meas = magnitude;
1266  if (fabs(magnitude) > cs_math_zero_threshold) {
1267 
1268  const cs_real_t inv = 1/magnitude;
1269  qv->unitv[0] = inv * v[0];
1270  qv->unitv[1] = inv * v[1];
1271  qv->unitv[2] = inv * v[2];
1272 
1273  }
1274  else
1275  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1276 }
1277 
1278 /*=============================================================================
1279  * Public function prototypes
1280  *============================================================================*/
1281 
1282 /*----------------------------------------------------------------------------*/
1286 /*----------------------------------------------------------------------------*/
1287 
1288 void
1289 cs_math_set_machine_epsilon(void);
1290 
1291 /*----------------------------------------------------------------------------*/
1295 /*----------------------------------------------------------------------------*/
1296 
1297 double
1298 cs_math_get_machine_epsilon(void);
1299 
1300 /*----------------------------------------------------------------------------*/
1310 /*----------------------------------------------------------------------------*/
1311 
1312 void
1313 cs_math_3_length_unitv(const cs_real_t xa[3],
1314  const cs_real_t xb[3],
1315  cs_real_t *len,
1316  cs_real_3_t unitv);
1317 
1318 /*----------------------------------------------------------------------------*/
1330 /*----------------------------------------------------------------------------*/
1331 
1332 void
1333 cs_math_sym_33_eigen(const cs_real_t m[6],
1334  cs_real_t eig_vals[3]);
1335 
1336 /*----------------------------------------------------------------------------*/
1349 /*----------------------------------------------------------------------------*/
1350 
1351 void
1352 cs_math_33_eigen(const cs_real_t m[3][3],
1353  cs_real_t *eig_ratio,
1354  cs_real_t *eig_max);
1355 
1356 /*----------------------------------------------------------------------------*/
1367 /*----------------------------------------------------------------------------*/
1368 
1369 double
1370 cs_math_surftri(const cs_real_t xv[3],
1371  const cs_real_t xe[3],
1372  const cs_real_t xf[3]);
1373 
1374 /*----------------------------------------------------------------------------*/
1386 /*----------------------------------------------------------------------------*/
1387 
1388 double
1389 cs_math_voltet(const cs_real_t xv[3],
1390  const cs_real_t xe[3],
1391  const cs_real_t xf[3],
1392  const cs_real_t xc[3]);
1393 
1394 /*----------------------------------------------------------------------------*/
1407 /*----------------------------------------------------------------------------*/
1408 
1409 void
1410 cs_math_33_eig_val_vec(const cs_real_t m_in[3][3],
1411  const cs_real_t tol_err,
1412  cs_real_t eig_val[restrict 3],
1413  cs_real_t eig_vec[restrict 3][3]);
1414 
1415 /*----------------------------------------------------------------------------*/
1425 /*----------------------------------------------------------------------------*/
1426 
1427 void
1428 cs_math_fact_lu(cs_lnum_t n_blocks,
1429  const int b_size,
1430  const cs_real_t *a,
1431  cs_real_t *a_lu);
1432 
1433 /*----------------------------------------------------------------------------*/
1443 /*----------------------------------------------------------------------------*/
1444 
1445 void
1446 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1447  const int n,
1448  cs_real_t x[],
1449  const cs_real_t b[]);
1450 
1451 /*----------------------------------------------------------------------------*/
1452 
1453 END_C_DECLS
1454 
1455 #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:142
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:337
cs_math_big_r
const cs_real_t cs_math_big_r
cs_math_sym_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:510
cs_math_fw_and_bw_lu
void cs_math_fw_and_bw_lu(const cs_real_t a_lu[], const int n, cs_real_t x[], const cs_real_t b[])
Block Jacobi utilities. Compute forward and backward to solve an LU P*P system.
Definition: cs_math.c:715
cs_math_epzero
const cs_real_t cs_math_epzero
cs_math_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:351
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:322
cs_nvec3_t
Definition: cs_defs.h:368
cs_math_surftri
double cs_math_surftri(const cs_real_t xv[3], const cs_real_t xe[3], const cs_real_t xf[3])
Compute the area of the convex_hull generated by 3 points. This corresponds to the computation of the...
Definition: cs_math.c:471
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
cs_math_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:245
cs_math_set_machine_epsilon
void cs_math_set_machine_epsilon(void)
Compute the value related to the machine precision.
Definition: cs_math.c:224
bft_error.h
cs_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:501
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:370
cs_math_sym_33_eigen
void cs_math_sym_33_eigen(const cs_real_t m[6], cs_real_t eig_vals[3])
Compute all eigenvalues of a 3x3 symmetric matrix with symmetric storage.
Definition: cs_math.c:265
cs_math_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:533
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:335
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:371
XY
Definition: cs_math.h:66
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:316
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:511
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:342
t
Definition: cs_field_pointer.h:92
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 (Euclidean norm) between two points xa and xb in a Cartesian coordinate system of ...
Definition: cs_math.c:438
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:657
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:146