floatarrayf.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This program 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
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
 
#ifndef floatarrayf_h
#define floatarrayf_h
 
#include "oofemenv.h"
#include "contextioresulttype.h"
#include "contextmode.h"
#include "datastream.h"
#include "error.h"
#include "floatarray.h"
 
#include <array>
#include <initializer_list>
#include <cmath>
#include <algorithm>
#include <numeric>
 
namespace oofem {

template<std::size_t N>
class OOFEM_EXPORT FloatArrayF
#ifdef _USE_EIGEN
    : public VectorNd<N>
#endif
{
#ifdef _USE_EIGEN
public:
    OOFEM_EIGEN_DERIVED(FloatArrayF,VectorNd<N>);
    Index size() const { return VectorNd<N>::size(); }
    const double *data() const { return VectorNd<N>::data(); }
    double *data() {  return VectorNd<N>::data(); }
 
    template<
        typename... Args,
        class = typename std::enable_if_t<sizeof...(Args) == N>,
        class = std::enable_if_t<(std::conjunction_v<std::is_same<double, Args>...>)>
    >
    FloatArrayF(Args... args){
        std::initializer_list<double> ini{ args ... };
        int i=0; for(const double& x: ini){ (*this)[i++]=x; }
    }
 
 
    // since we re-declare operator[ {...} ] below, we need to include those as well for overload resolution
    double &operator[] (Index i){ return VectorNd<N>::operator[](i); }
    const double &operator[] (Index i) const { return VectorNd<N>::operator[](i); }
#else
protected:
    std::array< double, N > values;
 
public:

    auto begin() { return this->values.begin(); }
    auto end() { return this->values.end(); }
    auto begin() const { return this->values.begin(); }
    auto end() const { return this->values.end(); }

    FloatArrayF(const FloatArrayF<N> &x) : values{x.values} { }
    FloatArrayF(const FloatArray &x)
    {
#ifndef NDEBUG
        if ( x.giveSize() != N ) {
            OOFEM_ERROR("Can't convert dynamic float array of size %d to fixed size %ld\n", x.giveSize(), N);
        }
#endif
        std::copy_n(x.begin(), N, values.begin());
    }

    template<typename... V, class = typename std::enable_if_t<sizeof...(V) == N>>
    FloatArrayF(V... x) : values{x...} { }
    FloatArrayF() : values{} { }

    void operator = (const FloatArrayF<N> &src) { values = src.values; }
 

    inline double &operator[] (std::size_t i)
    {
        #ifndef NDEBUG
        return values.at( i );
        #else
        return values [ i ];
        #endif
    }

    inline const double &operator[] (std::size_t i) const {
        #ifndef NDEBUG
        return values.at( i );
        #else
        return values [ i ];
        #endif
    }

    int size() const { return N; }
    const double *data() const { return values.data(); }
    double *data() { return values.data(); }
 
#endif

    inline double &at(std::size_t i) { return (*this)[i-1]; }
    inline double at(std::size_t i) const { return (*this)[i-1]; }
 
    #ifdef _USE_EIGEN
        template<std::size_t M> inline FloatArrayF<M> operator[] (std::size_t const (&c)[M]) const { return (*this)(c).eval(); }
        template<size_t M> inline void assign(const FloatArrayF<M> &x, int const (&c)[M] ){ (*this)(c)=x; }
        template<size_t M> inline void assemble(const FloatArrayF<M> &x, int const (&c)[M] ){ (*this)(c)+=x; }
    #else
        template<std::size_t M>
        inline FloatArrayF<M> operator[] (std::size_t const (&c)[M]) const
        {
            FloatArrayF<M> x;
            for ( std::size_t i = 0; i < M; ++i ) {
                x[i] = (*this)[ c[i] ];
            }
            return x;
        }

        template<size_t M>
        inline void assign(const FloatArrayF<M> &x, int const (&c)[M] )
        {
            for ( std::size_t i = 0; i < M; ++i ) {
                (*this)[c[i]] = x[i];
            }
        }

        template<size_t M>
        inline void assemble(const FloatArrayF<M> &x, int const (&c)[M] )
        {
            for ( std::size_t i = 0; i < M; ++i ) {
                (*this)[c[i]] += x[i];
            }
        }
    #endif
    int giveSize() const { return size(); }
    void printYourself(const std::string &name="FloatArrayF") const
    {
        printf("%s (%d): \n", name.c_str(), this->size());
        for ( double x: *this ) {
            printf( "%10.3e  ", x );
        }
        printf("\n");
    }

    const double *givePointer() const { return data(); }
    double *givePointer() { return data(); }
 
    contextIOResultType storeYourself(DataStream &stream) const
    {
        if ( !stream.write(this->givePointer(), this->size()) ) {
            return CIO_IOERR;
        }
        return CIO_OK;
    }
 
    contextIOResultType restoreYourself(DataStream &stream)
    {
        if ( !stream.read(this->givePointer(), this->size()) ) {
            return CIO_IOERR;
        }
        return CIO_OK;
    }
 
    int givePackSize(DataStream &buff) const
    {
        return buff.givePackSizeOfDouble(this->size());
    }
 
    //friend class FloatMatrixF;
};
 

template<std::size_t N>
std::ostream & operator << ( std::ostream & out, const FloatArrayF<N> & x )
{
    out << x.size();
    for ( double xi : x ) {
        out << " " << xi;
    }
    return out;
}
 
 
#ifdef _USE_EIGEN
// cast args to VectorNd<N> to avoid the operator calling itself; it also prevents Eigen from returning expression template, which would fail where the result is fed into a function arg (is this a limitation of derived class in Eigen?)
template<size_t N, size_t M> inline FloatArrayF<N> assemble(const FloatArrayF<M> &x, int const (&c)[M] ){ FloatArrayF<N> out=VectorNd<N>::Zero(); out(c)=x; return out; }
template<std::size_t N> FloatArrayF<N> operator * ( double a, const FloatArrayF<N> & x ){ return a*(VectorNd<N>)x; }
template<std::size_t N> FloatArrayF<N> operator * ( const FloatArrayF<N> & x, double a ){ return a*(VectorNd<N>)x; }
template<std::size_t N> FloatArrayF<N> mult ( const FloatArrayF<N> & x, const FloatArrayF<N> & y ){ return (x.array()*y.array()).matrix(); }
template<std::size_t N> FloatArrayF<N> operator / ( const FloatArrayF<N> & x, double a ){ return ((VectorNd<N>)x)/a; }
template<std::size_t N> FloatArrayF<N> operator ^ ( const FloatArrayF<N> & x, double a){ return x.array().pow(a); }
template<std::size_t N> FloatArrayF<N> operator + ( const FloatArrayF<N> & x, const FloatArrayF<N> & y ){ return ((VectorNd<N>)x)+((VectorNd<N>)y); }
template<std::size_t N> FloatArrayF<N> operator - ( const FloatArrayF<N> & x, const FloatArrayF<N> & y ){ return ((VectorNd<N>)x)-((VectorNd<N>)y); }
template<std::size_t N> FloatArrayF<N> operator - ( const FloatArrayF<N> & x ){ return -((VectorNd<N>)x); };
//template<std::size_t N> FloatArrayF<N> &operator += ( FloatArrayF<N> & x, const FloatArrayF<N> & y ){ (VectorNd<N>)x+=y; return x; }
//template<std::size_t N> FloatArrayF<N> &operator -= ( FloatArrayF<N> & x, const FloatArrayF<N> & y ){ (VectorNd<N>)x-=y; return x; }
//template<std::size_t N> FloatArrayF<N> &operator *= ( FloatArrayF<N> & x, double a ){ (VectorNd<N>)x*=a; return x; }
//template<std::size_t N> FloatArrayF<N> &operator /= ( FloatArrayF<N> & x, double a ){ (VectorNd<N>)x/=a; return x; }
//template<std::size_t N> FloatArrayF<N> operator ^= ( FloatArrayF<N> & x, double a){ x=x.array().pow(a); return x; }
template<std::size_t N> bool iszero(const FloatArrayF<N> &x){ return (x.array()==0).all(); }
template<std::size_t N> bool isfinite( const FloatArrayF<N> &x ){ return x.isfinite().all(); }
template<std::size_t N> double norm_squared( const FloatArrayF<N> & x ){ return x.squaredNorm(); }
template<std::size_t N> double norm( const FloatArrayF<N> & x ){ return x.norm(); }
template<std::size_t N> FloatArrayF<N> normalize( const FloatArrayF<N> & x ){ return x/x.norm(); }
template<std::size_t N> double sum( const FloatArrayF<N> & x ){ return x.array().sum(); }
template<std::size_t N> double product( const FloatArrayF<N> & x ){ return x.array().prod(); }
inline FloatArrayF<3> cross( const FloatArrayF<3> & x, const FloatArrayF<3> & y ){ return x.cross(y); }
template<std::size_t N> double dot( const FloatArrayF<N> & x, const FloatArrayF<N> & y ){ return x.dot(y); }
template<std::size_t N> FloatArrayF<N> max(const FloatArrayF<N> &a, const FloatArrayF<N> &b){ return a.array().max(b); }
template<std::size_t N> FloatArrayF<N> min(const FloatArrayF<N> &a, const FloatArrayF<N> &b){ return a.array().min(b); }
template<std::size_t N> FloatArrayF<N> zeros() { return FloatArrayF<N>::Zero(); }
#else

template<size_t N, size_t M>
inline FloatArrayF<N> assemble(const FloatArrayF<M> &x, int const (&c)[M] )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < M; ++i ) {
        out[c[i]] = x[i];
    }
    return out;
}
 
 

