My Project
programmer's documentation
Loading...
Searching...
No Matches
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
49BEGIN_C_DECLS
50
51/*=============================================================================
52 * Local Macro definitions
53 *============================================================================*/
54
55/*============================================================================
56 * Type definition
57 *============================================================================*/
58
59/* Symetric tensor component name */
60
71
72/*============================================================================
73 * Global variables
74 *============================================================================*/
75
76/* Numerical constants */
77
78extern const cs_real_t cs_math_zero_threshold;
79extern const cs_real_t cs_math_1ov3;
80extern const cs_real_t cs_math_2ov3;
81extern const cs_real_t cs_math_4ov3;
82extern const cs_real_t cs_math_5ov3;
83extern const cs_real_t cs_math_1ov6;
84extern const cs_real_t cs_math_1ov12;
85extern const cs_real_t cs_math_1ov24;
86extern const cs_real_t cs_math_epzero;
87extern const cs_real_t cs_math_infinite_r;
88extern const cs_real_t cs_math_big_r;
89extern const cs_real_t cs_math_pi;
90
91/*=============================================================================
92 * Inline static function prototypes
93 *============================================================================*/
94
95/*----------------------------------------------------------------------------*/
104/*----------------------------------------------------------------------------*/
105
106static inline int
107cs_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
136static inline cs_real_t
137cs_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
154static inline cs_real_t
155cs_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
173static inline cs_real_t
174cs_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
192static inline cs_real_t
193cs_math_sq(cs_real_t x)
194{
195 return x*x;
196}
197
198/*----------------------------------------------------------------------------*/
206/*----------------------------------------------------------------------------*/
207
208static inline cs_real_t
209cs_math_pow2(cs_real_t x)
210{
211 return x*x;
212}
213
214/*----------------------------------------------------------------------------*/
222/*----------------------------------------------------------------------------*/
223
224static inline cs_real_t
225cs_math_pow3(cs_real_t x)
226{
227 return x*x*x;
228}
229
230/*----------------------------------------------------------------------------*/
238/*----------------------------------------------------------------------------*/
239
240static inline cs_real_t
241cs_math_pow4(cs_real_t x)
242{
243 return (x*x)*(x*x);
244}
245
246/*----------------------------------------------------------------------------*/
256/*----------------------------------------------------------------------------*/
257
258static inline cs_real_t
259cs_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
283static inline cs_real_t
284cs_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
303static inline cs_real_t
304cs_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
325static inline cs_real_t
326cs_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
348static inline cs_real_t
349cs_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
371static inline cs_real_t
372cs_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
387static inline cs_real_t
388cs_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
404static inline void
405cs_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
428static inline void
429cs_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
449static inline void
450cs_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
471static inline void
472cs_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
496static inline void
497cs_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
517static inline void
518cs_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
538static inline void
539cs_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
560static inline void
561cs_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
582static inline void
583cs_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
603static inline void
604cs_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
625static inline void
626cs_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
646static inline cs_real_t
647cs_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
666static inline cs_real_t
667cs_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
690static inline void
691cs_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
716static inline cs_real_t
717cs_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
735static inline void
736cs_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
767static inline void
768cs_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
802static inline void
803cs_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
837static inline void
838cs_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
871static inline void
872cs_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
900static inline void
901cs_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
934static inline void
935cs_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
963static inline void
964cs_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
1012static inline void
1013cs_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
1067static inline void
1068cs_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
1097void
1098cs_math_set_machine_epsilon(void);
1099
1100/*----------------------------------------------------------------------------*/
1104/*----------------------------------------------------------------------------*/
1105
1106double
1107cs_math_get_machine_epsilon(void);
1108
1109/*----------------------------------------------------------------------------*/
1119/*----------------------------------------------------------------------------*/
1120
1121void
1122cs_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
1141void
1142cs_math_sym_33_eigen(const cs_real_t m[6],
1143 cs_real_t eig_vals[3]);
1144
1145/*----------------------------------------------------------------------------*/
1158/*----------------------------------------------------------------------------*/
1159
1160void
1161cs_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
1178double
1179cs_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
1197double
1198cs_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
1215void
1216cs_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
1233void
1234cs_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
1241END_C_DECLS
1242
1243#endif /* __CS_MATH_H__ */
bft_error.h
cs_defs.h
restrict
#define restrict
Definition cs_defs.h:127
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition cs_defs.h:467
cs_real_t
double cs_real_t
Floating-point value.
Definition cs_defs.h:302
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition cs_defs.h:315
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition cs_defs.h:317
END_C_DECLS
#define END_C_DECLS
Definition cs_defs.h:468
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_lnum_t
int cs_lnum_t
local mesh entity id
Definition cs_defs.h:298
t
@ t
Definition cs_field_pointer.h:98
k
@ k
Definition cs_field_pointer.h:70
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_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_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_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_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_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_1ov6
const cs_real_t cs_math_1ov6
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
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
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_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_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
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_infinite_r
const cs_real_t cs_math_infinite_r
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_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_4ov3
const cs_real_t cs_math_4ov3
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
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_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_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
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_2ov3
const cs_real_t cs_math_2ov3
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_1ov12
const cs_real_t cs_math_1ov12
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_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_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_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_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
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
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_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_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_1ov24
const cs_real_t cs_math_1ov24
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_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_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
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_sym_tensor_component_t
cs_math_sym_tensor_component_t
Definition cs_math.h:61
ZZ
@ ZZ
Definition cs_math.h:65
XY
@ XY
Definition cs_math.h:66
XZ
@ XZ
Definition cs_math.h:68
YZ
@ YZ
Definition cs_math.h:67
YY
@ YY
Definition cs_math.h:64
XX
@ XX
Definition cs_math.h:63
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_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_1ov3
const cs_real_t cs_math_1ov3
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_5ov3
const cs_real_t cs_math_5ov3
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
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_epzero
const cs_real_t cs_math_epzero
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_big_r
const cs_real_t cs_math_big_r
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
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_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_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_pi
const cs_real_t cs_math_pi
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_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_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_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
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_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
cs_math_zero_threshold
const cs_real_t cs_math_zero_threshold
len
size_t len
Definition mei_scanner.c:560
cs_nvec3_t
Definition cs_defs.h:343
cs_nvec3_t::meas
double meas
Definition cs_defs.h:345
cs_nvec3_t::unitv
double unitv[3]
Definition cs_defs.h:346
Generated on Wed Mar 13 2024 18:14:52 for My Project by doxygen 1.9.8