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-2015 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 
39 
40 /*=============================================================================
41  * Local Macro definitions
42  *============================================================================*/
43 
44 /*============================================================================
45  * Type definition
46  *============================================================================*/
47 
48 /*============================================================================
49  * Global variables
50  *============================================================================*/
51 
52 /*============================================================================
53  * Public function prototypes for Fortran API
54  *============================================================================*/
55 
56 /*----------------------------------------------------------------------------
57  * Wrapper to cs_math_sym_33_inv_cramer
58  *----------------------------------------------------------------------------*/
59 
60 void CS_PROCF (symmetric_matrix_inverse, SYMMETRIC_MATRIX_INVERSE)
61 (
62  const cs_real_6_t s,
63  cs_real_6_t sout
64 );
65 
66 /*----------------------------------------------------------------------------
67  * Wrapper to cs_math_sym_33_product
68  *----------------------------------------------------------------------------*/
69 
70 void CS_PROCF (symmetric_matrix_product, SYMMETRIC_MATRIX_PRODUCT)
71 (
72  const cs_real_6_t s1,
73  const cs_real_6_t s2,
74  cs_real_6_t sout
75 );
76 
77 /*=============================================================================
78  * Public function prototypes
79  *============================================================================*/
80 
81 /*----------------------------------------------------------------------------*/
91 /*----------------------------------------------------------------------------*/
92 
93 static inline void
95  const cs_real_t v[3],
96  cs_real_3_t mv)
97 {
98  for (int ii = 0; ii < 3; ii++)
99  mv[ii] = m[ii][0] * v[0] + m[ii][1] * v[1] + m[ii][2] * v[2];
100 }
101 
102 /*----------------------------------------------------------------------------*/
113 /*----------------------------------------------------------------------------*/
114 
115 static inline void
117  const cs_real_t v[3],
118  cs_real_t mv[restrict 3])
119 {
120  mv[0] = m[0] * v[0] + m[3] * v[1] + m[5] * v[2];
121  mv[1] = m[3] * v[0] + m[1] * v[1] + m[4] * v[2];
122  mv[2] = m[5] * v[0] + m[4] * v[1] + m[2] * v[2];
123 }
124 
125 /*----------------------------------------------------------------------------*/
134 /*----------------------------------------------------------------------------*/
135 
136 static inline cs_real_t
138  const cs_real_t v[3])
139 {
140  cs_real_t uv = u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
141 
142  return uv;
143 }
144 
145 /*----------------------------------------------------------------------------*/
153 /*----------------------------------------------------------------------------*/
154 
155 static inline cs_real_t
157 {
158  cs_real_t v2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
159 
160  return v2;
161 }
162 
163 /*----------------------------------------------------------------------------*/
173 /*----------------------------------------------------------------------------*/
174 
175 static inline void
177  cs_real_t sout[restrict 6])
178 {
179  double detinv;
180 
181  sout[0] = s[1]*s[2] - s[4]*s[4];
182  sout[1] = s[0]*s[2] - s[5]*s[5];
183  sout[2] = s[0]*s[1] - s[3]*s[3];
184  sout[3] = s[4]*s[5] - s[3]*s[2];
185  sout[4] = s[3]*s[5] - s[0]*s[4];
186  sout[5] = s[3]*s[4] - s[1]*s[5];
187 
188  detinv = 1. / (s[0]*sout[0] + s[3]*sout[3] + s[5]*sout[5]);
189 
190  sout[0] *= detinv;
191  sout[1] *= detinv;
192  sout[2] *= detinv;
193  sout[3] *= detinv;
194  sout[4] *= detinv;
195  sout[5] *= detinv;
196 }
197 
198 /*----------------------------------------------------------------------------*/
209 /*----------------------------------------------------------------------------*/
210 
211 static inline void
213  const cs_real_t s2[6],
214  cs_real_t sout[restrict 6])
215 {
216  /* S11 */
217  sout[0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
218  /* S22 */
219  sout[1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
220  /* S33 */
221  sout[2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
222  /* S12 = S21 */
223  sout[3] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
224  /* S23 = S32 */
225  sout[4] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
226  /* S13 = S31 */
227  sout[5] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
228 }
229 
230 /*----------------------------------------------------------------------------*/
242 /*----------------------------------------------------------------------------*/
243 
244 static inline void
246  const cs_real_t s2[6],
247  const cs_real_t s3[6],
248  cs_real_t sout[restrict 3][3])
249 {
250  cs_real_33_t _sout;
251 
252  /* S11 */
253  _sout[0][0] = s1[0]*s2[0] + s1[3]*s2[3] + s1[5]*s2[5];
254  /* S22 */
255  _sout[1][1] = s1[3]*s2[3] + s1[1]*s2[1] + s1[4]*s2[4];
256  /* S33 */
257  _sout[2][2] = s1[5]*s2[5] + s1[4]*s2[4] + s1[2]*s2[2];
258  /* S12 */
259  _sout[0][1] = s1[0]*s2[3] + s1[3]*s2[1] + s1[5]*s2[4];
260  /* S21 */
261  _sout[1][0] = s2[0]*s1[3] + s2[3]*s1[1] + s2[5]*s1[4];
262  /* S23 */
263  _sout[1][2] = s1[3]*s2[5] + s1[1]*s2[4] + s1[4]*s2[2];
264  /* S32 */
265  _sout[2][1] = s2[3]*s1[5] + s2[1]*s1[4] + s2[4]*s1[2];
266  /* S13 */
267  _sout[0][2] = s1[0]*s2[5] + s1[3]*s2[4] + s1[5]*s2[2];
268  /* S31 */
269  _sout[2][0] = s2[0]*s1[5] + s2[3]*s1[4] + s2[5]*s1[2];
270 
271  sout[0][0] = _sout[0][0]*s3[0] + _sout[0][1]*s3[3] + _sout[0][2]*s3[5];
272  /* S22 */
273  sout[1][1] = _sout[1][0]*s3[3] + _sout[1][1]*s3[1] + _sout[1][2]*s3[4];
274  /* S33 */
275  sout[2][2] = _sout[2][0]*s3[5] + _sout[2][1]*s3[4] + _sout[2][2]*s3[2];
276  /* S12 */
277  sout[0][1] = _sout[0][0]*s3[3] + _sout[0][1]*s3[1] + _sout[0][2]*s3[4];
278  /* S21 */
279  sout[1][0] = s3[0]*_sout[1][0] + s3[3]*_sout[1][1] + s3[5]*_sout[1][2];
280  /* S23 */
281  sout[1][2] = _sout[1][0]*s3[5] + _sout[1][1]*s3[4] + _sout[1][2]*s3[2];
282  /* S32 */
283  sout[2][1] = s3[2]*_sout[2][0] + s3[1]*_sout[2][1] + s3[4]*_sout[2][2];
284  /* S13 */
285  sout[0][2] = _sout[0][0]*s3[5] + _sout[0][1]*s3[4] + _sout[0][2]*s3[2];
286  /* S31 */
287  sout[2][0] = s3[0]*_sout[2][0] + s3[3]*_sout[2][1] + s3[5]*_sout[2][2];
288 }
289 
290 /*----------------------------------------------------------------------------*/
291 
293 
294 #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's rule.
Definition: cs_math.h:176
void symmetric_matrix_inverse(const cs_real_6_t s, cs_real_6_t sout)
Definition: cs_math.c:92
#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:137
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:309
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:106
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:94
#define BEGIN_C_DECLS
Definition: cs_defs.h:429
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:245
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:156
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:212
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:307
#define END_C_DECLS
Definition: cs_defs.h:430
double cs_real_t
Definition: cs_defs.h:296
#define CS_PROCF(x, y)
Definition: cs_defs.h:453
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:311
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:116