Cartesian3d.h

Go to the documentation of this file.

 
#ifndef GMATELASTIC_CARTESIAN3D_H
#define GMATELASTIC_CARTESIAN3D_H
 
#include <GMatTensor/Cartesian3d.h>
 
#include "config.h"
#include "version.h"
 
namespace GMatElastic {
 
namespace Cartesian3d {
 
template <class T>
inline auto Epseq(const T& A) -> typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type
{
    GMATELASTIC_ASSERT(A.dimension() >= 2);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 1) == 3);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 2) == 3);
 
    return xt::eval(std::sqrt(2.0 / 3.0) * GMatTensor::Cartesian3d::Norm_deviatoric(A));
}
 
template <class T, class U>
inline void epseq(const T& A, U& ret)
{
    GMATELASTIC_ASSERT(A.dimension() >= 2);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 1) == 3);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 2) == 3);
    GMATELASTIC_ASSERT(xt::has_shape(A, ret.shape()));
 
    GMatTensor::Cartesian3d::norm_deviatoric(A, ret);
    ret *= std::sqrt(2.0 / 3.0);
}
 
template <class T>
inline auto Sigeq(const T& A) -> typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type
{
    GMATELASTIC_ASSERT(A.dimension() >= 2);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 1) == 3);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 2) == 3);
 
    return xt::eval(std::sqrt(1.5) * GMatTensor::Cartesian3d::Norm_deviatoric(A));
}
 
template <class T, class U>
inline void sigeq(const T& A, U& ret)
{
    GMATELASTIC_ASSERT(A.dimension() >= 2);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 1) == 3);
    GMATELASTIC_ASSERT(A.shape(A.dimension() - 2) == 3);
    GMATELASTIC_ASSERT(xt::has_shape(A, ret.shape()));
 
    GMatTensor::Cartesian3d::norm_deviatoric(A, ret);
    ret *= std::sqrt(1.5);
}
 
template <size_t N>
class Elastic : public GMatTensor::Cartesian3d::Array<N> {
protected:
    array_type::tensor<double, N> m_K; 
    array_type::tensor<double, N> m_G; 
    array_type::tensor<double, N + 2> m_Eps; 
    array_type::tensor<double, N + 2> m_Sig; 
    array_type::tensor<double, N + 4> m_C; 
 
    using GMatTensor::Cartesian3d::Array<N>::m_ndim;
    using GMatTensor::Cartesian3d::Array<N>::m_stride_tensor2;
    using GMatTensor::Cartesian3d::Array<N>::m_stride_tensor4;
    using GMatTensor::Cartesian3d::Array<N>::m_size;
    using GMatTensor::Cartesian3d::Array<N>::m_shape;
    using GMatTensor::Cartesian3d::Array<N>::m_shape_tensor2;
    using GMatTensor::Cartesian3d::Array<N>::m_shape_tensor4;
 
public:
    using GMatTensor::Cartesian3d::Array<N>::rank;
 
    Elastic() = default;
 
    template <class T>
    Elastic(const T& K, const T& G)
    {
        this->init_Elastic(K, G);
    }
 
protected:
    template <class T>
    void init_Elastic(const T& K, const T& G)
    {
        GMATELASTIC_ASSERT(K.dimension() == N);
        GMATELASTIC_ASSERT(xt::has_shape(K, G.shape()));
        std::copy(K.shape().cbegin(), K.shape().cend(), m_shape.begin());
        this->init(m_shape);
 
        m_K = K;
        m_G = G;
        m_Eps = xt::zeros<double>(m_shape_tensor2);
        m_Sig = xt::zeros<double>(m_shape_tensor2);
        m_C = xt::empty<double>(m_shape_tensor4);
 
#pragma omp parallel
        {
            auto C = xt::adapt(m_C.data(), {m_ndim, m_ndim, m_ndim, m_ndim});
            double K;
            double G;
            auto II = GMatTensor::Cartesian3d::II();
            auto I4d = GMatTensor::Cartesian3d::I4d();
 
#pragma omp for
            for (size_t i = 0; i < m_size; ++i) {
                C.reset_buffer(&m_C.flat(i * m_stride_tensor4), m_stride_tensor4);
                K = m_K.flat(i);
                G = m_G.flat(i);
                C = K * II + 2.0 * G * I4d;
            }
        }
    }
 
public:
    const array_type::tensor<double, N>& K() const
    {
        return m_K;
    }
 
    const array_type::tensor<double, N>& G() const
    {
        return m_G;
    }
 
    template <class T>
    void set_Eps(const T& arg)
    {
        GMATELASTIC_ASSERT(xt::has_shape(arg, m_shape_tensor2));
        std::copy(arg.cbegin(), arg.cend(), m_Eps.begin());
        this->refresh();
    }
 
    template <class T>
    void set_Eps(const T& arg, bool compute_tangent)
    {
        GMATELASTIC_ASSERT(xt::has_shape(arg, m_shape_tensor2));
        std::copy(arg.cbegin(), arg.cend(), m_Eps.begin());
        this->refresh(compute_tangent);
    }
 
    virtual void refresh(bool compute_tangent = true)
    {
        (void)(compute_tangent);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
 
            double K = m_K.flat(i);
            double G = m_G.flat(i);
 
            const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
            double* Sig = &m_Sig.flat(i * m_stride_tensor2);
 
            double epsm = GMatTensor::Cartesian3d::pointer::Hydrostatic(Eps);
 
            Sig[0] = (3.0 * K - 2.0 * G) * epsm + 2.0 * G * Eps[0];
            Sig[1] = 2.0 * G * Eps[1];
            Sig[2] = 2.0 * G * Eps[2];
            Sig[3] = 2.0 * G * Eps[3];
            Sig[4] = (3.0 * K - 2.0 * G) * epsm + 2.0 * G * Eps[4];
            Sig[5] = 2.0 * G * Eps[5];
            Sig[6] = 2.0 * G * Eps[6];
            Sig[7] = 2.0 * G * Eps[7];
            Sig[8] = (3.0 * K - 2.0 * G) * epsm + 2.0 * G * Eps[8];
        }
    }
 
    const array_type::tensor<double, N + 2>& Eps() const
    {
        return m_Eps;
    }
 
    array_type::tensor<double, N + 2>& Eps()
    {
        return m_Eps;
    }
 
    const array_type::tensor<double, N + 2>& Sig() const
    {
        return m_Sig;
    }
 
    const array_type::tensor<double, N + 4>& C() const
    {
        return m_C;
    }
 
    virtual array_type::tensor<double, N> energy() const
    {
        array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
        namespace GT = GMatTensor::Cartesian3d::pointer;
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
 
            double K = m_K.flat(i);
            double G = m_G.flat(i);
 
            const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
 
            std::array<double, m_stride_tensor2> Epsd;
            double epsm = GT::Hydrostatic_deviatoric(Eps, &Epsd[0]);
            double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&Epsd[0], &Epsd[0]));
 
            double U = 3.0 * K * std::pow(epsm, 2.0);
            double V = 2.0 * G * std::pow(epsd, 2.0);
 
            ret.flat(i) = U + V;
        }
 
        return ret;
    }
};
 
} // namespace Cartesian3d
} // namespace GMatElastic
 
#endif