AQUAgpusph 5.0.4
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
22
23#define vec float4
24#define dvec double4
25#define ivec int4
26#define lvec long4
27#define uivec uint4
28#define ulvec ulong4
29#define svec usize4
30#define ssvec ssize4
31#define matrix float16
32
35#define VEC_ZERO ((float4)(0.f,0.f,0.f,0.f))
38#define VEC_ONE ((float4)(1.f, 1.f, 1.f, 0.f))
41#define VEC_ALL_ONE ((float4)(1.f, 1.f, 1.f, 1.f))
45#define VEC_INFINITY ((float4)(INFINITY, INFINITY, INFINITY, 0.f))
48#define VEC_ALL_INFINITY ((float4)(INFINITY, INFINITY, INFINITY, INFINITY))
52#define VEC_NEG_INFINITY (-VEC_INFINITY)
55#define VEC_ALL_NEG_INFINITY (-VEC_ALL_INFINITY)
56
59#define MAT_ZERO ((float16)(0.f, 0.f, 0.f, 0.f, \
60 0.f, 0.f, 0.f, 0.f, \
61 0.f, 0.f, 0.f, 0.f, \
62 0.f, 0.f, 0.f, 0.f))
63
66#define MAT_ONE ((float16)(1.f, 1.f, 1.f, 0.f, \
67 1.f, 1.f, 1.f, 0.f, \
68 1.f, 1.f, 1.f, 0.f, \
69 0.f, 0.f, 0.f, 0.f))
70
72#define MAT_ALL_ONE ((float16)(1.f, 1.f, 1.f, 1.f, \
73 1.f, 1.f, 1.f, 1.f, \
74 1.f, 1.f, 1.f, 1.f, \
75 1.f, 1.f, 1.f, 1.f))
76
83#define MAT_EYE ((float16)(1.f, 0.f, 0.f, 0.f, \
84 0.f, 1.f, 0.f, 0.f, \
85 0.f, 0.f, 1.f, 0.f, \
86 0.f, 0.f, 0.f, 0.f))
87
91#define MAT_ALL_EYE ((float16)(1.f, 0.f, 0.f, 0.f, \
92 0.f, 1.f, 0.f, 0.f, \
93 0.f, 0.f, 1.f, 0.f, \
94 0.f, 0.f, 0.f, 1.f))
95
106#define vec_xyz vec3
107#define dvec_xyz dvec3
108#define ivec_xyz ivec3
109#define lvec_xyz lvec3
110#define uivec_xyz uivec3
111#define ulvec_xyz ulvec3
112#define svec_xyz svec3
113#define ssvec_xyz ssvec3
115
122#define XYZ xyz
123
146#define BEGIN_LOOP_OVER_NEIGHS() \
147 C_I(); \
148 for(int ci = -1; ci <= 1; ci++) { \
149 for(int cj = -1; cj <= 1; cj++) { \
150 for(int ck = -1; ck <= 1; ck++) { \
151 const usize c_j = c_i + \
152 ci + \
153 cj * n_cells.x + \
154 ck * n_cells.x * n_cells.y; \
155 usize j = ihoc[c_j]; \
156 while((j < N) && (icell[j] == c_j)) {
157
162#define END_LOOP_OVER_NEIGHS() \
163 j++; \
164 } \
165 } \
166 } \
167 }
168
197#define BEGIN_NEIGHS(CELL, NPARTS, NCELLS, ICELL, IHOC) \
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 = CELL + \
202 __ci + \
203 __cj * NCELLS.x + \
204 __ck * NCELLS.x * NCELLS.y; \
205 uint j = IHOC[__c_j]; \
206 while((j < NPARTS) && (ICELL[j] == __c_j)) {
207
214#define END_NEIGHS() \
215 j++; \
216 } \
217 } \
218 } \
219 }
220
221
226#define MATRIX_DOT(_M, _V) \
227 ((float4)(dot(_M.s012, _V.xyz), \
228 dot(_M.s456, _V.xyz), \
229 dot(_M.s89A, _V.xyz), \
230 0.f))
231
234#define MATRIX_DOT_ALL(_M, _V) \
235 ((float4)(dot(_M.s0123, _V), \
236 dot(_M.s4567, _V), \
237 dot(_M.s89AB, _V), \
238 dot(_M.sCDEF, _V)))
239
245#define MATRIX_MUL(_M1, _M2) ((float16)( \
246 dot(_M1.s012, _M2.s048), dot(_M1.s012, _M2.s159), dot(_M1.s012, _M2.s26A), 0.f, \
247 dot(_M1.s456, _M2.s048), dot(_M1.s456, _M2.s159), dot(_M1.s456, _M2.s26A), 0.f, \
248 dot(_M1.s89A, _M2.s048), dot(_M1.s89A, _M2.s159), dot(_M1.s89A, _M2.s26A), 0.f, \
249 0.f, 0.f, 0.f, 0.f))
250
256#define MATRIX_MUL_ALL(_M1, _M2) ((float16)( \
257 dot(_M1.s0123, _M2.s048C), dot(_M1.s0123, _M2.s159D), dot(_M1.s0123, _M2.s26AE), dot(_M1.s0123, _M2.s37BF), \
258 dot(_M1.s4567, _M2.s048C), dot(_M1.s4567, _M2.s159D), dot(_M1.s4567, _M2.s26AE), dot(_M1.s4567, _M2.s37BF), \
259 dot(_M1.s89AB, _M2.s048C), dot(_M1.s89AB, _M2.s159D), dot(_M1.s89AB, _M2.s26AE), dot(_M1.s89AB, _M2.s37BF), \
260 dot(_M1.sCDEF, _M2.s048C), dot(_M1.sCDEF, _M2.s159D), dot(_M1.sCDEF, _M2.s26AE), dot(_M1.sCDEF, _M2.s37BF)))
261
264#define TRANSPOSE s048C159D26AE37BF
265
268#define DIAG s05A
269
273#define MATRIX_FROM_DIAG(_V) \
274 ((float16)(_V.x, 0.f, 0.f, 0.f, \
275 0.f, _V.y, 0.f, 0.f, \
276 0.f, 0.f, _V.z, 0.f, \
277 0.f, 0.f, 0.f, 0.f))
278
283#define MATRIX_TRACE(_M) (_M.s0 + _M.s5 + _M.sA)
284
290matrix outer(const vec3 v1, const vec3 v2)
291{
292 matrix m = MAT_ZERO;
293 m.s012 = v1.x * v2;
294 m.s456 = v1.y * v2;
295 m.s89A = v1.z * v2;
296 return m;
297}
298
304float det(const matrix m)
305{
306 return m.s0 * (m.s5 * m.sA - m.s6 * m.s9) +
307 m.s1 * (m.s6 * m.s8 - m.s4 * m.sA) +
308 m.s2 * (m.s4 * m.s9 - m.s5 * m.s8);
309}
310
311
322matrix inv(const matrix m)
323{
324 const float d = 1.f / det(m);
325 if(fabs(d) > 1.e16f) {
326 return MAT_ALL_EYE;
327 }
328 return ((matrix)(
329 (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,
330 (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,
331 (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,
332 0.f, 0.f, 0.f, 1.f));
333}
334
340#define MATRIX_INV(_M) \
341 MATRIX_MUL(inv(MATRIX_MUL(_M.TRANSPOSE, _M)), _M.TRANSPOSE)
MAT_ALL_EYE
#define MAT_ALL_EYE
MAT_EYE.
Definition 2D.h:74
det
float det(const matrix m)
Determinant of a matrix.
Definition 3D.h:304
inv
matrix inv(const matrix m)
Inverse of a matrix.
Definition 3D.h:322
outer
matrix outer(const vec3 v1, const vec3 v2)
Perform the outer product of two vectors.
Definition 3D.h:290
MAT_ZERO
#define MAT_ZERO
Definition Reduction.hcl.in:75
Generated by doxygen 1.14.0