add calculate_derivatives.py

"""

Symbolic differentiation of the spontaneous strain and piezoelectric tensor wrt P 

"""

import sympy as sym

import numpy as np

import itertools

def Calculate_Derivatives_2D():

    ndim = 2

    Px = sym.Symbol('Px')

    Py = sym.Symbol('Py')

    P0 = sym.Symbol('P0')

    e33 = sym.Symbol('e33')

    e31 = sym.Symbol('e31')

    e15 = sym.Symbol('e15')

    e00 = sym.Symbol('e00')

    norm = sym.sqrt(Px**2 + Py**2)

    P = np.array([Px, Py])

    I = np.eye(ndim)

    p = P/norm

    # Piezoelectric tensor

    eP = np.zeros((ndim,ndim,ndim),object)

    for k,i,j in itertools.product(range(ndim), repeat=3):

        eP[k,i,j] = (norm/P0)**3*(e33*p[i]*p[j]*p[k] + e31*(I[i,j]-p[i]*p[j])*p[k]

                                      + e15*0.5*( (I[k,i]-p[k]*p[i])*p[j] + (I[k,j]-p[k]*p[j])*p[i] ))

    # derivative of eP wrt Px

    deP_dPx = np.zeros_like(eP, object)

    for k,i,j in itertools.product(range(ndim), repeat=3):

        deP_dPx[k,i,j] = sym.diff(eP[k,i,j], Px)

    # derivative of eP wrt Py

    deP_dPy = np.zeros_like(eP, object)

    for k,i,j in itertools.product(range(ndim), repeat=3):

        deP_dPy[k,i,j] = sym.diff(eP[k,i,j], Py)

    # calculate spon. strain

    eps0_P = np.zeros((ndim,ndim), object)

    for i,j in itertools.product(range(ndim), repeat=2):

        eps0_P[i,j] = 3/2*e00*(norm/P0)**2*(p[i]*p[j] - 1/3*I[i,j])

    # get its deivative wrt Px

    deps0_P_dPx = np.zeros_like(eps0_P, object)

    for i,j in itertools.product(range(ndim), repeat=2):

        deps0_P_dPx[i,j] = sym.diff(eps0_P[i,j], Px)

    # get its deivative wrt Py

    deps0_P_dPy = np.zeros_like(eps0_P, object)

    for i,j in itertools.product(range(ndim), repeat=2):

        deps0_P_dPy[i,j] = sym.diff(eps0_P[i,j], Py)

    eP_val = sym.lambdify([Px, Py, P0, e33, e31, e15], eP, 'numpy')

    deP_dPx_val = sym.lambdify([Px, Py, P0, e33, e31, e15], deP_dPx, 'numpy')

    deP_dPy_val = sym.lambdify([Px, Py, P0, e33, e31, e15], deP_dPy, 'numpy')

    eps0_val = sym.lambdify([Px, Py, P0, e00], eps0_P, 'numpy')

    deps0_P_dPx_val = sym.lambdify([Px, Py, P0, e00], deps0_P_dPx, 'numpy')

    deps0_P_dPy_val = sym.lambdify([Px, Py, P0, e00], deps0_P_dPy, 'numpy')

    del eP, deP_dPx, deP_dPy, eps0_P, deps0_P_dPy, deps0_P_dPx

    return eP_val, deP_dPx_val, deP_dPy_val, eps0_val, deps0_P_dPx_val, deps0_P_dPy_val

def Calculate_Derivatives_Land_Potential_2D():

    Px = sym.Symbol('Px')

    Py = sym.Symbol('Py')

    P0 = sym.Symbol('P0')

    a0 = sym.Symbol('a0')

    a1 = sym.Symbol('a1')

    a2 = sym.Symbol('a2')

    a3 = sym.Symbol('a3')

    a4 = sym.Symbol('a4')

    Psi = a0 + a1/P0**2 * (Px**2 + Py**2) + a2/P0**4 * (Px**4 + Py**4) \

             + a3/P0**4 * Px**2 * Py**2 + a4/P0**6 * (Px**6 + Py**6)

    dPsi_dPx = sym.diff(Psi, Px)

    dPsi_dPy = sym.diff(Psi, Py)

    Psi_val = sym.lambdify([Px, Py, P0, a0, a1, a2, a3, a4], Psi, 'numpy')

    dPsi_dPx_val = sym.lambdify([Px, Py, P0, a1, a2, a3, a4], dPsi_dPx, 'numpy')

    dPsi_dPy_val = sym.lambdify([Px, Py, P0, a1, a2, a3, a4], dPsi_dPy, 'numpy')

    del Psi, dPsi_dPx, dPsi_dPy

    return Psi_val, dPsi_dPx_val, dPsi_dPy_val