In [69]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [70]:

    

import sys
import numpy as np
import matplotlib.pyplot as plt

import codim1.core.tools as tools
from codim1.core import *
from codim1.assembly import *
from codim1.fast_lib import *


    

In [71]:

    

# Elastic parameters
shear_modulus = 1.0
poisson_ratio = 0.25

# Quadrature points for the various circumstances
quad_min = 4
quad_max = 10
quad_logr = 10
quad_oneoverr = 10
n_elements = 50
degree = 3


    

In [72]:

    

bf = BasisFunctions.from_degree(degree)
mesh = Mesh.circular_mesh(n_elements, 1.0)
qs = QuadStrategy(mesh, quad_max, quad_max, quad_logr, quad_oneoverr)
dh = DOFHandler(mesh, bf, range(n_elements))


    

In [73]:

    

k_d = DisplacementKernel(shear_modulus, poisson_ratio)
k_t = TractionKernel(shear_modulus, poisson_ratio)
k_ta = AdjointTractionKernel(shear_modulus, poisson_ratio)
k_rh = RegularizedHypersingularKernel(shear_modulus, poisson_ratio)
k_h = HypersingularKernel(shear_modulus, poisson_ratio)


    

In [74]:

    

# Assemble the matrix of displacements induced by displacements at
# another location.
# Gup multiplies displacements
# Guu mutliplies tractions
matrix_assembler = MatrixAssembler(mesh, bf, dh, qs)
Gup = matrix_assembler.assemble_matrix(k_t)

# This mass matrix term arises from considering the cauchy singular form
# of the Gup matrix.
mass_matrix = MassMatrix(mesh, bf, bf, dh, QuadGauss(degree + 1), compute_on_init = True)
Gup -= 0.5 * mass_matrix.M


    

In [75]:

    

# Make the input function behave like a basis -- for internal reasons,
# this makes assembly easier.
# The theta width over which to apply the load to our cylinder.
alpha = (1 / 50.) * np.pi
alpha = np.cos((np.pi / 2) - alpha)

def section_traction(x, d):
    # Only apply tractions over the arcs near y=1, y=-1
    if np.abs(x[0]) < alpha:
        x_length = np.sqrt(x[0] ** 2 + x[1] ** 2)
        return -x[d] / x_length
    return 0.0
traction_function = BasisFunctions.from_function(section_traction)

# Assemble the rhs, composed of the displacements induced by the
# traction inputs.
rhs_assembler = RHSAssembler(mesh, bf, dh, qs)
rhs = rhs_assembler.assemble_rhs(traction_function, k_d)


    

In [76]:

    

soln_coeffs = np.linalg.solve(Gup, rhs)
# Create a solution object that pairs the coefficients with the basis
soln = Solution(bf, dh, soln_coeffs)


    

In [77]:

    

# Evaluate that solution at 400 points around the circle
x, u = tools.evaluate_boundary_solution(400 / n_elements, soln, mesh)


    

In [78]:

    

plt.plot(x[:, 0], u[:, 0])
plt.xlabel(r'X')
plt.ylabel(r'$u_x$', fontsize = 18)


    

    
Out[78]:





<matplotlib.text.Text at 0x749ef50>






    

In [79]:

    

plt.plot(x[:, 0], u[:, 1])
plt.xlabel(r'X')
plt.ylabel(r'$u_y$', fontsize = 18)


    

    
Out[79]:





<matplotlib.text.Text at 0x74af5d0>






    

In [80]:

    

ip = InteriorPoint(mesh, dh, qs)
x_pts = 20
y_pts = 22
x = np.linspace(-0.95, 0.95, x_pts)
y = np.linspace(-0.95, 0.95, y_pts)
int_ux = np.zeros((y_pts, x_pts))
int_uy = np.zeros((y_pts, x_pts))
for i in range(x_pts):
    for j in range(y_pts):
        x_val = x[i]
        y_val = y[j]
        if ((x_val ** 2) + (y_val ** 2)) > 0.99:
            int_ux[j, i] = 0
            int_uy[j, i] = 0
            continue
        traction_effect = ip.compute((x_val, y_val), np.zeros(2),
                                     k_d, traction_function)
        displacement_effect = ip.compute((x_val, y_val),
                                          np.zeros(2),
                                          k_t, soln)
        int_ux[j, i] = traction_effect[0] + displacement_effect[0]
        int_uy[j, i] = traction_effect[1] + displacement_effect[1]
int_u = np.array([int_ux, int_uy])
x_pts = 20
y_pts = 22
x = np.linspace(-0.95, 0.95, x_pts)
y = np.linspace(-0.95, 0.95, y_pts)
X, Y = np.meshgrid(x, y)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
quiver_plot = plt.quiver(X, Y, int_u[0, :, :], int_u[1, :, :])
plt.quiverkey(quiver_plot, 0.60, 0.95, 0.05, r'0.05', labelpos='W')
plt.xlabel('X')
plt.ylabel('Y')
plt.ylim([-1.2, 1.2])
plt.xlim([-1.2, 1.2])
circ = plt.Circle((0, 0), radius = 1.0, color = 'g', fill = False)
ax.add_patch(circ)
plt.title(r'Displacement vectors for a disk compressed in plane strain')
plt.show()


    

In [80]: