"%pylab inline" sets matplotlib to show plots in the browser rather than opening them on the server's desktop.



In [67]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

This line sets the default matplotlib figure size to be larger.



In [68]:

    
matplotlib.rcParams['savefig.dpi'] = 120

Import numpy, matplotlib and the necessary parts of the codim1 library.



In [69]:

    
import numpy as np
import matplotlib.pyplot as plt

from codim1.assembly import *
from codim1.core import *
import codim1.core.tools as tools
from codim1.fast_lib import AdjointTractionKernel, RegularizedHypersingularKernel, DisplacementKernel, TractionKernel

Set the elastic material parameters, the range of possible quadrature point counts and the degree/number of elements for the approximation



In [70]:

    
shear_modulus = 1.0
poisson_ratio = 0.25

# Quadrature points for the various circumstances
quad_min = 2
quad_max = 10
quad_logr = 10
quad_oneoverr = 10
interior_quad_pts = 8

degree = 3
n_elements = 50

Create the appropriate kernel objects.



In [71]:

    
k_d = DisplacementKernel(shear_modulus, poisson_ratio)
k_t = TractionKernel(shear_modulus, poisson_ratio)
k_tp = AdjointTractionKernel(shear_modulus, poisson_ratio)
k_h = RegularizedHypersingularKernel(shear_modulus, poisson_ratio)

Create a simple mesh along the x axis.



In [72]:

    
mesh = Mesh.simple_line_mesh(n_elements, (-1.0, 0.0), (1.0, -0.0))
tools.plot_mesh(mesh)

Define the approximation using basis functions of the chosen degree. And create a quadrature strategy using the chosen quadrature orders. And then create the local to global degree of freedom mapping.



In [73]:

    
bf = BasisFunctions.from_degree(degree)
qs = QuadStrategy(mesh, quad_min, quad_max, quad_logr, quad_oneoverr)
dh = DOFHandler(mesh, bf)

(Further explanation needed) The integral equation I am solving is: ($\int_{S_x} \int_{S_y} \phi_i(x) G_{pp}(x, y) \phi_j(y) dS_x dS_y) \Delta u_j = \int_{S_x} \phi_i(x) f(x) dS_x$ , where $\phi_i(x)$ is a basis function, $G_{pp}$ is the hypersingular kernel from the traction boundary integral equation, and $f(x)$ is the applied traction boundary condition. Once discretized, this becomes a linear system like: $\hat{G}_{pp} \Delta u = m$ So, the right hand side will consist of the "Mass Matrix" where the source basis is the approximation basis, and the solution basis is the boundary condition. The left hand side will consist of the matrix formed by each double integral over any pair of two elements.

Here, I form the RHS.



In [74]:

    
def fnc(x, d):
    if d == 1:
        return -0.0
    return 1.0
traction_func = BasisFunctions.from_function(fnc)
mass_matrix = MassMatrix(mesh, bf, traction_func, dh,
                         QuadGauss(degree + 1), compute_on_init = True)
rhs = -mass_matrix.for_rhs()

Now, I form the $\hat{G}_{pp}$ matrix.



In [75]:

    
derivs_assembler = MatrixAssembler(mesh, bf.get_gradient_basis(mesh), dh, qs)
Gpp = derivs_assembler.assemble_matrix(k_h)

Because I only applied traction boundary conditions, any rigid body motion could be added to the solution while not violating the boundary conditions. In terms of the discretized equations, this means that the matrix $\hat{G}_{pp}$ has a big old null space. In order to make the problem well-posed, I fix the displacement to be zero at the first and last degrees of freedom.



In [76]:

    
first_x_dof = dh.dof_map[0, 0, 0]
last_x_dof = dh.dof_map[0, -1, -1]
Gpp[first_x_dof, :] = 0.0
Gpp[first_x_dof, first_x_dof] = 1.0
rhs[first_x_dof] = 0.0
Gpp[last_x_dof, :] = 0.0
Gpp[last_x_dof, last_x_dof] = 1.0
rhs[last_x_dof] = 0.0

Finally, solve the linear system.



In [77]:

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

The solution of the linear system consists of the coefficients of the polynomial approximations on each element. To get the actual value of the solution at any point, I evaluate these polynomials. Here, I evaluate the solution at 200 points along the x axis.



In [78]:

    
x, s = tools.evaluate_boundary_solution(200 / n_elements, soln, mesh)

Plot the estimated solution along with the exact solution:

$\frac{3}{2}\sqrt(1 - x^2)$



In [79]:

    
correct = 1.5 * np.sqrt(1.0 - x[:, 0] ** 2)
resid = np.abs(correct - s[:, 0])/np.abs(correct) * 100
#plt.plot(x[:, 0], s[:, 0], label = 'estimated')
#plt.plot(x[:, 0], correct)
plt.plot(x[:, 0], resid, label = 'residual')
plt.xlabel(r'X', fontsize = 18)
plt.ylabel(r'$u_x$', fontsize = 18)
plt.show()

In order to get the displacement at interior points, we need to use the displacement boundary integral equation:

$u(y) = \int_{S} G_{uu}(x, y) t(x) dS - \int_{S} G_{up}(x, y) u(x) dS$

But, because there two crack surfaces, the equation reduces to:

$u(y) = -\int_{S} G_{up}(x, y) \Delta u(x) dS$

Here, I do this for a 20x20 grid of points.



In [80]:

    
x_pts = 20
y_pts = 20
x = np.linspace(-5, 5, x_pts)
# Doesn't sample 0.0!
y = np.linspace(-5, 5, y_pts)
int_ux = np.zeros((x_pts, y_pts))
int_uy = np.zeros((x_pts, y_pts))
ip = InteriorPoint(mesh, dh, qs)
for i in range(x_pts):
    for j in range(y_pts):
        displacement_effect = -ip.compute((x[i], y[j]), np.array([0.0, 0.0]), k_t, soln)
        int_ux[j, i] = displacement_effect[0]
        int_uy[j, i] = displacement_effect[1]

Make some pretty plots!



In [81]:

    
X, Y = np.meshgrid(x, y)
plt.figure(3)
plt.title('Displacement field for a constant stress drop crack.')
plt.quiver(X, Y, int_ux, int_uy)
plt.xlabel('x')
plt.ylabel('y')
plt.show()



In [82]:

    
plt.streamplot(X, Y, int_ux, int_uy)
plt.title('Disp streamlines for const. stress drop crack.')
plt.xlabel('x')
plt.ylabel('y')
plt.show()