In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import sys
import numpy as np
import matplotlib.pyplot as plt
from codim1.core import *
from codim1.assembly import *
from codim1.fast_lib import *
import codim1.core.tools as tools


    

In [3]:

    

shear_modulus = 1.0
poisson_ratio = 0.25
offset = 1.0
x_pts = 30
y_pts = 30
n_elements = 100
degree = 7
quad_min = degree + 1
quad_max = 3 * degree
quad_logr = 3 * degree + 1
quad_oneoverr = 3 * degree + 1
interior_quad_pts = 13


    

In [4]:

    

k_d = DisplacementKernel(shear_modulus, poisson_ratio)
k_t = TractionKernel(shear_modulus, poisson_ratio)
k_tp = AdjointTractionKernel(shear_modulus, poisson_ratio)
k_h = HypersingularKernel(shear_modulus, poisson_ratio)
k_sh = SemiRegularizedHypersingularKernel(shear_modulus, poisson_ratio)
k_rh = RegularizedHypersingularKernel(shear_modulus, poisson_ratio)


    

In [6]:

    

left_end = np.array((-1.0, 0.0))
right_end = np.array((1.0, -0.0))
mesh = Mesh.simple_line_mesh(n_elements, left_end, right_end)
tools.plot_mesh(mesh)
plt.show()


    

In [7]:

    

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


    

In [8]:

    

def point_src(pt, normal):
    src_strength = (right_end - left_end) /\
                   np.linalg.norm(right_end - left_end)
    stress = np.array(k_sh.call(pt - left_end, normal, normal))
    src_strength2 = -src_strength
    stress2 = np.array(k_sh.call(pt - right_end, normal, normal))
    traction = -2 * stress.dot(src_strength)
    traction -= 2 * stress2.dot(src_strength2)
    return traction

soln_coeffs = tools.interpolate(point_src, dh, bf, mesh)


    

In [10]:

    

soln = Solution(bf, dh, soln_coeffs)
x, t = tools.evaluate_boundary_solution(8, soln, mesh)
tx = t[:, 0]
ty = t[:, 1]
distance_to_left = np.sqrt((x[:, 0] - left_end[0]) ** 2 +
                           (x[:, 1] - left_end[1]) ** 2)
plt.figure(1)
plt.plot(distance_to_left, tx, label = 'tx', linewidth = '3')
plt.plot(distance_to_left, ty, label = 'ty', linewidth = '3')
plt.show()


    

In [11]:

    

def exact_edge_dislocation_trac(X, Y, nx, ny):
    # Analytical traction field due to an edge dislocation on a surface with
    # normal n
    # Swap X and Y from the normally given solution.
    nu = poisson_ratio
    factor = 2 * shear_modulus * offset / (2 * np.pi * (1 - nu))
    denom = (Y ** 2 + X ** 2) ** 2
    sxx = -Y * (3 * X ** 2 + Y ** 2) * (factor / denom)
    syy = Y * (X ** 2 - Y ** 2) * (factor / denom)
    sxy = X * (X ** 2 - Y ** 2) * (factor / denom)
    tx = sxx * nx + sxy * ny
    ty = sxy * nx + syy * ny
    return tx, ty


    

In [12]:

    

exact_tx, exact_ty = \
    exact_edge_dislocation_trac(x[:, 0] + 1, x[:, 1], 0.0, 1.0)
exact_tx2, exact_ty2 = \
    exact_edge_dislocation_trac(x[:, 0] - 1, x[:, 1], 0.0, 1.0)
exact_tx -= exact_tx2
exact_ty -= exact_ty2


    

    




-c:8: RuntimeWarning: divide by zero encountered in divide
-c:8: RuntimeWarning: invalid value encountered in multiply
-c:9: RuntimeWarning: divide by zero encountered in divide
-c:9: RuntimeWarning: invalid value encountered in multiply
-c:10: RuntimeWarning: divide by zero encountered in divide
-c:10: RuntimeWarning: invalid value encountered in multiply

In [13]:

    

plt.figure(2)
error_x = np.abs((tx - exact_tx) / exact_tx)
plt.plot(error_x * 100)
plt.xlabel('x')
plt.ylabel(r'$|t - t_{exact}|$')
plt.title('Percent error')
plt.show()


    

In [ ]: