My Project
programmer's 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-2019 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 /*----------------------------------------------------------------------------*/
346 /*----------------------------------------------------------------------------*/
347 
348 static inline cs_real_t
349 cs_math_3_33_3_dot_product(const cs_real_t n1[3],
350  const cs_real_t t[3][3],
351  const cs_real_t n2[3])
352 {
353  cs_real_t n_t_n =
354  ( n1[0] * t[0][0] * n2[0] + n1[1] * t[1][0] * n2[0] + n1[2] * t[2][0] * n2[0]
355  + n1[0] * t[0][1] * n2[1] + n1[1] * t[1][1] * n2[1] + n1[2] * t[2][1] * n2[1]
356  + n1[0] * t[0][2] * n2[2] + n1[1] * t[1][2] * n2[2] + n1[2] * t[2][2] * n2[2]);
357  return n_t_n;
358 }
359 
360 
361 /*----------------------------------------------------------------------------*/
369 /*----------------------------------------------------------------------------*/
370 
371 static inline cs_real_t
372 cs_math_3_norm(const cs_real_t v[3])
373 {
374  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
375 }
376 
377 /*----------------------------------------------------------------------------*/
385 /*----------------------------------------------------------------------------*/
386 
387 static inline cs_real_t
388 cs_math_3_square_norm(const cs_real_t v[3])
389 {
390  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
391 
392  return v2;
393 }
394 
395 /*----------------------------------------------------------------------------*/
402 /*----------------------------------------------------------------------------*/
403 
404 static inline void
405 cs_math_3_normalise(const cs_real_t vin[3],
406  cs_real_t vout[restrict 3])
407 {
408  cs_real_t norm = cs_math_3_norm(vin);
409 
410  cs_real_t inv_norm = ((norm > cs_math_zero_threshold) ? 1. / norm : 0);
411 
412  vout[0] = inv_norm * vin[0];
413  vout[1] = inv_norm * vin[1];
414  vout[2] = inv_norm * vin[2];
415 }
416 
417 /*----------------------------------------------------------------------------*/
426 /*----------------------------------------------------------------------------*/
427 
428 static inline void
429 cs_math_3_orthogonal_projection(const cs_real_t n[3],
430  const cs_real_t v[3],
431  cs_real_t vout[restrict 3])
432 {
433  vout[0] = v[0]*(1.-n[0]*n[0])- v[1]* n[1]*n[0] - v[2]* n[2]*n[0];
434  vout[1] = -v[0]* n[0]*n[1] + v[1]*(1.-n[1]*n[1])- v[2]* n[2]*n[1];
435  vout[2] = -v[0]* n[0]*n[2] - v[1]* n[1]*n[2] + v[2]*(1.-n[2]*n[2]);
436 }
437 
438 /*----------------------------------------------------------------------------*/
447 /*----------------------------------------------------------------------------*/
448 
449 static inline void
450 cs_math_3_normal_scaling(const cs_real_t n[3],
451  cs_real_t factor,
452  cs_real_3_t v)
453 {
454  cs_real_t v_dot_n = (factor -1.) * cs_math_3_dot_product(v, n);
455  for (int i = 0; i < 3; i++)
456  v[i] += v_dot_n * n[i];
457 }
458 
459 /*----------------------------------------------------------------------------*/
469 /*----------------------------------------------------------------------------*/
470 
471 static inline void
472 cs_math_33_normal_scaling_add(const cs_real_t n[3],
473  cs_real_t factor,
474  cs_real_33_t t)
475 {
476  cs_real_t n_t_n = (factor -1.) *
477  ( n[0] * t[0][0] * n[0] + n[1] * t[1][0] * n[0] + n[2] * t[2][0] * n[0]
478  + n[0] * t[0][1] * n[1] + n[1] * t[1][1] * n[1] + n[2] * t[2][1] * n[1]
479  + n[0] * t[0][2] * n[2] + n[1] * t[1][2] * n[2] + n[2] * t[2][2] * n[2]);
480  for (int i = 0; i < 3; i++) {
481  for (int j = 0; j < 3; j++)
482  t[i][j] += n_t_n * n[i] * n[j];
483  }
484 }
485 /*----------------------------------------------------------------------------*/
494 /*----------------------------------------------------------------------------*/
495 
496 static inline void
497 cs_math_33_3_product(const cs_real_t m[3][3],
498  const cs_real_t v[3],
499  cs_real_3_t mv)
500 {
501  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
502  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
503  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
504 }
505 
506 /*----------------------------------------------------------------------------*/
515 /*----------------------------------------------------------------------------*/
516 
517 static inline void
518 cs_math_33_3_product_add(const cs_real_t m[3][3],
519  const cs_real_t v[3],
520  cs_real_3_t mv)
521 {
522  mv[0] += m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
523  mv[1] += m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
524  mv[2] += m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
525 }
526 
527 /*----------------------------------------------------------------------------*/
536 /*----------------------------------------------------------------------------*/
537 
538 static inline void
539 cs_math_33t_3_product(const cs_real_t m[3][3],
540  const cs_real_t v[3],
541  cs_real_3_t mv)
542 {
543  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
544  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
545  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
546 }
547 
548 /*----------------------------------------------------------------------------*/
558 /*----------------------------------------------------------------------------*/
559 
560 static inline void
561 cs_math_sym_33_3_product(const cs_real_t m[6],
562  const cs_real_t v[3],
563  cs_real_t mv[restrict 3])
564 {
565  mv[0] = m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
566  mv[1] = m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
567  mv[2] = m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
568 }
569 
570 /*----------------------------------------------------------------------------*/
580 /*----------------------------------------------------------------------------*/
581 
582 static inline void
583 cs_math_sym_33_3_product_add(const cs_real_t m[6],
584  const cs_real_t v[3],
585  cs_real_t mv[restrict 3])
586 {
587  mv[0] += m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
588  mv[1] += m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
589  mv[2] += m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
590 }
591 
592 /*----------------------------------------------------------------------------*/
601 /*----------------------------------------------------------------------------*/
602 
603 static inline void
604 cs_math_66_6_product(const cs_real_t m[6][6],
605  const cs_real_t v[6],
606  cs_real_t mv[restrict 6])
607 {
608  for (int i = 0; i < 6; i++) {
609  for (int j = 0; j < 6; j++)
610  mv[i] = m[i][j] * v[j];
611  }
612 }
613 
614 /*----------------------------------------------------------------------------*/
623 /*----------------------------------------------------------------------------*/
624 
625 static inline void
626 cs_math_66_6_product_add(const cs_real_t m[6][6],
627  const cs_real_t v[6],
628  cs_real_t mv[restrict 6])
629 {
630  for (int i = 0; i < 6; i++) {
631  for (int j = 0; j < 6; j++)
632  mv[i] += m[i][j] * v[j];
633  }
634 }
635 
636 /*----------------------------------------------------------------------------*/
644 /*----------------------------------------------------------------------------*/
645 
646 static inline cs_real_t
647 cs_math_33_determinant(const cs_real_t m[3][3])
648 {
649  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
650  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
651  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
652 
653  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
654 }
655 
656 /*----------------------------------------------------------------------------*/
664 /*----------------------------------------------------------------------------*/
665 
666 static inline cs_real_t
667 cs_math_sym_33_determinant(const cs_real_6_t m)
668 {
669  const cs_real_t com0 = m[1]*m[2] - m[4]*m[4];
670  const cs_real_t com1 = m[4]*m[5] - m[3]*m[2];
671  const cs_real_t com2 = m[3]*m[4] - m[1]*m[5];
672 
673  return m[0]*com0 + m[3]*com1 + m[5]*com2;
674 }
675 
676 /*----------------------------------------------------------------------------*/
684 /*----------------------------------------------------------------------------*/
685 
686 #if defined(__INTEL_COMPILER)
687 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
688 #endif
689 
690 static inline void
691 cs_math_3_cross_product(const cs_real_t u[3],
692  const cs_real_t v[3],
693  cs_real_t uv[restrict 3])
694 {
695  uv[0] = u[1]*v[2] - u[2]*v[1];
696  uv[1] = u[2]*v[0] - u[0]*v[2];
697  uv[2] = u[0]*v[1] - u[1]*v[0];
698 }
699 
700 /*----------------------------------------------------------------------------*/
710 /*----------------------------------------------------------------------------*/
711 
712 #if defined(__INTEL_COMPILER)
713 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
714 #endif
715 
716 static inline cs_real_t
717 cs_math_3_triple_product(const cs_real_t u[3],
718  const cs_real_t v[3],
719  const cs_real_t w[3])
720 {
721  return (u[1]*v[2] - u[2]*v[1]) * w[0]
722  + (u[2]*v[0] - u[0]*v[2]) * w[1]
723  + (u[0]*v[1] - u[1]*v[0]) * w[2];
724 }
725 
726 /*----------------------------------------------------------------------------*/
733 /*----------------------------------------------------------------------------*/
734 
735 static inline void
736 cs_math_33_inv_cramer(const cs_real_t in[3][3],
737  cs_real_t out[3][3])
738 {
739  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
740  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
741  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
742 
743  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
744  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
745  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
746 
747  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
748  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
749  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
750 
751  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
752  const double invdet = 1/det;
753 
754  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
755  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
756  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
757 }
758 
759 /*----------------------------------------------------------------------------*/
765 /*----------------------------------------------------------------------------*/
766 
767 static inline void
768 cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
769 {
770  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
771  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
772  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
773  cs_real_t a10 = a[2][0]*a[1][2] - a[1][0]*a[2][2];
774  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
775  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
776  cs_real_t a20 = a[1][0]*a[2][1] - a[2][0]*a[1][1];
777  cs_real_t a21 = a[2][0]*a[0][1] - a[0][0]*a[2][1];
778  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
779 
780  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
781 
782  a[0][0] = a00 * det_inv;
783  a[0][1] = a01 * det_inv;
784  a[0][2] = a02 * det_inv;
785  a[1][0] = a10 * det_inv;
786  a[1][1] = a11 * det_inv;
787  a[1][2] = a12 * det_inv;
788  a[2][0] = a20 * det_inv;
789  a[2][1] = a21 * det_inv;
790  a[2][2] = a22 * det_inv;
791 }
792 
793 /*----------------------------------------------------------------------------*/
800 /*----------------------------------------------------------------------------*/
801 
802 static inline void
803 cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
804 {
805  cs_real_t a00 = a[1][1]*a[2][2] - a[2][1]*a[1][2];
806  cs_real_t a01 = a[2][1]*a[0][2] - a[0][1]*a[2][2];
807  cs_real_t a02 = a[0][1]*a[1][2] - a[1][1]*a[0][2];
808  cs_real_t a11 = a[0][0]*a[2][2] - a[2][0]*a[0][2];
809  cs_real_t a12 = a[1][0]*a[0][2] - a[0][0]*a[1][2];
810  cs_real_t a22 = a[0][0]*a[1][1] - a[1][0]*a[0][1];
811 
812  double det_inv = 1. / (a[0][0]*a00 + a[1][0]*a01 + a[2][0]*a02);
813 
814  a[0][0] = a00 * det_inv;
815  a[0][1] = a01 * det_inv;
816  a[0][2] = a02 * det_inv;
817  a[1][0] = a01 * det_inv;
818  a[1][1] = a11 * det_inv;
819  a[1][2] = a12 * det_inv;
820  a[2][0] = a02 * det_inv;
821  a[2][1] = a12 * det_inv;
822  a[2][2] = a22 * det_inv;
823 }
824 
825 /*----------------------------------------------------------------------------*/
835 /*----------------------------------------------------------------------------*/
836 
837 static inline void
838 cs_math_sym_33_inv_cramer(const cs_real_t s[6],
839  cs_real_t sout[restrict 6])
840 {
841  double detinv;
842 
843  sout[0] = s[1]*s[2] - s[4]*s[4];
844  sout[1] = s[0]*s[2] - s[5]*s[5];
845  sout[2] = s[0]*s[1] - s[3]*s[3];
846  sout[3] = s[4]*s[5] - s[3]*s[2];
847  sout[4] = s[3]*s[5] - s[0]*s[4];
848  sout[5] = s[3]*s[4] - s[1]*s[5];
849 
850  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
851 
852  sout[0] *= detinv;
853  sout[1] *= detinv;
854  sout[2] *= detinv;
855  sout[3] *= detinv;
856  sout[4] *= detinv;
857  sout[5] *= detinv;
858 }
859 
860 /*----------------------------------------------------------------------------*/
869 /*----------------------------------------------------------------------------*/
870 
871 static inline void
872 cs_math_33_product(const cs_real_t m1[3][3],
873  const cs_real_t m2[3][3],
874  cs_real_33_t mout)
875 {
876  mout[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
877  mout[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
878  mout[0][2] = m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
879 
880  mout[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
881  mout[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
882  mout[1][2] = m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
883 
884  mout[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
885  mout[2][1] = m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
886  mout[2][2] = m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
887 }
888 
889 /*----------------------------------------------------------------------------*/
898 /*----------------------------------------------------------------------------*/
899 
900 static inline void
901 cs_math_33_product_add(const cs_real_t m1[3][3],
902  const cs_real_t m2[3][3],
903  cs_real_33_t mout)
904 {
905  mout[0][0] += m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
906  mout[0][1] += m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1];
907  mout[0][2] += m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2];
908 
909  mout[1][0] += m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
910  mout[1][1] += m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1];
911  mout[1][2] += m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2];
912 
913  mout[2][0] += m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
914  mout[2][1] += m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1];
915  mout[2][2] += m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2];
916 }
917 
918 
919 /*----------------------------------------------------------------------------*/
932 /*----------------------------------------------------------------------------*/
933 
934 static inline void
935 cs_math_sym_33_product(const cs_real_t s1[6],
936  const cs_real_t s2[6],
937  cs_real_t sout[restrict 6])
938 {
939  /* S11 */
940  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
941  /* S22 */
942  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
943  /* S33 */
944  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
945  /* S12 = S21 */
946  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
947  /* S23 = S32 */
948  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
949  /* S13 = S31 */
950  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
951 }
952 
953 /*----------------------------------------------------------------------------*/
961 /*----------------------------------------------------------------------------*/
962 
963 static inline void
964 cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3],
965  cs_real_t sout[restrict 6][6])
966 {
967  int tens2vect[3][3];
968  int iindex[6], jindex[6];
969 
970  tens2vect[0][0] = 0; tens2vect[0][1] = 3; tens2vect[0][2] = 5;
971  tens2vect[1][0] = 3; tens2vect[1][1] = 1; tens2vect[1][2] = 4;
972  tens2vect[2][0] = 5; tens2vect[2][1] = 4; tens2vect[2][2] = 2;
973 
974  iindex[0] = 0; iindex[1] = 1; iindex[2] = 2;
975  iindex[3] = 0; iindex[4] = 1; iindex[5] = 0;
976 
977  jindex[0] = 0; jindex[1] = 1; jindex[2] = 2;
978  jindex[3] = 1; jindex[4] = 2; jindex[5] = 2;
979 
980  /* Consider : W = R*s^t + s*R.
981  * W_ij = Sum_{k<3} [s_jk*r_ik + s_ik*r_jk]
982  * We look for A such as A*R = W
983  */
984  for (int i = 0; i < 6; i++) {
985  int ii = iindex[i];
986  int jj = jindex[i];
987  for (int k = 0; k < 3; k++) {
988  int ik = tens2vect[k][ii];
989  int jk = tens2vect[k][jj];
990 
991  sout[ik][i] += s[k][jj];
992 
993  sout[jk][i] += s[k][ii];
994  }
995  }
996 }
997 
998 /*----------------------------------------------------------------------------*/
1010 /*----------------------------------------------------------------------------*/
1011 
1012 static inline void
1013 cs_math_sym_33_double_product(const cs_real_t s1[6],
1014  const cs_real_t s2[6],
1015  const cs_real_t s3[6],
1016  cs_real_t sout[restrict 3][3])
1017 {
1018  cs_real_33_t _sout;
1019 
1020  /* S11 */
1021  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
1022  /* S22 */
1023  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
1024  /* S33 */
1025  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
1026  /* S12 */
1027  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
1028  /* S21 */
1029  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
1030  /* S23 */
1031  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
1032  /* S32 */
1033  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
1034  /* S13 */
1035  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
1036  /* S31 */
1037  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
1038 
1039  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
1040  /* S22 */
1041  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
1042  /* S33 */
1043  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
1044  /* S12 */
1045  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
1046  /* S21 */
1047  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
1048  /* S23 */
1049  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
1050  /* S32 */
1051  sout[2][1] = s3[3]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
1052  /* S13 */
1053  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
1054  /* S31 */
1055  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
1056 }
1057 
1058 /*----------------------------------------------------------------------------*/
1065 /*----------------------------------------------------------------------------*/
1066 
1067 static inline void
1068 cs_nvec3(const cs_real_3_t v,
1069  cs_nvec3_t *qv)
1070 {
1071  cs_real_t magnitude = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1072 
1073  qv->meas = magnitude;
1074  if (fabs(magnitude) > cs_math_zero_threshold) {
1075 
1076  const cs_real_t inv = 1/magnitude;
1077  qv->unitv[0] = inv * v[0];
1078  qv->unitv[1] = inv * v[1];
1079  qv->unitv[2] = inv * v[2];
1080 
1081  }
1082  else
1083  qv->unitv[0] = qv->unitv[1] = qv->unitv[2] = 0;
1084 
1085 }
1086 
1087 /*=============================================================================
1088  * Public function prototypes
1089  *============================================================================*/
1090 
1091 /*----------------------------------------------------------------------------*/
1095 /*----------------------------------------------------------------------------*/
1096 
1097 void
1098 cs_math_set_machine_epsilon(void);
1099 
1100 /*----------------------------------------------------------------------------*/
1104 /*----------------------------------------------------------------------------*/
1105 
1106 double
1107 cs_math_get_machine_epsilon(void);
1108 
1109 /*----------------------------------------------------------------------------*/
1119 /*----------------------------------------------------------------------------*/
1120 
1121 void
1122 cs_math_3_length_unitv(const cs_real_t xa[3],
1123  const cs_real_t xb[3],
1124  cs_real_t *len,
1125  cs_real_3_t unitv);
1126 
1127 /*----------------------------------------------------------------------------*/
1139 /*----------------------------------------------------------------------------*/
1140 
1141 void
1142 cs_math_sym_33_eigen(const cs_real_t m[6],
1143  cs_real_t eig_vals[3]);
1144 
1145 /*----------------------------------------------------------------------------*/
1158 /*----------------------------------------------------------------------------*/
1159 
1160 void
1161 cs_math_33_eigen(const cs_real_t m[3][3],
1162  cs_real_t *eig_ratio,
1163  cs_real_t *eig_max);
1164 
1165 /*----------------------------------------------------------------------------*/
1176 /*----------------------------------------------------------------------------*/
1177 
1178 double
1179 cs_math_surftri(const cs_real_t xv[3],
1180  const cs_real_t xe[3],
1181  const cs_real_t xf[3]);
1182 
1183 /*----------------------------------------------------------------------------*/
1195 /*----------------------------------------------------------------------------*/
1196 
1197 double
1198 cs_math_voltet(const cs_real_t xv[3],
1199  const cs_real_t xe[3],
1200  const cs_real_t xf[3],
1201  const cs_real_t xc[3]);
1202 
1203 /*----------------------------------------------------------------------------*/
1213 /*----------------------------------------------------------------------------*/
1214 
1215 void
1216 cs_math_fact_lu(cs_lnum_t n_blocks,
1217  const int b_size,
1218  const cs_real_t *a,
1219  cs_real_t *a_lu);
1220 
1221 /*----------------------------------------------------------------------------*/
1231 /*----------------------------------------------------------------------------*/
1232 
1233 void
1234 cs_math_fw_and_bw_lu(const cs_real_t a_lu[],
1235  const int n,
1236  cs_real_t x[],
1237  const cs_real_t b[]);
1238 
1239 /*----------------------------------------------------------------------------*/
1240 
1241 END_C_DECLS
1242 
1243 #endif /* __CS_MATH_H__ */
cs_math_pow4
static cs_real_t cs_math_pow4(cs_real_t x)
Compute the 4-th power of a real value.
Definition: cs_math.h:241
cs_math_pi
const cs_real_t cs_math_pi
cs_math_3_norm
static cs_real_t cs_math_3_norm(const cs_real_t v[3])
Compute the euclidean norm of a vector of dimension 3.
Definition: cs_math.h:372
cs_defs.h
cs_math_1ov6
const cs_real_t cs_math_1ov6
cs_math_33t_3_product
static void cs_math_33t_3_product(const cs_real_t m[3][3], const cs_real_t v[3], cs_real_3_t 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:539
bft_error.h
cs_math_33_3_product_add
static void cs_math_33_3_product_add(const cs_real_t m[3][3], const cs_real_t v[3], cs_real_3_t mv)
Compute the product of a matrix of 3x3 real values by a vector of 3 real values add.
Definition: cs_math.h:518
len
size_t len
Definition: mei_scanner.c:560
cs_fuel_incl::a
double precision, save a
Definition: cs_fuel_incl.f90:146
cs_math_3_square_norm
static cs_real_t cs_math_3_square_norm(const cs_real_t v[3])
Compute the square norm of a vector of 3 real values.
Definition: cs_math.h:388
cs_math_sym_33_product
static void cs_math_sym_33_product(const cs_real_t s1[6], const cs_real_t s2[6], cs_real_t sout[restrict 6])
Compute the product of two symmetric matrices. Warning: this is valid if and only if s1 and s2 commut...
Definition: cs_math.h:935
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:218
restrict
#define restrict
Definition: cs_defs.h:127
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:315
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:482
cs_math_binom
static int cs_math_binom(int n, int k)
Computes the binomial coefficient of n and k.
Definition: cs_math.h:107
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:146
XX
Definition: cs_math.h:63
cs_math_3_cross_product
static void cs_math_3_cross_product(const cs_real_t u[3], const cs_real_t v[3], cs_real_t uv[restrict 3])
Compute the cross product of two vectors of 3 real values.
Definition: cs_math.h:691
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:468
XZ
Definition: cs_math.h:68
atimbr::v
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:302
cs_math_3_distance
static cs_real_t cs_math_3_distance(const cs_real_t xa[3], const cs_real_t xb[3])
Compute the (euclidean) distance between two points xa and xb in a cartesian coordinate system of dim...
Definition: cs_math.h:259
cs_math_sq
static cs_real_t cs_math_sq(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:193
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:467
cs_math_33_inv_cramer
static void cs_math_33_inv_cramer(const cs_real_t in[3][3], cs_real_t out[3][3])
Inverse a 3x3 matrix.
Definition: cs_math.h:736
cs_math_big_r
const cs_real_t cs_math_big_r
cs_math_3_distance_dot_product
static cs_real_t cs_math_3_distance_dot_product(const cs_real_t xa[3], const cs_real_t xb[3], const cs_real_t xc[3])
Compute .
Definition: cs_math.h:284
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:423
cs_math_get_machine_epsilon
double cs_math_get_machine_epsilon(void)
Get the value related to the machine precision.
Definition: cs_math.c:198
cs_math_3_normal_scaling
static void cs_math_3_normal_scaling(const cs_real_t n[3], cs_real_t factor, cs_real_3_t v)
Add the dot product with a normal vector to the normal direction to a vector.
Definition: cs_math.h:450
cs_math_33_determinant
static cs_real_t cs_math_33_determinant(const cs_real_t m[3][3])
Compute the determinant of a 3x3 matrix.
Definition: cs_math.h:647
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
cs_math_fmin
static cs_real_t cs_math_fmin(cs_real_t x, cs_real_t y)
Compute the min value of two real values.
Definition: cs_math.h:155
XY
Definition: cs_math.h:66
cs_math_33_product_add
static void cs_math_33_product_add(const cs_real_t m1[3][3], const cs_real_t m2[3][3], cs_real_33_t mout)
Compute the product of a matrix of 3x3 real values by a matrix of 3x3 real values and add.
Definition: cs_math.h:901
cs_math_reduce_sym_prod_33_to_66
static void cs_math_reduce_sym_prod_33_to_66(const cs_real_t s[3][3], cs_real_t sout[restrict 6][6])
Compute a 6x6 matrix A, equivalent to a 3x3 matrix s, such as: A*R_6 = R*s^t + s*R.
Definition: cs_math.h:964
cs_math_epzero
const cs_real_t cs_math_epzero
cs_math_sym_33_3_product
static void cs_math_sym_33_3_product(const cs_real_t m[6], const cs_real_t v[3], cs_real_t mv[restrict 3])
Compute the product of a symmetric matrix of 3x3 real values by a vector of 3 real values....
Definition: cs_math.h:561
y
void const cs_lnum_t *const const cs_real_t *const const cs_real_t *const const cs_real_t *const const cs_real_t *const y
Definition: cs_wall_functions.h:1147
cs_nvec3_t::unitv
double unitv[3]
Definition: cs_defs.h:346
cs_math_infinite_r
const cs_real_t cs_math_infinite_r
cs_math_sym_33_determinant
static cs_real_t cs_math_sym_33_determinant(const cs_real_6_t m)
Compute the determinant of a 3x3 symmetric matrix.
Definition: cs_math.h:667
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:317
cs_math_66_6_product
static void cs_math_66_6_product(const cs_real_t m[6][6], const cs_real_t v[6], cs_real_t mv[restrict 6])
Compute the product of a matrix of 6x6 real values by a vector of 6 real values.
Definition: cs_math.h:604
cs_nvec3
static void cs_nvec3(const cs_real_3_t v, cs_nvec3_t *qv)
Define a cs_nvec3_t structure from a cs_real_3_t.
Definition: cs_math.h:1068
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:298
numvar::ik
integer, save ik
Definition: numvar.f90:75
cs_math_3_orthogonal_projection
static void cs_math_3_orthogonal_projection(const cs_real_t n[3], const cs_real_t v[3], cs_real_t vout[restrict 3])
Orthogonal projection of a vector with respect to a normalised vector.
Definition: cs_math.h:429
cs_math_66_6_product_add
static void cs_math_66_6_product_add(const cs_real_t m[6][6], const cs_real_t v[6], cs_real_t mv[restrict 6])
Compute the product of a matrix of 6x6 real values by a vector of 6 real values and add it to the vec...
Definition: cs_math.h:626
cs_math_3_square_distance
static cs_real_t cs_math_3_square_distance(const cs_real_t xa[3], const cs_real_t xb[3])
Compute the squared distance between two points xa and xb in a cartesian coordinate system of dimensi...
Definition: cs_math.h:304
cs_math_1ov12
const cs_real_t cs_math_1ov12
cs_math_sym_33_double_product
static void cs_math_sym_33_double_product(const cs_real_t s1[6], const cs_real_t s2[6], const cs_real_t s3[6], cs_real_t sout[restrict 3][3])
Compute the product of three symmetric matrices.
Definition: cs_math.h:1013
cs_math_1ov3
const cs_real_t cs_math_1ov3
cs_math_pow3
static cs_real_t cs_math_pow3(cs_real_t x)
Compute the cube of a real value.
Definition: cs_math.h:225
cs_math_33_inv_cramer_in_place
static void cs_math_33_inv_cramer_in_place(cs_real_t a[3][3])
Inverse a 3x3 matrix in place, using Cramer's rule.
Definition: cs_math.h:768
cs_math_3_dot_product
static cs_real_t cs_math_3_dot_product(const cs_real_t u[3], const cs_real_t v[3])
Compute the dot product of two vectors of 3 real values.
Definition: cs_math.h:326
cs_math_fmax
static cs_real_t cs_math_fmax(cs_real_t x, cs_real_t y)
Compute the max value of two real values.
Definition: cs_math.h:174
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:321
cs_math_3_normalise
static void cs_math_3_normalise(const cs_real_t vin[3], cs_real_t vout[restrict 3])
Normalise a vector of 3 real values.
Definition: cs_math.h:405
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:303
cs_math_5ov3
const cs_real_t cs_math_5ov3
t
Definition: cs_field_pointer.h:98
cs_math_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition: cs_math.h:61
cs_math_sym_33_3_product_add
static void cs_math_sym_33_3_product_add(const cs_real_t m[6], const cs_real_t v[3], cs_real_t mv[restrict 3])
Compute the product of a symmetric matrix of 3x3 real values by a vector of 3 real values and add it ...
Definition: cs_math.h:583
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:540
cs_math_33_3_product
static void cs_math_33_3_product(const cs_real_t m[3][3], const cs_real_t v[3], cs_real_3_t mv)
Compute the product of a matrix of 3x3 real values by a vector of 3 real values.
Definition: cs_math.h:497
cs_math_fabs
static cs_real_t cs_math_fabs(cs_real_t x)
Compute the absolute value of a real value.
Definition: cs_math.h:137
cs_math_4ov3
const cs_real_t cs_math_4ov3
atimbr::u
double precision, dimension(:,:,:), allocatable u
Definition: atimbr.f90:113
cs_math_sym_33_inv_cramer
static void cs_math_sym_33_inv_cramer(const cs_real_t s[6], cs_real_t sout[restrict 6])
Compute the inverse of a symmetric matrix using Cramer's rule.
Definition: cs_math.h:838
cs_math_33_inv_cramer_sym_in_place
static void cs_math_33_inv_cramer_sym_in_place(cs_real_t a[3][3])
Inverse a 3x3 symmetric matrix (with non-symmetric storage) in place, using Cramer's rule.
Definition: cs_math.h:803
cs_math_1ov24
const cs_real_t cs_math_1ov24
cs_nvec3_t::meas
double meas
Definition: cs_defs.h:345
cs_math_33_normal_scaling_add
static void cs_math_33_normal_scaling_add(const cs_real_t n[3], cs_real_t factor, cs_real_33_t t)
Add the dot product with a normal vector to the normal,normal component of a tensor: t += factor * n....
Definition: cs_math.h:472
YY
Definition: cs_math.h:64
cs_nvec3_t
Definition: cs_defs.h:343
cs_math_2ov3
const cs_real_t cs_math_2ov3
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:390
cs_math_3_triple_product
static cs_real_t cs_math_3_triple_product(const cs_real_t u[3], const cs_real_t v[3], const cs_real_t w[3])
Compute the triple product.
Definition: cs_math.h:717
cs_math_set_machine_epsilon
void cs_math_set_machine_epsilon(void)
Compute the value related to the machine precision.
Definition: cs_math.c:177
ZZ
Definition: cs_math.h:65
cs_math_pow2
static cs_real_t cs_math_pow2(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:209
YZ
Definition: cs_math.h:67
cs_math_3_33_3_dot_product
static cs_real_t cs_math_3_33_3_dot_product(const cs_real_t n1[3], const cs_real_t t[3][3], const cs_real_t n2[3])
Compute the dot product of a tensor t with two vectors n1, and n2 n1 t n2.
Definition: cs_math.h:349
xa
void const int const int const int const cs_real_t const int const cs_real_t const cs_real_t const cs_real_t const cs_real_t const cs_real_t const cs_real_t const cs_real_t const cs_real_t cs_real_t cs_real_t xa[]
Definition: cs_matrix_building.h:65
cs_math_33_product
static void cs_math_33_product(const cs_real_t m1[3][3], const cs_real_t m2[3][3], cs_real_33_t mout)
Compute the product of a matrix of 3x3 real values by a matrix of 3x3 real values.
Definition: cs_math.h:872
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:453
k
Definition: cs_field_pointer.h:70
Generated on Thu Apr 9 2020 17:26:56 for My Project by   doxygen 1.8.16