#!/usr/bin/env python
# coding: utf-8

# # PGE383 - Advanced Geomechanics - Midterm Exam
# 
# **Instructions:** This exam will be assigned at 11:50 AM Friday, 10/2 via [Canvas](https://utexas.instructure.com/courses/1150226) and is due at 11:00 AM Monday, 10/5.
# 
# You are allowed any resource (e.g. books, notes, web) to complete the exam except consultation with another student or professor. You are reminded of the [UT Student Honor Code](http://www.utexas.edu/about-ut/mission-core-purpose-honor-code), violations of academic integrity with be [reported](http://deanofstudents.utexas.edu/sjs/acadint_reportingprocedures.php).

# ## Problem 1
# 
# Valid constitutive models must obey the *Principle of Material Reference Frame Indifference*, i.e. they should be invarient to rigid rotations.  Mathematically, for a Cauchy stress $\boldsymbold{\sigma}$ and a orthoganal time-dependent rotation tensor $\mathbf{R}$ we might say that the stress is invariant to a rotation if
# 
# \begin{equation}
# \boldsymbol{\sigma} = \mathbf{R} \boldsymbol{\sigma} \mathbf{R}^{\intercal}.
# \label{eqn:1}
# \end{equation}
# 
# However, often in plasticity modeling, we use rate-forms of constitutive models.  Is the rate-of-Cauchy stress material reference frame indifferent?  Explain.

# **Solution**
# 
# Take the time derivative of both sides of the equation, i.e.
# 
# \begin{align*}
#     \dot{\boldsymbol{\sigma}} = \dot{\mathbf{R}} \boldsymbol{\sigma} \mathbf{R}^\intercal + \mathbf{R} \dot{\boldsymbol{\sigma}} \mathbf{R}^\intercal + \mathbf{R} \boldsymbol{\sigma} \dot{\mathbf{R}}^\intercal
# \end{align*}
# 
# This clearly is not equivalent to $\ref{eqn:1}$, therefore the rate-of-Cauchy stress **is not** material reference frame indifferent.

# ## Problem 2
# 
# If the Lagrangian strain components at a point are as shown
# 
# $$
# \mathbf{E} = 
# \begin{bmatrix}
# -1 & 0 & 0 \\
#  0 & 1.84 & -0.28 \\
#  0 & -0.28 & 1.84
# \end{bmatrix}
# $$
# 
# Compute:
# 
# 1. the principle strains
# 
# 2. the maximum shear strain component

# **Solution**
# 
# 1. We need to simply solve an eigenvalue problem.

# In[18]:


import numpy as np
import numpy.linalg


w, v = np.linalg.eig(np.array([[-1., 0., 0.], [0., 1.84, -0.28], [0., -0.28, 1.84]]));
w


# The maximum shear strain is one-half the difference between the maximum and minumum principle strains (from Mohr's circles)

# In[9]:


Emax = 0.5 * (w[0] - w[2]); Emax


# ## Problem 3
# 
# For the stress tensor
# 
# $$
# \boldsymbol{\sigma} =
# \begin{bmatrix}
# 2 & 0 & 4 \\
# 0 & 3 & 6 \\
# 4 & 6 & 0
# \end{bmatrix}
# $$
# 
# find a direction $\mathbf{\hat{n}}$, such that the traction vector on a plane normal to $\mathbf{\hat{n}}$ has components $t_1 = t_2 = 0$, and determine $t_3$ on that plane.

# **Solution**
# 
# First we will solve for the unknown normal vector using the Cauchy stress equation, i.e. $\vec{t} = \sigma^T \hat{n}$.  Where $\vec{t}^T$ = [0 0 1], we can use 1 here in place of $t_3$ because any arbrary $t_3$ that appears in the solution, will cancel when dividing the solution by its norm to create the normal vector.

# In[21]:


sigma = np.array([[2., 0, 4],[0, 3, 6], [4, 6, 0]])
#First compute the solution
n = np.linalg.solve(sigma.T, np.array([0,0,1]))
#Now divide by the norm to get a normal vector
n_hat = n / np.linalg.norm(n); n_hat


# Now can use this normal vector to determine $t_3$.

# In[22]:


np.dot(sigma.T, n_hat)


# So the final answers in rational form are
# 
# $$\hat{n}^T = \left[ \frac{2}{3}, \frac{2}{3}, -\frac{1}{3} \right]$$
# 
# and
# 
# $$t_3 = \frac{20}{3}$$

# ## Problem 4
# 
# Evaluate the second invariant of the *deviatoric* stress for the stress tensor
# 
# $$
# \boldsymbol{\sigma} =
# \begin{bmatrix}
#   2 & -1 & 1 \\
#   -1 & 0 & 1 \\
#   1 & 1 & 2
# \end{bmatrix}
# $$

