cartesian_linear-elasticity_B-matrix.py
c3116612864cneighbour

cartesian_linear-elasticity_B-matrix.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

			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

			# ==================================================================================================
			# B-matrix at each integration point
			# ==================================================================================================

			# allocate
			B = np.empty((nelem,nip,nne,3,3,3))
			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)

			# compute B-matrix and integration-point volume and store for later use
			for m in range(nne):

			Be = np.zeros((3,3,3))

			Be[0,0,0] = dNdxe[m,0]
			Be[0,1,1] = dNdxe[m,0]
			Be[1,0,0] = dNdxe[m,1]
			Be[1,1,1] = dNdxe[m,1]

			B[e,q,m,:,:,:] = Be

			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
			# ==================================================================================================

			# allocate strain tensor per integration point
			Eps = np.empty((nelem,nip,3,3))

			# loop over nodes
			for e, nodes in enumerate(conn):

			# nodal displacements in 3-d
			# u_z = 0
			ue = np.zeros((nne,3))
			ue[:,0] = disp[nodes,0]
			ue[:,1] = disp[nodes,1]

			# loop over integration points
			for q, (w, xi) in enumerate(zip(W, Xi)):

			# alias integration point values
			Be = B[e,q,:,:,:,:]

			# displacement gradient
			gradu = np.einsum('mijk,mk->ij', Be, ue)

			# compute strain tensor, and store per integration point
			Eps[e,q] = .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
			sig = Sig[e,q,:,:]
			Be = B [e,q,:,:,:,:]
			dV = vol[e,q]

			# assemble to element internal force
			# (plane strain: force in z-direction irrelevant)
			# Bm = B[e,q,m,:,:,:]
			# fm = ddot32(transpose3(Bm), sig) * dV
			# fe[[mndim+0, mndim+1]] = [fm[m,2], fm[m,0]]
			fq = np.einsum('mijk,ij->mk', Be, sig) * dV
			fe += fq[:,:2].reshape(nne*ndim)

			# assemble to global stiffness matrix
			iie = dofs[nodes,:].ravel()
			fint[iie] += fe

			# ==================================================================================================
			# stiffness matrix from 'tangent'
			# ==================================================================================================

			# allocate stiffness matrix
			K = np.zeros((ndof, ndof))

			# loop over nodes
			for e, nodes in enumerate(conn):

			# allocate element stiffness matrix
			Ke = np.zeros((nnendim, nnendim))

			# loop over integration points
			for q, (w, xi) in enumerate(zip(W, Xi)):

			# alias integration point values
			c4 = C4 [e,q,:,:,:,:]
			Be = B [e,q,:,:,:,:]
			dV = vol[e,q]

			# assemble to element stiffness matrix
			# Bm = B[e,q,m,:,:,:]
			# Bn = B[e,q,n,:,:,:]
			# Kmn = ddot33(transpose3(Bm),ddot43(C4,Bn)) * dV
			# Ke[[mndim+0,mndim+1], [nndim+0,nndim+1]] += Kmn[[2,0], [2,0]]
			Kq = np.einsum('mabc,abde,nedf->mcnf', Be, c4, Be) * dV
			Ke += Kq[:,:2,:,:2].reshape(nnendim, nnendim)

			# 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[np.ix_(iiu, iiu)]
			Kup = K[np.ix_(iiu, iip)]
			Kpu = K[np.ix_(iip, iiu)]
			Kpp = K[np.ix_(iip, iip)]
			# - residual force
			ru = r[iiu]

			# solve for unknown displacement DOFs
			uu = np.linalg.solve(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
			# ==================================================================================================

			# allocate strain tensor per integration point
			Eps = np.empty((nelem,nip,3,3))

			# loop over nodes
			for e, nodes in enumerate(conn):

			# nodal displacements in 3-d
			# u_z = 0
			ue = np.zeros((nne,3))
			ue[:,0] = disp[nodes,0]
			ue[:,1] = disp[nodes,1]

			# loop over integration points
			for q, (w, xi) in enumerate(zip(W, Xi)):

			# alias integration point values
			Be = B[e,q,:,:,:,:]

			# displacement gradient
			gradu = np.einsum('mijk,mk->ij', Be, ue)

			# compute strain tensor, and store per integration point
			Eps[e,q] = .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
			sig = Sig[e,q,:,:]
			Be = B [e,q,:,:,:,:]
			dV = vol[e,q]

			# assemble to element internal force
			# (plane strain: force in z-direction irrelevant)
			# Bm = B[e,q,m,:,:,:]
			# fm = ddot32(transpose3(Bm), sig) * dV
			# fe[[mndim+0, mndim+1]] = [fm[m,2], fm[m,0]]
			fq = np.einsum('mijk,ij->mk', Be, sig) * dV
			fe += fq[:,:2].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()

rGOOSEFEM/docs/theory/examples/cartesian_linear-elasticity_B-matrix.pyc3116612864cneighbour

cartesian_linear-elasticity_B-matrix.pyView Options

rGOOSEFEM/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py
c3116612864cneighbour

cartesian_linear-elasticity_B-matrix.py
View Options