%matplotlib notebook
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from solidspy import *
from IPython.display import Image

"""
Explanation: Stiffness matrix for a general bilinear element.
Summary
This notebook describes the computational steps required in the computation of the displacement based finite element stiffness matrix using an isoparametric formulation. In order to show the detailed steps we first describe the computations at a single integration point. These steps are later written in condensed form using available SOLIDSPy subroutines.
Stiffness matrix
The stiffness matrix in the displacement based finite elment method is defined like:
$${K^{QP}} = \int\limits_V {B_{ij}^Q{C_{ijkl}}B_{kl}^PdV} $$
where:
$$C = \frac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu )}}\left[ {\begin{array}{*{20}{c}}
1&{\frac{\nu }{{1 - \nu }}}&0\
{\frac{\nu }{{1 - \nu }}}&1&0\
0&0&{\frac{{1 - 2\nu }}{{2(1 - \nu )}}}
\end{array}} \right]$$
is the elastic tensor and $B_{ij}^Q$ is the strain-displacement interpolator associated to the $Q$-th degree of freedom.
Mathematical aspects
The isoparametric formulation can be easily understood if one assumes that the physical and the natural spaces are regarded as the result of deforming the element (in the physical space) into a perfect square (in the natural space). This pseudo-deformation process is ilustrated in the figure below:
<center><img src="img/Picture1.png" alt="Distorted quad" style="width:300px"></center>
and where $t$ is a pseudo-time and ${\vec r}$ and ${\vec x}$ denote the natural and physical space respectively. Mathematically, the relation between these spaces can be written like:
$${x_i} = {x_i}(r)$$
$${r_I} = {r_I}(x).$$
Relating differential elements in both spaces we have:
$$d{x_i} = \frac{{\partial {x_i}}}{{\partial {r_J}}}d{r_J} \equiv {J_{iJ}}d{r_J}.$$
where $J_{iJ}$ is the Jacobian of the transformation. We can establish the relation between the gradient operators in both spaces like:
$$\frac{\partial }{{\partial {x_i}}} = \frac{{\partial {r_J}}}{{\partial {x_i}}}\frac{\partial }{{\partial {r_J}}}$$
or in compact notation:
$$\nabla i^x = J{iJ}^{ - 1}\nabla _J^r$$
and where
$$J_{iJ}^{ - 1} = \frac{{\partial {r_J}}}{{\partial {x_i}}}.$$
Finite element approach
In the isoparametric finite element implementation the space variable (i.e., $x_i$) is approximatted via interpolation using the same shape functions $N_i^K(r)$ as in the approximation of the primary variable $u_i$:
$${x_i} = N_i^Q(r){{\hat x}^Q}.$$
In the above ${\hat x}^Q$ corresponds to the physical coordinates of node $Q$. The contribution to the numerical-Jacobian from the Q-th degree of freedom can now be written like:
$${J_{iJ}} \equiv \frac{{\partial {x_i}}}{{\partial {r_J}}} = \frac{{\partial N_i^Q(r)}}{{\partial {r_J}}}{{\hat x}^Q}.$$
For the computation of the term $B_{ij}$ appearing in the stiffness matrix $K$, and relating nodal displacements to displacement gradients according to a particular definition of the strain tensor, we require the fundamental interpolator $L_{ij}^K$. This will allows us to obtain displacement gradients ${u_{i,j}}$ out of the nodal displacements ${\hat u}^K$.
Using the expression for the displacement vector ${u_i}$ in terms of nodal displacements and the definition of the gradient operator derived pteviously, we have for the displacement gradients:
$${u_{i,j}} = J_{jQ}^{ - 1}\nabla Q^rN_i^K(r){{\hat u}^K} \equiv L{ij}^K{{\hat u}^K}$$
from which we identify the fundamental interpolator $L_{ij}^K(r)$ according to :
$$L_{ij}^K(r) = J_{jQ}^{ - 1}N_{i,Q}^K(r).$$
The strain-displacement operator for a given strain definition can now be assembled using combinations of the fundamental interpolator. Also, notice that in the computation of $L_{ij}^K(r)$ we require the inverse of the Jacobian of the transformation. In the actual implementation we evaluate directly the Jacobian using the interpolated version of $x_i$ and perform the invesrion numerically. For our particular definition of the strain tensor the strain-displacements operator takes the form:
$$B_{ij}^Q = \frac{1}{2}\left( {L_{ij}^Q + L_{ji}^Q} \right).$$
Computational aspects
End of explanation
"""


coord = np.array([[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]])
kl = np.zeros([8, 8])
B = np.zeros((3, 2*4))
enu  = 1.0/3.0
Emod = 8.0/3.0
C = femutil.umat(enu, Emod)

