Cartesian2d.h

Go to the documentation of this file.

 
#ifndef GMATTENSOR_CARTESIAN2D_H
#define GMATTENSOR_CARTESIAN2D_H
 
#include <xtensor/xadapt.hpp>
#include <xtensor/xnoalias.hpp>
#include <xtensor/xrandom.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
 
#include "config.h"
#include "detail.hpp"
#include "version.h"
 
namespace GMatTensor {
 
namespace Cartesian2d {
 
namespace pointer {
 
template <class T>
inline void O2(T* ret)
{
    std::fill(ret, ret + 4, T(0));
}
 
template <class T>
inline void O4(T* ret)
{
    std::fill(ret, ret + 16, T(0));
}
 
template <class T>
inline void I2(T* ret)
{
    ret[0] = 1.0;
    ret[1] = 0.0;
    ret[2] = 0.0;
    ret[3] = 1.0;
}
 
template <class T>
inline void II(T* ret)
{
    std::fill(ret, ret + 16, T(0));
 
    for (size_t i = 0; i < 2; ++i) {
        for (size_t j = 0; j < 2; ++j) {
            for (size_t k = 0; k < 2; ++k) {
                for (size_t l = 0; l < 2; ++l) {
                    if (i == j && k == l) {
                        ret[i * 8 + j * 4 + k * 2 + l] = 1.0;
                    }
                }
            }
        }
    }
}
 
template <class T>
inline void I4(T* ret)
{
    std::fill(ret, ret + 16, T(0));
 
    for (size_t i = 0; i < 2; ++i) {
        for (size_t j = 0; j < 2; ++j) {
            for (size_t k = 0; k < 2; ++k) {
                for (size_t l = 0; l < 2; ++l) {
                    if (i == l && j == k) {
                        ret[i * 8 + j * 4 + k * 2 + l] = 1.0;
                    }
                }
            }
        }
    }
}
 
template <class T>
inline void I4rt(T* ret)
{
    std::fill(ret, ret + 16, T(0));
 
    for (size_t i = 0; i < 2; ++i) {
        for (size_t j = 0; j < 2; ++j) {
            for (size_t k = 0; k < 2; ++k) {
                for (size_t l = 0; l < 2; ++l) {
                    if (i == k && j == l) {
                        ret[i * 8 + j * 4 + k * 2 + l] = 1.0;
                    }
                }
            }
        }
    }
}
 
template <class T>
inline void I4s(T* ret)
{
    I4(ret);
 
    std::array<double, 16> i4rt;
    I4rt(&i4rt[0]);
 
    std::transform(ret, ret + 16, &i4rt[0], ret, std::plus<T>());
 
    std::transform(ret, ret + 16, ret, std::bind(std::multiplies<T>(), std::placeholders::_1, 0.5));
}
 
template <class T>
inline void I4d(T* ret)
{
    I4s(ret);
 
    std::array<double, 16> ii;
    II(&ii[0]);
 
    std::transform(
        &ii[0], &ii[0] + 16, &ii[0], std::bind(std::multiplies<T>(), std::placeholders::_1, 0.5));
 
    std::transform(ret, ret + 16, &ii[0], ret, std::minus<T>());
}
 
template <class T>
inline T Trace(const T* A)
{
    return A[0] + A[3];
}
 
template <class T>
inline T Hydrostatic(const T* A)
{
    return T(0.5) * Trace(A);
}
 
template <class T>
inline void sym(const T* A, T* ret)
{
    ret[0] = A[0];
    ret[1] = 0.5 * (A[1] + A[2]);
    ret[2] = ret[1];
    ret[3] = A[3];
}
 
template <class T>
inline T Hydrostatic_deviatoric(const T* A, T* ret)
{
    T m = Hydrostatic(A);
    ret[0] = A[0] - m;
    ret[1] = A[1];
    ret[2] = A[2];
    ret[3] = A[3] - m;
    return m;
}
 
template <class T>
inline T Deviatoric_ddot_deviatoric(const T* A)
{
    T m = Hydrostatic(A);
    return (A[0] - m) * (A[0] - m) + (A[3] - m) * (A[3] - m) + T(2) * A[1] * A[2];
}
 
template <class T>
inline T Norm_deviatoric(const T* A)
{
    return std::sqrt(Deviatoric_ddot_deviatoric(A));
}
 
template <class T>
inline T A2_ddot_B2(const T* A, const T* B)
{
    return A[0] * B[0] + A[3] * B[3] + A[1] * B[2] + A[2] * B[1];
}
 
template <class T>
inline T A2s_ddot_B2s(const T* A, const T* B)
{
    return A[0] * B[0] + A[3] * B[3] + T(2) * A[1] * B[1];
}
 
template <class T>
inline void A2_dyadic_B2(const T* A, const T* B, T* ret)
{
    for (size_t i = 0; i < 2; ++i) {
        for (size_t j = 0; j < 2; ++j) {
            for (size_t k = 0; k < 2; ++k) {
                for (size_t l = 0; l < 2; ++l) {
                    ret[i * 8 + j * 4 + k * 2 + l] = A[i * 2 + j] * B[k * 2 + l];
                }
            }
        }
    }
}
 
template <class T>
inline void A2_dot_B2(const T* A, const T* B, T* ret)
{
    ret[0] = A[1] * B[2] + A[0] * B[0];
    ret[1] = A[0] * B[1] + A[1] * B[3];
    ret[2] = A[2] * B[0] + A[3] * B[2];
    ret[3] = A[2] * B[1] + A[3] * B[3];
}
 
template <class T>
inline void A4_ddot_B2(const T* A, const T* B, T* ret)
{
    std::fill(ret, ret + 4, T(0));
 
    for (size_t i = 0; i < 2; i++) {
        for (size_t j = 0; j < 2; j++) {
            for (size_t k = 0; k < 2; k++) {
                for (size_t l = 0; l < 2; l++) {
                    ret[i * 2 + j] += A[i * 8 + j * 4 + k * 2 + l] * B[l * 2 + k];
                }
            }
        }
    }
}
 
} // namespace pointer
 
inline array_type::tensor<double, 2> Random2()
{
    array_type::tensor<double, 2> ret = xt::random::randn<double>({2, 2});
    return ret;
}
 
inline array_type::tensor<double, 4> Random4()
{
    array_type::tensor<double, 4> ret = xt::random::randn<double>({2, 2, 2, 2});
    return ret;
}
 
inline array_type::tensor<double, 2> O2()
{
    return xt::zeros<double>({2, 2});
}
 
inline array_type::tensor<double, 4> O4()
{
    return xt::zeros<double>({2, 2, 2, 2});
}
 
inline array_type::tensor<double, 2> I2()
{
    array_type::tensor<double, 2> ret = xt::empty<double>({2, 2});
    pointer::I2(ret.data());
    return ret;
}
 
inline array_type::tensor<double, 4> II()
{
    array_type::tensor<double, 4> ret = xt::empty<double>({2, 2, 2, 2});
    pointer::II(ret.data());
    return ret;
}
 
inline array_type::tensor<double, 4> I4()
{
    array_type::tensor<double, 4> ret = xt::empty<double>({2, 2, 2, 2});
    pointer::I4(ret.data());
    return ret;
}
 
inline array_type::tensor<double, 4> I4rt()
{
    array_type::tensor<double, 4> ret = xt::empty<double>({2, 2, 2, 2});
    pointer::I4rt(ret.data());
    return ret;
}
 
inline array_type::tensor<double, 4> I4s()
{
    array_type::tensor<double, 4> ret = xt::empty<double>({2, 2, 2, 2});
    pointer::I4s(ret.data());
    return ret;
}
 
inline array_type::tensor<double, 4> I4d()
{
    array_type::tensor<double, 4> ret = xt::empty<double>({2, 2, 2, 2});
    pointer::I4d(ret.data());
    return ret;
}
 
template <class T>
inline auto Trace(const T& A)
{
    return detail::impl_A2<T, 2>::ret0(A, [](const auto& a) { return pointer::Trace(a); });
}
 
template <class T, class R>
inline void trace(const T& A, R& ret)
{
    detail::impl_A2<T, 2>::ret0(A, ret, [](const auto& a) { return pointer::Trace(a); });
}
 
template <class T>
inline auto Hydrostatic(const T& A)
{
    return detail::impl_A2<T, 2>::ret0(A, [](const auto& a) { return pointer::Hydrostatic(a); });
}
 
template <class T, class R>
inline void hydrostatic(const T& A, R& ret)
{
    detail::impl_A2<T, 2>::ret0(A, ret, [](const auto& a) { return pointer::Hydrostatic(a); });
}
 
template <class T>
inline auto A2_ddot_B2(const T& A, const T& B)
{
    return detail::impl_A2<T, 2>::B2_ret0(
        A, B, [](const auto& a, const auto& b) { return pointer::A2_ddot_B2(a, b); });
}
 
template <class T, class R>
inline void A2_ddot_B2(const T& A, const T& B, R& ret)
{
    detail::impl_A2<T, 2>::B2_ret0(
        A, B, ret, [](const auto& a, const auto& b) { return pointer::A2_ddot_B2(a, b); });
}
 
template <class T>
inline auto A2s_ddot_B2s(const T& A, const T& B)
{
    return detail::impl_A2<T, 2>::B2_ret0(
        A, B, [](const auto& a, const auto& b) { return pointer::A2s_ddot_B2s(a, b); });
}
 
template <class T, class R>
inline void A2s_ddot_B2s(const T& A, const T& B, R& ret)
{
    detail::impl_A2<T, 2>::B2_ret0(
        A, B, ret, [](const auto& a, const auto& b) { return pointer::A2s_ddot_B2s(a, b); });
}
 
template <class T>
inline auto Norm_deviatoric(const T& A)
{
    return detail::impl_A2<T, 2>::ret0(
        A, [](const auto& a) { return pointer::Norm_deviatoric(a); });
}
 
template <class T, class R>
inline void norm_deviatoric(const T& A, R& ret)
{
    detail::impl_A2<T, 2>::ret0(A, ret, [](const auto& a) { return pointer::Norm_deviatoric(a); });
}
 
template <class T>
inline auto Deviatoric(const T& A)
{
    return detail::impl_A2<T, 2>::ret2(
        A, [](const auto& a, const auto& r) { return pointer::Hydrostatic_deviatoric(a, r); });
}
 
template <class T, class R>
inline void deviatoric(const T& A, R& ret)
{
    detail::impl_A2<T, 2>::ret2(
        A, ret, [](const auto& a, const auto& r) { return pointer::Hydrostatic_deviatoric(a, r); });
}
 
template <class T>
inline auto Sym(const T& A)
{
    return detail::impl_A2<T, 2>::ret2(
        A, [](const auto& a, const auto& r) { return pointer::sym(a, r); });
}
 
template <class T, class R>
inline void sym(const T& A, R& ret)
{
    detail::impl_A2<T, 2>::ret2(
        A, ret, [](const auto& a, const auto& r) { return pointer::sym(a, r); });
}
 
template <class T>
inline auto A2_dot_B2(const T& A, const T& B)
{
    return detail::impl_A2<T, 2>::B2_ret2(A, B, [](const auto& a, const auto& b, const auto& r) {
        return pointer::A2_dot_B2(a, b, r);
    });
}
 
template <class T, class R>
inline void A2_dot_B2(const T& A, const T& B, R& ret)
{
    detail::impl_A2<T, 2>::B2_ret2(A, B, ret, [](const auto& a, const auto& b, const auto& r) {
        return pointer::A2_dot_B2(a, b, r);
    });
}
 
template <class T>
inline auto A2_dyadic_B2(const T& A, const T& B)
{
    return detail::impl_A2<T, 2>::B2_ret4(A, B, [](const auto& a, const auto& b, const auto& r) {
        return pointer::A2_dyadic_B2(a, b, r);
    });
}
 
template <class T, class R>
inline void A2_dyadic_B2(const T& A, const T& B, R& ret)
{
    detail::impl_A2<T, 2>::B2_ret4(A, B, ret, [](const auto& a, const auto& b, const auto& r) {
        return pointer::A2_dyadic_B2(a, b, r);
    });
}
 
template <class T, class U>
inline auto A4_ddot_B2(const T& A, const U& B)
{
    return detail::impl_A4<T, 2>::B2_ret2(A, B, [](const auto& a, const auto& b, const auto& r) {
        return pointer::A4_ddot_B2(a, b, r);
    });
}
 
template <class T, class U, class R>
inline void A4_ddot_B2(const T& A, const U& B, R& ret)
{
    detail::impl_A4<T, 2>::B2_ret2(A, B, ret, [](const auto& a, const auto& b, const auto& r) {
        return pointer::A4_ddot_B2(a, b, r);
    });
}
 
template <class T>
inline size_t underlying_size_A2(const T& A)
{
    return detail::impl_A2<T, 2>::toSizeT0(A.shape());
}
 
template <class T>
inline size_t underlying_size_A4(const T& A)
{
    return detail::impl_A4<T, 2>::toSizeT0(A.shape());
}
 
template <class T>
inline auto underlying_shape_A2(const T& A) -> std::array<size_t, detail::impl_A2<T, 2>::rank>
{
    return detail::impl_A2<T, 2>::toShapeT0(A.shape());
}
 
template <class T>
inline auto underlying_shape_A4(const T& A) -> std::array<size_t, detail::impl_A4<T, 2>::rank>
{
    return detail::impl_A4<T, 2>::toShapeT0(A.shape());
}
 
template <size_t N>
class Array {
public:
    constexpr static std::size_t rank = N;
 
