Page MenuHomec4science
Files F91103438

cartesian_linear-elasticity_sparse.py
No OneTemporary
Actions

File Metadata

Created
Thu, Nov 7, 22:43

cartesian_linear-elasticity_sparse.py
View Options

# ==================================================================================================
#
# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
#
# ==================================================================================================
import numpy as np
from scipy.sparse import dok_matrix
from scipy.sparse.linalg import spsolve
np.set_printoptions(linewidth=200)
# ==================================================================================================
# tensor products
# ==================================================================================================
ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
# ==================================================================================================
# mesh definition
# ==================================================================================================
# number of elements
nx = 20
ny = 20
# mesh dimensions
nelem = nx * ny # number of elements
nnode = (nx+1) * (ny+1) # number of nodes
nne = 4 # number of nodes per element
ndim = 2 # number of dimensions
ndof = nnode * ndim # total number of degrees of freedom
# out-of-plane thickness
thick = 1.
# zero-initialise coordinates, displacements, and connectivity
coor = np.zeros((nnode,ndim), dtype='float')
disp = np.zeros((nnode,ndim), dtype='float')
conn = np.zeros((nelem,nne ), dtype='int' )
# coordinates
# - grid of points
x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
# - store from grid of points
coor[:,0] = x.ravel()
coor[:,1] = y.ravel()
# connectivity
# - grid of node numbers
inode = np.arange(nnode).reshape(ny+1, nx+1)
# - store from grid of node numbers
conn[:,0] = inode[ :-1, :-1].ravel()
conn[:,1] = inode[ :-1,1: ].ravel()
conn[:,2] = inode[1: ,1: ].ravel()
conn[:,3] = inode[1: , :-1].ravel()
# - node sets
nodesLeft = inode[ :, 0]
nodesRight = inode[ :,-1]
nodesBottom = inode[ 0, :]
nodesTop = inode[-1, :]
# DOF-numbers per node
dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
# ==================================================================================================
# quadrature definition
# ==================================================================================================
# integration point coordinates (local element coordinates)
Xi = np.array([
[-1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), +1./np.sqrt(3.)],
[-1./np.sqrt(3.), +1./np.sqrt(3.)],
])
# integration point weights
W = np.array([
[1.],
[1.],
[1.],
[1.],
])
# number of integration points
nip = 4
# ==================================================================================================
# gradient of the shape functions at each integration point
# ==================================================================================================
# allocate
dNx = np.empty((nelem,nip,nne,ndim))
vol = np.empty((nelem,nip))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal coordinates
xe = coor[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# shape function gradients (w.r.t. the local element coordinates)
dNdxi = np.array([
[-.25*(1.-xi[1]), -.25*(1.-xi[0])],
[+.25*(1.-xi[1]), -.25*(1.+xi[0])],
[+.25*(1.+xi[1]), +.25*(1.+xi[0])],
[-.25*(1.+xi[1]), +.25*(1.-xi[0])],
])
# Jacobian
Je = np.einsum('mi,mj->ij', dNdxi, xe)
invJe = np.linalg.inv(Je)
detJe = np.linalg.det(Je)
# shape function gradients (w.r.t. the global coordinates)
dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
# store for later use
dNx[e,q,:,:] = dNdxe
vol[e,q] = w * detJe * thick
# ==================================================================================================
# stiffness tensor at each integration point (provides constitutive response and 'tangent')
# ==================================================================================================
# identity tensors (per integration point)
i = np.eye(3)
I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
I4s = (I4+I4rt)/2.
II = dyad22(I,I)
I4d = I4s-II/3.
# bulk and shear modulus (homogeneous)
K = 0.8333333333333333 * np.ones((nelem, nip))
G = 0.3846153846153846 * np.ones((nelem, nip))
# elasticity tensor (per integration point)
C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
# ==================================================================================================
# strain from nodal displacements, stress from constitutive response
# ==================================================================================================
# zero-initialise strain tensor per integration point
# (plain strain -> all z-components are not written and should remain zero at all times)
Eps = np.zeros((nelem,nip,3,3))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal displacements
ue = disp[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
dNdxe = dNx[e,q,:,:]
# displacement gradient
gradu = np.einsum('mi,mj->ij', dNdxe, ue)
# compute strain tensor, and store per integration point
# (use plain strain to convert 2-d to 3-d tensor)
Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
# constitutive response: strain tensor -> stress tensor (per integration point)
Sig = ddot42(C4, Eps)
# ==================================================================================================
# internal force from stress
# ==================================================================================================
# allocate internal force
fint = np.zeros((ndof))
# loop over nodes
for e, nodes in enumerate(conn):
# allocate element internal force
fe = np.zeros((nne*ndim))
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
# (plane strain: stress in z-direction irrelevant)
sig = Sig[e,q,:2,:2]
dNdxe = dNx[e,q,:,:]
dV = vol[e,q]
# assemble to element internal force
# fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
# assemble to global stiffness matrix
iie = dofs[nodes,:].ravel()
fint[iie] += fe
# ==================================================================================================
# stiffness matrix from 'tangent'
# ==================================================================================================
# allocate stiffness matrix
K = dok_matrix((ndof, ndof), dtype=np.float)
# loop over nodes
for e, nodes in enumerate(conn):
# allocate element stiffness matrix
Ke = np.zeros((nne*ndim, nne*ndim))
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
c4 = C4 [e,q,:2,:2,:2,:2]
dNdxe = dNx[e,q,:,:]
dV = vol[e,q]
# assemble to element stiffness matrix
# Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * dV
Ke += ( np.einsum('mi,ijkl,nl->mjnk', dNdxe, c4, dNdxe) * dV ).reshape(nne*ndim, nne*ndim)
# assemble to global stiffness matrix
iie = dofs[nodes,:].ravel()
K[np.ix_(iie,iie)] += Ke
# ==================================================================================================
# partition and solve
# ==================================================================================================
# prescribed external force: zero on all free DOFs
# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
fext = np.zeros((ndof))
# DOF-partitioning: ['u'nknown, 'p'rescribed]
# - prescribed
iip = np.hstack((
dofs[nodesBottom,1],
dofs[nodesLeft ,0],
dofs[nodesRight ,0],
))
# - unknown
iiu = np.setdiff1d(dofs.ravel(), iip)
# fixed displacements
up = np.hstack((
0.0 * np.ones((len(nodesBottom))),
0.0 * np.ones((len(nodesLeft ))),
0.1 * np.ones((len(nodesRight ))),
))
# residual force
r = fext - fint
# partition
# - stiffness matrix
Kuu = K.tocsr()[iiu,:].tocsc()[:,iiu]
Kup = K.tocsr()[iiu,:].tocsc()[:,iip]
Kpu = K.tocsr()[iip,:].tocsc()[:,iiu]
Kpp = K.tocsr()[iip,:].tocsc()[:,iip]
# - residual force
ru = r[iiu]
# solve for unknown displacement DOFs
uu = spsolve(Kuu, ru - Kup.dot(up))
# assemble
u = np.empty((ndof))
u[iiu] = uu
u[iip] = up
# convert to nodal displacements
disp = u[dofs]
# ==================================================================================================
# strain from nodal displacements, stress from constitutive response
# ==================================================================================================
# zero-initialise strain tensor per integration point
# (plain strain -> all z-components are not written and should remain zero at all times)
Eps = np.zeros((nelem,nip,3,3))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal displacements
ue = disp[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
dNdxe = dNx[e,q,:,:]
# displacement gradient
gradu = np.einsum('mi,mj->ij', dNdxe, ue)
# compute strain tensor, and store per integration point
# (use plain strain to convert 2-d to 3-d tensor)
Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
# constitutive response: strain tensor -> stress tensor (per integration point)
Sig = ddot42(C4, Eps)
# ==================================================================================================
# internal force from stress, reaction (external) force from internal force on fixed nodes
# ==================================================================================================
# allocate internal force
fint = np.zeros((ndof))
# loop over nodes
for e, nodes in enumerate(conn):
# allocate element internal force
fe = np.zeros((nne*ndim))
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
# (plane strain: stress in z-direction irrelevant)
sig = Sig[e,q,:2,:2]
dNdxe = dNx[e,q,:,:]
dV = vol[e,q]
# assemble to element internal force
# fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
# assemble to global stiffness matrix
iie = dofs[nodes,:].ravel()
fint[iie] += fe
# reaction force
fext[iip] = fint[iip]
# ==================================================================================================
# plot
# ==================================================================================================
import matplotlib.pyplot as plt
fig, axes = plt.subplots(ncols=2, figsize=(16,8))
# reconstruct external force as nodal vectors
fext = fext[dofs]
# plot original nodal positions + displacement field as arrows
axes[0].scatter(coor[:,0], coor[:,1])
axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
# plot new nodal positions + external (reaction) force field as arrows
axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1)
# fix axes limits
for axis in axes:
axis.set_xlim([-0.4, 1.5])
axis.set_ylim([-0.4, 1.5])
plt.show()

Event Timeline

c4science · Help