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-2016 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 /*----------------------------------------------------------------------------
31  * Local headers
32  *----------------------------------------------------------------------------*/
33 
34 #include "cs_defs.h"
35 
36 /*----------------------------------------------------------------------------
37  * Standard C library headers
38  *----------------------------------------------------------------------------*/
39 
40 #include <math.h>
41 
42 /*----------------------------------------------------------------------------*/
43 
45 
46 /*=============================================================================
47  * Local Macro definitions
48  *============================================================================*/
49 
50 /*============================================================================
51  * Type definition
52  *============================================================================*/
53 
54 /*============================================================================
55  * Global variables
56  *============================================================================*/
57 
58 /* Numerical constants */
59 
61 extern const cs_real_t cs_math_onethird;
62 extern const cs_real_t cs_math_onesix;
63 extern const cs_real_t cs_math_onetwelve;
64 extern const cs_real_t cs_math_epzero;
65 extern const cs_real_t cs_math_infinite_r;
66 extern const cs_real_t cs_math_big_r;
67 extern const cs_real_t cs_math_pi;
68 
69 /*============================================================================
70  * Public function prototypes for Fortran API
71  *============================================================================*/
72 
73 /*----------------------------------------------------------------------------
74  * Wrapper to cs_math_sym_33_inv_cramer
75  *----------------------------------------------------------------------------*/
76 
77 void CS_PROCF (symmetric_matrix_inverse, SYMMETRIC_MATRIX_INVERSE)
78 (
79  const cs_real_6_t s,
80  cs_real_6_t sout
81 );
82 
83 /*----------------------------------------------------------------------------
84  * Wrapper to cs_math_sym_33_product
85  *----------------------------------------------------------------------------*/
86 
87 void CS_PROCF (symmetric_matrix_product, SYMMETRIC_MATRIX_PRODUCT)
88 (
89  const cs_real_6_t s1,
90  const cs_real_6_t s2,
91  cs_real_6_t sout
92 );
93 
94 /*=============================================================================
95  * Inline static function prototypes
96  *============================================================================*/
97 
98 /*----------------------------------------------------------------------------*/
106 /*----------------------------------------------------------------------------*/
107 
108 static inline cs_real_t
110 {
111  return x*x;
112 }
113 
114 /*----------------------------------------------------------------------------*/
124 /*----------------------------------------------------------------------------*/
125 
126 static inline cs_real_t
128  const cs_real_t xb[3])
129 {
130  cs_real_3_t v;
131 
132  v[0] = xb[0] - xa[0];
133  v[1] = xb[1] - xa[1];
134  v[2] = xb[2] - xa[2];
135 
136  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
137 }
138 
139 /*----------------------------------------------------------------------------*/
148 /*----------------------------------------------------------------------------*/
149 
150 static inline cs_real_t
152  const cs_real_t v[3])
153 {
154  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
155 
156  return uv;
157 }
158 
159 /*----------------------------------------------------------------------------*/
167 /*----------------------------------------------------------------------------*/
168 
169 static inline cs_real_t
171 {
172  return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
173 }
174 
175 /*----------------------------------------------------------------------------*/
183 /*----------------------------------------------------------------------------*/
184 
185 static inline cs_real_t
187 {
188  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
189 
190  return v2;
191 }
192 
193 /*----------------------------------------------------------------------------*/
202 /*----------------------------------------------------------------------------*/
203 
204 static inline void
206  const cs_real_t v[3],
207  cs_real_3_t mv)
208 {
209  mv[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
210  mv[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
211  mv[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
212 }
213 
214 /*----------------------------------------------------------------------------*/
223 /*----------------------------------------------------------------------------*/
224 
225 static inline void
227  const cs_real_t v[3],
228  cs_real_3_t mv)
229 {
230  mv[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
231  mv[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
232  mv[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
233 }
234 
235 /*----------------------------------------------------------------------------*/
245 /*----------------------------------------------------------------------------*/
246 
247 static inline void
249  const cs_real_t v[3],
250  cs_real_t mv[restrict 3])
251 {
252  mv[0] = m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
253  mv[1] = m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
254  mv[2] = m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
255 }
256 
257 /*----------------------------------------------------------------------------*/
265 /*----------------------------------------------------------------------------*/
266 
267 static inline cs_real_t
269 {
270  const cs_real_t com0 = m[1][1]*m[2][2] - m[2][1]*m[1][2];
271  const cs_real_t com1 = m[2][1]*m[0][2] - m[0][1]*m[2][2];
272  const cs_real_t com2 = m[0][1]*m[1][2] - m[1][1]*m[0][2];
273 
274  return m[0][0]*com0 + m[1][0]*com1 + m[2][0]*com2;
275 }
276 
277 /*----------------------------------------------------------------------------*/
285 /*----------------------------------------------------------------------------*/
286 
287 #if defined(__INTEL_COMPILER)
288 #pragma optimization_level 0 /* Bug with O1 or above with icc 15.0.1 20141023 */
289 #endif
290 
291 static inline void
293  const cs_real_t v[3],
294  cs_real_t uv[restrict 3])
295 {
296  uv[0] = u[1]*v[2] - u[2]*v[1];
297  uv[1] = u[2]*v[0] - u[0]*v[2];
298  uv[2] = u[0]*v[1] - u[1]*v[0];
299 }
300 
301 /*----------------------------------------------------------------------------*/
308 /*----------------------------------------------------------------------------*/
309 
310 static inline void
311 cs_math_33_inv(const cs_real_t in[3][3],
312  cs_real_t out[3][3])
313 {
314  out[0][0] = in[1][1]*in[2][2] - in[2][1]*in[1][2];
315  out[0][1] = in[2][1]*in[0][2] - in[0][1]*in[2][2];
316  out[0][2] = in[0][1]*in[1][2] - in[1][1]*in[0][2];
317 
318  out[1][0] = in[2][0]*in[1][2] - in[1][0]*in[2][2];
319  out[1][1] = in[0][0]*in[2][2] - in[2][0]*in[0][2];
320  out[1][2] = in[1][0]*in[0][2] - in[0][0]*in[1][2];
321 
322  out[2][0] = in[1][0]*in[2][1] - in[2][0]*in[1][1];
323  out[2][1] = in[2][0]*in[0][1] - in[0][0]*in[2][1];
324  out[2][2] = in[0][0]*in[1][1] - in[1][0]*in[0][1];
325 
326  const double det = in[0][0]*out[0][0]+in[1][0]*out[0][1]+in[2][0]*out[0][2];
327  const double invdet = 1/det;
328 
329  out[0][0] *= invdet, out[0][1] *= invdet, out[0][2] *= invdet;
330  out[1][0] *= invdet, out[1][1] *= invdet, out[1][2] *= invdet;
331  out[2][0] *= invdet, out[2][1] *= invdet, out[2][2] *= invdet;
332 }
333 
334 /*----------------------------------------------------------------------------*/
344 /*----------------------------------------------------------------------------*/
345 
346 static inline void
348  cs_real_t sout[restrict 6])
349 {
350  double detinv;
351 
352  sout[0] = s[1]*s[2] - s[4]*s[4];
353  sout[1] = s[0]*s[2] - s[5]*s[5];
354  sout[2] = s[0]*s[1] - s[3]*s[3];
355  sout[3] = s[4]*s[5] - s[3]*s[2];
356  sout[4] = s[3]*s[5] - s[0]*s[4];
357  sout[5] = s[3]*s[4] - s[1]*s[5];
358 
359  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
360 
361  sout[0] *= detinv;
362  sout[1] *= detinv;
363  sout[2] *= detinv;
364  sout[3] *= detinv;
365  sout[4] *= detinv;
366  sout[5] *= detinv;
367 }
368 
369 /*----------------------------------------------------------------------------*/
380 /*----------------------------------------------------------------------------*/
381 
382 static inline void
384  const cs_real_t s2[6],
385  cs_real_t sout[restrict 6])
386 {
387  /* S11 */
388  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
389  /* S22 */
390  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
391  /* S33 */
392  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
393  /* S12 = S21 */
394  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
395  /* S23 = S32 */
396  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
397  /* S13 = S31 */
398  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
399 }
400 
401 /*----------------------------------------------------------------------------*/
413 /*----------------------------------------------------------------------------*/
414 
415 static inline void
417  const cs_real_t s2[6],
418  const cs_real_t s3[6],
419  cs_real_t sout[restrict 3][3])
420 {
421  cs_real_33_t _sout;
422 
423  /* S11 */
424  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
425  /* S22 */
426  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
427  /* S33 */
428  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
429  /* S12 */
430  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
431  /* S21 */
432  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
433  /* S23 */
434  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
435  /* S32 */
436  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
437  /* S13 */
438  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
439  /* S31 */
440  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
441 
442  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
443  /* S22 */
444  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
445  /* S33 */
446  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
447  /* S12 */
448  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
449  /* S21 */
450  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
451  /* S23 */
452  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
453  /* S32 */
454  sout[2][1] = s3[2]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
455  /* S13 */
456  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
457  /* S31 */
458  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
459 }
460 
461 /*=============================================================================
462  * Public function prototypes
463  *============================================================================*/
464 
465 /*----------------------------------------------------------------------------*/
469 /*----------------------------------------------------------------------------*/
470 
471 void
473 
474 /*----------------------------------------------------------------------------*/
478 /*----------------------------------------------------------------------------*/
479 
480 double
482 
483 /*----------------------------------------------------------------------------*/
493 /*----------------------------------------------------------------------------*/
494 
495 void
497  const cs_real_t xb[3],
498  cs_real_t *len,
499  cs_real_3_t unitv);
500 
501 /*----------------------------------------------------------------------------*/
512 /*----------------------------------------------------------------------------*/
513 
514 void
515 cs_math_33_eigen(const cs_real_t m[3][3],
516  cs_real_t *eig_ratio,
517  cs_real_t *eig_max);
518 
519 /*----------------------------------------------------------------------------*/
530 /*----------------------------------------------------------------------------*/
531 
532 double
533 cs_math_surftri(const cs_real_t xv[3],
534  const cs_real_t xe[3],
535  const cs_real_t xf[3]);
536 
537 /*----------------------------------------------------------------------------*/
549 /*----------------------------------------------------------------------------*/
550 
551 double
552 cs_math_voltet(const cs_real_t xv[3],
553  const cs_real_t xe[3],
554  const cs_real_t xf[3],
555  const cs_real_t xc[3]);
556 
557 /*----------------------------------------------------------------------------*/
558 
560 
561 #endif /* __CS_MATH_H__ */
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&#39;s rule.
Definition: cs_math.h:347
void symmetric_matrix_inverse(const cs_real_6_t s, cs_real_6_t sout)
Definition: cs_math.c:114
#define restrict
Definition: cs_defs.h:122
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:151
const cs_real_t cs_math_onesix
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:310
size_t len
Definition: mei_scanner.c:600
const cs_real_t cs_math_big_r
void symmetric_matrix_product(const cs_real_6_t s1, const cs_real_6_t s2, cs_real_6_t sout)
Definition: cs_math.c:128
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:205
const cs_real_t cs_math_pi
static void cs_math_33_inv(const cs_real_t in[3][3], cs_real_t out[3][3])
Inverse a 3x3 matrix.
Definition: cs_math.h:311
#define BEGIN_C_DECLS
Definition: cs_defs.h:448
const cs_real_t cs_math_epzero
void cs_math_33_eigen(const cs_real_t m[3][3], cs_real_t *eig_ratio, cs_real_t *eig_max)
Compute the eigenvalues of a 3x3 matrix which is symmetric and real -> Oliver K. Smith "eigenvalues o...
Definition: cs_math.c:190
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:416
double cs_real_t
Floating-point value.
Definition: cs_defs.h:296
static cs_real_t cs_math_3_length(const cs_real_t xa[3], const cs_real_t xb[3])
Compute the (euclidean) length between two points xa and xb in a cartesian coordinate system of dimen...
Definition: cs_math.h:127
Definition: cs_field_pointer.h:67
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:186
static cs_real_t cs_math_sq(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:109
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:308
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
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.
Definition: cs_math.h:383
double cs_math_get_machine_epsilon(void)
Get the value related to the machine precision.
Definition: cs_math.c:171
const cs_real_t cs_math_onetwelve
void cs_math_set_machine_epsilon(void)
Compute the value related to the machine precision.
Definition: cs_math.c:150
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:338
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:268
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:170
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:308
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:292
const cs_real_t cs_math_zero_threshold
const cs_real_t cs_math_onethird
#define END_C_DECLS
Definition: cs_defs.h:449
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:226
#define CS_PROCF(x, y)
Definition: cs_defs.h:472
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:314
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:248
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:275
const cs_real_t cs_math_infinite_r