MEK2500 Fall 2015: Source code associated with mandatory assignment 1

oblig1.py

# Copyright (C) 2015 Marie E. Rognes
#
# Run this file with python oblig1.py

import pylab
import numpy
from numpy import linalg

# Define parameters
L = 10
m = 10
n = 10*m
k = 0.1

# Set-up set of points in reference domain
X1s = numpy.linspace(0, L, n)
X2s = numpy.linspace(0, 1, m)
X1, X2 = numpy.meshgrid(X1s, X2s)

# Define deformation
x_1 = L*k*X1*X2
x_2 = - k*X1**2

# Plot deformation as vector field
pylab.figure()
pylab.title("Deformation")
Q = pylab.quiver(X1, X2, x_1, x_2)
pylab.xlabel("X_1")
pylab.ylabel("X_2")
pylab.xlim(-1, L+1)
pylab.ylim(-1, 2)
pylab.grid(True)

# Compute strain tensor e and its eigenvalues
epsilon = numpy.zeros(shape=(n, m))
q = 0
for (i, X_1) in enumerate(X1s):
    for (j, X_2) in enumerate(X2s):
        e = numpy.matrix([[L*k*X_2 - 1, 0.5*L*k*X_1 - k*X_1],
                          [0.5*L*k*X_1 - k*X_1, -1]])
        w, v = linalg.eig(e)

        epsilon[i][j] = w[q]

# Plot largest principal strains (first eigenvalue)
pylab.figure()
pylab.title("Largest principal strain")
pylab.contourf(X1, X2, epsilon.T)
pylab.xlabel("X_1")
pylab.ylabel("X_2")
pylab.xlim(-1, L+1)
pylab.ylim(-0.5, 1.5)
pylab.colorbar()
pylab.grid(True)
pylab.show()