"""
Explanation: Let us define ${\hat x}^Q$ as the vector of nodal points coordinates and $E$ and $\nu$ as the elastic modulus and Poisson's ratio respectively for the current material.
End of explanation
"""


XW, XP = gaussutil.gpoints2x2()

"""
Explanation: Computations at a single integration point
Retrieve the Gauss points and weights in vectors $XW$ and $XP$ respectively;
End of explanation
"""


r = XP[0, 0]
s = XP[0, 1]
alf = XW[0]

"""
Explanation: and let us consider the natural coordinates at a single point:
End of explanation
"""


dhdx = 0.25*np.array([
        [s - 1, -s + 1, s + 1, -s - 1],
        [r - 1, -r - 1, r + 1, -r + 1]])
det, jaco_inv = femutil.jacoper(dhdx, coord)

"""
Explanation: The first step involves the computation of the Jacobian operator and its inverse at the Gauss point:
$$J = \left[ {\begin{array}{{20}{c}}
{\frac{{\partial x}}{{\partial r}}}&{\frac{{\partial y}}{{\partial r}}}\
{\frac{{\partial x}}{{\partial s}}}&{\frac{{\partial s}}{{\partial s}}}
\end{array}} \right] \equiv \left[ {\begin{array}{{20}{c}}
{}&{\frac{{\partial {N^Q}}}{{\partial r}}}&{}\
{...}&{}&{...}\
{}&{\frac{{\partial {N^Q}}}{{\partial s}}}&{}
\end{array}} \right]\left{ {\begin{array}{*{20}{c}}
.&.\
.&.\
{{x^Q}}&{{y^Q}}\
.&.\
.&.
\end{array}} \right}$$
End of explanation
"""


dhdx = np.dot(jaco_inv, dhdx)

"""
Explanation: Denoting the terms of the gradient operator in the physical space like ${L_x^Q}$ and ${L_Y^Q}$ respectivly these can be obtained like:
$$\left{ {\begin{array}{{20}{c}}
{L_x^Q}\
{}\
{L_y^Q}
\end{array}} \right} = {\left[ J \right]^{ - 1}}\left{ {\begin{array}{{20}{c}}
{}&{\frac{{\partial {N^Q}}}{{\partial r}}}&{}\
{...}&{}&{...}\
{}&{\frac{{\partial {N^Q}}}{{\partial s}}}&{}
\end{array}} \right}$$
End of explanation
"""


B[0, ::2]  = dhdx[0, :]
B[1, 1::2] = dhdx[1, :]
B[2, ::2]  = dhdx[1, :]
B[2, 1::2] = dhdx[0, :]

"""
Explanation: allowing us to assemble the strain-displacement operator $B$ and to compute the strain tensor at the Gauss point as:
$$\left{ {\begin{array}{{20}{c}}
{\frac{{\partial u}}{{\partial x}}}\
{\frac{{\partial v}}{{\partial y}}}\
{\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}}
\end{array}} \right} = \left[ {\begin{array}{{20}{c}}
{}&{L_x^Q}&0&{}\
{...}&0&{L_y^Q}&{...}\
{}&{L_y^Q}&{L_x^Q}&{}
\end{array}} \right]\left{ {\begin{array}{*{20}{c}}
.\
.\
.\
{{u^Q}}\
{{v^Q}}\
.\
.\
.
\end{array}} \right}$$
End of explanation
"""


kl = np.dot(np.dot(B.T,C), B)*alf*det

np.round(kl, 3)

"""
Explanation: The stiffness matrix at the current integration point is computed according to:
$$k = {B^T}(r,s)CB(r,s)\left\| {J(r,s)} \right\|\alpha $$
End of explanation
"""


kl = np.zeros([8, 8])
B  = np.zeros([3, 8])
XW, XP = gaussutil.gpoints2x2()
ngpts = 4
for i in range(ngpts):
    ri = XP[i, 0]
    si = XP[i, 1]
    alf = XW[i]
    ddet, B = femutil.stdm4NQ(ri , si , coord)
    kl = kl + np.dot(np.dot(B.T,C), B)*alf*ddet

"""
Explanation: Complete loop over the Gauss points. UEL() subroutine.
We now compute the complete stiffness matrix after looping along all the Gauss points. In the actual implementation Python uses the strain-displacement matrix subroutine stdm4n() contained within the femutils module. In SOLIDSPy the resulting block of code is implemented in the form of a single subroutine termed UEL().
End of explanation
"""


np.round(kl , 3)

from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()

"""
Explanation: Final stiffness matrix
End of explanation
"""