    Array() = default;
 
    virtual ~Array() = default;
 
    Array(const std::array<size_t, N>& shape)
    {
        this->init(shape);
    }
 
    const std::array<size_t, N>& shape() const
    {
        return m_shape;
    }
 
    const std::array<size_t, N + 2>& shape_tensor2() const
    {
        return m_shape_tensor2;
    }
 
    const std::array<size_t, N + 4>& shape_tensor4() const
    {
        return m_shape_tensor4;
    }
 
    array_type::tensor<double, N + 2> O2() const
    {
        return xt::zeros<double>(m_shape_tensor2);
    }
 
    array_type::tensor<double, N + 4> O4() const
    {
        return xt::zeros<double>(m_shape_tensor4);
    }
 
    array_type::tensor<double, N + 2> I2() const
    {
        array_type::tensor<double, N + 2> ret = xt::empty<double>(m_shape_tensor2);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
            Cartesian2d::pointer::I2(&ret.flat(i * m_stride_tensor2));
        }
 
        return ret;
    }
 
    array_type::tensor<double, N + 4> II() const
    {
        array_type::tensor<double, N + 4> ret = xt::empty<double>(m_shape_tensor4);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
            Cartesian2d::pointer::II(&ret.flat(i * m_stride_tensor4));
        }
 
