AQUAgpusph 4.1.2
Loading...
Searching...
No Matches
3D.h
Go to the documentation of this file.
1/*
2 * This file is part of AQUAgpusph, a free CFD program based on SPH.
3 * Copyright (C) 2012 Jose Luis Cercos Pita <jl.cercos@upm.es>
4 *
5 * AQUAgpusph is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * AQUAgpusph is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with AQUAgpusph. If not, see <http://www.gnu.org/licenses/>.
17 */
18
23#define unit unsigned int
24#define vec2 float2
25#define vec3 float3
26#define vec4 float4
27#define ivec2 int2
28#define ivec3 int3
29#define ivec4 int4
30#define uivec2 uint2
31#define uivec3 uint3
32#define uivec4 uint4
33
34#define vec float4
35#define ivec int4
36#define uivec uint4
37#define matrix float16
38
46#define _CONVERT(TYPE) convert_ ## TYPE
47
59#define CONVERT(TYPE, v) _CONVERT(TYPE)(v)
60
63#define VEC_ZERO ((float4)(0.f,0.f,0.f,0.f))
66#define VEC_ONE ((float4)(1.f, 1.f, 1.f, 0.f))
69#define VEC_ALL_ONE ((float4)(1.f, 1.f, 1.f, 1.f))
73#define VEC_INFINITY ((float4)(INFINITY, INFINITY, INFINITY, 0.f))
76#define VEC_ALL_INFINITY ((float4)(INFINITY, INFINITY, INFINITY, INFINITY))
80#define VEC_NEG_INFINITY (-VEC_INFINITY)
83#define VEC_ALL_NEG_INFINITY (-VEC_ALL_INFINITY)
84
87#define MAT_ZERO ((float16)(0.f, 0.f, 0.f, 0.f, \
88 0.f, 0.f, 0.f, 0.f, \
89 0.f, 0.f, 0.f, 0.f, \
90 0.f, 0.f, 0.f, 0.f))
94#define MAT_ONE ((float16)(1.f, 1.f, 1.f, 0.f, \
95 1.f, 1.f, 1.f, 0.f, \
96 1.f, 1.f, 1.f, 0.f, \
97 0.f, 0.f, 0.f, 0.f))
100#define MAT_ALL_ONE ((float16)(1.f, 1.f, 1.f, 1.f, \
101 1.f, 1.f, 1.f, 1.f, \
102 1.f, 1.f, 1.f, 1.f, \
103 1.f, 1.f, 1.f, 1.f))
111#define MAT_EYE ((float16)(1.f, 0.f, 0.f, 0.f, \
112 0.f, 1.f, 0.f, 0.f, \
113 0.f, 0.f, 1.f, 0.f, \
114 0.f, 0.f, 0.f, 0.f))
119#define MAT_ALL_EYE ((float16)(1.f, 0.f, 0.f, 0.f, \
120 0.f, 1.f, 0.f, 0.f, \
121 0.f, 0.f, 1.f, 0.f, \
122 0.f, 0.f, 0.f, 1.f))
123
133#define vec_xyz vec3
134
144#define ivec_xyz ivec3
145
155#define uivec_xyz uivec3
156
163#define XYZ xyz
164
172#define C_I() const uint c_i = icell[i]
173
196#define BEGIN_LOOP_OVER_NEIGHS() \
197 C_I(); \
198 for(int ci = -1; ci <= 1; ci++) { \
199 for(int cj = -1; cj <= 1; cj++) { \
200 for(int ck = -1; ck <= 1; ck++) { \
201 const uint c_j = c_i + \
202 ci + \
203 cj * n_cells.x + \
204 ck * n_cells.x * n_cells.y; \
205 uint j = ihoc[c_j]; \
206 while((j < N) && (icell[j] == c_j)) {
207
212#define END_LOOP_OVER_NEIGHS() \
213 j++; \
214 } \
215 } \
216 } \
217 }
218
223#define MATRIX_DOT(_M, _V) \
224 ((float4)(dot(_M.s012, _V.xyz), \
225 dot(_M.s456, _V.xyz), \
226 dot(_M.s89A, _V.xyz), \
227 0.f))
228
231#define MATRIX_DOT_ALL(_M, _V) \
232 ((float4)(dot(_M.s0123, _V), \
233 dot(_M.s4567, _V), \
234 dot(_M.s89AB, _V), \
235 dot(_M.sCDEF, _V)))
236
242#define MATRIX_MUL(_M1, _M2) ((float16)( \
243 dot(_M1.s012, _M2.s048), dot(_M1.s012, _M2.s159), dot(_M1.s012, _M2.s26A), 0.f, \
244 dot(_M1.s456, _M2.s048), dot(_M1.s456, _M2.s159), dot(_M1.s456, _M2.s26A), 0.f, \
245 dot(_M1.s89A, _M2.s048), dot(_M1.s89A, _M2.s159), dot(_M1.s89A, _M2.s26A), 0.f, \
246 0.f, 0.f, 0.f, 0.f))
247
253#define MATRIX_MUL_ALL(_M1, _M2) ((float16)( \
254 dot(_M1.s0123, _M2.s048C), dot(_M1.s0123, _M2.s159D), dot(_M1.s0123, _M2.s26AE), dot(_M1.s0123, _M2.s37BF), \
255 dot(_M1.s4567, _M2.s048C), dot(_M1.s4567, _M2.s159D), dot(_M1.s4567, _M2.s26AE), dot(_M1.s4567, _M2.s37BF), \
256 dot(_M1.s89AB, _M2.s048C), dot(_M1.s89AB, _M2.s159D), dot(_M1.s89AB, _M2.s26AE), dot(_M1.s89AB, _M2.s37BF), \
257 dot(_M1.sCDEF, _M2.s048C), dot(_M1.sCDEF, _M2.s159D), dot(_M1.sCDEF, _M2.s26AE), dot(_M1.sCDEF, _M2.s37BF)))
258
261#define TRANSPOSE s048C159D26AE37BF
262
265#define DIAG s05A
266
270#define MATRIX_FROM_DIAG(_V) \
271 ((float16)(_V.x, 0.f, 0.f, 0.f, \
272 0.f, _V.y, 0.f, 0.f, \
273 0.f, 0.f, _V.z, 0.f, \
274 0.f, 0.f, 0.f, 0.f))
275
280#define MATRIX_TRACE(_M) (_M.s0 + _M.s5 + _M.sA)
281
287matrix outer(const vec3 v1, const vec3 v2)
288{
289 matrix m = MAT_ZERO;
290 m.s012 = v1.x * v2;
291 m.s456 = v1.y * v2;
292 m.s89A = v1.z * v2;
293 return m;
294}
295
301float det(const matrix m)
302{
303 return m.s0 * (m.s5 * m.sA - m.s6 * m.s9) +
304 m.s1 * (m.s6 * m.s8 - m.s4 * m.sA) +
305 m.s2 * (m.s4 * m.s9 - m.s5 * m.s8);
306}
307
308
319matrix inv(const matrix m)
320{
321 const float d = 1.f / det(m);
322 if(fabs(d) > 1.e16f) {
323 return MAT_ALL_EYE;
324 }
325 return ((matrix)(
326 (m.s5 * m.sA - m.s6 * m.s9) * d, (m.s2 * m.s9 - m.s1 * m.sA) * d, (m.s1 * m.s6 - m.s2 * m.s5) * d, 0.f,
327 (m.s6 * m.s8 - m.s4 * m.sA) * d, (m.s0 * m.sA - m.s2 * m.s8) * d, (m.s2 * m.s4 - m.s0 * m.s6) * d, 0.f,
328 (m.s4 * m.s9 - m.s5 * m.s8) * d, (m.s1 * m.s8 - m.s0 * m.s9) * d, (m.s0 * m.s5 - m.s1 * m.s4) * d, 0.f,
329 0.f, 0.f, 0.f, 1.f));
330}
331
337#define MATRIX_INV(_M) \
338 MATRIX_MUL(inv(MATRIX_MUL(_M.TRANSPOSE, _M)), _M.TRANSPOSE)
det
float det(const matrix m)
Determinant of a matrix.
Definition: 3D.h:301
MAT_ALL_EYE
#define MAT_ALL_EYE
Eye matrix.
Definition: 3D.h:119
inv
matrix inv(const matrix m)
Inverse of a matrix.
Definition: 3D.h:319
outer
matrix outer(const vec3 v1, const vec3 v2)
Perform the outer product of two vectors.
Definition: 3D.h:287
matrix
#define matrix
Definition: 3D.h:37
MAT_ZERO
#define MAT_ZERO
Null matrix, i.e. filled with zero components.
Definition: 3D.h:87
vec3
#define vec3
Definition: 3D.h:25
Generated by doxygen 1.9.6