## Computer Implementation Algorithm

107 The computer implementation of a numerically integrated isoparametric element is summarized as follows.

But first, it is assumed that this operation is to be performed in a function called stiff and it takes as input arguments elcod, young, poiss, type, ndime, ndofn, ngaus. In turn it will compute the stiffness matrix KELEM of element ielem.

1. Retrieve element geometry and material properties for the current element

2. Zero the stiffness matrix

3. Call function dmat to set the constitutive matrix De of the element

4. Enter (nested) loop covering all integration points

(a) Look up the sampling position of the current point (£p, nq) (s, t) and their weights (weigp)

(b) Call shape function routine sfr given £p, nq which will return the shape function Ni

(sfr) and the derivatives -gf- and (deriv) at the point £p,r]q

(c) Call another subroutine (jacob), given Nf, and at point £p,r)q will return dN? dN?

cartesian shape function derivatives -¿^f- and (cartd), the Jacobian matrix Je (jacm), its inverse Je_1 (jaci), its determinant det Je (djac) and the x and y coordinates all at the point £p,nq

(d) Call strain matrix (bmatps) routine, given Nf, , q^ , at £p, nq will return the strain matrix B^ (bmat)

(e) Call a routine (dbmat) to evaluate DeB^ (dbmat)

(f) Evaluate BeTDeBedetJe x integration weights and assemble them into the element stiffness matrix Kj

(g) Assemble DeBe (smat into a stress array for later evaluation of stresses from the nodal displacements.

5. Write Stiffness matrix

Suggested list of variables:

idime, ndime idofn,ndofn ielem,nelem igaus,jgaus,ngaus inode, nnode kgasp,ngasp type poiss young

Index, Number of dimensions (2 for 2D) Index, Number of degree of freedom per node Index, number of elements Index, Index, Number of Gauss rule adopted Index, number of nodes per element Kounter, number of Gauss points used 1 for plane stress; 2 for plane strain Poisson's ratio Young's modulus elcod(ndime,nnode)

gpcod(ndime,ngasp)

posgp(mgaus) weigp(mgaus) shape(nnode)

deriv(ndime,nnode)

cartd(ndime,nnode)

djacb jacm(ndime,ndime) jaci(ndime,ndime) bmat(nstre,nevab)

eration

Local array of nodal cartesian coordinates of the element currently under consid-

x(£i,Vi) ••• x(£s,ns) y(£i,m) ••• y(£.8,m)

^ coordinate of sampling point n coordinate of sampling point

Local array of cartesian coordinates of the Gauss points for element currently under consideration x(£ci ,VGi y(S,Gi ,nGi £ coordinates of Gauss point Weight factor for Gauss point Shape function associated with

nG5 VG5

each node of current element at (S,P,nP)

Shape function

Cartesian shape current element

SN-, Oy derivative

function sampled

Ox y^p onh

Oy at sampling point (£p ,np) within the element np) np) j derivatives associated with the nodes of the at any point (£,p,np) within the element Vp)

Determinant of the Jacobian matrix sampled at any point (S,P,nP) within the element

Jacobian matrix at sampling point Inverse of Jacobian matrix at sampling point