        return ret;
    }
 
    array_type::tensor<double, N + 4> I4() const
    {
        array_type::tensor<double, N + 4> ret = xt::empty<double>(m_shape_tensor4);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
            Cartesian2d::pointer::I4(&ret.flat(i * m_stride_tensor4));
        }
 
        return ret;
    }
 
    array_type::tensor<double, N + 4> I4rt() const
    {
        array_type::tensor<double, N + 4> ret = xt::empty<double>(m_shape_tensor4);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
            Cartesian2d::pointer::I4rt(&ret.flat(i * m_stride_tensor4));
        }
 
        return ret;
    }
 
    array_type::tensor<double, N + 4> I4s() const
    {
        array_type::tensor<double, N + 4> ret = xt::empty<double>(m_shape_tensor4);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
            Cartesian2d::pointer::I4s(&ret.flat(i * m_stride_tensor4));
        }
 
        return ret;
    }
 
    array_type::tensor<double, N + 4> I4d() const
    {
        array_type::tensor<double, N + 4> ret = xt::empty<double>(m_shape_tensor4);
 
#pragma omp parallel for
        for (size_t i = 0; i < m_size; ++i) {
            Cartesian2d::pointer::I4d(&ret.flat(i * m_stride_tensor4));
        }
 
        return ret;
    }
 
protected:
    void init(const std::array<size_t, N>& shape)
    {
        m_shape = shape;
        size_t nd = m_ndim;
        std::copy(m_shape.begin(), m_shape.end(), m_shape_tensor2.begin());
        std::copy(m_shape.begin(), m_shape.end(), m_shape_tensor4.begin());
        std::fill(m_shape_tensor2.begin() + N, m_shape_tensor2.end(), nd);
        std::fill(m_shape_tensor4.begin() + N, m_shape_tensor4.end(), nd);
        m_size = std::accumulate(m_shape.begin(), m_shape.end(), 1, std::multiplies<size_t>());
    }
 
    static constexpr size_t m_ndim = 2;
 
    static constexpr size_t m_stride_tensor2 = 4;
 
    static constexpr size_t m_stride_tensor4 = 16;
 
    size_t m_size;
 
    std::array<size_t, N> m_shape;
 
    std::array<size_t, N + 2> m_shape_tensor2;
 
    std::array<size_t, N + 4> m_shape_tensor4;
};
 
} // namespace Cartesian2d
} // namespace GMatTensor
 
#endif