vec.h

Go to the documentation of this file.

/*
**  vec.h
**
**  This is part of Shiokaze, a research-oriented fluid solver for computer graphics.
**  Created by Ryoichi Ando <rand@nii.ac.jp> on Dec 29, 2017.
**
**  Permission is hereby granted, free of charge, to any person obtaining a copy of
**  this software and associated documentation files (the "Software"), to deal in
**  the Software without restriction, including without limitation the rights to use,
**  copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the
**  Software, and to permit persons to whom the Software is furnished to do so,
**  subject to the following conditions:
**
**  The above copyright notice and this permission notice shall be included in all copies
**  or substantial portions of the Software.
**
**  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
**  INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
**  PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
**  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
**  CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
**  OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
//
#ifndef SHKZ_VEC_H
#define SHKZ_VEC_H
//
#include <cmath>
#include <ostream>
#include <cassert>
#include <shiokaze/core/common.h>
//
SHKZ_BEGIN_NAMESPACE
//
template <class T, unsigned D> struct vec {
    //
    T v[D];
    vec() {
        for( unsigned dim=0; dim<D; ++dim ) this->v[dim] = T();
    }
    template<class Y>  vec( const Y *v ) {
        for( unsigned dim=0; dim<D; ++dim ) this->v[dim] = v[dim];
    }
    template<class Y> vec( const vec<Y,D> &v) {
        for( unsigned dim=0; dim<D; ++dim ) this->v[dim] = v.v[dim];
    }
    vec( T v ) {
        for( unsigned dim=0; dim<D; ++dim ) this->v[dim] = v;
    }
    vec( T x, T y ) {
        assert( D == 2 );
        v[0] = x; v[1] = y;
    }
    vec( T x, T y, T z ) {
        assert( D == 3 );
        v[0] = x; v[1] = y; v[2] = z;
    }
    const T &operator[](unsigned idx) const {
        return v[idx];
    }
    T &operator[](unsigned idx) {
        return v[idx];
    }
    bool empty() const {
        for( unsigned dim=0; dim<D; ++dim ) if( v[dim] ) return false;
        return true;
    }
    bool operator==(const vec &v) const {
        for( unsigned dim=0; dim<D; ++dim ) if( this->v[dim] != v[dim]) return false;
        return true;
    }
    bool operator!=(const vec &v) const {
        return ! (*this == v);
    }
    template<class Y> vec<T,D> operator+(const vec<Y,D> &v) const {
        struct vec result(*this);
        for( unsigned dim=0; dim<D; ++dim ) result[dim] += v[dim];
        return result;
    }
    template<class Y> vec<T,D> operator-(const vec<Y,D> &v) const {
        struct vec result(*this);
        for( unsigned dim=0; dim<D; ++dim ) result[dim] -= v[dim];
        return result;
    }
    template<class Y> void operator+=(const vec<Y,D> &v) {
        for( unsigned dim=0; dim<D; ++dim ) this->v[dim] += v[dim];
    }
    template<class Y> void operator-=(const vec<Y,D> &v) {
        for( unsigned dim=0; dim<D; ++dim ) this->v[dim] -= v[dim];
    }
    template<class Y> vec operator*(Y s) const {
        struct vec result(*this);
        for( unsigned dim=0; dim<D; ++dim ) result[dim] *= s;
        return result;
    }
    template<class Y> vec operator/(Y s) const {
        struct vec result(*this);
        for( unsigned dim=0; dim<D; ++dim ) result[dim] /= s;
        return result;
    }
    template<class Y> void operator*=(Y s) {
        for( unsigned dim=0; dim<D; ++dim ) v[dim] *= s;
    }
    template<class Y> void operator/=(Y s) {
        for( unsigned dim=0; dim<D; ++dim ) v[dim] /= s;
    }
    template<class Y> T operator*(vec<Y,D> v) const {
        T result = T();
        for( unsigned dim=0; dim<D; ++dim ) result += this->v[dim]*v[dim];
        return result;
    }
    double norm2() const {
        double result = T();
        for( unsigned dim=0; dim<D; ++dim ) result += this->v[dim]*this->v[dim];
        return result;
    }
    double norm_inf() const {
        double result = T();
        for( unsigned dim=0; dim<D; ++dim ) result = std::max(result,(double)std::abs(this->v[dim]));
        return result;
    }
    double len() const {
        return sqrt(norm2());
    }
    vec normal() const {
        vec copy = *this;
        copy.normalize();
        return copy;
    }
    bool normalize() {
        T length = len();
        if( length ) {
            for( unsigned dim=0; dim<D; ++dim ) v[dim] /= length;
        }
        return length > 0.0;
    }
    vec rotate90() const {
        assert( D == 2 );
        struct vec result;
        result[0] = -v[1];
        result[1] = v[0];
        return result;
    }
    template<class Y=double> T cross2(const vec<Y,2> &r) const {
        assert( D == 2 );
        return v[0]*r[1]-v[1]*r[0];
    }
    template<class Y=double> vec<T,3> cross(const vec<Y,3> &r) const {
        assert( D == 3 );
        return vec<T,3>(v[1]*r[2]-v[2]*r[1],v[2]*r[0]-v[0]*r[2],v[0]*r[1]-v[1]*r[0]);
    }
    template<class Y=double> vec<Y,D> cell() const {
        vec<Y,D> result;
        for( int dim=0; dim<D; ++dim ) result[dim] = v[dim]+0.5;
        return result;
    }
    template<class Y=double> vec<Y,D> nodal() const {
        vec<Y,D> result;
        for( int dim=0; dim<D; ++dim ) result[dim] = v[dim];
        return result;
    }
    template<class Y=double> vec<Y,D> face( int dim ) const {
        vec<Y,D> result;
        for( int _dim=0; _dim<D; ++_dim ) result[_dim] = v[_dim]+0.5*(dim!=_dim);
        return result;
    }
    template<class Y=double> vec<Y,D> edge( int dim ) const {
        assert( D == 3 );
        return vec<Y,3>(v[0]+0.5*(dim==0),v[1]+0.5*(dim==1),v[2]+0.5*(dim==2));
    }
};
template <class T, class Y> static inline T operator^(const vec<T,2> &l, const vec<Y,2> &r) {
    return l.cross2(r);
}
template <class T, class Y> static inline vec<T,3> operator^(const vec<T,3> &l, const vec<Y,3> &r) {
    return l.cross(r);
}
template <class Y, class T, unsigned D> static inline vec<T,D> operator*(Y s, const vec<T,D> &vec) {
    return vec*s;
}
template <class T, unsigned D> static inline vec<T,D> operator-(const vec<T,D> &v) {
    return -1.0 * v;
}
//
template <class T> using vec2 = vec<T,2>;
template <class T> using vec3 = vec<T,3>;
//
using vec2r = vec2<Real>;
using vec2f = vec2<float>;
using vec2d = vec2<double>;
using vec2i = vec2<int>;
using vec3r = vec3<Real>;
using vec3f = vec3<float>;
using vec3d = vec3<double>;
using vec3i = vec3<int>;
//
SHKZ_END_NAMESPACE
#endif