template<std::size_t N>
FloatArrayF<N> operator * ( double a, const FloatArrayF<N> & x )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = x[i] * a;
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> operator * ( const FloatArrayF<N> & x, double a )
{
    return a*x;
}

template<std::size_t N>
FloatArrayF<N> mult ( const FloatArrayF<N> & x, const FloatArrayF<N> & y )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = x[i] * y[i];
    }
    return out;
}
 
 
template<std::size_t N>
FloatArrayF<N> operator / ( const FloatArrayF<N> & x, double a )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = x[i] / a;
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> operator ^ ( const FloatArrayF<N> & x, double a)
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = std::pow(x[i], a);
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> operator + ( const FloatArrayF<N> & x, const FloatArrayF<N> & y )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = x[i] + y[i];
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> operator - ( const FloatArrayF<N> & x, const FloatArrayF<N> & y )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = x[i] - y[i];
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> operator - ( const FloatArrayF<N> & x )
{
    FloatArrayF<N> out;
    for ( std::size_t i = 0; i < N; ++i ) {
        out[i] = - x[i];
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> &operator += ( FloatArrayF<N> & x, const FloatArrayF<N> & y )
{
    for ( std::size_t i = 0; i < N; ++i ) {
        x[i] += y[i];
    }
    return x;
}
 
template<std::size_t N>
FloatArrayF<N> &operator -= ( FloatArrayF<N> & x, const FloatArrayF<N> & y )
{
    for ( std::size_t i = 0; i < N; ++i ) {
        x[i] -= y[i];
    }
    return x;
}
 
template<std::size_t N>
FloatArrayF<N> &operator *= ( FloatArrayF<N> & x, double a )
{
    for ( auto &v : x ) {
        v *= a;
    }
    return x;
}
 
template<std::size_t N>
FloatArrayF<N> &operator /= ( FloatArrayF<N> & x, double a )
{
    for ( auto &v : x ) {
        v /= a;
    }
    return x;
}
 
template<std::size_t N>
FloatArrayF<N> operator ^= ( FloatArrayF<N> & x, double a)
{
    for ( auto &v : x ) {
        v = std::pow(v, a);
    }
}
 

template<std::size_t N>
bool iszero(const FloatArrayF<N> &x)
{
    for ( auto &v : x ) {
        if ( v != 0. ) {
            return false;
        }
    }
    return true;
}

template<std::size_t N>
bool isfinite( const FloatArrayF<N> &x )
{
    for ( auto &val : x ) {
        if ( !std::isfinite(val) ) {
            return false;
        }
    }
    return true;
}

template<std::size_t N>
double norm_squared( const FloatArrayF<N> & x )
{
    double ans = 0.;
    for ( auto &val : x ) {
        ans += val * val;
    }
    return ans;
}

template<std::size_t N>
double norm( const FloatArrayF<N> & x )
{
    return std::sqrt(norm_squared(x));
}

template<std::size_t N>
FloatArrayF<N> normalize( const FloatArrayF<N> & x )
{
    return x / norm(x);
}

template<std::size_t N>
double sum( const FloatArrayF<N> & x )
{
    return std::accumulate(x.begin(), x.end(), 0.);
}

template<std::size_t N>
double product( const FloatArrayF<N> & x )
{
    return std::accumulate(x.begin(), x.end(), 1.0, [](double a, double b) { return a*b; });
}

inline FloatArrayF<3> cross( const FloatArrayF<3> & x, const FloatArrayF<3> & y )
{
    return {
        x[1] * y[2] - x[2] * y[1],
        x[2] * y[0] - x[0] * y[2],
        x[0] * y[1] - x[1] * y[0],
    };
}

template<std::size_t N>
double dot( const FloatArrayF<N> & x, const FloatArrayF<N> & y )
{
    double ans = 0.;
    for ( std::size_t i = 0; i < N; ++i ) {
        ans += x[i] * y[i];
    }
    return ans;
}
 
template<std::size_t N>
FloatArrayF<N> max(const FloatArrayF<N> &a, const FloatArrayF<N> &b)
{
    FloatArrayF<N> out;
    for (std::size_t i = 0; i < N; ++i) {
        out[i] = std::max(a[i], b[i]);
    }
    return out;
}
 
template<std::size_t N>
FloatArrayF<N> min(const FloatArrayF<N> &a, const FloatArrayF<N> &b)
{
    FloatArrayF<N> out;
    for (std::size_t i = 0; i < N; ++i) {
        out[i] = std::min(a[i], b[i]);
    }
    return out;
}

template<std::size_t N>
FloatArrayF<N> zeros() {
    return FloatArrayF<N>();
}
 
 
 
#endif
 

template<std::size_t N>
FloatArrayF<N> distance_squared(const FloatArrayF<N> &a, const FloatArrayF<N> &b)
{
    return norm_squared(a-b);
}

template<std::size_t N>
FloatArrayF<N> distance(const FloatArrayF<N> &a, const FloatArrayF<N> &b)
{
    return norm(a-b);
}
 
 

inline void swap_46(FloatArrayF<6> &t)
{
    std::swap( t[3], t[5] );
}
 
inline void swap_46(FloatArrayF<9> &t)
{
    // OOFEM: 11, 22, 33, 23, 13, 12, 32, 31, 21
    // UMAT:  11, 22, 33, 12, 13, 23, 32, 21, 31
    std::swap( t[3], t[5] );
    std::swap( t[7], t[8] );
}
 

const FloatArrayF<6> I6 {1., 1., 1., 0., 0., 0.};

inline FloatArrayF<6> to_voigt_strain_6(const FloatArrayF<6> &s)
{
    return {s[0], s[1], s[2], 0.5*s[3], 0.5*s[4], 0.5*s[5]};
}

inline FloatArrayF<6> to_voigt_stress_6(const FloatArrayF<6> &e)
{
    return {e[0], e[1], e[2], 2*e[3], 2*e[4], 2*e[5]};
}
 
 
} // end namespace oofem
#endif // floatarrayf_h