Code Collection for AER 1403 Project 2
Jian Wu 1001040517
This project uses 19 M files and they are arranged in a sequence: Main entry, Meshing, Stiffness
Matrix, Solution and Post-process.
Main Entry
% project2.m
% 2D problem: thin plate with a circular hole
% jian wu 2014
% clear memory
clc;close all;clear all;colordef white;clf
% delete previous diary files (if there is)
delete('physicalCoord [Link]','[Link]')
tol=1e-6; % tolerance for find function
% materials
E = 75e9; % Pa
poisson = 0.30;
thickness = 0.001; % meter
% matriz C
C = E/(1-poisson^2)*[1 poisson 0;poisson 1 0;0 0 (1-poisson)/2];
% load
P = 200e6; % Pa
% Mesh generation
Lx = 0.8; % meter
Ly = 0.8; % meter
global R;
R = 0.05; % meter
numberElementsX = 12;
numberElementsY = 12;
numberElements = numberElementsX*numberElementsY;
[nodeCoordinates, elementNodes] = ...
quadPlateMesh(Lx,Ly,numberElementsX,numberElementsY); % Mesh for a fourth
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
drawingMesh(nodeCoordinates,elementNodes,'Q4','k-');
numberNodes=size(xx,1);
% create a table of physical coordinates
createTable(nodeCoordinates,numberElements,elementNodes,numberNodes,5)
% GDof: global number of degrees of freedom
GDof=2*numberNodes;
% calculation of the system stiffness matrix
stiffness=formStiffness2D(GDof,numberElements,...
elementNodes,numberNodes,nodeCoordinates,C,1,thickness);
% boundary conditions
fixedNodeX=find(abs(nodeCoordinates(:,1))<tol); % fixed in XX
fixedNodeY=find(abs(nodeCoordinates(:,2))<tol); % fixed in YY
prescribedDof=[fixedNodeX; fixedNodeY+numberNodes];
% force vector or element load vector (distributed load applied at yy=Ly/2)
force=zeros(GDof,1); % second numberNodes elements for y direction
% upperBord=find(nodeCoordinates(:,2)==Ly/2) % miss node. use tolerance
% instead
upperBord=find(abs(nodeCoordinates(:,2)-Ly/2)<tol);
middle=upperBord(2:end-1);
force(middle+numberNodes)=P*thickness*(xx(middle-1)-xx(middle+1))/2;
force(upperBord(1)+numberNodes)=P*thickness*(xx(upperBord(1))-xx(upperBord(2)))/2;
force(upperBord(end)+numberNodes)=P*thickness*(xx(upperBord(end-1))-xx(upperBord(end)))/
2;
% solution
displacements=solution(GDof,prescribedDof,stiffness,force);
% displacements
disp('Displacements')
jj=1:GDof; format
f=[jj; displacements'];
fprintf('node U\n');
fprintf('%3d %12.8f\n',f)
UX=displacements(1:numberNodes);
UY=displacements(numberNodes+1:GDof);
% create a table of displacement
createTable([UX UY],numberElements,elementNodes,numberNodes,1)
% scaleFactor depends on displacements
scaleFactor=computeScaleFactor(numberNodes,nodeCoordinates,UX,UY);
% deformed shape
figure
drawingField(nodeCoordinates+scaleFactor*[UX UY],elementNodes,'xReflection',UX);%U XX
hold on
drawingMesh(nodeCoordinates+scaleFactor*[UX UY],elementNodes,'Q4','k-');
drawingMesh(nodeCoordinates,elementNodes,'Q4','k--');
colorbar
title('U XX (one deformed shape)')
axis off
% stresses at nodes
[strain,vonMisesStress]=stresses2D(GDof,numberElements,elementNodes,numberNodes,...
nodeCoordinates,displacements,UX,UY,C);
% drawing strain fields
% on top of the deformed shape
figure
drawingField(nodeCoordinates+scaleFactor*[UX
UY],elementNodes,'fullRotation',strain(:,:,1));%epsilon XX
hold on
drawingMesh(nodeCoordinates+scaleFactor*[UX UY],elementNodes,'Q4','k-');
drawingMesh(nodeCoordinates,elementNodes,'Q4','k--');
colorbar
title('Epsilon X strain (on deformed shape)')
axis off
figure
drawingField(nodeCoordinates+scaleFactor*[UX
UY],elementNodes,'fullRotation',strain(:,:,2));%epsilon YY
hold on
drawingMesh(nodeCoordinates+scaleFactor*[UX UY],elementNodes,'Q4','k-');
drawingMesh(nodeCoordinates,elementNodes,'Q4','k--');
colorbar
title('Epsilon Y strain (on deformed shape)')
axis off
figure
drawingField(nodeCoordinates+scaleFactor*[UX
UY],elementNodes,'fullRotation',strain(:,:,3));%epsilon XY
hold on
drawingMesh(nodeCoordinates+scaleFactor*[UX UY],elementNodes,'Q4','k-');
drawingMesh(nodeCoordinates,elementNodes,'Q4','k--');
colorbar
title('Epsilon XY strain (on deformed shape)')
axis off
% drawing stress fields
% on top of the deformed shape
figure
drawingField(nodeCoordinates+scaleFactor*[UX
UY],elementNodes,'fullRotation',vonMisesStress);%sigma XX
hold on
drawingMesh(nodeCoordinates+scaleFactor*[UX UY],elementNodes,'Q4','k-');
drawingMesh(nodeCoordinates,elementNodes,'Q4','k--');
colorbar
title('von Mises stress (on deformed shape)')
axis off
% Node average displacement
devDisp=sqrt(sum(UX.^2+UY.^2)/numberNodes)
% Elements average stress
devStressX=mean(mean(vonMisesStress.^2,2))
Meshing
% quadPlateMesh.m
% To mesh a plate with hole at center using Transfinite Interpolation (TFI)
function [nodeCoordinates, elementNodes] = ...
quadPlateMesh(Lx,Ly,numberElementsX,numberElementsY)
if any(mod(numberElementsY,2)~=0)
error('Number of elments in Y must be even');
end
% Dimensions of the plate
L = Lx ; % Length of the plate
B = Ly ; % Breadth of the plate
% Number of discretizations along xi and eta axis
m = numberElementsY/2+1 ; % eta axis
n = numberElementsX+1 ; % xi axis
%
% Model plate as two regions which lie in first quadrant
% global R theta;
global R theta;
% R = 0.08 ; % Radius of the hole at center
%%%%%%%%%%%%%%%%%%%%%%%Dont change from
here%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta = pi/2 ; % Quarter angle of the hole
global O P1 P2 P3 P4 P5 CMP ;
O = [0. 0.] ; % Centre of plate and hole
P1 = [R 0.] ; % Edge of the hole and plate
P2 = [L/2 0.] ; % Edge of the plate
P3 = [L/2 B/2] ; % Edge of the plate
P4 = [0. B/2] ; % Edge of the plate
P5 = [0. R] ; % Edge of the hole and plate
CMP = [R*cos(theta/2.) R*sin(theta/2.)] ;
% discretize along xi and eta axis
xi = linspace(0.,1,m) ;
eta = linspace(0.,1.,n) ;
% Number of Domains
Domain = 2 ;
DX = cell(1,Domain) ;
DY = cell(1,Domain) ;
for d = 1:Domain % Loop for two domains lying in first coordinate
% Initialize matrices in x and y axis
X = zeros(m,n) ;
Y = zeros(m,n) ;
for i = 1:m
Xi = xi(i) ;
for j = 1:n
Eta = eta(j) ;
% Transfinite Interpolation
XY = (1-Eta)*Xb(Xi,d)+Eta*Xt(Xi,d)+(1-Xi)*Xl(Eta,d)+Xi*Xr(Eta,d)......
-(Xi*Eta*Xt(1,d)+Xi*(1-Eta)*Xb(1,d)+Eta*(1-Xi)*Xt(0,d)+(1-Xi)*(1-Eta)*Xb(0,d)) ;
X(i,j) = XY(1) ;
Y(i,j) = XY(2) ;
end
end
DX{d} = X ;
DY{d} = Y ;
end
% Arrange the coordinates for each domain
X1 = DX{1} ; Y1 = DY{1} ; % Grid for first domain
X2 = DX{2} ; Y2 = DY{2} ; % Grid for second domain
X = [X1 ;X2(m-1:-1:1,:)] ; % Merge both the domains
Y = [Y1 ;Y2(m-1:-1:1,:)] ;
nodeCoordinates = [X(:) Y(:)];
numberElements = numberElementsX*numberElementsY;
elementNodes = zeros(numberElements, 4);
% First node for every element
node1 = 0;
for i=1:numberElements
node1 = floor((i-1)/numberElementsY)+i;
% order nodes counterclockwise
elementNodes(i,:) = [node1 node1+numberElementsY+1 node1+numberElementsY+2
node1+1];
end
end % end of function quadPlateMesh
% Xb.m
function xyb = Xb(s,Domain)
global R ;
r=R;
global O P1 P2 P3 P4 P5 CMP ;
switch Domain
case 1
x = O(1)+r*cos(pi/4*s) ;
y = O(2)+r*sin(pi/4*s) ;
case 2
x = O(1)+r*cos(pi/4*s) ;
y = O(2)+r*sin(pi/4*s) ;
end
xyb = [x ; y] ;
% Xl.m
function xyl = Xl(s,Domain)
global O P1 P2 P3 P4 P5 CMP ;
switch Domain
case 1
x = P1(1)+(P2(1)-P1(1))*s ;
y = P1(2)+(P2(2)-P1(2))*s ;
case 2
x = P5(1)+(P4(1)-P5(1))*s ;
y = P5(2)+(P4(2)-P5(2))*s ;
end
xyl = [x ; y] ;
% Xr.m
function xyr = Xr(s,Domain)
global O P1 P2 P3 P4 P5 CMP ;
switch Domain
case 1
x = CMP(1)+(P3(1)-CMP(1))*s ;
y = CMP(2)+(P3(2)-CMP(2))*s ;
case 2
x = CMP(1)+(P3(1)-CMP(1))*s ;
y = CMP(2)+(P3(2)-CMP(2))*s ;
end
xyr = [x ; y] ;
% Xt.m
function xyt = Xt(s,Domain)
global O P1 P2 P3 P4 P5 CMP ;
switch Domain
case 1
x = P2(1)+(P3(1)-P2(1))*s ;
y = P2(2)+(P3(2)-P2(2))*s ;
case 2
x = P3(1)+(P4(1)-P3(1))*s ;
y = P3(2)+(P4(2)-P3(2))*s ;
end
xyt = [x ; y] ;
Stiffness Matrix
% formStiffness2D.m
%................................................................
function [stiffness,mass]=formStiffness2D(GDof,numberElements,...
elementNodes,numberNodes,nodeCoordinates,C,rho,thickness)
% compute stiffness matrix (and mass matrix)
% for plane stress Q4 elements
stiffness=zeros(GDof);
mass=zeros(GDof);
% 2 by 2 quadrature
[gaussWeights,gaussLocations]=gaussQuadrature('complete');
for e=1:numberElements
indice=elementNodes(e,:);
elementDof=[ indice indice+numberNodes ];
ndof=length(indice);
% cycle for Gauss point
for q=1:size(gaussWeights,1)
GaussPoint=gaussLocations(q,:);
xi=GaussPoint(1);
eta=GaussPoint(2);
% shape functions and derivatives
[shapeFunction,naturalDerivatives]=shapeFunctionQ4(xi,eta);
% Jacobian matrix, inverse of Jacobian,
% derivatives w.r.t. x,y
[Jacob,invJacobian,XYderivatives]=...
Jacobian(nodeCoordinates(indice,:),naturalDerivatives);
% B matrix
B=zeros(3,2*ndof);
B(1,1:ndof) = XYderivatives(:,1)';
B(2,ndof+1:2*ndof) = XYderivatives(:,2)';
B(3,1:ndof) = XYderivatives(:,2)';
B(3,ndof+1:2*ndof) = XYderivatives(:,1)';
% stiffness matrix
stiffness(elementDof,elementDof)=...
stiffness(elementDof,elementDof)+...
B'*C*thickness*B*gaussWeights(q)*det(Jacob);
% mass matrix
mass(indice,indice)=mass(indice,indice)+...
shapeFunction*shapeFunction'*...
rho*thickness*gaussWeights(q)*det(Jacob);
mass(indice+numberNodes,indice+numberNodes)=...
mass(indice+numberNodes,indice+numberNodes)+...
shapeFunction*shapeFunction'*...
rho*thickness*gaussWeights(q)*det(Jacob);
end
end
% shapeFunctionQ4.m
% .............................................................
function [shape,naturalDerivatives]=shapeFunctionQ4(xi,eta)
% shape function and derivatives for Q4 elements
% shape : Shape functions
% naturalDerivatives: derivatives w.r.t. xi and eta
% xi, eta: natural coordinates (-1 ... +1)
shape=1/4*[ (1-xi)*(1-eta);(1+xi)*(1-eta);
(1+xi)*(1+eta);(1-xi)*(1+eta)];
naturalDerivatives=...
1/4*[-(1-eta), -(1-xi);1-eta, -(1+xi);
1+eta, 1+xi;-(1+eta), 1-xi];
end % end function shapeFunctionQ4
% gaussQuadrature.m
% .............................................................
function [weights,locations]=gaussQuadrature(option)
% Gauss quadrature for Q4 elements
% option 'complete' (2x2)
% option 'reduced' (1x1)
% locations: Gauss point locations
% weights: Gauss point weights
switch option
case 'complete'
locations=...
[ -0.577350269189626 -0.577350269189626;
0.577350269189626 -0.577350269189626;
0.577350269189626 0.577350269189626;
-0.577350269189626 0.577350269189626];
weights=[ 1;1;1;1];
case 'reduced'
locations=[0 0];
weights=[4];
end
end % end function gaussQuadrature
% Jacobian.m
% .............................................................
function [JacobianMatrix,invJacobian,XYDerivatives]=...
Jacobian(nodeCoordinates,naturalDerivatives)
% JacobianMatrix : Jacobian matrix
% invJacobian : inverse of Jacobian Matrix
% XYDerivatives : derivatives w.r.t. x and y
% naturalDerivatives : derivatives w.r.t. xi and eta
% nodeCoordinates : nodal coordinates at element level
JacobianMatrix=nodeCoordinates'*naturalDerivatives;
invJacobian=inv(JacobianMatrix);
XYDerivatives=naturalDerivatives*invJacobian;
end % end function Jacobian
Solution
% solution.m
%................................................................
function displacements=solution(GDof,prescribedDof,stiffness,force)
% function to find solution in terms of global displacements
activeDof=setdiff([1:GDof]', ...
[prescribedDof]);
U=stiffness(activeDof,activeDof)\force(activeDof);
displacements=zeros(GDof,1);
displacements(activeDof)=U;
% stresses2D.m
%................................................................
function [avgStrain,vonMisesStress] = stresses2D(GDof,numberElements,...
elementNodes,numberNodes,nodeCoordinates,...
displacements,UX,UY,C)
% 2 by 2 quadrature
[gaussWeights,gaussLocations]=gaussQuadrature('complete');
% Nodal stresses
nodalStrain=zeros(numberElements,size(elementNodes,2),3);
avgStrain=zeros(numberElements,size(elementNodes,2),3);
vonMisesStress=zeros(numberElements,size(elementNodes,2));
% strains and stresses at Gauss nodes
gaussStrain=zeros(numberElements,size(elementNodes,2),3);
gaussStress=zeros(numberElements,size(elementNodes,2),3);
stressPoints=[-1 -1;1 -1;1 1;-1 1];
for e=1:numberElements
indice=elementNodes(e,:);
elementDof=[ indice indice+numberNodes ];
nn=length(indice);
for q=1:size(gaussWeights,1)
pt=gaussLocations(q,:);
wt=gaussWeights(q);
xi=pt(1);
eta=pt(2);
% shape functions and derivatives
[shapeFunction,naturalDerivatives]=shapeFunctionQ4(xi,eta);
% Jacobian matrix, inverse of Jacobian,
% derivatives w.r.t. x,y
[Jacob,invJacobian,XYderivatives]=...
Jacobian(nodeCoordinates(indice,:),naturalDerivatives);
% B matrix
B=zeros(3,2*nn);
B(1,1:nn) = XYderivatives(:,1)';
B(2,nn+1:2*nn) = XYderivatives(:,2)';
B(3,1:nn) = XYderivatives(:,2)';
B(3,nn+1:2*nn) = XYderivatives(:,1)';
% element deformation
dofStrain=B*displacements(elementDof);
gaussStrain(e,q,:)=dofStrain;
gaussStress(e,q,:)=C*dofStrain;
end % loop end of q
% Calculating the Von Mises stress
vonMisesStress(e,:)=VonmisesStresses(squeeze(gaussStress(e,:,:)))';
% Extrapolating stress at Gaussian points to get nodal stress values
nodalStrain(e,:,:)=Extrapolation(squeeze(gaussStrain(e,:,:)));
end
% nodal strain averaging
type = 'average' ; % type = 'sum' ;
avgStrain=NodalAveraging(nodalStrain,numberElements,elementNodes,numberNodes,type);
end % end function stresses2D
% VonmisesStresses.m
function [vonmises] = VonmisesStresses(stressGP)
% To find the Principal Stresses
toumax = sqrt((0.5*(stressGP(:,1)-stressGP(:,2))).^2+stressGP(:,3).^2);
sigma1 = 0.5*(stressGP(:,1)+stressGP(:,2))+toumax ;
toumin = -sqrt((0.5*(stressGP(:,1)-stressGP(:,2))).^2+stressGP(:,3).^2);
sigma2 = 0.5*(stressGP(:,1)+stressGP(:,2))+toumin ;
%shear = [toumax toumin] ;
%principal = [sigma1 sigma2] ;
% To find the Von-Mises stresses
vonmises = sqrt(0.5*((sigma1-sigma2).^2+sigma2.^2+sigma1.^2));
end % end function VonmisesStresses
% NodalAveraging.m
function [avgData] =
NodalAveraging(rawData,numberElements,elementNodes,numberNodes,type)
%--------------------------------------------------------------------------
% Purpose:
% Nodal averaging of elemental nodal data
% Synopsis :
% avgData =
NodalAveraging(rawData,numberElements,elementNodes,numberNodes,type)
% Variable Description:
% rawData - elemental nodal data, 16*4*3
% type = average does the averaging
% = sum adds the data
% avgData - averaged nodal data
%--------------------------------------------------------------------------
nnode = numberNodes;
avg1=zeros(nnode,1);
avg2=zeros(nnode,1);
avg3=zeros(nnode,1);
avgData = zeros(numberElements,4,3);
for ind = 1:nnode
[r c] = find(elementNodes==ind);
switch type
case 'average'
share = length(r) ;
case 'sum'
share = 1 ;
end
sum1=0;
sum2=0;
sum3=0;
for i=1:length(r)
sum1=rawData(r(i),c(i),1)+sum1;
sum2=rawData(r(i),c(i),2)+sum2;
sum3=rawData(r(i),c(i),3)+sum3;
end
avg1(ind)=sum1/share;
avg2(ind)=sum2/share;
avg3(ind)=sum3/share;
for j=1:length(r)
avgData(r(j),c(j),:)=[avg1(ind) avg2(ind) avg3(ind)];
end
end
% Elements average strain xx
devStrainX=mean(avg1.^2)
end % end function NodalAveraging
% Extrapolation.m
% Procedure to get nodal stresses using Extrapolation method
function [sigma] = Extrapolation(stressGP)
%--------------------------------------------------------------------------
% Purpose : Extrapolating stress at Gaussian points to get nodal stress
% values
%
% Synopsis : sigma = extrapolation(stressGP)
%
% Variable Description :
% sigma - Stress at nodal points
% stressGP - Stress at gaussian points
%--------------------------------------------------------------------------
% Extrapolation matrix obtained from Shape Functions
% Stress points are at (-1,-1),(-1,1),(1,1),1,-1)
explmt = [1+sqrt(3)/2 -1/2 1-sqrt(3)/2 -1/2;
-1/2 1-sqrt(3)/2 -1/2 1+sqrt(3)/2 ;
1-sqrt(3)/2 -1/2 1+sqrt(3)/2 -1/2;
-1/2 1+sqrt(3)/2 -1/2 1-sqrt(3)/2 ] ;
sigmax = explmt*stressGP(:,1) ;
sigmay = explmt*stressGP(:,2) ;
sigmaxy = explmt*stressGP(:,3) ;
sigma = [sigmax sigmay sigmaxy];
Post-process
% drawingMesh.m
function drawingMesh(nodeCoordinates,elementNodes,option,lineStyle)
% drawingMesh: To plot structured grid.
switch option
case 'Q4'
% Plot grid for quadrilateral element Q4
set(gcf,'color','w') ;
axis equal
axis off
box on
hold on
% Plot other domains of plate by imaging coordinates
vec = [1 1 ; -1 1 ; -1 -1 ; 1 -1] ;
numberElements=size(elementNodes,1);
for e=1:numberElements
indice=elementNodes(e,:);
XX=[nodeCoordinates(indice,1);nodeCoordinates(indice(1),1)];
YY=[nodeCoordinates(indice,2);nodeCoordinates(indice(1),2)];
% Plot internal grid lines
for quadrant=1:4
plot(vec(quadrant,1)*XX,vec(quadrant,2)*YY,lineStyle,'linewidth',1);
hold on
end
end
case 'L2'
% future extension for 2D truss problem
end
hold off
end % end function drawingMesh
% drawingField.m
function drawingField(nodeCoordinates,elementNodes,symOption,component)
% drawingField: drawing stress fields on top of the deformed shape
% symOption: specify the type of symmetry for the component
nel = length(elementNodes) ; % number of elements
nnode = length(nodeCoordinates) ; % total number of nodes in system
% nnel = size(nodes,2)-1; % number of nodes per element 4
%
% Initialization of the required matrices
X = zeros(4,nel) ;
Y = zeros(4,nel) ;
profile = zeros(4,nel) ;
% Plot other domains of plate by imaging coordinates
vec = [1 1 ; -1 1 ; -1 -1 ; 1 -1] ;
for iel=1:nel
for i=1:4
nd(i)=elementNodes(iel,i); % extract connected node for (iel)-th element
X(i,iel)=nodeCoordinates(nd(i),1); % extract x value of the node
Y(i,iel)=nodeCoordinates(nd(i),2); % extract y value of the node
end
profile(:,iel) = component(nd') ; % extract component value of the node
end
% Plot grid for quadrilateral element Q4
set(gcf,'color','w') ;
axis equal
axis off
box on
hold on
switch symOption
case {'xReflection'}
% Plotting the FEM mesh and profile of the given component
for quadrant=1:2
plot(vec(quadrant,1)*X,vec(quadrant,2)*Y,'k-')
fill(vec(quadrant,1)*X,vec(quadrant,2)*Y,profile)
hold on
end
for quadrant=3:4
plot(vec(quadrant,1)*X,vec(quadrant,2)*Y,'k-')
fill(vec(quadrant,1)*X,vec(quadrant,2)*Y,-profile)
hold on
end
case {'yReflection'}
% Plotting the FEM mesh and profile of the given component
for quadrant=1:3:4
plot(vec(quadrant,1)*X,vec(quadrant,2)*Y,'k-')
fill(vec(quadrant,1)*X,vec(quadrant,2)*Y,profile)
hold on
end
for quadrant=2:3
plot(vec(quadrant,1)*X,vec(quadrant,2)*Y,'k-')
fill(vec(quadrant,1)*X,vec(quadrant,2)*Y,-profile)
hold on
end
case {'fullRotation'}
% Plotting the FEM mesh and profile of the given component
for quadrant=1:4
plot(vec(quadrant,1)*X,vec(quadrant,2)*Y,'k-')
fill(vec(quadrant,1)*X,vec(quadrant,2)*Y,profile)
hold on
end
end
title('Profile of component on Mesh') ;
axis off ;
hold off
end % end function drawingField
% createTable.m
function createTable(input,numberElements,elementNodes,numberNodes,type)
%--------------------------------------------------------------------------
% Purpose:
% Print outputs in tabular form
% Synopsis :
% createTable(input,numberElements,numberNodes,type)
% Variable Description:
% input - displacement or Stress matrix
% type = 1 displacement
% =2 Von Mises stress
% =3 Stress matrix
% =4 Averaged stress matrix
% =5 physical coordinates
%--------------------------------------------------------------------------
nel = numberElements ; % number of elements
nnel=4; % number of nodes per element
ndof=2; % number of dofs per node (UX,UY)
nnode = numberNodes ; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
if type == 1
diary('[Link]')
UX = input(1:2:sdof) ;
UY = input(2:2:sdof) ;
display('======PER NODE DISPLACEMENTS====== ');
fprintf('================================================\n');
fprintf(' node UX UY \n')
fprintf('================================================\n');
for i = 1:nnode
node = i;
ux = UX(i);
uy = UY(i);
fprintf(' %5d %+6.5e %+6.5e \n',node,ux,uy);
end
fprintf('=================================================\n');
diary off
elseif type == 2
display('=====VON MISES STRESSES PER NODE=======')
fprintf('=======================================\n');
fprintf('element node Von Mises \n')
fprintf('======================================\n');
for ielp = 1:nel
pos = 4*(ielp-1)+(1:4);
for i = 1:nnel ;
nodeno = nodes(ielp,i);
vmis = input(pos(i));
fprintf('%5d %5d %+6.5e \n',ielp,nodeno,vmis);
end
end
fprintf('=======================================\n');
elseif type == 3
fprintf('===============================================================\n');
fprintf('element node sigmax sigmay sigmaxy \n')
fprintf('===============================================================\n');
for ielp = 1:nel
pos = 4*(ielp-1)+(1:4);
for i = 1:nnel ;
nodeno = nodes(ielp,i);
sigmaX = input(pos(i),1);
sigmaY = input(pos(i),2);
sigmaXY = input(pos(i),3);
fprintf('%5d %5d %+6.5e %+6.5e %+6.5e
\n',ielp,nodeno,sigmaX,sigmaY,sigmaXY);
end
end
fprintf('===============================================================\n');
elseif type == 4
fprintf('==========================================================\n');
fprintf(' node sigmax sigmay sigmaxy \n')
fprintf('==========================================================\n');
for i = 1:nnode
nodeno = i;
sigmaX = input(i,1) ;
sigmaY = input(i,2) ;
sigmaXY = input(i,3) ;
fprintf('%5d %+6.5e %+6.5e %+6.5e \n',nodeno,sigmaX,sigmaY,sigmaXY);
end
fprintf('===========================================================\n');
elseif type == 5
diary('physicalCoord [Link]')
UX = input(1:2:sdof) ;
UY = input(2:2:sdof) ;
display('======MESH PHYSICAL COORDINATES====== ');
fprintf('================================================\n');
fprintf(' element node UX UY \n')
fprintf('================================================\n');
for j=1:nel
for i = 1:nnel
node = elementNodes(j,i);
ux = input(node,1);
uy = input(node,2);
fprintf(' %5d %5d %+6.5e %+6.5e \n',j,node,ux,uy);
end
end
fprintf('=================================================\n');
diary off
else
fprintf('undefined print')
end
% computeScaleFactor.m
function scaleFactor = computeScaleFactor( numberNodes,nodeCoordinates,UX,UY )
% number of total nodes
nno=numberNodes;
% number of spatial dimensions
nsd=2;
a = zeros(nsd); % Find the scale factor for deformation shape
for i=1:nsd
a(i) = max(nodeCoordinates(:,i)) - min(nodeCoordinates(:,i));
end
aa = max(a);
b = zeros(nsd);
scaleFactor = 0;
b(1) = max(UX) - min(UX);
b(2) = max(UY) - min(UY);
bb = max(b);
scaleFactor = 0.1 * aa / bb;
end