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/>.
 */

 
#define vec float2
#define dvec double2
#define ivec int2
#define lvec long2
#define uivec uint2
#define ulvec ulong2
#define svec usize2
#define ssvec ssize2
#define matrix float4

#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 dvec_xyz dvec2
#define ivec_xyz ivec2
#define lvec_xyz lvec2
#define uivec_xyz uivec2
#define ulvec_xyz ulvec2
#define svec_xyz svec2
#define ssvec_xyz ssvec2

#define XYZ xy

#define BEGIN_LOOP_OVER_NEIGHS()                                               \
    C_I();                                                                     \
    for(int ci = -1; ci <= 1; ci++) {                                          \
        for(int cj = -1; cj <= 1; cj++) {                                      \
            const usize c_j = c_i +                                            \
                              ci +                                             \
                              cj * n_cells.x;                                  \
            usize j = ihoc[c_j];                                               \
            while((j < N) && (icell[j] == c_j)) {

#define END_LOOP_OVER_NEIGHS()                                                 \
                j++;                                                           \
            }                                                                  \
        }                                                                      \
    }

#define BEGIN_NEIGHS(CELL, NPARTS, NCELLS, ICELL, IHOC)                        \
    for(int __ci = -1; __ci <= 1; __ci++) {                                    \
        for(int __cj = -1; __cj <= 1; __cj++) {                                \
            const usize __c_j = CELL +                                         \
                                __ci +                                         \
                                __cj * NCELLS.x;                               \
            usize j = IHOC[__c_j];                                             \
            while((j < NPARTS) && (ICELL[j] == __c_j)) {

#define END_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)