Anisotropic Elasticity for Inversion-Safety and Element Rehabilitation Eq. 47

An example from Anisotropic Elasticity for Inversion-Safety and Element Rehabilitation Eq. 47

The original equation:

placeholder

I❤️LA implementation:

from linearalgebra: tr

`J₃` = 1₃,₃
`κ_angle(Dₘ)` = 3(√2 v)^(2/3)(7/4‖`Dₘ`‖_F^2-1/4tr(`J₃``Dₘ`ᵀ`Dₘ`))⁻¹

where

`Dₘ` ∈ ℝ^(3×3)  
v ∈ ℝ

I❤️LA compiled to C++/Eigen:

/*
from linearalgebra: tr

`J₃` = 1₃,₃
`κ_angle(Dₘ)` = 3(√2 v)^(2/3)(7/4‖`Dₘ`‖_F^2-1/4tr(`J₃``Dₘ`ᵀ`Dₘ`))⁻¹

where

`Dₘ` ∈ ℝ^(3×3)  
v ∈ ℝ
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct anisotropic_elasticity_47ResultType {
    Eigen::Matrix<double, 3, 3> J₃;
    double κ_angle_left_parenthesis_Dₘ_right_parenthesis;
    anisotropic_elasticity_47ResultType(const Eigen::Matrix<double, 3, 3> & J₃,
               const double & κ_angle_left_parenthesis_Dₘ_right_parenthesis)
    : J₃(J₃),
    κ_angle_left_parenthesis_Dₘ_right_parenthesis(κ_angle_left_parenthesis_Dₘ_right_parenthesis)
    {}
};

anisotropic_elasticity_47ResultType anisotropic_elasticity_47(
    const Eigen::Matrix<double, 3, 3> & Dₘ,
    const double & v)
{
    Eigen::Matrix<double, 3, 3> J₃ = Eigen::MatrixXd::Ones(3, 3);

    double κ_angle_left_parenthesis_Dₘ_right_parenthesis = 3 * pow((sqrt(2) * v), (2 / double(3))) * 1 / ((7 / double(4) * pow((Dₘ).norm(), 2) - 1 / double(4) * (J₃ * Dₘ.transpose() * Dₘ).trace()));

    return anisotropic_elasticity_47ResultType(J₃, κ_angle_left_parenthesis_Dₘ_right_parenthesis);
}


void generateRandomData(Eigen::Matrix<double, 3, 3> & Dₘ,
    double & v)
{
    v = rand() % 10;
    Dₘ = Eigen::MatrixXd::Random(3, 3);
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    Eigen::Matrix<double, 3, 3> Dₘ;
    double v;
    generateRandomData(Dₘ, v);
    anisotropic_elasticity_47ResultType func_value = anisotropic_elasticity_47(Dₘ, v);
    std::cout<<"return value:\n"<<func_value.κ_angle_left_parenthesis_Dₘ_right_parenthesis<<std::endl;
    return 0;
}

I❤️LA compiled to Python/NumPy/SciPy:

"""
from linearalgebra: tr

`J₃` = 1₃,₃
`κ_angle(Dₘ)` = 3(√2 v)^(2/3)(7/4‖`Dₘ`‖_F^2-1/4tr(`J₃``Dₘ`ᵀ`Dₘ`))⁻¹

where

`Dₘ` ∈ ℝ^(3×3)  
v ∈ ℝ
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize


class anisotropic_elasticity_47ResultType:
    def __init__( self, J3, κ_angle_left_parenthesis_Dₘ_right_parenthesis):
        self.J3 = J3
        self.κ_angle_left_parenthesis_Dₘ_right_parenthesis = κ_angle_left_parenthesis_Dₘ_right_parenthesis


def anisotropic_elasticity_47(Dₘ, v):
    Dₘ = np.asarray(Dₘ, dtype=np.float64)

    assert Dₘ.shape == (3, 3)
    assert np.ndim(v) == 0

    J3 = np.ones((3, 3))
    κ_angle_left_parenthesis_Dₘ_right_parenthesis = 3 * np.power((np.sqrt(2) * v), (2 / 3)) * 1 / ((7 / 4 * np.power(np.linalg.norm(Dₘ, 'fro'), 2) - 1 / 4 * np.trace(J3 @ Dₘ.T @ Dₘ)))
    return anisotropic_elasticity_47ResultType(J3, κ_angle_left_parenthesis_Dₘ_right_parenthesis)


def generateRandomData():
    v = np.random.randn()
    Dₘ = np.random.randn(3, 3)
    return Dₘ, v


if __name__ == '__main__':
    Dₘ, v = generateRandomData()
    print("Dₘ:", Dₘ)
    print("v:", v)
    func_value = anisotropic_elasticity_47(Dₘ, v)
    print("return value: ", func_value.κ_angle_left_parenthesis_Dₘ_right_parenthesis)

I❤️LA compiled to MATLAB:

function output = anisotropic_elasticity_47(D_m, v)
% output = anisotropic_elasticity_47(`Dₘ`, v)
%
%    from linearalgebra: tr
%    
%    `J₃` = 1₃,₃
%    `κ_angle(Dₘ)` = 3(√2 v)^(2/3)(7/4‖`Dₘ`‖_F^2-1/4tr(`J₃``Dₘ`ᵀ`Dₘ`))⁻¹
%    
%    where
%    
%    `Dₘ` ∈ ℝ^(3×3)  
%    v ∈ ℝ
    if nargin==0
        warning('generating random input data');
        [D_m, v] = generateRandomData();
    end
    function [D_m, v] = generateRandomData()
        v = randn();
        D_m = randn(3, 3);
    end

    assert( isequal(size(D_m), [3, 3]) );
    assert(numel(v) == 1);

    J3 = ones(3, 3);
    kappa_angle_D_m = 3 * (sqrt(2) * v).^(2 / 3) * 1 / ((7 / 4 * norm(D_m, 'fro').^2 - 1 / 4 * trace(J3 * D_m' * D_m)));
    output.J3 = J3;

    output.kappa_angle_D_m = kappa_angle_D_m;
end

I❤️LA compiled to LaTeX:

\documentclass[12pt]{article}
\usepackage{mathdots}
\usepackage[bb=boondox]{mathalfa}
\usepackage{mathtools}
\usepackage{amssymb}
\usepackage{libertine}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\usepackage[paperheight=8in,paperwidth=4in,margin=.3in,heightrounded]{geometry}
\let\originalleft\left
\let\originalright\right
\renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft}
\renewcommand{\right}{\aftergroup\egroup\originalright}
\begin{document}

\begin{center}
\resizebox{\textwidth}{!} 
{
\begin{minipage}[c]{\textwidth}
\begin{align*}
\text{from linearalgebra import tr}\\
\textit{J₃} & = \mathbb{ 1 }_{ 3,3 } \\
\textit{κ\_angle(Dₘ)} & = 3{\left( \sqrt{2}\mathit{v} \right)}^{\frac{2}{3}}\left( \frac{7}{4}\left\|\textit{Dₘ}\right\|_F^{2} - \frac{1}{4}tr\left( \textit{J₃}{\textit{Dₘ}}^T\textit{Dₘ} \right) \right)^{-1} \\
\intertext{where} 
\textit{Dₘ} & \in \mathbb{R}^{ 3 \times 3 } \\
\mathit{v} & \in \mathbb{R} \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}

I❤️LA LaTeX output: