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

`∂²I₅/∂f²` = 2[A₁,₁I₃  A₁,₂I₃  A₁,₃I₃
               A₂,₁I₃  A₂,₂I₃  A₂,₃I₃
               A₃,₁I₃  A₃,₂I₃  A₃,₃I₃] 

where

A ∈ ℝ^(3×3)

/*
`∂²I₅/∂f²` = 2[A₁,₁I₃  A₁,₂I₃  A₁,₃I₃
               A₂,₁I₃  A₂,₂I₃  A₂,₃I₃
               A₃,₁I₃  A₃,₂I₃  A₃,₃I₃] 

where

A ∈ ℝ^(3×3)
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct anisotropic_elasticity_7ResultType {
    Eigen::Matrix<double, 9, 9> partial_differential_²I₅_solidus_partial_differential_f²;
    anisotropic_elasticity_7ResultType(const Eigen::Matrix<double, 9, 9> & partial_differential_²I₅_solidus_partial_differential_f²)
    : partial_differential_²I₅_solidus_partial_differential_f²(partial_differential_²I₅_solidus_partial_differential_f²)
    {}
};

anisotropic_elasticity_7ResultType anisotropic_elasticity_7(const Eigen::Matrix<double, 3, 3> & A)
{
    Eigen::Matrix<double, 9, 9> partial_differential_²I₅_solidus_partial_differential_f²_0;
    partial_differential_²I₅_solidus_partial_differential_f²_0 << A(1-1, 1-1) * Eigen::MatrixXd::Identity(3, 3), A(1-1, 2-1) * Eigen::MatrixXd::Identity(3, 3), A(1-1, 3-1) * Eigen::MatrixXd::Identity(3, 3),
    A(2-1, 1-1) * Eigen::MatrixXd::Identity(3, 3), A(2-1, 2-1) * Eigen::MatrixXd::Identity(3, 3), A(2-1, 3-1) * Eigen::MatrixXd::Identity(3, 3),
    A(3-1, 1-1) * Eigen::MatrixXd::Identity(3, 3), A(3-1, 2-1) * Eigen::MatrixXd::Identity(3, 3), A(3-1, 3-1) * Eigen::MatrixXd::Identity(3, 3);
    Eigen::Matrix<double, 9, 9> partial_differential_²I₅_solidus_partial_differential_f² = 2 * partial_differential_²I₅_solidus_partial_differential_f²_0;

    return anisotropic_elasticity_7ResultType(partial_differential_²I₅_solidus_partial_differential_f²);
}


void generateRandomData(Eigen::Matrix<double, 3, 3> & A)
{
    A = Eigen::MatrixXd::Random(3, 3);
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    Eigen::Matrix<double, 3, 3> A;
    generateRandomData(A);
    anisotropic_elasticity_7ResultType func_value = anisotropic_elasticity_7(A);
    std::cout<<"return value:\n"<<func_value.partial_differential_²I₅_solidus_partial_differential_f²<<std::endl;
    return 0;
}

"""
`∂²I₅/∂f²` = 2[A₁,₁I₃  A₁,₂I₃  A₁,₃I₃
               A₂,₁I₃  A₂,₂I₃  A₂,₃I₃
               A₃,₁I₃  A₃,₂I₃  A₃,₃I₃] 

where

A ∈ ℝ^(3×3)
"""
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_7ResultType:
    def __init__( self, partial_differential_2I5_solidus_partial_differential_f2):
        self.partial_differential_2I5_solidus_partial_differential_f2 = partial_differential_2I5_solidus_partial_differential_f2


def anisotropic_elasticity_7(A):
    A = np.asarray(A, dtype=np.float64)

    assert A.shape == (3, 3)

    partial_differential_2I5_solidus_partial_differential_f2_0 = np.block([[A[1-1, 1-1] * np.identity(3), A[1-1, 2-1] * np.identity(3), A[1-1, 3-1] * np.identity(3)], [A[2-1, 1-1] * np.identity(3), A[2-1, 2-1] * np.identity(3), A[2-1, 3-1] * np.identity(3)], [A[3-1, 1-1] * np.identity(3), A[3-1, 2-1] * np.identity(3), A[3-1, 3-1] * np.identity(3)]])
    partial_differential_2I5_solidus_partial_differential_f2 = 2 * partial_differential_2I5_solidus_partial_differential_f2_0
    return anisotropic_elasticity_7ResultType(partial_differential_2I5_solidus_partial_differential_f2)


def generateRandomData():
    A = np.random.randn(3, 3)
    return A


if __name__ == '__main__':
    A = generateRandomData()
    print("A:", A)
    func_value = anisotropic_elasticity_7(A)
    print("return value: ", func_value.partial_differential_2I5_solidus_partial_differential_f2)

function output = anisotropic_elasticity_7(A)
% output = anisotropic_elasticity_7(A)
%
%    `∂²I₅/∂f²` = 2[A₁,₁I₃  A₁,₂I₃  A₁,₃I₃
%                   A₂,₁I₃  A₂,₂I₃  A₂,₃I₃
%                   A₃,₁I₃  A₃,₂I₃  A₃,₃I₃] 
%    
%    where
%    
%    A ∈ ℝ^(3×3)
    if nargin==0
        warning('generating random input data');
        [A] = generateRandomData();
    end
    function [A] = generateRandomData()
        A = randn(3, 3);
    end

    assert( isequal(size(A), [3, 3]) );

    partial_differential_2I5_solidus_partial_differential_f2_0 = [[A(1, 1) * speye(3), A(1, 2) * speye(3), A(1, 3) * speye(3)]; [A(2, 1) * speye(3), A(2, 2) * speye(3), A(2, 3) * speye(3)]; [A(3, 1) * speye(3), A(3, 2) * speye(3), A(3, 3) * speye(3)]];
    partial_differential_2I5_solidus_partial_differential_f2 = 2 * partial_differential_2I5_solidus_partial_differential_f2_0;
    output.partial_differential_2I5_solidus_partial_differential_f2 = partial_differential_2I5_solidus_partial_differential_f2;
end

\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*}
\textit{∂²I₅/∂f²} & = 2\begin{bmatrix}
\mathit{A}_{1, 1}I_{ 3 } & \mathit{A}_{1, 2}I_{ 3 } & \mathit{A}_{1, 3}I_{ 3 }\\
\mathit{A}_{2, 1}I_{ 3 } & \mathit{A}_{2, 2}I_{ 3 } & \mathit{A}_{2, 3}I_{ 3 }\\
\mathit{A}_{3, 1}I_{ 3 } & \mathit{A}_{3, 2}I_{ 3 } & \mathit{A}_{3, 3}I_{ 3 }\\
\end{bmatrix} \\
\intertext{where} 
\mathit{A} & \in \mathbb{R}^{ 3 \times 3 } \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}