ReadElem[fileName_String]:= Module[{file,expr,expr1,expr2}, file = OpenRead[fileName]; If [ file === $Failed, Return[] ]; expr = ReadList[file,Table[Number,{9}]]; expr1 = Transpose[expr]; expr2 = Transpose[Take[expr1,-4]]; elem = expr2; ielem = Length[elem]; Close[file] ] ReadCo[fileName_String]:= Module[{file,expr4,expr5,expr6}, file = OpenRead[fileName]; If [ file=== $Failed, Return[] ]; expr4 = ReadList[file,Real,RecordLists -> True]; expr5 = Transpose [expr4]; expr6 = Transpose [ Take [expr5, {2,3} ] ]; co = expr6; xmin = Min[ Take[ Transpose[ expr6], {1} ] ]; xmax = Max[ Take[ Transpose[ expr6], {1} ] ]; ymin = Min[ Take[ Transpose[ expr6], {-1} ] ]; ymax = Max[ Take[ Transpose[ expr6], {-1} ] ]; inode = Length[co]; Close[file] ] ModelPlot2D[e_,n_,ie_,in_,nnode_]:= Module[{gpoly,xypoly,p1,p2,p3,p4list1,list2,coxy,coave}, gpoly=Table[0,{i,ie}]; xypoly=Table[0,{i,4},{j,2}]; coave=Table[0,{i,ie},{j,2}]; Do[ Do[ k=e[[j,i]];xypoly[[i,1]]=n[[k,1]]; xypoly[[i,2]]=n[[k,2]]; coave[[j,1]]+=xypoly[[i,1]]; coave[[j,2]]+=xypoly[[i,2]]; ,{i,nnode}]; gpoly[[j]]=Polygon[xypoly]; coave[[j,1]]/=4;coave[[j,2]]/=4; ,{j,ie}]; p1=Graphics[ Table[{Hue[j/ie],gpoly[[j]]},{j,ie}]]; coxy=Table[Point[{n[[i,1]],n[[i,2]]}],{i,in}]; p2=Graphics[{PointSize[0.03],coxy}]; p3=Graphics[Table[ Text[i,{n[[i,1]],n[[i,2]]},{0,1}],{i,in}]]; p4=Graphics[Table[ Text[i,{coave[[i,1]],coave[[i,2]]},{0,0}],{i,ie}]]; elemplot=Show[p1,p2,p4,disbc,loabc,PlotRange->All, AspectRatio-> Automatic , PlotLabel->"Element and Nodal Plot"]; ] ModelMesh2D[e_,n_,ie_,in_,nnode_]:= Module[ {lpoly,xypoly,p1,p2,p4,list1,list2,coxy,coave}, lpoly=Table[0,{i,nnode*ie}]; xy=Table[0,{i,4},{j,2}]; coave=Table[0,{i,ie},{j,2}]; Do[ Do[ l=e[[k,i]]; xy[[i,1]]=n[[l,1]]; xy[[i,2]]=n[[l,2]]; coave[[k,1]]+=xy[[i,1]]; coave[[k,2]]+=xy[[i,2]]; , {i,nnode}]; Do[ If[j===nnode, m=1,m=j+1]; lpoly[[(k-1)*4+j]]= Line[{{xy[[j,1]],xy[[j,2]]},{xy[[m,1]],xy[[m,2]]}}] ,{j,nnode}]; coave[[k,1]]/=4;coave[[k,2]]/=4; ,{k,ie}]; p1=Graphics[lpoly]; coxy=Table[Point[{n[[i,1]],n[[i,2]]}],{i,in}]; p2=Graphics[{PointSize[0.03],coxy}]; p3=Graphics[Table[ Text[i,{n[[i,1]],n[[i,2]]},{1,1}],{i,in}]]; p4=Graphics[Table[ Text[i,{coave[[i,1]],coave[[i,2]]},{0,0}],{i,ie}]]; elemplotmesh=Show[p1,p2,p4,disbc,loabc,PlotRange->All, AspectRatio-> Automatic, PlotLabel->"Element and Nodal Plot - Mesh"]; wiremesh={p1,p2,disbc,loabc}; ] ReadDispAnsys[fileName_String]:= Module[{file,expr7,expr8,expr9}, file = OpenRead[fileName]; If [ file === $Failed, Return [] ]; (* Skip[file,String,4]; *) expr7= ReadList[file, {Number,Word,Real,Real}]; expr8=Transpose[expr7]; expr9=Transpose[ Take[expr8, {1,3} ] ]; disp=expr9; Close[file] ] ReadLoadAnsys[fileName_String]:= Module[{file,expr10,expr11,expr12}, file = OpenRead[fileName]; If [ file === $Failed, Return [] ]; (* Skip[file,String,3]; *) expr10= ReadList[file, {Number,Word,Real,Real}]; expr11=Transpose[expr10]; expr12=Transpose[ Take[expr11, {1,3} ] ]; load=expr12; Close[file] ] ReadDisp[fileName_String]:= Module[{file,expr7}, file = OpenRead[fileName]; If [ file === $Failed, Return [] ]; expr7= ReadList[file, {Number,Number,Real}]; disp=expr7; Close[file] ] ReadLoad[fileName_String]:= Module[{file,expr8}, file = OpenRead[fileName]; If [ file === $Failed, Return [] ]; expr8= ReadList[file, {Number,Number,Real}]; load=expr8; Close[file] ] FormGlobalStifMat[]:= Module[{s3,tj,tji,tm,dm,c,t,tc,ct,cc,tmi}, (* Declare zero matrices *) xy=Table[0,{i,nnode},{j,ndofn}];c=Table[0,{i,3},{j,ndof}]; ct=Table[0,{i,3},{j,4}]; totdof=inode*ndofn; s=Table[0,{i,totdof},{j,totdof}]; fileout1=OpenWrite["stress.m",FormatType -> OutputForm]; fileout2=OpenWrite["stress1.m",FormatType -> OutputForm]; fileout3=OpenWrite["stressre.m"]; (* Do Loop for number of elem -- till ielem *) Do[ em=Table[0,{i,8},{j,8}];kee=Table[0,{i,4},{j,4}]; ker=Table[0,{i,4},{j,8}];kre=Table[0,{i,8},{j,4}]; (* Do Loop for numerical integration *) Do[ s3:=0.577350269; If[ii===1,x=-s3;y=-s3]; If[ii===2,x=s3;y=-s3]; If[ii===3,x=s3;y=s3]; If[ii===4,x=-s3;y=s3]; (* Assigning of coordinates to appropriate nodes within an element *) Do[k=elem[[n,i]];xy[[i,1]]=co[[k,1]]; xy[[i,2]]=co[[k,2]]; ,{i,4}]; tj=dkl.xy; tm=dkm.xy; dj=tj[[1,1]]*tj[[2,2]]-tj[[1,2]]*tj[[2,1]]; dm=tm[[1,1]]*tm[[2,2]]-tm[[1,2]]*tm[[2,1]]; tji=tj[[1,1]]; tj[[1,1]]=tj[[2,2]]/dj; tj[[2,2]]=tji/dj; tj[[1,2]]=-tj[[1,2]]/dj; tj[[2,1]]=-tj[[2,1]]/dj; tjm=tm[[1,1]]; tm[[1,1]]=tm[[2,2]]/dm; tm[[2,2]]=tjm/dm; tm[[1,2]]=-tm[[1,2]]/dm; tm[[2,1]]=-tm[[2,1]]/dm; t=tj.dkl; ttm=tm.dkj *dm/dj; tjj=sha.xy; dj=Abs[dj]; (* Do Loop to form strain matrix *) Do[ c[[1,2*j-1]] = t[[1,j]]; c[[3,2*j]] = t[[1,j]]; c[[2,2*j]] = t[[2,j]]; c[[3,2*j-1]] = t[[2,j]]; ,{j,4}]; Do[ ct[[1,2*j-1]] = ttm[[1,j]]; ct[[3,2*j]] = ttm[[1,j]]; ct[[2,2*j]] = ttm[[2,j]]; ct[[3,2*j-1]] = ttm[[2,j]]; ,{j,2}]; (*Put[dpss,"trace.m"];PutAppend[c,"trace.m"]; PutAppend[ct,"trace.m"];PutAppend[ii,"trace.m"]; PutAppend[n,"trace.m"];*) t=dpss.c; tc=dpss.ct; Do[ Write[fileout1,FortranForm[t[[i,1]]]," ",FortranForm[t[[i,2]]]," ", FortranForm[t[[i,3]]]," ",FortranForm[t[[i,4]]]," ", FortranForm[t[[i,5]]]," ",FortranForm[t[[i,6]]]," ", FortranForm[t[[i,7]]]," ",FortranForm[t[[i,8]]] ], {i,3} ]; (* Write x and y coordinates of gauss point *) Write[fileout1,FortranForm[tjj[[1,1]]]," ",FortranForm[tjj[[1,2]]] ]; (*Print[ "x and y coordinates of gauss point are ", tjj[[1,1]], " " ,tjj[[1,2]] ];*) Do[ Write[fileout2,FortranForm[tc[[i,1]]]," ",FortranForm[tc[[i,2]]]," ", FortranForm[tc[[i,3]]]," ",FortranForm[tc[[i,4]]] ], {i,3} ]; em=em + Transpose[c].t *dj*thick; kee=kee + Transpose[ct].tc *dj*thick; kre=kre + Transpose[c].tc *dj*thick; ,{ii,4}]; (* closing loop for numerical integration *) ker=Transpose[kre]; cc=N[Inverse[kee].ker]; Do[ Do[ Write[fileout3,FortranForm[cc[[i,j]]] ], {j,8}], {i,4} ]; em=em-kre.cc; (*If[n===1, Put[em,"em.m"], PutAppend[em,"em.m"] ];*) (* Assembly of Element Stiffness Matrix to Global Stiffness Matrix *) Do[ k=elem[[n,i]]; ii1=2*i-1; ii2=2*i; ki1=2*k-1; ki2=2*k; Do[ l=elem[[n,j]]; jj1=2*j-1; jj2=2*j; lj1=2*l-1; lj2=2*l; s[[ki1,lj1]] = s[[ki1,lj1]] + em[[ii1,jj1]]; s[[ki1,lj2]] = s[[ki1,lj2]] + em[[ii1,jj2]]; s[[ki2,lj1]] = s[[ki2,lj1]] + em[[ii2,jj1]]; s[[ki2,lj2]] = s[[ki2,lj2]] + em[[ii2,jj2]]; ,{j,4}] ,{i,4}]; ,{n,ielem}]; (* closing loop for total number of element *) Close[fileout1]; Close[fileout2]; Close[fileout3]; ] SolutionLinear[]:= Module[ {cs,n1,k1,pro }, (* Solution of simultaneous equations by Gauss Elimination *) n1=totdof-1; Do[ cs=s[[k,k]]; k1=k+1; (* Try to interchange rows to get non zero diagonal coefficient *) If [ Abs[cs] < 0.000000001, Do[ If[ Abs[s[[j,k]]] > 0.00000001, Do[ cs=s[[k,l]]; s[[k,l]]=s[[j,l]]; s[[j,l]]=cs; ,{l,k,totdof} ]; cs=paload[[k,1]]; paload[[k,1]]=paload[[j,1]]; paload[[j,1]]=cs; cs=s[[k,k]]; Break[]; ]; Print[{"A loop on J ",j," Row Number ",k}]; ,{j,k1,totdof}]; If[j===totdof,Print[{"****** Singularity in Row ", k}]; Print[{" No solution "}] ]; ]; (* Divide row by diagonal coefficient *) Do[ s[[k,j]]=s[[k,j]]/cs; , {j,k1,totdof} ]; paload[[k,1]]=paload[[k,1]]/cs; (* Elimination of left-most column from each row *) Do[ cs=s[[i,k]]; Do[ s[[i,j]]=s[[i,j]]-cs*s[[k,j]]; ,{j,k1,totdof}]; paload[[i,1]]=paload[[i,1]]-cs*paload[[k,1]]; , {i,k1,totdof}]; (*pro=N[k*100/n1];*) ,{k,n1}]; If[ Abs[s[[totdof,totdof]]] > 0.00000001, (* Compute last unknown *) dadisp[[totdof,1]]=paload[[totdof,1]]/s[[totdof,totdof]]; (* Back substitution process to compute remaining unknowns *) Do[ k=totdof-l; k1=k+1; Do[ dadisp[[k,1]]=dadisp[[k,1]] + paload[[k,1]] - s[[k,j]]*dadisp[[j,1]]; paload[[k,1]]=0.; ,{j,k1,totdof}]; (*pro=N[l*100/n1];*) (*Print[{" -- Back Sub. Process ", pro, " % finished --"}];*) ,{l,n1}]; Put[dadisp,"dadisp.m"] ,Print[{"****** Singularity in Row ",totdof}]; Print[{" last unknown equal to ",s[[totdof,totdof]] }] ]; ] ApplyDispAnsys[]:= Module[{huge,in,jm,j1,length,r,list,numdis,a,b}, huge=10000000. ; length=Length[disp]; list=Take[Transpose[co],1]; r=(Max[list]-Min[list])/30.; numdis=Table[0,{i,length}]; Do[ x1=co[[disp[[i,1]],1]]; y1=co[[disp[[i,1]],2]]; If[ disp[[i,2]]==="UX", in=1; jm=2*disp[[i,1]]-in; s[[jm,jm]]=s[[jm,jm]]*huge; paload[[jm,1]]=s[[jm,jm]]*disp[[i,3]]; numdis[[i]]=Polygon[{{x1,y1},{x1-2 r,y1-r},{x1-2 r,y1+r}}]; ]; If[ disp[[i,2]]==="UY", jm=2*disp[[i,1]]; s[[jm,jm]]=s[[jm,jm]]*huge; paload[[jm,1]]=s[[jm,jm]]*disp[[i,3]]; numdis[[i]]=Polygon[{{x1,y1},{x1-r,y1-2 r},{x1+r,y1-2 r}}]; ]; ,{i,length}]; disbc=Graphics[Table[{GrayLevel[.3],numdis[[i]]}, {i,length}]]; ] ApplyLoadAnsys[]:= Module[{length,r,numloa,rat,list1,list2,in}, length=Length[load]; list1=Take[Transpose[co],1]; list2=Take[Transpose[load],-1]; rat=(Max[list1]-Min[list1])/(15.*Max[list2]); numloa=Table[0,{i,length}]; Do[ r=Abs[rat*load[[i,3]]]; x1=co[[load[[i,1]],1]]; y1=co[[load[[i,1]],2]]; If[ load[[i,2]]==="FX", If[ load[[i,3]] > 0., numloa[[i]]={Hue[.6],Polygon[{ {x1,y1},{x1-r,y1-(r/2)},{x1-r,y1-(r/4)},{x1-3 r,y1-(r/4)}, {x1-3 r,y1+(r/4)},{x1-r,y1+(r/4)},{x1-r,y1+(r/2)}}]} ]; If[ load[[i,3]] < 0., numloa[[i]]={Hue[.1],Polygon[{ {x1,y1}, {r+x1,-r/2+y1},{r+x1,-r/4+y1}, {3 r+x1,-r/4+y1}, {3 r+x1,r/4+y1},{r+x1,r/4+y1},{r+x1,r/2+y1}}]} ]; in=1 ]; If[ load[[i,2]]==="FY", If[ load[[i,3]] > 0., numloa[[i]]={Hue[.6],Polygon[{ {x1,y1},{x1-(r/2),y1-r},{x1-(r/4),y1-r},{x1-(r/4),y1-3 r}, {x1+(r/4),y1-3 r},{x1+(r/4),y1-r},{x1+(r/2),y1-r}}]} ]; If[ load[[i,3]] < 0., numloa[[i]]={Hue[.1],Polygon[{ {x1,y1},{x1-(r/2),y1+r},{x1-(r/4),y1+r},{x1-(r/4),y1+3 r}, {x1+(r/4),y1+3 r},{x1+(r/4),y1+r},{x1+(r/2),y1+r}}]} ]; in=0 ]; jm=2*load[[i,1]]-in ; paload[[jm,1]]=load[[i,3]]; ,{i,length}]; loabc=Graphics[Table[numloa[[i]],{i,length}]]; ] ApplyDisp[]:= Module[{huge,len }, huge=10000000000. ; len=Length[disp]; Do[ If[ disp[[i,2]]===1, inum=1; jnum=2*disp[[i,1]]-inum; s[[jnum,jnum]]=s[[jnum,jnum]]*huge ; paload[[jnum,1]]=s[[jnum,jnum]]*disp[[i,3]] ]; If[ disp[[i,2]]===2, inum=1; jnum=2*disp[[i,1]]; s[[jnum,jnum]]=s[[jnum,jnum]]*huge; paload[[jnum,1]]=s[[jnum,jnum]]*disp[[i,3]] ]; If[ disp[[i,2]]===12, Do[ inum=j-1; jnum=2*disp[[i,1]]-inum; s[[jnum,jnum]]=s[[jnum,jnum]]*huge; paload[[jnum,1]]=s[[jnum,jnum]]*disp[[i,3]]; ,{j,2}] ]; ,{i,len}]; ] ApplyLoad[]:= Module[{len}, len=Length[load]; Do[ If[ load[[i,2]]===1, inum=1 ]; If[ load[[i,2]]===2, inum=0 ]; jnum=2*load[[i,1]]-inum; paload[[jnum,1]]=load[[i,3]]; ,{i,len}]; ] StressCompu[]:= Module[ {filein1,filein2,filein3,fileout4,ed,est,cc,tc,tt,t}, filein1=OpenRead["stress.m"]; filein2=OpenRead["stress1.m"]; filein3=OpenRead["stressre.m"]; fileout4=OpenWrite["results.m",FormatType -> OutputForm]; ed=Table[0,{i,8},{j,1}];cc=Table[0,{i,4},{j,8}]; num=4*ielem; sx=Table[0,{i,num}];sy=Table[0,{i,num}];sxy=Table[0,{i,num}]; xco=Table[0,{i,num}];yco=Table[0,{i,num}]; gaussx=Table[0,{i,4}];gaussy=Table[0,{i,4}]; Do[ (* Do loop for number of elements *) est=Table[0,{i,4},{j,3}]; Write[fileout4," Stress values for element no : ",n]; Do[ k=elem[[n,i]]; ed[[2*i-1,1]]=dadisp[[2*k-1,1]]; ed[[2*i,1]]=dadisp[[2*k,1]]; If[k===8,index=i]; ,{i,4} ]; (* cc=Read[filein3,{{Real,Real,Real,Real,Real,Real, Real,Real},{Real,Real,Real,Real,Real,Real,Real,Real}, {Real,Real,Real,Real,Real,Real,Real,Real}, {Real,Real,Real,Real,Real,Real,Real,Real}}];*) Do[Do[ cc[[i,j]]=Read[filein3,Real] ,{j,8}] ,{i,4}]; Do[ (* Do loop for 4 stress values at gauss points *) tt=Read[filein1,{{Real,Real,Real,Real,Real,Real, Real,Real},{Real,Real,Real,Real,Real,Real,Real,Real}, {Real,Real,Real,Real,Real,Real,Real,Real}}]; gaussx[[ii]]=Read[filein1,Real]; gaussy[[ii]]=Read[filein1,Real]; tc=Read[filein2,{{Real,Real,Real,Real}, {Real,Real,Real,Real},{Real,Real,Real,Real}}]; t=tt-tc.cc; estress=t.ed; Do[est[[ii,jj]]=estress[[jj,1]],{jj,3}] ,{ii,4}]; Do[ sx[[n*4-4+ii]]=est[[ii,1]]; sy[[n*4-4+ii]]=est[[ii,2]]; sxy[[n*4-4+ii]]=est[[ii,3]]; xco[[n*4-4+ii]]=gaussx[[ii]]; yco[[n*4-4+ii]]=gaussy[[ii]]; , {ii,4}]; extramat={{1.866,-0.5,0.133975,-0.5},{-0.5,1.866,-0.5,0.133975}, {0.133975,-0.5,1.866,-0.5},{-0.5,0.133975,-0.5,1.866}}; estext=extramat.est; If[n===1 , calstress=estext[[index,1]] ]; Do[ Write[fileout4, ii1, " ", estext[[ii1,1]]," ",estext[[ii1,2]], " ", estext[[ii1,3]] ]; ,{ii1,4} ]; (*Do[ Write[fileout4,cc[[i,1]]," ",cc[[i,2]]," ", cc[[i,3]]," ",cc[[i,4]]," ", cc[[i,5]]," ",cc[[i,6]]," ",cc[[i,7]]," ",cc[[i,8]] ], {i,4} ];*) ,{n,ielem}]; Close[filein1]; Close[filein2]; Close[filein3]; Close[fileout4]; ] DispPlot2D[e_,n_,ie_,in_,nnode_]:= Module[ {gpoly,xypoly,p1,p2,p3,list1,list2,coxy}, gpoly=Table[0,{i,ie}]; xypoly=Table[0,{i,4},{j,2}]; Do[ Do[ k=e[[j,i]];xypoly[[i,1]]=n[[k,1]]; xypoly[[i,2]]=n[[k,2]]; ,{i,nnode}]; gpoly[[j]]=Polygon[xypoly]; ,{j,ie}]; p1=Graphics[ Table[{Hue[j/ie],gpoly[[j]]},{j,ie}]]; coxy=Table[Point[{n[[i,1]],n[[i,2]]}],{i,in}]; p2=Graphics[{PointSize[0.03],coxy}]; p3=Graphics[Table[ Text[i,{n[[i,1]],n[[i,2]]},{0,1}],{i,in}]]; pldisp=Show[p1,p2,PlotRange->All, AspectRatio-> Automatic, PlotLabel->"Displacement Plot"]; ] StressFunction[ ]:= Module[ {num,list,u,v}, num=nnode*ielem; list=Table[ {xco[[i]],yco[[i]],sx[[i]]} ,{i,num}]; eqnsx=Fit[list,{1,x,y,x y},{x,y}]; xpoly[x_,y_]=Fit[list,{1,x,y,x y},{x,y}]; Print["polynomial of stress x ",eqnsx]; list=Table[ {xco[[i]],yco[[i]],sy[[i]]} ,{i,num}]; eqnsy=Fit[list,{1},{x,y}]; ypoly[x_,y_]=Fit[list,{1,x,y,x y},{x,y}]; Print["polynomial of stress y ",eqnsy]; list=Table[ {xco[[i]],yco[[i]],sxy[[i]]} ,{i,num}]; eqnsxy=Fit[list,{1},{x,y}]; xypoly[x_,y_]=Fit[list,{1,x,y,x y},{x,y}]; Print["polynomial of stress xy ",eqnsxy]; ] Plnssx[ ]:= Module[ { }, ContourPlot[eqnsx,{x,xmin,xmax},{y,ymin,ymax}, ColorFunction->Hue,AspectRatio->Automatic, Contours->10] ] Plnssy[ ]:= Module[ { }, ContourPlot[eqnsy,{x,xmin,xmax},{y,ymin,ymax}, ColorFunction->Hue,AspectRatio->Automatic, Contours->10] ] Plnssxy[ ]:= Module[ { }, ContourPlot[eqnsxy,{x,xmin,xmax},{y,ymin,ymax}, ColorFunction->Hue,AspectRatio->Automatic, Contours->10] ] StressAtLocation[a_,b_]:= Module[ {}, Print[" Stress at coordinates x = ",a," and y = ",b]; Print[" "]; Print[" SX SY SXY "]; Print[" "]; Print[xpoly[a,b]," ",ypoly[a,b]," ",xypoly[a,b] ]; ]