In [38]:
from numpy import *
from matplotlib import pyplot as plt
%matplotlib inline


We will solve a 1-d finite element problem for a tensioned cable using the method of Freytag. This derivation comes directly from the paper "Finite Element Analysis in situ", and is reproduced merely for learning purposes.

A cable with tension $\lambda$ runs from $x=a$ to $x=b$ and have an applied load of $q$. The governing equation is $\lambda \frac{\partial^2 u}{\partial x^2} + q = 0$, where $u$ gives the height of the cable. We have fixed boundary conditions $u(a)=u_1$ and $u(b)=u_2$.

In [392]:
a = .3; b = 3.5
lamb = 1.
u1 = 1.; u2 = 2.5


We make use of an approximate distance function $\omega$ from the boundary. In this case, we can construct it using the distances from each end of the interval, $\omega_1=x-a$ and $\omega_2=b-x$. We take a suitable $\omega = \omega_1^2 + \omega_2^2 - \sqrt{\omega_1^2 + \omega_2^2}$

In [357]:
t=linspace(a,b,100); w1 = t-a; w2 = b-t; w = w1 + w2 - sqrt(w1*w1+w2*w2)
plt.plot(t,w1,c='r',label='$\omega_1$')
plt.plot(t,w2,c='r',label='$\omega_2$')
plt.plot(t,w,c='b',label='$\omega$')
plt.legend()

Out[357]:
<matplotlib.legend.Legend at 0x1137574d0>


We assume the solution has a form $u=u_0+u^*=\sum_{i=1}^n C_i\eta_i + u^*$, where $u_0$ satisfies homogeneous Dirichlet boundary conditions, and $u^*$ is a function satisfying the nonhomogeneous boundary conditions. We use the distance functions to construct these fuctions.

$u^* = \frac{\omega_1 u_2 + \omega_2 u_1}{\omega_1 + \omega_2}$.

$\eta_i = \omega \chi_i$ where $\chi_i$ is any shape function (that doesn't know about the boundary).

In [408]:
ustar = (w1*u2 + w2*u1)/(w1+w2)
plt.plot(t,w1,c='r',label='$\omega_1$');
plt.plot(t,w2,c='r',label='$\omega_2$');
plt.plot(t,ustar,c='b',label='$u^*$');


To define the shape functions, we use a fill the interval from $x_i$ to $x_f$, which are independent of $a$ and $b$, but should be chosen so $[x_i,x_f]$ contains $[a,b]$. We use $n$ shape functions on this interval, resulting in a grid size of $h=(x_f-x_i)/n$. Using linear hat shape functions

In [410]:
x0 = a-.1; x1 = b+.1; n=20; h=(x1-x0)/n
res_per_sf = 200
res = res_per_sf*n #numerical resolution
t=linspace(x0,x1,res); w1 = t-a; w2 = b-t; w = w1 + w2 - sqrt(w1*w1+w2*w2); ustar = (w1*u2 + w2*u1)/(w1+w2) #redefine distance functions in terms of res
x=linspace(x0,x1,n)
chi = zeros((n,res))
for i in range(n):
if i!=0:
interval = where( (x[i-1]<=t)&(t<=x[i]) )[0]
chi[i,interval] = linspace(0,1,shape(interval)[0])
if i!=n-1:
interval = where( (x[i]<=t)&(t<=x[i+1]) )[0]
chi[i,interval] = linspace(1,0,shape(interval)[0])
eta = w*chi;
for i in range(n):
plt.plot(t,chi[i],c='r',label='$\chi_i$')
plt.plot(t,eta[i],c='b',label='$\eta_i$')
plt.xlim(x0,x1);
#plt.savefig('eta.png',dpi=200)

Out[410]:
(0.19999999999999998, 3.6)


Using the assumed form of the solution, we plug this into the governing equation.

$\int_a^b \left( \lambda\left(\frac{\partial^2 u_0}{\partial x^2} + \frac{\partial^2 u^*}{\partial x^2}\right) + q(x) \right)\eta_j(x) dx = 0, \hspace{10pt} j=1,...,n$

$-\sum_{i=1}^n C_i \int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial \eta_i}{\partial x}dx = - \int_a^b q(x) \eta_j(x) dx + \int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial u^*}{\partial x}dx \hspace{10pt} j=1,...,n$

That is, we have a linear system $Ax=b$ where

$A_{ij} = -\int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial \eta_i}{\partial x}dx$

$b_i = \int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial u^*}{\partial x}dx - \int_a^b q(x) \eta_j(x) dx$

$x_i = C_i$.

Differentiating, we have $\frac{\partial \eta_i}{\partial x} = \frac{\partial \omega}{\partial x}\chi_i + \omega \frac{\partial \chi_i}{\partial x}$.

Using $\omega$ as above, we have $\frac{\partial \omega}{\partial x} = \frac{\omega_1 - \omega_2}{\sqrt{\omega_1^2 +\omega_2^2}}$.

Similarly, $\frac{\partial u^*}{\partial x} = \frac{u_2-u_1}{\omega_1+\omega_2} = \frac{u_2-u_1}{b-a}$

In [387]:
dw = -(w1-w2)/sqrt(w1*w1+w2*w2)
dchi = zeros((n,res))
for i in range(n):
if i!=0:
dchi[i,where( (x[i-1]<t)&(t<x[i]) )[0]] = 1./h
if i!=n-1:
dchi[i,where( (x[i]<t)&(t<x[i+1]) )[0]] = -1./h
deta = dw*chi + w*dchi
for i in range(n):
plt.plot(t,dw,c='m',label='dw')
plt.plot(t,chi[i],c='c')
plt.plot(t,dw*chi[i],c='g',label='$\dw*chi_i$')
plt.plot(t,w*dchi[i],c='y',label='$\dw*chi_i$')
plt.plot(t,dchi[i],c='r',label='$\dchi_i$')
plt.plot(t,deta[i],c='b',label='$\eta_i$')
plt.xlim(x0,x1);
plt.ylim(-1/h*1.1,1/h*1.1);

In [402]:
#Ready to set up linear system.
#We need to discount the shape functions which lie completely outside the integration interval!
#For now, assume x0 and x1 are chosen close enough to a and b so all shape functions affect the integral.
interval = where((a<=t)&(t<=b))[0] #integration interval
dx = (x1-x0)/res
A = -lamb*sum((deta.reshape(n,1,res)*deta)[...,interval],axis=-1)*dx
dustar = (u2-u1)/(b-a)
B = sum( (lamb*dustar*deta - q*eta)[...,interval] ,axis=-1)*dx

In [403]:
C = linalg.solve(A,B)

In [412]:
u0 = dot(C,eta)
plt.plot(t,u0+ustar,c='g',label='solution');
plt.plot([a,b],[u1,u2],c='black',marker='o',linestyle='');
plt.ylim(-.1,3.5);