The Pre-Processing involve the building up of finite element model by nodal coordinates, user have to prepare nodal and element listing in order to be read in by program.

Table of nodal and elemental listing will be produced for user to counter check.  Nodal plot in colour will be displayed as well, shape functions, element properties could be read in for pre-processing.

<<readqm6.m;

ReadCo["node"];
OutputForm[TableForm[co,TableHeadings->{Automatic,
{"X - Coord","Y - Coord"}}]] >> resfile.m;

ReadElem["elem"];
OutputForm[TableForm[elem,TableHeadings->{Automatic,
{"Node - 1","Node -2","Node - 3","Node - 4"}}]] >>>resfile.m

(*  Reading in Shape Functions  -- dkl  *)
<<shqm6.m;

(*  Reading in Elasticity Matrix  --  dpss *)
e=1500;v=0.3; (* e is the young's modulus and v poisson ratio*)
<<dpss.m;

(* Declare thickness for Plane Stress *)
thick=1.  ; (* default thickness  *)

"Elasticity Matrix">>>resfile.m
OutputForm[MatrixForm[dpss]]>>>resfile.m

FormGlobalStifMat[]

(*  Reading in Displacement and Mechanical  
    Boundary conditions   *)

ReadDispAnsys["disp"];  ReadLoadAnsys["load"];
paload=Table[0,{i,totdof},{j,1}];
dadisp=Table[0,{i,totdof},{j,1}];

ApplyDispAnsys[];  ApplyLoadAnsys[];

(* Routine to do element plot  *)
ModelPlot2D[elem, co, ielem, inode, nnode];

ModelMesh2D[elem, co, ielem, inode, nnode];

In this phase, the assembled stiffness matrix with the known boundary conditions are called into user programmed routine to solve for the displacement unknowns.  Gauss Elimination Method is used for the solution of simultaneous equations. Pivoting and row interchange procedures are used to treat ill-conditioned matrix.

SolutionLinear[];

Post-Processing involved the print out of displacement and stress values of the analysis.  Displacement plots are produced to give a better picture of final shape with magnification factor adjustable.

StressCompu[];x=.;y=.;

StressFunction[];

newco=Table[0,{i,inode},{j,2}];

(*  Set displacement manification *)
mag=1./100.;

Do[newco[[i,1]]=co[[i,1]]+dadisp[[i*2-1,1]]*mag;
   newco[[i,2]]=co[[i,2]]+dadisp[[i*2,1]]*mag;
,{i,inode}];


OutputForm[TableForm[newco-co,TableHeadings->
{Automatic,{"X - Disp",
 "Y - Disp"}}]]  >>> resfile.m

DispPlot2D[elem,newco,ielem,inode,nnode];

Show[pldisp,wiremesh,AspectRatio->Automatic];

Plnssx[]

Show[%,wiremesh,Frame->False]

Plnssy[]

Show[%,wiremesh,Frame->False]

Plnssxy[]

Show[%,wiremesh,Frame->False]

StressAtLocation[0.,2.]

From the theory of beam bending, the stress at node 8 could be calculated through the following formulae :

Moment at node 8  = M = force x distance
 	                     = 300 x 10 Nm
	                     = 3000 Nm

Stress = ( M * y) / I          where   y = 1 m    and    I = (b d^3) / 12  =  0.6666 m^4

Theoretical Stress = ( 3000 * 1 ) / 0.6666  N/m^2   =  -4500  N/m^2  (Compressive)

theorystress=-4500;polystress=xpoly[0.,2.];
errorstress1=(-4500 - calstress )*100/(-4500);
errorstress2=(-4500 - polystress )*100/(-4500);

TableForm[{{theorystress,calstress,errorstress1,
polystress,errorstress2}},
TableHeadings -> {{"Node 8"}, {"T. Stress",
"C. Stress","% Error","P. Stress","% Error"}}]

Note :

T. Stress --  Theoretical Stress
C. Stress --  Calculated Stress
P. Stress --  Smoothed Polynomial Stress using Mathematica function Fit