Element strain matrix at any point within the element B = [ B1 B2

ONj Ox 0

where Bi =

ON, dbmat(istre,ievab) stores DB

ONj Oy

Oy ONj

9.6.2 MATLAB Code

9.6.2.1 stiff.m function KELEM = stiff(ngaus,posgp,weigp,type,nelem,lnods,coord,nnode) %---------------------------------------------------------------------------

% The purpose of this function is calculate element stiffness matricies for % bilinear and biquadratic isoparemtric elements using Gaussian integration.

% Functions called by this function are: dmat,sfr,jacob,bmatps %---------------------------------------------------------------------------

% VARIABLES

%---------------------------------------------------------------------------

% nelem Global variable NELEM

% nnode Global variable NNODES

% posgp Global variable POSGP

% weigp Global variable WEIGP

% ngaus Global variable NGAUS

% type Global variable TYPE

Number of columns in the element stiffness matrix Position indicator for element stiffness matrix Current element for formation of stiffness matrix Modulus of elasticity for current element Poisson's ration for current element Constitutive matix Matrix of element coordinates Counter

Element stiffness matrix

Element stiffness matricies for all elements returned by function

Current integration position

Current integration position

Shape function at current point

Derivative of shape function at current point

Cartesian shape function derivatives

Jacobian matrix

Jacobian matrix inverse

Determinant of Jacobian matrix x and y coordinates at the current point in the element Strain matrix [B]

Strain matrix * constitutive matrix [B]*[D]

Global variable LNODS Global variable COORD

%---------------------------------------------------

tic fprintf('Calculating ELEMENT STIFFNESS matricies\n') stifsize = 2*nnode; nrowcount = stifsize-1; for ielem = 1:nelem elmt = lnods(ielem,1); young = lnods(ielem,2);

poiss = lnods(ielem,3); %-----------------------------------------

% Extract element coordinates %

% [ Y1 Y2 Y3 . . . Yn] %-----------------------------------------

elcod = zeros(2,nnode); for inode = 1:nnode row = find(coord(:,1:1)==lnods(ielem,inode+3)); elcod(:,inode:inode) = coord(row:row,2:3)';

% Constitutive matrix %-----------------------------------------

D = dmat(young,poiss,type); %-----------------------------------------

% Element stiffness matrix - element ielem %-----------------------------------------

kelem = zeros(stifsize); for igaus = 1:ngaus

% Extract material constants from lnods %---------------------------------------

for jgaus = 1:ngaus s = posgp(igaus); t = posgp(jgaus); W = weigp(igaus)*weigp(jgaus); [shape,deriv] = sfr(s,t,nnode);

[cartd,jacm,jaci,djac,xy] = jacob(shape,deriv,elcod); [bmat,dbmat] = bmatps(shape,cartd,D); kelem = kelem + bmat'*dbmat*djac;

end end

Store element stiffness matricies in as a stack in a single matrix :

kelem(nelem)

startrow = stifsize*ielem - nrowcount; endrow = ielem*stifsize; KELEM(startrow:endrow,:) = kelem;

fprintf(1,'Elapsed time for this operation =/3.4fsec\n\n',t)

9.6.2.2 dmat.m function D = dmat(young,poiss,type) /-----------------------------------------------------------------

/ This function calculates the constitutive matrix for an element /-----------------------------------------------------------------

/ VARIABLES

/-----------------------------------------------------------------

/ young Young's modulus/Modulus of elasticity

/ poiss Poisson's ratio

/ E Modulus of elasticity

/ type type of problem - plain stress = 1, plain strain = 2

/ D Constitutive matrix returned by function /-----------------------------------------------------------------

/fprintf('Calculating element constitutive matrix\n')

else

 1 v 0; v 1 0; 0 0 (1-v)/2]

°/fprintf(1,'Elapsed time for this operation =°/3.4fsec\n\n',t)

9.6.2.3 sfr.m function [shape,deriv] = sfr(s,t,nnode) /---------------------------------------------------------------------------------

/ This function calculates the shape function and derivative for the current node /---------------------------------------------------------------------------------

% VARIABLES /---------------------------------------------------------------------------------

Shape function returned by function

Derivative of shape function returned by function

Number of nodes per element

Natural coordinate (xi) of sampling point - horizontal Natural coordinate (eta) of sampling point - vertical

/---------------------------------------------------------

/fprintf('Calculating shape functions and derivatives\n') /-------------------

 mode == 9 N9 = (1-s"2)*(1-t"2); N8 = ,5*(1-s)*(1-t"2) -,5*N9; N7 = ,5*(1-s"2)*(1+t) -,5*N9; N6 = ,5*(1+s)*(1-t"2) -,5*N9; N5 = ,5*(1-s"2)*(1-t) -,5*N9; N4 = ,25*(1-s)*(1+t) -,5*N7 - .5*N8 - .25*N9 N3 = ,25*(1+s)*(1+t) -,5*N6 - .5*N7 - .25*N9 N2 = ,25*(1+s)*(1-t) -,5*N5 - .5*N6 - .25*N9 N1 = ,25*(1-s)*(1-t) -,5*N5 - .5*N8 - .25*N9
 dN9ds = -2*s*(1-t"2); dN9dt = -2*t*(1-s"2); dN8ds = -.5*(1-t"2) - 5*dN9ds dN8dt = -t*(1-s) - 5*dN9dt dN7ds = -s*(1+t) - 5*dN9ds dN7dt = .5*(1-s"2) - 5*dN9dt dN6ds = .5*(1-t"2) - 5*dN9ds dN6dt = -t*(1+s) - 5*dN9dt dN5ds = -s*(1-t) - 5*dN9ds dN5dt = -.5*(1-s"2) - 5*dN9dt dN4ds = -.25*(1+t) - 5*dN7ds - ,5*dN8ds - .25*dN9ds dN4dt = .25*(1-s) - 5*dN7dt - ,5*dN8dt - .25*dN9dt dN3ds = .25*(1+t) - 5*dN6ds - ,5*dN7ds - .25*dN9ds dN3dt = .25*(1+s) - 5*dN6dt - ,5*dN7dt - .25*dN9dt dN2ds = .25*(1-t) - 5*dN5ds - ,5*dN6ds - .25*dN9ds dN2dt = -.25*(1+s) - 5*dN5dt - ,5*dN6dt - .25*dN9dt dN1ds = -.25*(1-t) - 5*dN5ds - ,5*dN8ds - .25*dN9ds dN1dt = -.25*(1-s) - 5*dN5dt - ,5*dN8dt - .25*dN9dt