# **Solution**
# 
# First we compute the deviatoric stress

# In[25]:


sigma = np.array([[2., -1., 1.], [-1., 0., 1], [1., 1., 2.]])

S = sigma - 1./3. * np.trace(sigma) * np.eye(3); S


# Now we perform the summation
# 
# $$
# J_2 = \frac{1}{2} S_{ij} S_{ij}
# $$

# In[26]:


J_2 = 0.0
for i in range(3):
    for j in range(3):
        J_2 += 0.5 * S[i,j] * S[i, j]
        
J_2


# Or in rational form we have
# 
# $$J_2 = \frac{13}{3}$$

# ## Problem 5
# 
# Under what condition are these two statements equavalent
# 
# $$ \det({\mathbf{F}}) = 1 $$ 
# 
# $$ \frac{\partial v_i}{\partial x_i} = 0 $$
# 
# Explain your answer.

# **Solution**
# 
# The correct answer to both conditions would be the *incompressibility* condition, from conservation of mass
# 
# $$
# \det(\mathbf{F}) = 1 = \frac{\rho_o}{\rho}
# $$
# 
# which implies that density doesn't change or the motion is *volume preserving* when viewed from the reference configuration.  Likewise, in Eulerian form, conservation of mass is
# 
# $$
# \frac{D \rho}{Dt} + \rho \frac{\partial v_i}{\partial x_i} = 0
# $$
# 
# When viewed from the position of a moving particle, incompressibility implies $\frac{D \rho}{Dt} = 0$ therefore
# 
# $$
# \frac{\partial v_i}{\partial x_i} = 0
# $$

# ### Problem 6
# 
# Show that the *plane stress* equilibrium equations, in the absence of body forces, implies
# 
# $$
# 2 \frac{\partial\sigma_{12}}{\partial x_1 \partial x_2} = -\frac{\partial^2 \sigma_{11}}{\partial x_1^2} -\frac{\partial^2 \sigma_{22}}{\partial x_2^2}
# $$

# **Solution**
# 
# Starting with the momentum equations
# 
# \begin{align}
#     \rho \frac{\partial^2 u_1}{\partial t^2} &= \frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{12}}{\partial x_2} + \frac{\partial \sigma_{13}}{\partial x_3} + b_1 \\
#     \rho \frac{\partial^2 u_2}{\partial t^2} &= \frac{\partial \sigma_{12}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2} + \frac{\partial \sigma_{23}}{\partial x_3} + b_2 \\
#     \rho \frac{\partial^2 u_3}{\partial t^2} &= \frac{\partial \sigma_{31}}{\partial x_1} + \frac{\partial \sigma_{23}}{\partial x_2} + \frac{\partial \sigma_{33}}{\partial x_3} + b_3
#  \end{align}
#  
# Now we apply the equalibrium condition, i.e. $ \rho \frac{\partial^2 \vec{u}}{\partial t^2} = 0$, the plane strain assumptions, i.e. $\sigma_{33} = \sigma_{13} = \sigma_{23} = 0$ and exclude the body forces, we are left with.
# 
# \begin{align}
#     0 &= \frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{12}}{\partial x_2}  \\
#     0 &= \frac{\partial \sigma_{12}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2}
#  \end{align}
#  
# Now we differentiate the first equation w.r.t. $x_1$ and the second w.r.t $x_2$ we have 
# 
# 
# \begin{align}
#     0 &= \frac{\partial^2 \sigma_{11}}{\partial x_1^2} + \frac{\partial^2 \sigma_{12}}{\partial x_1 \partial x_2}  \\
#     0 &= \frac{\partial^2 \sigma_{12}}{\partial x_1 \partial x_2} + \frac{\partial^2 \sigma_{22}}{\partial x_2^2}
#  \end{align}
#  
#  Adding the two equations together and rearranging, we have the result
#  
# $$
# 2 \frac{\partial\sigma_{12}}{\partial x_1 \partial x_2} = -\frac{\partial^2 \sigma_{11}}{\partial x_1^2} -\frac{\partial^2 \sigma_{22}}{\partial x_2^2}
# $$

# ## Problem 7
# 
# If a ductile material undergoing uniaxial stress, exhibits linear isotropic hardening, with hardening modulus $H$ and yield stress $Y$, write a von Mises yield function for this material as a function of these material constants and equivalent plastic strain $\varepsilon^p$.
# 
# **Solution**
# 
# $$ f = \sqrt{3 J_2} - Y - H \varepsilon^p$$