"""Constructs linear elasticity problems for first-order elements in 2D and 3D
"""
__docformat__ = "restructuredtext en"
__all__ = [ 'linear_elasticity', 'linear_elasticity_p1' ]
from scipy import array, matrix, ones, zeros, arange, empty, \
hstack, vstack, tile, ravel, mgrid, concatenate, \
cumsum, dot, eye, asarray
from scipy.linalg import inv, det
from scipy.sparse import coo_matrix, bsr_matrix
[docs]def linear_elasticity(grid, spacing=None, E=1e5, nu=0.3, format=None):
"""Linear elasticity problem discretizes with Q1 finite elements
on a regular rectangular grid
Parameters
----------
grid : tuple
length 2 tuple of grid sizes, e.g. (10,10)
spacing : tuple
length 2 tuple of grid spacings, e.g. (1.0,0.1)
E : float
Young's modulus
nu : float
Poisson's ratio
format : string
Format of the returned sparse matrix (eg. 'csr', 'bsr', etc.)
Returns
-------
A : {csr_matrix}
FE Q1 stiffness matrix
See Also
--------
linear_elasticity_p1
Notes
-----
- only 2d for now
Examples
--------
>>> from pyamg.gallery import linear_elasticity
>>> A,B = linear_elasticity((4,4))
References
----------
.. [1] J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI
"Matlab implementation of the finite element method in elasticity"
Computing, Volume 69, Issue 3 (November 2002) Pages: 239 - 263
http://www.math.hu-berlin.de/~cc/
"""
if len(grid) == 2:
return q12d(grid, spacing=spacing, E=E, nu=nu, format=format)
else:
raise NotImplemented('no support for grid=%s' % str(grid))
def q12d(grid, spacing=None, E=1e5, nu=0.3, dirichlet_boundary=True, format=None):
"""Q1 elements in 2 dimensions
See Also
--------
linear_elasticity
"""
X,Y = tuple(grid)
if X < 1 or Y < 1:
raise ValueError('invalid grid shape')
if dirichlet_boundary:
X += 1
Y += 1
pts = mgrid[0:X+1,0:Y+1]
pts = hstack((pts[0].T.reshape(-1,1) - X/2.0, pts[1].T.reshape(-1,1) - Y/2.0))
if spacing is None:
DX,DY = 1,1
else:
DX,DY = tuple(spacing)
pts = [DX,DY]
#compute local stiffness matrix
lame = E * nu / ((1 + nu) * (1 - 2*nu)) # Lame's first parameter
mu = E / (2 + 2*nu) # shear modulus
vertices = array([[0,0],[DX,0],[DX,DY],[0,DY]])
K = q12d_local(vertices, lame, mu )
nodes = arange( (X+1)*(Y+1) ).reshape(X+1,Y+1)
LL = nodes[:-1,:-1]
I = (2*LL).repeat( K.size ).reshape(-1,8,8)
J = I.copy()
I += tile( [0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8,1) )
J += tile( [0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8,1) ).T
V = tile( K, (X*Y,1) )
I = ravel(I)
J = ravel(J)
V = ravel(V)
A = coo_matrix( (V,(I,J)), shape=(pts.size,pts.size) ).tocsr() #sum duplicates
A = A.tobsr(blocksize=(2,2))
del I,J,V,LL,nodes
B = zeros(( 2 * (X+1)*(Y+1), 3))
B[0::2,0] = 1
B[1::2,1] = 1
B[0::2,2] = -pts[:,1]
B[1::2,2] = pts[:,0]
if dirichlet_boundary:
mask = zeros((X+1, Y+1),dtype='bool')
mask[1:-1,1:-1] = True
mask = ravel(mask)
data = zeros( ((X-1)*(Y-1),2,2) )
data[:,0,0] = 1
data[:,1,1] = 1
indices = arange( (X-1)*(Y-1) )
indptr = concatenate((array([0]),cumsum(mask)))
P = bsr_matrix((data,indices,indptr),shape=(2*(X+1)*(Y+1),2*(X-1)*(Y-1)))
Pt = P.T
A = P.T * A * P
B = Pt * B
return A.asformat(format),B
def q12d_local(vertices, lame, mu):
"""local stiffness matrix for two dimensional elasticity on a square element
Parameters
----------
lame : Float
Lame's first parameter
mu : Float
shear modulus
See Also
--------
linear_elasticity
Notes
-----
Vertices should be listed in counter-clockwise order::
[3]----[2]
| |
| |
[0]----[1]
Degrees of freedom are enumerated as follows::
[x=6,y=7]----[x=4,y=5]
| |
| |
[x=0,y=1]----[x=2,y=3]
"""
M = lame + 2*mu # P-wave modulus
R_11 = matrix([[ 2, -2, -1, 1],
[ -2, 2, 1, -1],
[ -1, 1, 2, -2],
[ 1, -1, -2, 2]]) / 6.0
R_12 = matrix([[ 1, 1, -1, -1],
[ -1, -1, 1, 1],
[ -1, -1, 1, 1],
[ 1, 1, -1, -1]]) / 4.0
R_22 = matrix([[ 2, 1, -1, -2],
[ 1, 2, -2, -1],
[ -1, -2, 2, 1],
[ -2, -1, 1, 2]]) / 6.0
F = inv( vstack( (vertices[1] - vertices[0], vertices[3] - vertices[0]) ) )
K = zeros((8,8)) # stiffness matrix
E = F.T * matrix([[M, 0],[0, mu]]) * F
K[0::2,0::2] = E[0,0] * R_11 + E[0,1] * R_12 + E[1,0] * R_12.T + E[1,1] * R_22
E = F.T * matrix([[mu, 0],[0, M]]) * F
K[1::2,1::2] = E[0,0] * R_11 + E[0,1] * R_12 + E[1,0] * R_12.T + E[1,1] * R_22
E = F.T * matrix([[0, mu],[lame, 0]]) * F;
K[1::2,0::2] = E[0,0] * R_11 + E[0,1] * R_12 + E[1,0] * R_12.T + E[1,1] * R_22
K[0::2,1::2] = K[1::2,0::2].T
K /= det(F)
return K
[docs]def linear_elasticity_p1(vertices, elements, E=1e5, nu=0.3, format=None):
"""P1 elements in 2 or 3 dimensions
Parameters
----------
vertices : array_like
array of vertices of a triangle or tets
elements : array_like
array of vertex indices for tri or tet elements
E : float
Young's modulus
nu : float
Poisson's ratio
format : string
'csr', 'csc', 'coo', 'bsr'
Returns
-------
A : {csr_matrix}
FE Q1 stiffness matrix
Notes
-----
- works in both 2d and in 3d
Examples
--------
>>> from pyamg.gallery import linear_elasticity_p1
>>> E = array([[0,1,2],[1,3,2]])
>>> V = array([[0.0,0.0],[1.0,0.0],[0.0,1.0],[1.0,1.0]])
>>> A,B = linear_elasticity_p1(V,E)
References
----------
.. [1] J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI
"Matlab implementation of the finite element method in elasticity"
Computing, Volume 69, Issue 3 (November 2002) Pages: 239 - 263
http://www.math.hu-berlin.de/~cc/
"""
#compute local stiffness matrix
lame = E * nu / ((1 + nu) * (1 - 2*nu)) # Lame's first parameter
mu = E / (2 + 2*nu) # shear modulus
vertices = asarray(vertices)
elements = asarray(elements)
D = vertices.shape[1] # spatial dimension
DoF = D*vertices.shape[0] # number of degrees of freedom
NE = elements.shape[0] # number of elements
if elements.shape[1] != D + 1:
raise ValueError('dimension mismatch')
if D == 2:
local_K = p12d_local
elif D == 3:
local_K = p13d_local
else:
raise NotImplementedError('only dimension 2 and 3 are supported')
row = elements.repeat(D).reshape(-1,D)
row *= D
row += arange(D)
row = row.reshape(-1,D*(D+1)).repeat(D*(D+1), axis=0)
row = row.reshape(-1,D*(D+1),D*(D+1))
col = row.swapaxes(1,2)
data = empty( (NE, D*(D+1), D*(D+1)), dtype=float)
for i in xrange(NE):
element_indices = elements[i,:]
element_vertices = vertices[element_indices,:]
data[i] = local_K(element_vertices, lame, mu)
row = row.ravel()
col = col.ravel()
data = data.ravel()
A = coo_matrix( (data,(row,col)), shape=(DoF,DoF)).tocsr() #sum duplicates
A = A.tobsr(blocksize=(D,D))
# compute rigid body modes
if D == 2:
B = zeros( (DoF,3) )
B[0::2,0] = 1 # vector field in x direction
B[1::2,1] = 1 # vector field in y direction
B[0::2,2] = -vertices[:,1] # rotation vector field (-y,x)
B[1::2,2] = vertices[:,0]
else:
B = zeros( (DoF,6) )
B[0::3,0] = 1 # vector field in x direction
B[1::3,1] = 1 # vector field in y direction
B[2::3,2] = 1 # vector field in z direction
B[0::3,3] = -vertices[:,1] # rotation vector field (-y,x,0)
B[1::3,3] = vertices[:,0]
B[0::3,4] = -vertices[:,2] # rotation vector field (-z,0,x)
B[2::3,4] = vertices[:,0]
B[1::3,5] = -vertices[:,2] # rotation vector field (0,-z,y)
B[2::3,5] = vertices[:,1]
return A.asformat(format),B
def p12d_local(vertices, lame, mu):
"""local stiffness matrix for P1 elements in 2d"""
assert(vertices.shape == (3,2))
A = vstack( (ones((1,3)),vertices.T) )
PhiGrad = inv(A)[:,1:] #gradients of basis functions
R = zeros((3,6))
R[ [[0],[2]], [0,2,4] ] = PhiGrad.T
R[ [[2],[1]], [1,3,5] ] = PhiGrad.T
C = mu*array([[2,0,0],[0,2,0],[0,0,1]]) + lame*array([[1,1,0],[1,1,0],[0,0,0]])
K = det(A)/2.0*dot(dot(R.T,C),R)
return K
def p13d_local(vertices, lame, mu):
"""local stiffness matrix for P1 elements in 3d"""
assert(vertices.shape == (4,3))
A = vstack( (ones((1,4)),vertices.T) )
PhiGrad = inv(A)[:,1:] #gradients of basis functions
R = zeros((6,12))
R[[0,3,4], 0::3] = PhiGrad.T
R[[3,1,5], 1::3] = PhiGrad.T
R[[4,5,2], 2::3] = PhiGrad.T
C = zeros((6,6))
C[0:3,0:3] = lame + 2*mu*eye(3)
C[3:6,3:6] = mu*eye(3)
K = det(A)/6*dot(dot(R.T,C),R)
return K