deriv = [dNlds dN2ds dN3ds dN4ds dN5ds dN6ds dN7ds dN8ds dN9ds;

dNldt dN2dt dN3dt dN4dt dN5dt dN6dt dN7dt dN8dt dN9dt]; %-------------------

elseif nnode == 8

 N7 = 5*(1-s"2)*(1+t) ; N6 = 5*(1+s)*(1-t"2) ; N5 = 5*(1-s"2)*(1-t) ; N4 = 25*(1-s)*(1+t) -,5*N7 .5*N8; N3 = 25*(1+s)*(1+t) -,5*N6 .5*N7; N2 = 25*(1+s)*(1-t) -,5*N5 .5*N6; N1 = 25*(1-s)*(1-t) -,5*N5 .5*N8; shape = [N1 N2 N3 N4 N5 N6 N7 N8]'; dN8ds = -,5*(1-t"2); dN8dt = -t*(1-s); dN7ds = -s*(1+t); dN7dt = ,5*(1-s"2); dN6ds = ,5*(1-t"2); dN6dt = -t*(1+s); dN5ds = -s*(1-t); dN5dt = -,5*(1-s"2); dN4ds = -,25*(1+t) - .5*dN7ds - ,5*dN8ds dN4dt = ,25*(1-s) - .5*dN7dt - ,5*dN8dt dN3ds = ,25*(1+t) - .5*dN6ds - ,5*dN7ds dN3dt = ,25*(1+s) - .5*dN6dt - ,5*dN7dt dN2ds = ,25*(1-t) - .5*dN5ds - ,5*dN6ds dN2dt = -,25*(1+s) - .5*dN5dt - ,5*dN6dt dN1ds = -,25*(1-t) - .5*dN5ds - ,5*dN8ds dN1dt = -,25*(1-s) - .5*dN5dt - ,5*dN8dt

deriv = [dNlds dN2ds dN3ds dN4ds dN5ds dN6ds dN7ds dN8ds;

dNldt dN2dt dN3dt dN4dt dN5dt dN6dt dN7dt dN8dt]; %-------------------

else

 dN4ds = -.25*(1+t); dN4dt = .25*(1-s); dN3ds = .25*(1+t); dN3dt = .25*(1+s); dN2ds = .25*(1-t); dN2dt = -.25*(1+s); dN1ds = -.25*(1-t); dN1dt = -.25*(1-s);

deriv = [dN1ds dN2ds dN3ds dN4ds;

dN1dt dN2dt dN3dt dN4dt];

%fprintf(1,'Elapsed time for this operation =%3.4fsec\n\n',t)

9.6.2.4 jacob.m function [cartd,jacm,jaci,djac,xy] = jacob(shape,deriv,elcod)

% This function calculates the cartesian shape function derivatives % the Jacobian matrix, Jacobian inverse and Jacobian determinant

### VARIABLES

shape Shape function at current point deriv Derivative of shape function at current point cartd Cartesian shape function derivatives returned by function jacm Jacobian matrix returned by function jaci Jacobian matrix inverse returned by function djac Determinant of Jacobian matrix returned by function xy x and y coordinates at the current point in the element ic fprintf('Calculating Jacobian matrix\n')

The cartesian shape function derivatives {cartd} are given by:

{dN/dx} -1 {dN/ds} {cartd} = { } = [J] { } {dN/dy} {dN/dt}

Start by calculating Jacobian [J] = jacm:

jacm = deriv*elcod'; jaci = inv(jacm); djac = det(jacm); cartd = jaci*deriv; xy = elcod*shape; t = toc;

%fprintf(1,'Elapsed time for this operation =%3.4fsec\n\n',t)

9.6.2.5 bmatps.m function [bmat,dbmat] = bmatps(shape,cartd,D) %-------------------------------------------------

% This function calculates the strain matrix B %-------------------------------------------------

% VARIABLES %-------------------------------------------------

% shape Shape function at current point

% cartd Cartesian shape function derivatives

% bmat Strain matrix returned by function

% dbmat Strain matrix * constitutive matrix D

%--------------------------------------------------------------------------

%fprintf('Calculating strain matrix [B]\n') numcols = 2*length(cartd); bmat = zeros(3,numcols); cartdcol = 0;

for ibmatcol = 1:2:numcols cartdcol = cartdcol+1;

bmat(:,ibmatcol:ibmatcol+1) = [cartd(1,cartdcol) 0 ;

0 cartd(2,cartdcol);

cartd(2,cartdcol) cartd(1,cartdcol)];

%fprintf(1,'Elapsed time for this operation =%3.4fsec\n\n',t)

0 0