2D.h

Go to the documentation of this file.

/*
 *  This file is part of AQUAgpusph, a free CFD program based on SPH.
 *  Copyright (C) 2012  Jose Luis Cercos Pita <jl.cercos@upm.es>
 *
 *  AQUAgpusph is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  AQUAgpusph is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with AQUAgpusph.  If not, see <http://www.gnu.org/licenses/>.
 */
 
#ifndef INFINITY
    #define INFINITY FLT_MAX
#endif
 
#define unit unsigned int
#define vec2 float2
#define vec3 float3
#define vec4 float4
#define ivec2 int2
#define ivec3 int3
#define ivec4 int4
#define uivec2 uint2
#define uivec3 uint3
#define uivec4 uint4
 
#define vec float2
#define ivec int2
#define uivec uint2
#define matrix float4
 
#define _CONVERT(TYPE) convert_ ## TYPE
 
#define CONVERT(TYPE, v) _CONVERT(TYPE)(v)
 
#define VEC_ZERO ((float2)(0.f,0.f))
#define VEC_ONE ((float2)(1.f, 1.f))
#define VEC_ALL_ONE VEC_ONE
#define VEC_INFINITY ((float2)(INFINITY, INFINITY))
#define VEC_ALL_INFINITY VEC_INFINITY
#define VEC_NEG_INFINITY (-VEC_INFINITY)
#define VEC_ALL_NEG_INFINITY (-VEC_ALL_INFINITY)
 
#define MAT_ZERO ((float4)(0.f, 0.f,                                           \
                           0.f, 0.f))
#define MAT_ONE ((float4)(1.f, 1.f,                                            \
                          1.f, 1.f))
#define MAT_ALL_ONE MAT_ONE
#define MAT_EYE ((float4)(1.f, 0.f,                                            \
                             0.f, 1.f))
#define MAT_ALL_EYE MAT_EYE
 
#define vec_xyz vec2
 
#define ivec_xyz ivec2
 
#define uivec_xyz uivec2
 
#define XYZ xy
 
#define C_I() const uint c_i = icell[i]
 
#define BEGIN_LOOP_OVER_NEIGHS()                                               \
    C_I();                                                                     \
    for(int ci = -1; ci <= 1; ci++) {                                          \
        for(int cj = -1; cj <= 1; cj++) {                                      \
            const uint c_j = c_i +                                             \
                             ci +                                              \
                             cj * n_cells.x;                                   \
            uint j = ihoc[c_j];                                                \
            while((j < N) && (icell[j] == c_j)) {
 
#define END_LOOP_OVER_NEIGHS()                                                 \
                j++;                                                           \
            }                                                                  \
        }                                                                      \
    }
 
#define MATRIX_DOT(_M, _V)                                                     \
    ((float2)(dot(_M.s01, _V),                                                 \
              dot(_M.s23, _V)))
 
#define MATRIX_DOT_ALL MATRIX_DOT
 
#define MATRIX_MUL(_M1, _M2)                                                   \
    ((float4)(dot(_M1.s01, _M2.s02), dot(_M1.s01, _M2.s13),                    \
              dot(_M1.s23, _M2.s02), dot(_M1.s23, _M2.s13)))
 
#define MATRIX_MUL_ALL MATRIX_MUL
 
#define TRANSPOSE s0213
 
#define DIAG s03
 
#define MATRIX_FROM_DIAG(_V)                                                   \
    ((float4)(_V.x,  0.f,                                                      \
               0.f, _V.y))
 
#define MATRIX_TRACE(_M) (_M.s0 + _M.s3)
 
matrix outer(const vec v1, const vec v2)
{
    matrix m;
    m.s01 = v1.x * v2;
    m.s23 = v1.y * v2;
    return m;
}
 
float det(const matrix m)
{
    return m.s0 * m.s3 - m.s1 * m.s2;
}
 
matrix inv(const matrix m)
{
    const float d = 1.f / det(m);
    if(fabs(d) > 1.e16f) {
        return MAT_ALL_EYE;
    }
    return ((matrix)( m.s3, -m.s1,
                     -m.s2,  m.s0)) * d;
}
 
#define MATRIX_INV(_M)                                                         \
    MATRIX_MUL(inv(MATRIX_MUL(_M.TRANSPOSE, _M)), _M.TRANSPOSE)