(************** Content-type: application/mathematica **************
                     CreatedBy='Mathematica 4.2'

                    Mathematica-Compatible Notebook

This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do
one of the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the
  application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
*******************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[    118441,       2373]*)
(*NotebookOutlinePosition[    119540,       2408]*)
(*  CellTagsIndexPosition[    119496,       2404]*)
(*WindowFrame->Normal*)



Notebook[{
Cell["\<\
To use this program, here are the instructions:
1.  Execute Cells 1 through 10, in that order (twice each to  get rid of  
     annoying warning messages) to initialize all necessary modules. 
     Shortcut for Mathematica savvy users: select Kernel \
->Evaluation->Evaluate
     Initialization  (In Mathematica 2.2,Action->Evaluate Initialization).
2.  Cells 11 through 13 contain driver programs for doing test problems 
     Run them and check that the answers (displacements and stresses) are \
correct.
3.  It is recommended you place the main programs for the problems in \
following
     Cells.   Those in Cells 11 -- 13 may be used as \"templates\".
4.  If the plots produced by  running a driver program appear too small,
     click on the top one  (only) with the mouse, grab a corner \"handle\" 
     and enlarge it.  Then rerun the program.  \
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 0.727611, 0.140017]],

Cell[TextData[{
  StyleBox["  ",
    FontFamily->"Palatino"],
  "Cell 0.  Mathematica version-dependent settings (to get plots displayed \
correctly)."
}], "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
DisplayChannel=$DisplayFunction;
If [$VersionNumber>=6.0, DisplayChannel=Print]; 
  (* fix for Mathematica 6 & later *)  \
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  " ",
  StyleBox["Off[General::spell1]; Off[General::spell];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell  1.  Modules QuadGaussRuleInfo, TrigGaussRuleInfo and \
LineGaussRuleInfo provide 
quadrature information for  iso-P quadrilaterals, triangles, and lines, \
respectively.
 
For quadrilaterals the rules implemented are 1x1 through 5x5.  These are \
identified by 
the number of points in each direction, p=1, ... 5. See IFEM Notes Ch 23 for \
argument rules. 
For triangles the rules are p=1,3,-3,6,-6,7;  see IFEM Ch 24.

This code is used by all iso-P stiffness and stress modules. Hence it must be \
initialized 
before their appearance.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
QuadGaussRuleInfo[{rule_,numer_},point_]:= Module[
 {\[Xi],\[Eta],p1,p2,i,j,w1,w2,m,info={{Null,Null},0}},
  If [Length[rule]==2,  {p1,p2}=rule, p1=p2=rule];
  If [p1<0, Return[QuadNonProductGaussRuleInfo[
      {-p1,numer},point]]];
  If [Length[point]==2, {i,j}=point, m=point; 
      j=Floor[(m-1)/p1]+1; i=m-p1*(j-1) ];
  {\[Xi],w1}=  LineGaussRuleInfo[{p1,numer},i];
  {\[Eta],w2}=  LineGaussRuleInfo[{p2,numer},j];
  info={{\[Xi],\[Eta]},w1*w2};
  If [numer, Return[N[info]], Return[Simplify[info]]];
]; 

TrigGaussRuleInfo[{rule_,numer_},point_]:= Module[
 {zeta,p=rule,i=point,g1,g2,info={{Null,Null,Null},0}  },
  If [p== 1, info={{1/3,1/3,1/3},1}];
  If [p== 3, info={{1,1,1}/6,1/3}; info[[1,i]]=2/3];
  If [p==-3, info={{1,1,1}/2,1/3}; info[[1,i]]=0  ];
  If [p== 6,  g1=(8-Sqrt[10]+Sqrt[38-44*Sqrt[2/5]])/18;
              g2=(8-Sqrt[10]-Sqrt[38-44*Sqrt[2/5]])/18;
     If [i<4, info={{g1,g1,g1},(620+Sqrt[213125-
              53320*Sqrt[10]])/3720}; info[[1,i]]=1-2*g1];
     If [i>3, info={{g2,g2,g2},(620-Sqrt[213125-
              53320*Sqrt[10]])/3720}; info[[1,i-3]]=1-2*g2]];
  If [p== -6, 
     If [i<4, info={{1,1,1}/6,3/10}; info[[1,i]]=2/3];  
     If [i>3, info={{1,1,1}/2,1/30}; info[[1,i-3]]=0]]; 
  If [p== 7,  g1=(6-Sqrt[15])/21; g2=(6+Sqrt[15])/21;
     If [i<4,      info={{g1,g1,g1},(155-Sqrt[15])/1200};
                   info[[1,i]]=  1-2*g1];
     If [i>3&&i<7, info={{g2,g2,g2},(155+Sqrt[15])/1200};
                   info[[1,i-3]]=1-2*g2];
     If [i==7,     info={{1/3,1/3,1/3},9/40} ]];    
  If [numer, Return[N[info]], Return[Simplify[info]]];
];
 
LineGaussRuleInfo[{rule_,numer_},point_]:= Module[
  {g2={-1,1}/Sqrt[3],w3={5/9,8/9,5/9}, 
   g3={-Sqrt[3/5],0,Sqrt[3/5]}, 
   w4={(1/2)-Sqrt[5/6]/6, (1/2)+Sqrt[5/6]/6,
       (1/2)+Sqrt[5/6]/6, (1/2)-Sqrt[5/6]/6},
   g4={-Sqrt[(3+2*Sqrt[6/5])/7],-Sqrt[(3-2*Sqrt[6/5])/7],
        Sqrt[(3-2*Sqrt[6/5])/7], Sqrt[(3+2*Sqrt[6/5])/7]},
   g5={-Sqrt[5+2*Sqrt[10/7]],-Sqrt[5-2*Sqrt[10/7]],0, 
        Sqrt[5-2*Sqrt[10/7]], Sqrt[5+2*Sqrt[10/7]]}/3,
   w5={322-13*Sqrt[70],322+13*Sqrt[70],512,
       322+13*Sqrt[70],322-13*Sqrt[70]}/900,
   i=point,p=rule,info={Null,0}}, 
  If [p==1, info={0,2}];
  If [p==2, info={g2[[i]],1}];
  If [p==3, info={g3[[i]],w3[[i]]}]; 
  If [p==4, info={g4[[i]],w4[[i]]}];
  If [p==5, info={g5[[i]],w5[[i]]}];
  If [numer, Return[N[info]], Return[Simplify[info]]];
]; \
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2A:  Modules for 4-node iso-P bilinear quadrilateral ring \
element  
for axisymmetric solids. These are documented in AFEM Chapter 12.\
\>", "Text",\

  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad4IsoPRingStiffness[ncoor_,Emat_,options_]:=Module[\n  \
{p=2,numer=False,Jcons=False,Kfac=1,qcoor,k,\n   \
r1,r2,r3,r4,z1,z2,z3,z4,Nf,N1,N2,N3,N4,A0,Jdet,Be,\n   \
dNr1,dNr2,dNr3,dNr4,dNz1,dNz2,dNz3,dNz4,rk,w,c,Ke,\n  \
Ke0=Table[0,{8},{8}],modname=\"Quad4IsoPRingStiffness: \"},\n  If \
[Length[options]==1,{numer}=options];\n  If \
[Length[options]==2,{numer,p}=options];\n  If \
[Length[options]==3,{numer,p,Jcons}=options];\n  If \
[Length[options]==4,{numer,p,Jcons,Kfac}=options];\n  If [p<1||p>5, \
Print[modname,\"illegal p:\",p]; Return[Ke0]];\n  \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;\n  \
A0=((r3-r1)*(z4-z2)-(r4-r2)*(z3-z1))/2;\n  If [numer&&(A0<=0), Print[modname,\
\"Neg or zero area\"]; \n      Return[Ke0]];  Ke=Ke0;\n  For [k=1,k<=p*p,k++, \
\n      {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];\n     \
{Nf,{dNr1,dNr2,dNr3,dNr4},{dNz1,dNz2,dNz3,dNz4},\n      \
Jdet}=Quad4IsoPRingShapeFunDer[ncoor,qcoor,Jcons];\n      If \
[numer&&(Jdet<=0), Print[modname,\"Neg or zero\",\n         \" Gauss point \
Jacobian at k=\",k]; Return[Ke0]];\n      {N1,N2,N3,N4}=Nf; \
rk=r1*N1+r2*N2+r3*N3+r4*N4; \n      Be={{ dNr1,   0, dNr2,   0, dNr3,   0, \
dNr4,   0},\n          {    0,dNz1,    0,dNz2,    0,dNz3,    0,dNz4},\n       \
   {N1/rk,   0,N2/rk,   0,N3/rk,   0,N4/rk,   0},\n          { dNz1,dNr1, \
dNz2,dNr2, dNz3,dNr3, dNz4,dNr4}};\n      c=Kfac*w*rk*Jdet; If \
[numer,Be=N[Be]; c=N[c]]; \n      Ke+=c*Transpose[Be].(Emat.Be);   \n      ]; \
ClearAll[Ke0,Be]; Return[Ke] ];\n\n\
Quad4IsoPRingShapeFunDer[ncoor_,qcoor_,Jcons_]:= Module[\n  { \
r1,r2,r3,r4,z1,z2,z3,z4,\[Xi],\[Eta],Nf,dNr,dNz,A0,A1,A2,Jdet},\n  \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor; {\[Xi],\[Eta]}=qcoor;\n   Nf={(1-\
\[Xi])*(1-\[Eta]),(1+\[Xi])*(1-\[Eta]),(1+\[Xi])*(1+\[Eta]),(1-\[Xi])*(1+\
\[Eta])}/4; \n   A0=((r3-r1)*(z4-z2)-(r4-r2)*(z3-z1))/2;\n   \
A1=((r3-r4)*(z1-z2)-(r1-r2)*(z3-z4))/2;\n   \
A2=((r2-r3)*(z1-z4)-(r1-r4)*(z2-z3))/2;\n   Jdet=(A0+A1*\[Xi]+A2*\[Eta])/4; \
If [Jcons,Jdet=A0/4];\n   \
dNr={z2-z4+(z4-z3)*\[Xi]+(z3-z2)*\[Eta],z3-z1+(z3-z4)*\[Xi]+(z1-z4)*\[Eta],\n \
   z4-z2+(z1-z2)*\[Xi]+(z4-z1)*\[Eta],z1-z3+(z2-z1)*\[Xi]+(z2-z3)*\[Eta]}/(8*\
Jdet);\n   \
dNz={r4-r2+(r3-r4)*\[Xi]+(r2-r3)*\[Eta],r1-r3+(r4-r3)*\[Xi]+(r4-r1)*\[Eta],\n \
   r2-r4+(r2-r1)*\[Xi]+(r1-r4)*\[Eta],r3-r1+(r1-r2)*\[Xi]+(r3-r2)*\[Eta]}/(8*\
Jdet);\n   Return[{Nf,dNr,dNz,Jdet}]\n];\n\n\
Quad4IsoPRingBodyForces[ncoor_,bfor_,options_]:=Module[\n  \
{p=2,numer=False,Jcons=False,Kfac=1,qcoor,k,m,\n  \
r1,r2,r3,r4,z1,z2,z3,z4,N1,N2,N3,N4,dNr,dNz,Jdet,Be,\n  \
br1,bz1,br2,bz2,br3,bz3,br4,bz4,brc,bzc,bk,rk,w,c,fe,\n  \
fe0=Table[0,{8}],modname=\"Quad4IsoPRingBodyForces: \"},\n  If \
[Length[options]==1,{numer}=options];\n  If \
[Length[options]==2,{numer,p}=options];\n  If \
[Length[options]==3,{numer,p,Jcons}=options];\n  If \
[Length[options]==4,{numer,p,Jcons,Kfac}=options];\n  If [p<1||p>5, \
Print[modname,\"p out of range\"]; Return[fe0]];\n  \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;\n  \
A0=((r3-r1)*(z4-z2)-(r4-r2)*(z3-z1))/2;\n  If [numer&&(A0<=0), Print[modname,\
\"Neg or zero area\"]; \n      Return[fe0]]; fe=fe0; m=Length[bfor]; \n  If \
[m!=2&&m!=4, Print[modname,\" Illegal bfor\"]; Return[fe0]]; \n  If [m==2, \
br1=br2=br3=br4=bfor[[1]];bz1=bz2=bz3=bz4=bfor[[2]]];\n  If \
[m==4,{{br1,bz1},{br2,bz2},{br3,bz3},{br4,bz4}}=bfor];\n  For \
[k=1,k<=p*p,k++, \n      {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n     {{N1,N2,N3,N4},dNr,dNz,Jdet}=\n           \
Quad4IsoPRingShapeFunDer[ncoor,qcoor,Jcons];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      If [numer&&(Jdet<=0), Print[modname,\"Neg or zero\",\n         \" \
Gauss point Jacobian at k=\",k]; Return[fe0]];\n      \
rk=r1*N1+r2*N2+r3*N3+r4*N4;  c=Kfac*w*Jdet*rk; ",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      brk=br1*N1+br2*N2+br3*N3+br4*N4;\n      \
bzk=bz1*N1+bz2*N2+bz3*N3+bz4*N4;\n      bk={N1*brk,N1*bzk,N2*brk,N2*bzk,\n    \
      N3*brk,N3*bzk,N4*brk,N4*bzk};\n      If [numer,bk=N[bk]]; fe+=c*bk;   \n\
      ];  If[!numer, fe=Simplify[fe]];\n      Return[fe] ];\n",
  StyleBox[" ",
    FontColor->RGBColor[0, 0, 1]],
  "     \nQuad4IsoPRingStresses[ncoor_,Emat_,ue_,options_]:= \n  \
Module[{numer=False,g=1/Sqrt[3],Jcons=False,w0=0,\n     \
eps=10.^(-9),r1,r2,r3,r4,z1,z2,z3,z4,Nf,N1,N2,N3,N4,\n     \
dNr1,dNr2,dNr3,dNr4,dNz1,dNz2,dNz3,dNz4,\n     \
T1,T2,T3,T4,Td,Tg4,Jdet,qcoor,w,c,Be,\n     \
gctab={{0,0}},k,kg,rk,sigg,sige,udis=ue,\n     \
modname=\"Quad4IsoPRingStresses: \"},\n  If \
[Length[options]==1,{numer}=options];\n  If \
[Length[options]==2,{numer,g}=options];\n  If \
[Length[options]==3,{numer,g,w0}=options];\n  If [Head[g]==Symbol||g>0, \
Td=4*g^2*(4+w0); \n     T1=4*g^2*w0; T2=4+4*g^2+w0+2*g*(4+w0); \n     \
T3=-4+4*g^2-w0; T4=4+4*g^2+w0-2*g*(4+w0);\n     \
Tg4={{T1,T2,T3,T4,T3},{T1,T3,T2,T3,T4},\n          \
{T1,T4,T3,T2,T3},{T1,T3,T4,T3,T2}}/Td;",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "      \n     gctab={{0,0},{-1,-1},{1,-1},{1,1},{-1,1}}*g];\n  \
kg=Length[gctab]; sigg=Table[{0,0,0,0},{kg}];\n  If [numer, gctab=N[gctab]; \
Tg4=N[Tg4]; udis=N[ue] ]; \n  {{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;\n  For \
[k=1,k<=kg,k++, qcoor=gctab[[k]]; \n      \
{Nf,{dNr1,dNr2,dNr3,dNr4},{dNz1,dNz2,dNz3,dNz4},\n       \
Jdet}=Quad4IsoPRingShapeFunDer[ncoor,qcoor,Jcons];\n      {N1,N2,N3,N4}=Nf; \
rk=r1*N1+r2*N2+r3*N3+r4*N4;  \n      Be={{ dNr1,   0, dNr2,   0, dNr3,   0, \
dNr4,   0},\n          {    0,dNz1,    0, dNz2,   0,dNz3,    0,dNz4},\n       \
   {N1/rk,   0,N2/rk,   0,N3/rk,   0,N4/rk,   0},\n          { dNz1,dNr1, \
dNz2,dNr2, dNz3,dNr3, dNz4,dNr4}};\n      If [numer,Be=N[Be]]; \
sigg[[k]]=Emat.(Be.udis);\n      ] ; \n  If [kg==1, \
sige=Table[sigg[[1]],{4}], sige=Tg4.sigg];\n  If [numer, \
sige=Chop[sige,eps]]; \n  If [!numer,sige=Simplify[sige]]; Return[sige] ]; \n \
 \n",
  StyleBox["(* ClearAll[Em,\[Nu],a,b,d,h,p,numer]; \nEm=96; \[Nu]=1/3; a=4; \
b=2; d=1; e=3; Kfac=1;\nncoor={{d,0},{a+d,0},{a+d,b+e},{d,b}}; numer=True; \
Jcons=True;\n",
    FontColor->RGBColor[0, 0, 1]],
  "\[Sigma]rr",
  StyleBox["\nPrint[\"Emat=\",Emat//MatrixForm];\nFor [p=1,p<=5,p++, \
Print[\"Gauss rule p=\",p];   \n    \
Ke=Quad4IsoPRingStiffness[ncoor,Emat,{numer,p,Jcons,Kfac}];\n    \
Ke=Simplify[Ke]; Print[\"Ke=\",Ke//MatrixForm]; \n    Print[\"Eigenvalues of \
Ke=\",Chop[Eigenvalues[N[Ke]]]];\n    ]; *)\n    \n(* \
ClearAll[a,b,d,h,p,numer]; \na=6; b=2; d=1; br=3; bz=1;\nJcons=False; \
numer=True;\nncoor={{d,0},{a+d,0},{a+d,b},{d,b}}; \nFor [p=1,p<=2,p++",
    FontColor->RGBColor[0, 0, 1]],
  ",",
  StyleBox[" Print[\"Gauss rule p=\",p];   \n     \
fe=Quad4IsoPRingBodyForces[ncoor,{3,-1},{numer,p}];\n     Print[\"fe \
=\",Transpose[Partition[fe,2]]//MatrixForm];  \n     \
bfor={{1,0},{6,0},{6,0},{1,0}};\n     \
fe=Quad4IsoPRingBodyForces[ncoor,bfor,{numer,p}];\n     Print[\"fe \
=\",Transpose[Partition[fe,2]]//MatrixForm];\n     ]; *)\n     \n(* \
ClearAll[Em,\[Nu],a,b,d,err,ezz,grz,ur,uz,r,z]; \nEm=2500; \[Nu]=1/4; \n\
{err,ezz,err,grz}={3/80,-1/40,3/80,4/50};\n\
ncoor={{d,0},{a+d,0},{a+d,b},{d,b}}; num=False;\n\
Emat=Em/((1+\[Nu])*(1-2*\[Nu]))*{{1+\[Nu],\[Nu],\[Nu],0},{\[Nu],1+\[Nu],\[Nu],\
0},\n         {\[Nu],\[Nu],1+\[Nu],0},{0,0,0,1/2-\[Nu]}};\n\
{err,ezz,err,grz}={3/80,-1/40,3/80,4/50};\nur[r_,z_]:=err*r; \
uz[r_,z_]:=ezz*z+grz*r;\nue=Table[{0,0},{4}];\nFor \
[n=1,n<=4,n++,{rn,zn}=ncoor[[n]]; \n     ue[[n]]={ur[rn,zn],uz[rn,zn]} ];\n\
ue=Flatten[ue]; Print[\"ue=\",ue]; \n\
sige=Quad4IsoPRingStresses[ncoor,Emat,ue,{}];\nPrint[\"Corner \
stresses=\",sige//MatrixForm]; *)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2B:  Modules for 8-node iso-P serendipity quadrilateral ring \
element 
 for axisymmetric solids. These are documented in AFEM Chapter 12.\
\>", \
"Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad8IsoPRingStiffness[ncoor_,Emat_,options_]:=Module[\n  \
{p=2,numer=False,Jcons=False,Kfac=1,qcoor,k,\n   \
r1,r2,r3,r4,r5,r6,r7,r8,z1,z2,z3,z4,z5,z6,z7,z8,\n   \
Nf,N1,N2,N3,N4,N5,N6,N7,N8,\n   dNr1,dNr2,dNr3,dNr4,dNr5,dNr6,dNr7,dNr8,\n   \
dNz1,dNz2,dNz3,dNz4,dNz5,dNz6,dNz7,dNz8,\n   \
rk,w,c,A0,Jdet,Be,Ke,Ke0=Table[0,{16},{16}],\n   \
modname=\"Quad8IsoPRingStiffness: \"}, Ke=Ke0;\n  If \
[Length[options]==1,{numer}=options];\n  If \
[Length[options]==2,{numer,p}=options];\n  If \
[Length[options]==3,{numer,p,Jcons}=options];\n  If \
[Length[options]==4,{numer,p,Jcons,Kfac}=options];\n  If [p<1||p>5, \
Print[modname,\"illegal p:,p\"]; Return[Ke0]];\n  \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4},\n   {r5,z5},{r6,y6},{r7,z7},{r8,z8}}=ncoor;\
\n  A0=((r5-r7)*(z6-z8)-(r6-r8)*(z5-z7))/4;\n  If [numer&&(A0<=0), \
Print[modname,\"Neg or zero area\"]; \n      Return[Ke0]];  \n  For \
[k=1,k<=p*p,k++, \n      {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];\n     \
{{N1,N2,N3,N4,N5,N6,N7,N8},\n      {dNr1,dNr2,dNr3,dNr4,dNr5,dNr6,dNr7,dNr8},\
\n      {dNz1,dNz2,dNz3,dNz4,dNz5,dNz6,dNz7,dNz8},\n      \
Jdet}=Quad8IsoPRingShapeFunDer[ncoor,qcoor,Jcons];\n      If \
[numer&&(Jdet<=0), Print[modname,\"Neg or zero\",\n         \" Gauss point \
Jacobian at k=\",k]; Return[Ke0]];\n      \
rk=r1*N1+r2*N2+r3*N3+r4*N4+r5*N5+r6*N6+r7*N7+r8*N8; \n      Be={{ dNr1,   0, \
dNr2,   0, dNr3,   0, dNr4,   0,\n            dNr5,   0, dNr6,   0, dNr7,   \
0, dNr8,   0},\n          {    0,dNz1,    0,dNz2,    0,dNz3,    0,dNz4,\n     \
          0,dNz5,    0,dNz6,    0,dNz7,    0,dNz8},\n          {N1/rk,   \
0,N2/rk,   0,N3/rk,   0,N4/rk,   0,\n           N5/rk,   0,N6/rk,   0,N7/rk,  \
 0,N8/rk,   0},\n          { dNz1,dNr1, dNz2,dNr2, dNz3,dNr3, dNz4,dNr4,\n    \
        dNz5,dNr5, dNz6,dNr6, dNz7,dNr7, dNz8,dNr8}};\n      \
c=Kfac*w*rk*Jdet; If[!numer, Be=Simplify[Be]];\n      If [numer,Be=N[Be]; \
c=N[c]]; \n      Ke+=c*Transpose[Be].(Emat.Be); \n      ]; Return[Ke] ];\n\n\
Quad8IsoPRingShapeFunDer[ncoor_,qcoor_,Jcons_]:=Module[\n  \
{r1,r2,r3,r4,r5,r6,r7,r8,z1,z2,z3,z4,z5,z6,z7,z8,\n   \
\[Xi],\[Eta],rv,zv,A0,dN\[Xi],dN\[Eta],N1B,N2B,N3B,N4B,J11,J12,J21,J22,\n   \
Nf,dNr,dNz,Jdet}, {\[Xi],\[Eta]}=qcoor; \n  \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4},{r5,z5},{r6,z6},{r7,z7},\n  {r8,z8}}=ncoor; \
A0=((r5-r7)*(z6-z8)-(r6-r8)*(z5-z7))/4;  \n   N1B=(1-\[Xi])*(1-\[Eta])/4; \
N2B=(1+\[Xi])*(1-\[Eta])/4; \n   N3B=(1+\[Xi])*(1+\[Eta])/4; \
N4B=(1-\[Xi])*(1+\[Eta])/4;\n   Nf={-N1B*(1+\[Xi]+\[Eta]),  -N2B*(1-\[Xi]+\
\[Eta]),-N3B*(1-\[Xi]-\[Eta]),\n       -N4B*(1+\[Xi]-\[Eta]), \
2*N1B*(1+\[Xi]), 2*N3B*(1-\[Eta]),\n      2*N3B*(1-\[Xi]),   \
2*N4B*(1-\[Eta])};\n   dN\[Xi]={(1-\[Eta])*(2*\[Xi]+\[Eta]),(1-\[Eta])*(2*\
\[Xi]-\[Eta]),(1+\[Eta])*(2*\[Xi]+\[Eta]),\n         (1+\[Eta])*(2*\[Xi]-\
\[Eta]), 4*\[Xi]*(\[Eta]-1),2*(1-\[Eta]^2),\n         \
-4*\[Xi]*(1+\[Eta]),-2*(1-\[Eta]^2)}/4;\n   dN\[Eta]={(1-\[Xi])*(\[Xi]+2*\
\[Eta]),-(1+\[Xi])*(\[Xi]-2*\[Eta]),(1+\[Xi])*(\[Xi]+2*\[Eta]),\n         \
-(1-\[Xi])*(\[Xi]-2*\[Eta]), -2*(1-\[Xi]^2),-4*(1+\[Xi])*\[Eta],\n          \
2*(1-\[Xi]^2),-4*(1-\[Xi])*\[Eta]}/4;\n   rv={r1,r2,r3,r4,r5,r6,r7,r8}; \
zv={z1,z2,z3,z4,z5,z6,z7,z8}; \n   J11=dN\[Xi].rv; J12=dN\[Xi].zv; J21=dN\
\[Eta].rv; J22=dN\[Eta].zv;\n   Jdet=Simplify[J11*J22-J12*J21]; If \
[Jcons,Jdet=A0];\n   dNr=( J22*dN\[Xi]-J12*dN\[Eta])/Jdet; \n   dNz=(-J21*dN\
\[Xi]+J11*dN\[Eta])/Jdet;\n   Return[{Nf,dNr,dNz,Jdet}] ];\n   \n\
Quad8IsoPRingBodyForces[ncoor_,bfor_,options_]:=Module[\n  \
{p=2,numer=False,Jcons=False,Kfac=1,qcoor,k,m,mOK,\n   \
r1,r2,r3,r4,r5,r6,r7,r8,z1,z2,z3,z4,z5,z6,z7,z8,\n   \
Nf,N1,N2,N3,N4,N5,N6,N7,N8,dNr,dNz,\n   br1,br2,br3,br4,br5,br6,br7,br8,\n   \
bz1,bz2,bz3,bz4,bz5,bz6,bz7,bz8,\n   rk,w,c,A0,Jdet,fe,fe0=Table[0,{16}],\n   \
modname=\"Quad8IsoPRingBodyForces: \"}, fe=fe0;\n  If \
[Length[options]==1,{numer}=options];\n  If \
[Length[options]==2,{numer,p}=options];\n  If \
[Length[options]==3,{numer,p,Jcons}=options];\n  If \
[Length[options]==4,{numer,p,Jcons,Kfac}=options];\n  If [p<1||p>5, \
Print[modname,\"illegal p:\",p]; Return[fe0]];\n  \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4},\n   {r5,z5},{r6,z6},{r7,z7},{r8,z8}}=ncoor;\
\n  A0=((r5-r7)*(z6-z8)-(r6-r8)*(z5-z7))/4;\n  If [numer&&(A0<=0), \
Print[modname,\"Neg or zero area\"]; \n      Return[fe0]]; m=Length[bfor]; \
mOK=MemberQ[{2,4,8},m];\n  If [!mOK, Print[modname,\" Illegal bfor\"]; \
Return[fe0]]; \n  If [m==2, br1=br2=br3=br4=br5=br6=br7=br8=bfor[[1]]; \n     \
       bz1=bz2=bz3=bz4=bz5=bz6=bz7=bz8=bfor[[2]]];\n  If \
[m==4,{{br1,bz1},{br2,bz2},{br3,bz3},{br4,bz4}}=bfor;\n            \
{br5,bz5,br6,bz6,br7,bz7,br8,bz8}={br1+br2,\n             \
bz1+bz2,br2+br3,bz2+bz3,br3+br4,bz3+bz4,\n             br4+br1,bz4+bz1}/2];\n \
 If [m==8,{{br1,bz1},{br2,bz2},{br3,bz3},{br4,bz4},\n            \
{br5,bz5},{br6,bz6},{br7,bz7},{br8,bz8}}=bfor]; \n  For [k=1,k<=p*p,k++, \n   \
   {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n     {{N1,N2,N3,N4,N5,N6,N7,N8},dNr,dNz,Jdet}=\n           \
Quad8IsoPRingShapeFunDer[ncoor,qcoor,Jcons];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      If [numer&&(Jdet<=0), Print[modname,\"Neg or zero\",\n         \" \
Gauss point Jacobian at k=\",k]; Return[fe0]];\n      \
rk=r1*N1+r2*N2+r3*N3+r4*N4+r5*N5+r6*N6+r7*N7+r8*N8;",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      brk=br1*N1+br2*N2+br3*N3+br4*N4+br5*N5+br6*N6+br7*N7+br8*N8;\n     \
 bzk=bz1*N1+bz2*N2+bz3*N3+bz4*N4+bz5*N5+bz6*N6+bz7*N7+bz8*N8;\n      \
bk={N1*brk,N1*bzk,N2*brk,N2*bzk,N3*brk,N3*bzk,N4*brk,N4*bzk,\n          \
N5*brk,N5*bzk,N6*brk,N6*bzk,N7*brk,N7*bzk,N8*brk,N8*bzk};\n      \
c=Kfac*w*Jdet*rk; If [numer,bk=N[bk];c=N[c]]; fe+=c*bk;   \n      ];  \
Return[fe] ];\n      \nQuad8IsoPRingStresses[ncoor_,Emat_,ue_,options_]:= \n  \
Module[{numer=False,g=1/Sqrt[3],Jcons=False,w0=0,\n   \
eps=10.^(-9),r1,r2,r3,r4,r5,r6,r7,r8,z1,z2,z3,z4,\n   \
z5,z6,z7,z8,Nf,N1,N2,N3,N4,N5,N6,N7,N8,\n   \
dNr1,dNr2,dNr3,dNr4,dNr5,dNr6,dNr7,dNr8,\n   \
dNz1,dNz2,dNz3,dNz4,dNz5,dNz6,dNz7,dNz8,\n   \
T1,T2,T3,T4,T5,T6,Td,Tg8,Jdet,qcoor,w,c,Be,\n   \
gctab={{0,0}},k,kg,rk,sigg,sige,udis=ue,\n   modname=\"Quad8IsoPRingStresses: \
\"},\n  If [Length[options]==1,{numer}=options];\n  If \
[Length[options]==2,{numer,g}=options];\n  If \
[Length[options]==3,{numer,g,w0}=options];\n  If [Head[g]==Symbol||g>0, \
Td=4*g^2*(4+w0); \n     T1=4*g^2*w0; T2=4+4*g^2+w0+2*g*(4+w0); \n     \
T3=-4+4*g^2-w0; T4=4+4*g^2+w0-2*g*(4+w0);\n     T5=g*(4+4*g+w0); \
T6=g*(-4+4*g-w0);\n     Tg8={{T1,T2,T3,T4,T3},{T1,T3,T2,T3,T4},\n          \
{T1,T4,T3,T2,T3},{T1,T3,T4,T3,T2},\n          \
{T1,T5,T5,T6,T6},{T1,T6,T5,T5,T6},\n          \
{T1,T6,T6,T5,T5},{T1,T5,T6,T6,T5}}/Td;",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "     \n     gctab={{0,0},{-1,-1},{1,-1},{1,1},{-1,1}}*g];\n  \
kg=Length[gctab]; sigg=Table[{0,0,0,0},{kg}];\n  If [numer, gctab=N[gctab]; \
Tg8=N[Tg8]; udis=N[ue]]; \n  {{r1,z1},{r2,z2},{r3,z3},{r4,z4},\n   \
{r5,z5},{r6,z6},{r7,z7},{r8,z8}}=ncoor;\n  For [k=1,k<=kg,k++, \
qcoor=gctab[[k]];\n     {{N1,N2,N3,N4,N5,N6,N7,N8},\n      \
{dNr1,dNr2,dNr3,dNr4,dNr5,dNr6,dNr7,dNr8},\n      \
{dNz1,dNz2,dNz3,dNz4,dNz5,dNz6,dNz7,dNz8},\n      \
Jdet}=Quad8IsoPRingShapeFunDer[ncoor,qcoor,Jcons];\n      \
rk=r1*N1+r2*N2+r3*N3+r4*N4+r5*N5+r6*N6+r7*N7+r8*N8; \n      Be={{ dNr1,   0, \
dNr2,   0, dNr3,   0, dNr4,   0,\n            dNr5,   0, dNr6,   0, dNr7,   \
0, dNr8,   0},\n          {    0,dNz1,    0,dNz2,    0,dNz3,    0,dNz4,\n     \
          0,dNz5,    0,dNz6,    0,dNz7,    0,dNz8},\n          {N1/rk,   \
0,N2/rk,   0,N3/rk,   0,N4/rk,   0,\n           N5/rk,   0,N6/rk,   0,N7/rk,  \
 0,N8/rk,   0},\n          { dNz1,dNr1, dNz2,dNr2, dNz3,dNr3, dNz4,dNr4,\n    \
        dNz5,dNr5, dNz6,dNr6, dNz7,dNr7, dNz8,dNr8}};\n      If \
[numer,Be=N[Be]]; sigg[[k]]=Emat.(Be.udis)\n      ] ;\n  If [kg==1, \
sige=Table[sigg[[1]],{4}], sige=Tg8.sigg];\n  If [numer, \
sige=Chop[sige,eps]]; \n  If [!numer,sige=Simplify[sige]]; Return[sige] ];",
  StyleBox["\n      ",
    FontColor->RGBColor[0, 0, 1]],
  "  \n(*",
  StyleBox["ClearAll[Em,\[Nu],a,b,d,p,numer];\na=4; b=2; d=0; Em=96; \
\[Nu]=1/3; \nncoor={{d,0},{a+d,0},{a+d,b},{d,b},{a/2+d,0},{a+d,b/2},\n       \
{a/2+d,b},{d,b/2}};\n\
Emat=Em/((1+\[Nu])*(1-2*\[Nu]))*{{1-\[Nu],\[Nu],\[Nu],0},{\[Nu],1-\[Nu],\[Nu],\
0},\n                {\[Nu],\[Nu],1-\[Nu],0},{0,0,0,1/2-\[Nu]}};\n\
Print[\"Emat=\",Emat//MatrixForm];\nFor [p=1,p<=4,p++, \n     \
Ke=Quad8IsoPRingStiffness[ncoor,Emat,{True,p}];\n     Ke=Simplify[Chop[Ke]];  \
Print[\"Ke=\",Ke//MatrixForm]; \n     Print[\"Eigenvalues of \
Ke=\",Chop[Eigenvalues[N[Ke]],.000001]]; \n     ]; *)\n     \n\
(*ClearAll[a,b,d,p];\na=3; b=2; d=1;\n\
ncoor={{d,0},{a+d,0},{a+d,b},{d,b},{a/2+d,0},{a+d,b/2},\n       \
{a/2+d,b},{d,b/2}};\nFor [p=1,p<=3,p++, For [case=1,case<=2,case++,\n     If \
[case==1, bfor={36,-18}];\n     If [case==2, \
bfor=Table[{60*ncoor[[i,1]],0},{i,8}]];\n     \
fe=Quad8IsoPRingBodyForces[ncoor,bfor,{True,p}];\n     fe=Simplify[Chop[fe]]; \
 \n     Print[\"fe=\",Transpose[Partition[fe,2]]//MatrixForm];\n     \
frsum=Sum[fe[[2*i-1]],{i,8}]; fzsum=Sum[fe[[2*i]],{i,8}];\n     \
Print[\"frsum=\",frsum,\" fzsum=\",fzsum];\n     ]]; *)\n     \n\
(*ClearAll[Em,\[Nu],a,b,d,err,ezz,grz,ur,uz,r,z]; \nEm=2500; \[Nu]=1/4; d=1; \
a=3; b=2;\n{err,ezz,err,grz}={3/80,-1/40,3/80,4/50};\n\
ncoor={{d,0},{a+d,0},{a+d,b},{d,b},{a/2+d,0},{a+d,b/2},{a/2+d,b},{d,b/2}}; \n\
Emat=Em/((1+\[Nu])*(1-2*\[Nu]))*{{1+\[Nu],\[Nu],\[Nu],0},{\[Nu],1+\[Nu],\[Nu],\
0},\n         {\[Nu],\[Nu],1+\[Nu],0},{0,0,0,1/2-\[Nu]}};\n\
{err,ezz,err,grz}={3/80,-1/40,3/80,4/50};\nur[r_,z_]:=err*r; \
uz[r_,z_]:=ezz*z+grz*r;\nue=Table[{0,0},{8}];\nFor [n=1,n<=8,n++, \
{rn,zn}=ncoor[[n]]; \n     ue[[n]]={ur[rn,zn],uz[rn,zn]} ];\nue=Flatten[ue]; \
Print[\"ue=\",ue]; \nsige=Quad8IsoPRingStresses[ncoor,Emat,ue,{True}];\n\
Print[\"Corner stresses=\",sige//MatrixForm];*)",
    FontColor->RGBColor[0, 0, 1]],
  "\n"
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 3:  Assembly of master stiffness as a full symmetric matrix.  \
Axisymmetric,
2 dof per node.  Two iso-P quads implemented: 4-Node and 8-Node.
Uses element  stiffness modules,which must be initialized at this \
point.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "RingMasterStiffness[nodrzc_,eletyp_,elenod_,elemat_,defopt_]:=\n  \
Module[{numele=Length[elenod],numnod=Length[nodrzc],\n   \
numer=False,p=2,Jcons=False,Kfac=1,e,type,etype,enl,eftab,\n   \
i,ii,j,jj,n,nn,Em,\[Nu],Emat,ncoor,Ke,neldof,K,\n   \
modname=\"RingMasterStiffness: \"},\n   K=Table[0,{2*numnod},{2*numnod}];\n   \
If [Length[defopt]==1,{numer}=defopt];\n   If \
[Length[defopt]==2,{numer,p}=defopt];\n   If \
[Length[defopt]==3,{numer,p,Jcons}=defopt];\n   If \
[Length[defopt]==4,{numer,p,Jcons,Kfac}=defopt];\n   For [e=1,e<=numele,e++,\n\
       type=eletyp[[e]]; etype=StringTake[type,5]; \n       If \
[etype!=\"Quad4\"&&etype!=\"Quad8\", Print[modname,\n          \"Illegal \
element type, e=\",e]; Return[Null]];\n       enl=elenod[[e]]; \
nn=Length[enl];\n       \
eftab=Flatten[Table[{2*enl[[i]]-1,2*enl[[i]]},{i,nn}]];\n       \
ncoor=Table[nodrzc[[enl[[i]]]],{i,nn}];        \n       \
{Em,\[Nu]}=elemat[[e]]; \n        Emat=Em/((1+\[Nu])*(1-2*\[Nu]))*{{1-\[Nu],\
\[Nu],\[Nu],0},{\[Nu],1-\[Nu],\[Nu],0},\n                \
{\[Nu],\[Nu],1-\[Nu],0},{0,0,0,1/2-\[Nu]}};\n       If \
[numer,ncoor=N[ncoor];Emat=N[Emat]];\n       options={numer,p,Jcons,Kfac}; \n \
      If [etype==\"Quad4\", \n          If [type==\"Quad4.1\", \
options={numer,1,Jcons,Kfac}]; \n          \
Ke=Quad4IsoPRingStiffness[ncoor,Emat,options] ];\n       If \
[etype==\"Quad8\",\n          If [type==\"Quad8.3\", \
options={numer,3,Jcons,Kfac}];  \n          \
Ke=Quad8IsoPRingStiffness[ncoor,Emat,options] ];\n       If \
[numer,Ke=Chop[Ke]]; If [!numer,Ke=Simplify[Ke]];\n  ",
  StyleBox["     (*Print[\"Ke=\",Ke//MatrixForm];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n       neldof=Length[Ke];\n       For [i=1,i<=neldof,i++,ii=eftab[[i]];\n\
           For [j=i,j<=neldof,j++,jj=eftab[[j]];\n               \
K[[jj,ii]]=K[[ii,jj]]+=Ke[[i,j]] ]];\n       ]; \n  Return[K]; \n];\n\n(*",
  StyleBox["ClearAll[Em,\[Nu],a,b,e,p,numer,exx,eyy,gxy]; \na=2; b=1; Em=12; \
\[Nu]=1/3; \nnodrzc={{0,b},{0,0},{a/2,b},{a/2,0},{a,b},{a,0}};\neletyp=Table[\
\"Quad4\",{2}]; \nelenod={{1,2,4,3},{3,4,6,5}};\n\
elemat=Table[{Em,\[Nu]},{2}];\ndefopt={True,2,False};\n\
K=RingMasterStiffness[nodrzc,eletyp,elenod,elemat,defopt];\n\
Print[\"K=\",Chop[K]//MatrixForm];\nPrint[\"diag of \
K=\",Table[K[[i,i]],{i,12}]]; \nPrint[\"eigs of \
K=\",Chop[Eigenvalues[N[K]]]];",
    FontColor->RGBColor[0, 0, 1]],
  "\n",
  StyleBox["\nClearAll[Em,\[Nu],a,b,e,p,numer,exx,eyy,gxy]; \na=2; b=1; \
Em=12; \[Nu]=1/3; \n\
nodrzc={{0,b},{0,b/2},{0,0},{a/2,b},{a/2,0},{a,b},{a,b/2},{a,0}};\n\
eletyp=Table[\"Quad8\",{1}]; \nelenod={{1,3,8,6,2,5,7,4}};\nelemat=Table[{Em,\
\[Nu]},{1}];\ndefopt={True,2,False};\n\
K=RingMasterStiffness[nodrzc,eletyp,elenod,elemat,defopt];\n\
Print[\"K=\",Chop[K]//MatrixForm];\nPrint[\"diag of \
K=\",Table[K[[i,i]],{i,16}]]; \nPrint[\"eigs of \
K=\",Chop[Eigenvalues[N[K]]]];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 4:  Application of displacement boundary conditions on master \
stiffness equations.  
Code only takes single freedom constraints,but accepts nonzero prescribed \
displacements.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ModifiedMasterStiffness[pdof_,Km_]:= Module[{i,j,k,K},K=Km; \n      For \
[k=1,k<=Length[pdof],k++,i=pdof[[k]];\n          For [j=1,j<=Length[K],j++,\n \
             K[[i,j]]=K[[j,i]]=0]; K[[i,i]]=1;\n          ]; Return[K];\n];\n\
ModifiedNodeForces[pdof_,pdofv_,Km_,nfv_,bfor_]:= Module[\n    \
{i,j,k,l,d,kk=Length[pdof],n=Length[Km],fixed,rhs},\n    rhs=nfv+bfor; \
d=pdofv; fixed=Table[False,{n}]; \n    Do [i=pdof[[k]]; \
fixed[[i]]=True,{k,1,kk}];\n    For [k=1,k<=kk,k++,i=pdof[[k]]; \n        For \
[j=1,j<=n,j++,If [fixed[[j]],Continue[]];  \n             \
rhs[[j]]=rhs[[j]]-Km[[i,j]]*d[[k]]; \n            ]; rhs[[i]]=d[[k]]; \n      \
  ]; \n    Return[rhs];\n];\n        \n",
  StyleBox["(*K=Array[Ks,{6,6}];\nPrint[\"Master stiffness \
matrix:\"];Print[K//TableForm];\np=Array[ps,{6}];\nPrint[\"Node force vector:\
\"];Print[p];\np=ModifiedNodeForces[{1,2,4},{d1,d2,d3},K,p];\nPrint[\"Node \
Force vector modified for displacement B.C.:\"];\nPrint[p];\n\
K=ModifiedMasterStiffness[{1,2,4},K];\nPrint[\"Master stiffness modified for \
displacement B.C.:\"];\nPrint[K//TableForm];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 5:  Assembly of consistent force vector associated with body \
forces.
 Axisymmetric, 2 dof per node. Both 4-node and 8-node iso-P quads \
implemented.
Uses element  body force modules,which must be initialized at this \
point.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "RingBodyForces[nodrzc_,eletyp_,elenod_,elebfor_,defopt_]:=\n  \
Module[{numele=Length[elenod],numnod=Length[nodrzc],\n   \
numer=False,p=2,Jcons=False,Kfac=1,e,type,etype,enl,eftab,\n   i,ii,n,nn,Em,\
\[Nu],Emat,ncoor,neldof,fe,fe0,f,\n   modname=\"RingBodyForces: \"},\n   \
f=Table[0,{2*numnod}]; fe0=Table[0,{16}];\n   If \
[Length[defopt]==1,{numer}=defopt];\n   If \
[Length[defopt]==2,{numer,p}=defopt];\n   If \
[Length[defopt]==3,{numer,p,Jcons}=defopt];\n   If \
[Length[defopt]==4,{numer,p,Jcons,Kfac}=defopt];\n   \
options={numer,p,Jcons,Kfac}; \n   For [e=1,e<=numele,e++,\n       \
type=eletyp[[e]]; etype=StringTake[type,5]; \n       If \
[etype!=\"Quad4\"&&etype!=\"Quad8\", Print[modname,\n          \"Illegal \
element type,e=\",e]; Return[Null]];\n       enl=elenod[[e]]; nn=Length[enl];\
\n       eftab=Flatten[Table[{2*enl[[i]]-1,2*enl[[i]]},{i,nn}]];\n       \
ncoor=Table[nodrzc[[enl[[i]]]],{i,nn}]; \n       bfor=elebfor[[e]];\n       \
If [numer,ncoor=N[ncoor]; bfor=N[bfor]];\n       If [etype==\"Quad4\", \n     \
      fe=Quad4IsoPRingBodyForces[ncoor,bfor,options] ];\n       If \
[etype==\"Quad8\",       \n           \
fe=Quad8IsoPRingBodyForces[ncoor,bfor,options] ];\n       neldof=Length[fe];\n\
       For [i=1,i<=neldof,i++,ii=eftab[[i]];\n            f[[ii]]+=fe[[i]] ];\
\n       ];\n  Return[f]; \n];\n\n(*",
  StyleBox["ClearAll[a,b,e,p];  a=2; b=1;  \n\
nodrzc={{0,b},{0,0},{a/2,b},{a/2,0},{a,b},{a,0}};\nelenod=  \
{{1,2,4,3},{3,4,6,5}};\n(*Plot2DElementsAndNodes[nodrzc,elenod,1/2.5,\n       \
\"Mesh for body forces\",True,True];*)\neletyp=  {\"Quad4\",\"Quad4\"};\n\
elebfor={{1,-2},{1,-2}}; defopt={True};\n\
f=RingBodyForces[nodrzc,eletyp,elenod,elebfor,defopt];\nPrint[\"Body force \
vector:\",Transpose[Partition[f,2]]//MatrixForm]; \n\nClearAll[a,b,e,p]; a=2; \
b=1; \nnodrzc={{0,b},{0,b/2},{0,0},{a/2,b},{a/2,0},{a,b},{a,b/2},{a,0}};\n\
eletyp=Table[\"Quad8\",{1}]; \nelenod={{1,3,8,6,2,5,7,4}};\nelebfor={{1,-2}}; \
defopt={True,2};\nf=RingBodyForces[nodrzc,eletyp,elenod,elebfor,defopt];\n\
Print[\"Body force vector:\",Transpose[Partition[f,2]]//MatrixForm];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 6: Computation of nodal stresses (averaged across all elements \
meeting at a node)  
in Axisymmetric structure.  Quad 4 and Quad8 implemented. Stress modules are \
expected to be initialized.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "RingNodalStresses[nodrzc_,eletyp_,elenod_,elemat_,\n  \
noddis_,defopt_]:=Module[{numele=Length[elenod],\n  \
numnod=Length[nodrzc],numer=False,g=1/Sqrt[3],\n  \
w0=0,options,e,type,etype,enl,eftab,i,j,k,\n  \
n,nn,Em,\[Nu],Emat,ncoor,ue,ncount,esig,sige,sig,\n  \
modname=\"RingNodalStresses: \"}, \n  sige=Table[{0,0,0,0},{numele}]; \n  \
sig= Table[{0,0,0,0},{numnod}]; \n  If [Length[defopt]==1,{numer}=defopt];\n  \
If [Length[defopt]==2,{numer,g}=defopt];\n  If \
[Length[defopt]==3,{numer,g,w0}=defopt];\n  options={numer,g,w0}; \
ncount=Table[0,{numnod}]; ",
  StyleBox["\n  ",
    FontColor->RGBColor[1, 0, 0]],
  "For [e=1,e<=numele,e++, \n      type=eletyp[[e]]; \
etype=StringTake[type,5]; \n      enl=elenod[[e]]; nn=Length[enl]; \n      \
eftab=Flatten[Table[{2*enl[[i]]-1,2*enl[[i]]},{i,nn}]];\n      \
ncoor=Table[nodrzc[[enl[[i]]]],{i,nn}];\n      \
ue=Flatten[Table[noddis[[enl[[i]]]],{i,nn}]];\n      {Em,\[Nu]}=elemat[[e]]; \
\n      Emat=Em/((1+\[Nu])*(1-2*\[Nu]))*{{1-\[Nu],\[Nu],\[Nu],0},{\[Nu],1-\
\[Nu],\[Nu],0},\n                {\[Nu],\[Nu],1-\[Nu],0},{0,0,0,1/2-\[Nu]}}; \
\n      If [numer,ncoor=N[ncoor];ue=N[ue];Emat=N[Emat]];\n      If \
[etype==\"Quad4\",\n          \
esig=Quad4IsoPRingStresses[ncoor,Emat,ue,options]",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "];\n      If [etype==\"Quad8\",\n          \
esig=Quad8IsoPRingStresses[ncoor,Emat,ue,options]",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "];\n      For [i=1,i<=nn,i++, n=enl[[i]]; ncount[[n]]++; \n            \
sig[[n]]+=esig[[i]] ]; \n      ];\n   For [n=1,n<=numnod,n++,k=ncount[[n]];If \
[k>1,sig[[n]]/=k]];\n   Return[sig]];\n   \n(*",
  StyleBox["ClearAll[Em,\[Nu],a,b,exx,eyy,gxy]; a=2; b=1;\nEm=120; \[Nu]=1/3; \
exx=0.15; eyy=-0.125; gxy=-0.070;\n\
nodrzc={{0,b},{0,0},{a/2,b},{a/2,0},{a,b},{a,0}};\nelenod=  \
{{1,2,4,3},{3,4,6,5}};\nnumnod=Length[nodrzc]; numele=Length[elenod];\n\
eletyp=Table[\"Quad4\",{numele}]; \nelemat=Table[{Em,\[Nu]},{numele}]; \n\
defopt={True,2,1,1};\nr=Table[nodrzc[[i,1]],{i,numnod}]; \
z=Table[nodrzc[[i,2]],{i,numnod}];\n\
noddis=Table[{r[[i]]*exx,r[[i]]*gxy+z[[i]]*eyy},{i,numnod}];\nPrint[\"noddis=\
\",Transpose[noddis]//MatrixForm];\n\
sig=RingNodalStresses[nodrzc,eletyp,elenod,elemat,noddis,defopt];\n\
Print[\"Computed node stresses=\",Transpose[Chop[N[sig]]]//MatrixForm]; *)\n\n\
(*ClearAll[Em,\[Nu],a,b,e,p,numer,exx,eyy,gxy]; a=2; b=1;\nEm=120; \[Nu]=1/3; \
exx=0.15; eyy=-0.125; gxy=-0.070;\n\
nodrzc={{0,b},{0,b/2},{0,0},{a/2,b},{a/2,0},{a,b},{a,b/2},{a,0}};\n\
elenod={{1,3,8,6,2,5,7,4}};\nnumnod=Length[nodrzc]; numele=Length[elenod];\n\
eletyp=Table[\"Quad8\",{numele}];  elemat=Table[{Em,\[Nu]},{numele}];\n\
defopt={False};\nr=Table[nodrzc[[i,1]],{i,numnod}]; \
z=Table[nodrzc[[i,2]],{i,numnod}];\n\
noddis=Table[{r[[i]]*exx,r[[i]]*gxy+z[[i]]*eyy},{i,numnod}]; \n\
Print[\"noddis=\",Transpose[noddis]//MatrixForm];\n\
sig=RingNodalStresses[nodrzc,eletyp,elenod,elemat,noddis,defopt];\n\
Print[\"Computed node stresses=\",Transpose[Chop[sig]]//MatrixForm]; *)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 7:  Mesh plot modules.  Written by C. A. Felippa,1997/98. 
Mathematica adaptation of ancient (1966) Fortran code.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Plot2DElements[xycoor_,elenod_,aspect_,label_]:= \n\
Module[{enl,n,nc,xyc,poly=sides={},numele=Length[elenod]},\n Do \
[enl=elenod[[ne]]; nc=Length[enl]; \n    If [nc!=3 && nc!=4,Continue[]]; \
xyc=Table[0,{nc+1}];\n    enl=AppendTo[enl,enl[[1]]];    \n    Do \
[n=enl[[i]]; xyc[[i]]=xycoor[[n]],{i,1,nc+1}];\n    poly= \
AppendTo[poly,Graphics[Polygon[Take[xyc,nc]]]];\n    \
sides=AppendTo[sides,Graphics[Line[xyc]]],\n {ne,1,numele}];\n \
Show[Graphics[RGBColor[1,1,0]],poly,\n      \
Graphics[Thickness[0.008]],Graphics[RGBColor[0,0,0]],\n      \
sides,AspectRatio->aspect,PlotLabel->label ];\n      ];\n      \n\
Plot2DElementsAndNodes[xycoor_,elenod_,aspect_,label_,\n  \
elabels_,nlabels_]:= Module[{i,j,k,e,enl,elab,n,nlab,\n   \
x,y,delx,dely,xmin,xmax,ymin,ymax,xy0,eps,ex,ey,b,nb,\n   \
numele=Length[elenod],rade,gside,gpoly,\n   \
numnod=Length[xycoor],radn,pnod=pnlab={},\n   \
pefill=pesid=pelab=pecirc=pedisk={}},\n   \
x=Table[xycoor[[i,1]],{i,1,numnod}]; \n   \
y=Table[xycoor[[i,2]],{i,1,numnod}];\n   \
{xmin,xmax,ymin,ymax}=N[{Min[x],Max[x],Min[y],Max[y]}];\n   \
{delx,dely}={xmax-xmin,ymax-ymin}/Sqrt[N[numnod]];\n   \
radn=0.07*Min[delx,dely]; {ex,ey}={-2.3*radn,1.25*radn}; \n   \
rade=0.16*Min[delx,dely]; \n   \
(*Print[\"delx,dely,rade=\",{delx,dely,rade}];*) \n   Do [enl=elenod[[e]]; \
n=Length[enl];\n       If [n==3,b={1,2,3};nb=3];\n       If \
[n==4,b={1,2,3,4};nb=4];\n       If [n==6,b={1,4,2,5,3,6};nb=6];\n       If \
[n==8,b={1,5,2,6,3,7,4,8};nb=8];\n       If [n==9,b={1,5,2,6,3,7,4,8};nb=8];\n\
       gside=gpoly=Table[xycoor[[enl[[b[[i]]]]]],{i,nb}]; \n       \
AppendTo[gside,gpoly[[1]]];\n       xy0={0,0};Do \
[xy0+=gside[[i]],{i,nb}];xy0=xy0/nb; \n       \
pefill=AppendTo[pefill,Graphics[Polygon[gpoly]]];\n       pesid= \
AppendTo[pesid, Graphics[Line[gside]]];\n       If [n==9,xy0+={0,-2*rade}];\n \
      If [elabels,pedisk= AppendTo[pedisk,\n            \
Graphics[Disk[xy0,rade]]];\n          \
pecirc=AppendTo[pecirc,Graphics[Circle[xy0,rade]]];\n          \
elab=FontForm[e,{\"Times\",11}];\n          pelab= \
AppendTo[pelab,Graphics[Text[elab,xy0]]]],\n   {e,1,numele}];\n   Do \
[nlab=FontForm[n,{\"Times\",11}];\n       If [nlabels,\n       \
pnlab=AppendTo[pnlab,\n             \
Graphics[Text[nlab,xycoor[[n]]+{ex,ey}]]]];\n       pnod= AppendTo[pnod, \n   \
          Graphics[Disk[xycoor[[n]],radn]]],\n   {n,1,numnod}];  \
eps=0.1*(ymax-ymin);\n   Show[Graphics[RGBColor[1,1,0]],pefill,\n        \
Graphics[GrayLevel[0]],\n        Graphics[AbsoluteThickness[1.5]],pesid,\n    \
    Graphics[GrayLevel[0]],pnod,\n        Graphics[AbsoluteThickness[0.8]],\n \
       Graphics[GrayLevel[1]],pedisk,\n        \
Graphics[GrayLevel[0]],pecirc,pelab,\n        Graphics[GrayLevel[0]],pnlab,\n \
       PlotLabel->label,\n        \
PlotRange->{{xmin-eps,xmax+eps},{ymin-eps,ymax+eps}},\n        \
Axes->False,AspectRatio->aspect];\n];\n\n",
  StyleBox["(*xycoor=N[{{0,10},{-8.66,5},{-8.66,-5},{0,-10},{8.66,-5},{8.66,5}\
,\n     {-4.33,2.5},{-4.33,-2.5},{4.33,-2.5},{4.33,2.5},{0,0}}];\n\
elenod={{1,2,7},{2,3,8,7},{3,4,8},{1,7,11},{7,8,11},{8,4,11},\n    \
{1,11,10},{11,9,10},{11,4,9},{1,10,6},{10,9,5,6},{4,5,9}};\n           \n\
aspect=10/8.66;\nPlot2DElements[xycoor,elenod,aspect,\"hexagon test mesh\"];\n\
Plot2DElementsAndNodes[xycoor,elenod,aspect,\n   \"hexagon test mesh with \
elem & node labels\",True,True];*)",
    FontColor->RGBColor[0, 0, 1]],
  "\n   "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 8:  Contour plot modules.  Written by C. A. Felippa, 1997. 
Redone 2008. Mathematica adaptation of ancient (1966) Fortran  code.\
\>", \
"Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ContourBandColor[f_,fmin_,fmax_]:= Module[{fs,h,s,b=1},\n   If [fmax==0 \
||fmin>=fmax, (* White if wrong inputs *)\n       \
Return[Graphics[RGBColor[1,1,1]] ]];\n   If [f>fmax || f<fmin, (* Black if \
outside range *)\n       Return[Graphics[RGBColor[0,0,0]] ]]; \n   \
fs=(fmax-f)/(fmax-fmin); h=0.70*fs; s=1-1/(1+100*(fs-1/2)^2);",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   Return[Graphics[Hue[h,s,b]]]];\n \n\
TrigBandSegment[xs_,ys_,fs_,fv_,eps_]:=Module[{x1,x2,x3,\n  \
y1,y2,y3,f1,f2,f3,x21,y21,x32,y32,x13,y13,f21,f32,f31,\n  \
\[Zeta]211,\[Zeta]212,\[Zeta]311,\[Zeta]313,\[Zeta]322,\[Zeta]323,P1,P3,P21,\
P32,P13}, \n  {x1,x2,x3}=xs; {y1,y2,y3}=ys; {f1,f2,f3}=fs; P1={x1,y1}; \
P3={x3,y3};\n  If [fv<f1, Return[{{P1,P1},-1}]]; If [fv>f3, \
Return[{{P3,P3},0}]];\n  f21=f2-f1; If [f21<eps,f1=f1-eps;f21=eps];\n  \
f32=f3-f2; If [f32<eps,f3=f2+eps;f32=eps]; f31=f3-f1; \n  \
\[Zeta]212=(fv-f1)/f21; \[Zeta]212=Max[Min[\[Zeta]212,1],0]; \[Zeta]211=1-\
\[Zeta]212;\n  \[Zeta]313=(fv-f1)/f31; \[Zeta]313=Max[Min[\[Zeta]313,1],0]; \
\[Zeta]311=1-\[Zeta]313;",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "  P21={x2*\[Zeta]212+x1*\[Zeta]211,y2*\[Zeta]212+y1*\[Zeta]211};\n  \
P13={x1*\[Zeta]311+x3*\[Zeta]313,y1*\[Zeta]311+y3*\[Zeta]313};\n  If \
[fv<=f2,Return[{{P21,P13},1}]];\n  \[Zeta]323=(fv-f2)/f32; \
\[Zeta]323=Max[Min[\[Zeta]323,1],0]; \[Zeta]322=1-\[Zeta]323; \n  P32={x3*\
\[Zeta]323+x2*\[Zeta]322,y3*\[Zeta]323+y2*\[Zeta]322};\n  \
Return[{{P32,P13},2}]];\n \n\
ContourBandPlotFunctionOnTrig3[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\n  \
lines_]:=Module[{nv,iv,xc1,yc1,xc2,yc2,xc3,yc3,fc1,fc2,fc3,\n  \
xs,ys,fs,x1,x2,x3,y1,y2,y3,f1,f2,f3,firsthit,Pc1,Pc2,Pc3,\n  \
ftop,fbot,favg,Pb,Pt,tb,tt,P1,P2,P3,P4,Q1,Q2,Q3,Q4,ptab,\n  \
flast,Plast,tlast,kb,kl,p,pb,pl}, \n  If [fmax<fmin||finc<=0, \
Return[{{},{}}]];\n  nv=Floor[(fmax-fmin+eps)/finc]; nv=Min[nv,1000];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n  pb=Table[{},{nv}]; pl=Table[{},{nv+1}]; \n \
{{xc1,yc1},{xc2,yc2},{xc3,yc3}}=N[xyc]; {fc1,fc2,fc3}=N[fc];\n  \
{{x1,y1,f1},{x2,y2,f2},{x3,y3,f3}}=Sort[{{xc1,yc1,fc1},\n          \
{xc2,yc2,fc2},{xc3,yc3,fc3}},#1[[3]]<=#2[[3]]&]; \n  xs={x1,x2,x3}; \
ys={y1,y2,y3}; fs={f1,f2,f3}; \n  Pc1={x1,y1}; Pc2={x2,y2}; Pc3={x3,y3};\n  \
ptab={{{},{Pc1,Pc2,Pc3},{Pc1,Q4,Q3},{Pc1,Pc2,Q3,Q4}},\n     \
{{Null},{},{Null},{Null}},\n     \
{{Null},{Pc3,Pc2,Q1,Q2},{Q1,Q2,Q4,Q3},{Pc2,Q1,Q2,Q4,Q3}},\n     \
{{Null},{Pc3,Q1,Q2},{Null},{Q1,Q2,Q4,Q3}}};\n  flast=fmin-1000; \
firsthit=True; kb=kl=0;\n  For [iv=1,iv<=nv,iv++, \n       \
fbot=N[fmin+(iv-1)*finc]; ftop=N[fmin+iv*finc];\n       If [ftop<=f1 || \
fbot>=f3, Continue[]]; favg=(fbot+ftop)/2;\n       If [fbot==flast, \
{fbot,Pb,tb}={flast,Plast,tlast},\n          \
{Pb,tb}=TrigBandSegment[xs,ys,fs,fbot,eps],\n          \
{Pb,tb}=TrigBandSegment[xs,ys,fs,fbot,eps]]; \n       \
{Pt,tt}=TrigBandSegment[xs,ys,fs,ftop,eps];\n       {P1,P2}=Pb; {P3,P4}=Pt; \
{flast,Plast,tlast}={ftop,Pt,tt};\n       If [lines,\n          If \
[firsthit&&(tb>0), ",
  StyleBox[" \n              ",
    FontColor->RGBColor[1, 0, 0]],
  "pl[[++kl]]=Graphics[Line[Pb]]; firsthit=False];\n\t      If [tt>0, \
pl[[++kl]]=Graphics[Line[Pt]]]\n\t      ];\n       \
p=ptab[[tb+2,tt+2]]/.{Q1->P1,Q2->P2,Q3->P3,Q4->P4};\n       \
pb[[++kb]]={ContourBandColor[favg,fmin,fmax],Graphics[Polygon[p]]};\n\t   ];  \
",
  StyleBox["\n  ",
    FontColor->RGBColor[1, 0, 0]],
  "If [kb>0,pb=Take[pb,kb],{},{}]; If [kl>0,pl=Take[pl,kl],{},{}];\n  \
ClearAll[ptab]; Return[{pb,pl}];\n ];\n \n\
ContourBandPlotFunctionOnQuad4[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\n   \
lines_]:=Module[{x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,f1,f2,f3,f4,f5,\n   \
e,xytab,ftab,xyt,ft, pbq,plq,pbt,plt },\n   pbq=plq=Table[{},{4}];\n   \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4}}=xyc; {f1,f2,f3,f4}=fc;\n   \
x5=(x1+x2+x3+x4)/4; y5=(y1+y2+y3+y4)/4; f5=(f1+f2+f3+f4)/4;\n   \
xytab={{{x1,y1},{x2,y2},{x5,y5}},{{x2,y2},{x3,y3},{x5,y5}},\n          \
{{x3,y3},{x4,y4},{x5,y5}},{{x4,y4},{x1,y1},{x5,y5}}};\n   \
ftab={{f1,f2,f5},{f2,f3,f5},{f3,f4,f5},{f4,f1,f5}};\n   For [e=1,e<=4,e++, \
xyt=xytab[[e]]; ft=ftab[[e]];\n       \
{pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnTrig3[\n            \
xyt,ft,{fmin,fmax,finc,eps},lines];\n       ];       \n    Return[{pbq,plq}];\
\n\t];\n\t\nContourBandPlotFunctionOnQuad9[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\
\n   lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,y1,y2,y3,y4,\n   \
y5,y6,y7,y8,y9,f1,f2,f3,f4,f5,f6,f7,f8,f9,\n   e,xytab,ftab,xyt,ft, \
pbq,plq,pbt,plt },\n   pbq=plq=Table[{},{4}];\n   \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},\n    \
{x5,y5},{x6,y6},{x7,y7},{x8,y8},{x9,y9}}=xyc; \n    \
{f1,f2,f3,f4,f5,f6,f7,f8,f9}=fc;\n   \
xytab={{{x1,y1},{x5,y5},{x9,y9},{x8,y8}},\n          \
{{x2,y2},{x6,y6},{x9,y9},{x5,y5}},\n          \
{{x3,y3},{x7,y7},{x9,y9},{x6,y6}},\n          \
{{x4,y4},{x8,y8},{x9,y9},{x7,y7}}};\n   \
ftab={{f1,f5,f9,f8},{f2,f6,f9,f5},{f3,f7,f9,f6},{f4,f8,f9,f7}};\n   For \
[e=1,e<=4,e++, xye=xytab[[e]]; fe=ftab[[e]];\n       \
{pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnQuad4[\n             \
xye,fe,{fmin,fmax,finc,eps},lines];\n       ];       \n    Return[{pbq,plq}];\
\n\t];\n\t\nContourBandPlotFunctionOnQuad8[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\
\n   lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,y1,y2,y3,y4,\n   \
y5,y6,y7,y8,y9,f1,f2,f3,f4,f5,f6,f7,f8,f9,\n   e,xytab,ftab,xyt,ft, \
pbq,plq,pbt,plt },\n   pbq=plq=Table[{},{4}];\n   \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},\n    {x5,y5},{x6,y6},{x7,y7},{x8,y8}}=xyc; \
\n    {f1,f2,f3,f4,f5,f6,f7,f8}=fc;\n   x9=(2*(x5+x6+x7+x8)-(x1+x2+x3+x4))/4;\
\n   y9=(2*(y5+y6+y7+y8)-(y1+y2+y3+y4))/4;\n   \
f9=(2*(f5+f6+f7+f8)-(f1+f2+f3+f4))/4;\n   \
xytab={{{x1,y1},{x5,y5},{x9,y9},{x8,y8}},\n          \
{{x2,y2},{x6,y6},{x9,y9},{x5,y5}},\n          \
{{x3,y3},{x7,y7},{x9,y9},{x6,y6}},\n          \
{{x4,y4},{x8,y8},{x9,y9},{x7,y7}}};\n   \
ftab={{f1,f5,f9,f8},{f2,f6,f9,f5},{f3,f7,f9,f6},{f4,f8,f9,f7}};\n   For \
[e=1,e<=4,e++, xye=xytab[[e]]; fe=ftab[[e]];\n       \
{pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnQuad4[\n                   \
xye,fe,{fmin,fmax,finc,eps},lines];\n       ];       \n    Return[{pbq,plq}];\
\n\t];\n\t\n\
ContourBandPlotFunctionOnQuad16[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\n   \
lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,\n   \
x11,x12,x13,x14,x15,x16,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,\n   \
y11,y12,y13,y14,y15,y16,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,\n   \
f11,f12,f13,f14,f15,f16,\n   e,xytab,ftab,xyt,ft, pbq,plq,pbt,plt },\n   \
pbq=plq=Table[{},{9}];\n   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},\n    \
{x5,y5},{x6,y6},{x7,y7},{x8,y8},\n    {x9,y9},{x10,y10},{x11,y11},{x12,y12},\n\
    {x13,y13},{x14,y14},{x15,y15},{x16,y16}}=xyc; \n   \
{f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,\n    f13,f14,f15,f16}=fc;\n   \
xytab={{{x1,y1},  {x5,y5},  {x13,y13},{x12,y12}},\n          {{x2,y2},  \
{x7,y7},  {x14,y14},{x6,y6}},\n          {{x3,y3},  {x9,y9},  \
{x15,y15},{x8,y8}},\n          {{x4,y4},  {x11,y11},{x16,y16},{x10,y10}},\n   \
       {{x5,y5},  {x6,y6},  {x14,y14},{x13,y13}},\n          {{x7,y7},  \
{x8,y8},  {x15,y15},{x14,y14}},\n          {{x9,y9},  \
{x10,y10},{x16,y16},{x15,y15}},\n          \
{{x11,y11},{x12,y12},{x13,y13},{x16,y16}},\n          \
{{x13,y13},{x14,y14},{x15,y15},{x16,y16}}};\n   \
ftab={{f1,f5,f13,f12},{f2,f7,f14,f6},{f3,f9,f15,f8},\n        \
{f4,f11,f16,f10},{f5,f6,f14,f13},{f7,f8,f15,f14},\n        \
{f9,f10,f16,f15},{f11,f12,f13,f16},{f13,f14,f15,f16}};\n   For [e=1,e<=9,e++, \
xye=xytab[[e]]; fe=ftab[[e]];\n       \
{pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnQuad4[\n                \
xye,fe,{fmin,fmax,finc,eps},lines];\n       ];       \n    Return[{pbq,plq}];\
\n\t];\n\t\nContourBandPlotFunctionOnTrig6[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\
\n   lines_]:=Module[{x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,\n   \
f1,f2,f3,f4,f5,f6,e,xytab,ftab,xyt,ft,pbt,plt,pbe,ple},\n   \
pbt=plt=Table[{},{4}];\n   \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6}}=xyc; \n   \
{f1,f2,f3,f4,f5,f6}=fc;\n   \
xytab={{{x1,y1},{x4,y4},{x6,y6}},{{x2,y2},{x5,y5},{x4,y4}},\n          \
{{x3,y3},{x6,y6},{x5,y5}},{{x4,y4},{x5,y5},{x6,y6}}};\n   \
ftab={{f1,f4,f6},{f2,f5,f4},{f3,f6,f5},{f4,f5,f6}};\n   For [e=1,e<=4,e++, \
xyt=xytab[[e]]; ft=ftab[[e]];\n       \
{pbt[[e]],plt[[e]]}=ContourBandPlotFunctionOnTrig3[\n             \
xyt,ft,{fmin,fmax,finc,eps},lines];\n       ];       \n    Return[{pbt,plt}];\
\n\t];\n\t\n\
ContourBandPlotFunctionOnTrig10[xyc_,fc_,{fmin_,fmax_,finc_,eps_},\n   \
lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x0,y1,y2,y3,\n   \
y4,y5,y6,y7,y8,y9,y0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f0,\n   \
e,xytab,ftab,xyt,ft,pbt,plt},\n   pbt=plt=Table[{},{9}];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6},\n    \
{x7,y7},{x8,y8},{x9,y9},{x0,y0}}=xyc; \n   \
{f1,f2,f3,f4,f5,f6,f7,f8,f9,f0}=fc;\n   \
xytab={{{x1,y1},{x4,y4},{x9,y9}},{{x2,y2},{x6,y6},{x5,y5}},\n          \
{{x3,y3},{x8,y8},{x7,y7}},{{x9,y9},{x4,y4},{x0,y0}},\n          \
{{x4,y4},{x5,y5},{x0,y0}},{{x5,y5},{x6,y6},{x0,y0}},\n          \
{{x6,y6},{x7,y7},{x0,y0}},{{x7,y7},{x8,y8},{x0,y0}},\n          \
{{x8,y8},{x9,y9},{x0,y0}}};\n   \
ftab={{f1,f4,f9},{f2,f6,f5},{f3,f8,f7},{f9,f4,f0},\n         \
{f4,f5,f0},{f5,f6,f0},{f6,f7,f0},{f7,f8,f0},\n         {f8,f9,f0}};\n   For \
[e=1,e<=9,e++, xyt=xytab[[e]]; ft=ftab[[e]];\n       \
{pbt[[e]],plt[[e]]}=ContourBandPlotFunctionOnTrig3[\n            \
xyt,ft,{fmin,fmax,finc,eps},lines];\n       ];     \n    Return[{pbt,plt}];\n\
\t];\n\t\nContourPlotValueLegend[fmin_,fmax_,xyleg_,mv_]:=Module[\n  \
{f,fs,finc,fa=Abs[fmin],fb=Abs[fmax],xylinL,xylinR,xytext,\n   \
black=Graphics[RGBColor[0,0,0]],baselineskip=10,\n   \
white=Graphics[RGBColor[1,1,1]],preleg,plab,iv,\n   \
thick=Graphics[AbsoluteThickness[10]],xf1,xf2,xf3,xf4,\n   thin= \
Graphics[AbsoluteThickness[1]],color,padOK},\n   xf1=Offset[{-30,  \
10},xyleg]; xf2=Offset[{85,10},xyleg];\n   xf3=Offset[{ \
85,-10-baselineskip*(mv+1)},xyleg];\n   \
xf4=Offset[{-30,-10-baselineskip*(mv+1)},xyleg];\n   preleg={thin,white,\n    \
   Graphics[Polygon[{xf1,xf2,xf3,xf4}]], black,\n       \
Graphics[Line[{xf1,xf2,xf3,xf4,xf1}]],thick,\n       \
Graphics[Text[\"LEGEND\",Offset[{5,0},xyleg],{-1,0}]]};\n   \
finc=(fmax-fmin)/mv; plab=Table[{},{mv+1}];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   padOK=(Max[fa,fb]<=1000)&&(Min[fa,fb]>=1/100);\n   For \
[iv=1,iv<=mv+1,iv++, f=Chop[N[fmin+(iv-1)*finc]]; \n        \
xylinL=Offset[{-20,-baselineskip*iv},xyleg];\n        xylinR=Offset[{  \
5,-baselineskip*iv},xyleg]; \n        xytext=Offset[{ \
80,-baselineskip*iv},xyleg];\n        If [padOK, \
fs=PaddedForm[f,4,NumberSigns->{\"-\",\"+\"}],\n               \
fs=ScientificForm[f,4,NumberSigns->{\"-\",\"+\"}]];\n        \
color=ContourBandColor[f,fmin,fmax]; \n        \
plab[[iv]]={color,Graphics[Line[{xylinL,xylinR}]],\n                    \
black,Graphics[Text[fs,xytext,{1,0}]]}]; \n   Return[{preleg,plab}]];\n \n \
ContourBandPlotNodeFuncOver2DMesh[nodxyz_,elenod_,fn_,{fmin_,fmax_,finc_},\n  \
{colors_,mesh_,nodes_,bacg_,lines_,bound_},xyleg_,aspect_,label_]:= Module[\n \
 {enl,i,k,n,ne,xyc,fc,c,nn,pb,pl,eps,prebacg,pbacg,preband,pband,preline,\n  \
bnt,bn,nbn,enb,xyb,bseg,ebs,seg,rseg,nbseg,prebound,pbound,\n  \
pline,premesh,pmesh,prenode,preleg,pleg,pnode,numele=Length[elenod]},\n  \
bnt={{},{},{1,2,3,1},{1,2,3,4,1},{},{1,4,2,5,3,6,1},\n      \
{},{1,5,2,6,3,7,4,8,1},{1,5,2,6,3,7,4,8,1},\n      \
{1,4,5,2,6,7,3,8,9,1},{},{},{},{},{},{1,5,6,2,7,8,3,9,10,4,11,12,1}}; \n  \
eps=10.^(-8)*Max[Abs[fmax],Abs[fmin]]; pband=pmesh=pnode=pline=pbacg={}; \n  \
If [colors,pband=Table[{},{numele}]]; If [mesh, pmesh=Table[{},{numele}]]; \n \
 If [nodes, pnode=Table[{},{numele}]]; If [lines,pline=Table[{},{numele}]]; \n\
  If [bacg,  pbacg=Table[{},{numele}]]; preleg=pleg=bseg=pbound={};\n  \
prebacg= {Graphics[RGBColor[.7,.7,.7]]}; preband={};\n  preline= \
{Graphics[RGBColor[0,0,0]],Graphics[AbsoluteThickness[1]]};\n  premesh= \
{Graphics[RGBColor[1,0,0]],Graphics[AbsoluteThickness[2]],\n            \
Graphics[AbsoluteDashing[{1,10}]]};\n  prenode= \
{Graphics[RGBColor[0,0,0]],Graphics[AbsoluteDashing[{}]],\n            \
Graphics[AbsolutePointSize[5]]};\n  \
prebound={Graphics[RGBColor[0,0,1]],Graphics[AbsoluteDashing[{}]],\n          \
  Graphics[AbsoluteThickness[2]]};\n  For [e=1,e<=Length[elenod],e++, \n      \
enl=elenod[[e]]; ne=Length[enl];  \n      If [!MemberQ[{3,4,6,8,9,10,16},ne], \
Print[\"ContourBandPlot:\",\n           \"Elem \",e,\" with \",ne,\" nodes \
skipped\"]; Continue[]];\n      fc=Table[fn[[enl[[i]]]],{i,ne}]; \
xyc=Table[nodxyz[[enl[[i]]]],{i,ne}];\n      bn=bnt[[ne]];  nbn=Length[bn]; \
xyb=Table[xyc[[bn[[i]]]],{i,nbn}];\n      If [bound, \
enb=Table[enl[[bn[[i]]]],{i,nbn}]; \n          \
ebs=Table[{enb[[i]],enb[[i+1]]},{i,1,nbn-1}]; ",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "          For [i=1,i<=nbn-1,i++, seg=ebs[[i]]; rseg=Reverse[seg]; \n       \
       If [MemberQ[bseg,rseg], {{k}}=Position[bseg,rseg,1,1]; \n              \
bseg=Drop[bseg,{k}]; Continue[]]; \n              AppendTo[bseg,seg]];\n      \
   ];\n      If [ne==3, {pb,pl}=ContourBandPlotFunctionOnTrig3[\n             \
     xyc,fc,{fmin,fmax,finc,eps},lines]];\n      If [ne==4, \
{pb,pl}=ContourBandPlotFunctionOnQuad4[\n                  \
xyc,fc,{fmin,fmax,finc,eps},lines]];\n      If [ne==6, \
{pb,pl}=ContourBandPlotFunctionOnTrig6[\n                  \
xyc,fc,{fmin,fmax,finc,eps},lines]];\n      If [ne==8, \
{pb,pl}=ContourBandPlotFunctionOnQuad8[\n                  \
xyc,fc,{fmin,fmax,finc,eps},lines]];\n      If [ne==9, \
{pb,pl}=ContourBandPlotFunctionOnQuad9[\n                  \
xyc,fc,{fmin,fmax,finc,eps},lines]];\n      If \
[ne==10,{pb,pl}=ContourBandPlotFunctionOnTrig10[\n                  \
xyc,fc,{fmin,fmax,finc,eps},lines]];\n      If \
[ne==16,{pb,pl}=ContourBandPlotFunctionOnQuad16[\n                  \
xyc,fc,{fmin,fmax,finc,eps},lines]]; \n      If [colors, pband[[e]]=pb]; If \
[lines, pline[[e]]=pl];\n      If [mesh,   pmesh[[e]]=Graphics[Line[xyb]]]; \n\
      If [bacg,   pbacg[[e]]=Graphics[Polygon[xyb]]]; \n      If [nodes,  \
pnode[[e]]=Table[Graphics[Point[xyb[[i]] ]],{i,ne}]];  ",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      ];\n If [!colors,pband={}]; If [!mesh, pmesh={}]; If [!nodes, \
pnode={}]; \n If [!lines, pline={}]; If [!bacg, pbacg={}]; pleg={};",
  StyleBox["\n ",
    FontColor->RGBColor[1, 0, 0]],
  "nbseg=Length[bseg];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n If [nbseg>0, pbound=Table[0,{nbseg}];  \n     For [k=1,k<=nbseg,k++, \
{i,j}=bseg[[k]];\n         \
pbound[[k]]=Graphics[Line[{nodxyz[[i]],nodxyz[[j]]}]] ]];\n If \
[Length[xyleg]==0||!colors, Show[prebacg,pbacg,\n     \
preband,pband,preline,pline,premesh,pmesh,prenode,\n     \
pnode,prebound,pbound,TextStyle->{FontFamily->\"Courier\",\n     \
FontSlant->\"Plain\", FontSize->10},\n     \
AspectRatio->aspect,PlotLabel->label,      \n     \
DisplayFunction->DisplayChannel ]];\n If [Length[xyleg]>0&&colors,\n     \
{preleg,pleg}=ContourPlotValueLegend[fmin,fmax,xyleg,6];\n     \
Show[prebacg,pbacg,preband,pband,preline,pline,\n          \
premesh,pmesh,prenode,pnode,prebound,pbound,preleg,pleg, \n          \
TextStyle->{FontFamily->\"Courier\",\n                  FontSlant->\"Plain\", \
FontSize->9},\n          AspectRatio->aspect,PlotLabel->label,\n          \
DisplayFunction->DisplayChannel ]];\n \
ClearAll[pband,pline,pmesh,pnode,pbacg,pbound,pleg]; \n ];\n  \n",
  StyleBox["(*nodxyz=N[{{0,6},{0,3},{0,0},{5/2,6},{5/2,3},\n     \
{5/2,0},{5,6},{5,3},{5,0}}];        \nelenod=  {{1,3,9,7,2,6,8,4,5}};\n\
\[Sigma]rr=Table[8*nodxyz[[n,1]]/5-4,{n,1,9}]; \[Sigma]rrmax=4; \
\[Sigma]rrmin=-4;\n \[Sigma]rrinc=.5; aspect=6/5; Print[\"\[Sigma]rr=\",\
\[Sigma]rr];\nContourBandPlotNodeFuncOver2DMesh[nodxyz,elenod,\[Sigma]rr,{\
\[Sigma]rrmin,\[Sigma]rrmax,\[Sigma]rr\[Sigma]rrinc},\n  \
{True,True,True,True,True,True},{2,2.5},aspect,\"Sigma-xx (white: zero)\"];\n\
\[Sigma]zz=Table[4*nodxyz[[n,2]]/6,{n,1,9}]; \[Sigma]zzmax=4; \
\[Sigma]zzmin=-4;\n\[Sigma]zzinc=.5; aspect=6/5; Print[\"\[Sigma]zz=\",\
\[Sigma]zz];\nContourBandPlotNodeFuncOver2DMesh[nodxyz,elenod,\[Sigma]zz,{\
\[Sigma]zzmin,\[Sigma]zzmax,\[Sigma]zzinc},\n  \
{True,True,True,True,True,True},{2,2.5},aspect,\"Sigma-yy (white: \
zero)\"];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["Cell 9A: Print utility modules", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
PrintRingNodeCoordinates[nodrzc_,title_,digits_]:= Module[
  {numnod=Length[nodrzc],n,x,y,z,d=2,f=4,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++,{x,y}=nodrzc[[n]]; 
       tab[[n]]={ToString[n],PaddedForm[x,{d,f}],PaddedForm[y,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->{Right},   
      TableDirections->{Column,Row},TableSpacing->{0,2},
      TableHeadings ->{None,{\"node\", \"r-coor\", \"z-coor\"}}]];
   ClearAll[tab];
   ];
   
PrintRingElementTypeNodes[eletyp_,elenod_,title_,digits_]:=
  Module[{numele=Length[elenod],e,d=6,f=3,tab}, 
   tab=Table[0,{numele}]; If [Length[digits]==2,{d,f}=digits];
   For [e=1,e<=numele,e++, 
        tab[[e]]={ToString[e],eletyp[[e]],ToString[elenod[[e]]]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->Right,
         TableDirections->{Column,Row},TableSpacing->{0,2},
         TableHeadings->{None,{\"elem\", \"type\", \"node-list\"}}]];
   ClearAll[tab];
   ];
   
PrintRingElementMaterials[elemat_,title_,digits_]:= 
  Module[{numele=Length[elemat],e,mat,d=6,f=3,tab}, 
   tab=Table[0,{numele}]; If [Length[digits]==2,{d,f}=digits];
   For [e=1,e<=numele,e++, mat=elemat[[e]]; 
        tab[[e]]={ToString[e],ToString[mat]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->Right,
         TableDirections->{Column,Row},TableSpacing->{0,2},
         TableHeadings->{None,{\"elem\", \"material: {E,\[Nu]}\"}}]];
   ClearAll[tab];
   ];
   
PrintRingElemBodyForces[elebfor_,title_,digits_]:= 
  Module[{numele=Length[elebfor],e,br,bz,
   br1,bz1,br2,bz2,br3,bz3,br4,bz4,tab,d=8,f=4}, 
   If [numele<=0,Print[\"No body forces specified\"]; Return[]];
   tab=Table[0,{numele}]; If [Length[digits]==2,{d,f}=digits];
   If [StringLength[title]>0, Print[title]];
   If [Length[elebfor[[1]]]==2,
   For [e=1,e<=numele,e++, {br,bz}=elebfor[[e]]; 
       tab[[e]]={ToString[e],PaddedForm[br,{d,f}],
                             PaddedForm[bz,{d,f}]}];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings ->{None,{\"elem\", \"r-bforce\", \"z-bforce\"}}]]
       ];  
   If [Length[elebfor[[1]]]==4,
   For [e=1,e<=numele,e++, 
       {{br1,bz1},{br2,bz2},{br3,bz3},{br4,bz4}}=elebfor[[e]]; 
        tab[[e]]={ToString[e],
        PaddedForm[br1,{d,f}],PaddedForm[bz1,{d,f}],
        PaddedForm[br2,{d,f}],PaddedForm[bz2,{d,f}],
        PaddedForm[br3,{d,f}],PaddedForm[bz3,{d,f}],
        PaddedForm[br4,{d,f}],PaddedForm[bz4,{d,f}]}];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings ->{None,{\"elem\", \"br1\",\"bz1\",\"br2\",\"bz2\",
         \"br3\",\"bz3\",\"br4\",\"bz4\"}}]];
       ];      
   ClearAll[tab];
   ];
   
PrintRingElemTracForces[eletfor_,title_,digits_]:= 
  Module[{numele=Length[eletfor],e,tab,
  ne,side,tr1,tr2,tz1,tz2,d=8,f=4}, 
  If [numele<=0,Print[\"No surface traction forces specified\"]; 
      Return[]];
   tab=Table[0,{numele}]; If [Length[digits]==2,{d,f}=digits];
   For [e=1,e<=numele,e++, {ne,side,{tr1,tr2},{tz1,tz2}}=eletfor[[e]]; 
       tab[[e]]={ToString[ne],ToString[side],
                 PaddedForm[tr1,{d,f}],PaddedForm[tr2,{d,f}],
                 PaddedForm[tz1,{d,f}],PaddedForm[tz2,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings ->{None,{\"elem\", \"side\", \"r-tforce1\", \
\"r-tforce2\",
         \"z-tforce1\", \"z-tforce2\"}}]];         
   ClearAll[tab];
   ];
   
PrintRingFreedomActivity[nodtag_,nodval_,title_,digits_]:= Module[
  {numnod=Length[nodtag],n,t,v,d=6,f=2,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++, t=nodtag[[n]]; v=nodval[[n]];
        tab[[n]]={ToString[n],
        PaddedForm[t[[1]],d],PaddedForm[t[[2]],d],
        PaddedForm[v[[1]],{d,f}],PaddedForm[v[[2]],{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings->{None,{\"node\", \"r-tag\", \"z-tag\",
                                \"r-value\", \"z-value\"}} ]];
   ClearAll[tab];
   ];

PrintRingNodeDisplacements[noddis_,title_,digits_]:= Module[
  {numnod=Length[noddis],n,ur,uz,d=8,f=6,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits]; 
   For [n=1,n<=numnod,n++, {ur,uz}=noddis[[n]];
       tab[[n]]={ToString[n],PaddedForm[ur,{d,f}],
                  PaddedForm[uz,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings->{None,{\"node\", \"r-disp\", \"z-disp\"}}]];
   ClearAll[tab];
   ];
   
PrintRingNodeForces[nodfor_,title_,digits_]:= Module[
  {numnod=Length[nodfor],n,fr,fz,tab,d=8,f=4}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++, {fr,fz}=nodfor[[n]]; 
        tab[[n]]={ToString[n],
         PaddedForm[fr,{d,f}],PaddedForm[fz,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings ->{None,{\"node\", \"r-force\", \"z-force\"}}]];       \
  
   ClearAll[tab];
   ]; 
   
PrintRingNodeStresses[nodsig_,title_,digits_]:= 
  Module[{numnod=Length[nodsig],n,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\
\[Theta],\[Sigma]rz,d=8,f=4,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++, {\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\
\[Sigma]rz}=nodsig[[n]];
         tab[[n]]={ToString[n], 
            PaddedForm[\[Sigma]rr,{d,f}], PaddedForm[\[Sigma]zz,{d,f}],
            PaddedForm[\[Sigma]\[Theta]\[Theta],{d,f}], \
PaddedForm[\[Sigma]rz,{d,f}]} ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{\"node\", \"sigma-rr\", 
          \"sigma-zz\", \"sigma-\[Theta]\[Theta]\", \"sigma-xy\"}}]];
   ClearAll[tab];
   ]; 
   
PrintRingAnalysisSolution[noddis_,nodfor_,nodsig_,title_,digits_]:= 
  Module[{numnod=Length[noddis],n,ur,uz,fr,fz,
          \[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz,d=6,f=4,\
tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
  For [n=1,n<=numnod,n++, {ur,uz}=noddis[[n]]; {fr,fz}=nodfor[[n]]; 
        {\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}=nodsig[[n]\
];
        tab[[n]]={ToString[n], 
            PaddedForm[ur,{d,f}],  PaddedForm[uz,{d,f}],
            PaddedForm[\[Sigma]rr,{d,f}], PaddedForm[\[Sigma]zz,{d,f}],
            PaddedForm[\[Sigma]\[Theta]\[Theta],{d,f}], \
PaddedForm[\[Sigma]rz,{d,f}],
            PaddedForm[fr,{d,f}],  PaddedForm[fz,{d,f}]} ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{\"node\",\"r-disp\",\"z-disp\",\"sigma-rr\",
         \"sigma-zz\",\"sigma-\[Theta]\[Theta]\",\"sigma-xy\",\"r-force\",\"z-\
force\"}}]];
   ClearAll[tab];
   ];  
   
PrintRingNodeDispStresses[noddis_,nodsig_,title_,digits_]:= 
  Module[{numnod=Length[noddis],n,ur,uz,
          \[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz,d=6,f=4,\
tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
  For [n=1,n<=numnod,n++, {ur,uz}=noddis[[n]]; 
        {\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}=nodsig[[n]\
];
        tab[[n]]={ToString[n], 
            PaddedForm[ur,{d,f}],  PaddedForm[uz,{d,f}], 
            PaddedForm[\[Sigma]rr,{d,f}], PaddedForm[\[Sigma]zz,{d,f}],
            PaddedForm[\[Sigma]\[Theta]\[Theta],{d,f}], \
PaddedForm[\[Sigma]rz,{d,f}]} ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{\"node\",\"r-disp\",\"z-disp\",
        \"sigma-rr\",\"sigma-zz\",\"sigma-\[Theta]\[Theta]\",\"sigma-xy\"}}]];\

   ClearAll[tab];
   ]; \
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 9B: Mesh generation utility modules\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "GenQuad4NodeCoordinates[rzcorn_,Ner_,Nez_]:= Module[\n   \
{k=0,i,j,dr,dz,numnod,N1,N2,N3,N4,\[Xi],\[Eta],\n   \
rc1,rc2,rc3,rc4,yc1,yc2,yc3,yc4,r,z,nodrzc}, \n   numnod=(Ner+1)*(Nez+1); \
nodrzc=Table[{0,0},{numnod}];\n   \
{{rc1,zc1},{rc2,zc2},{rc3,zc3},{rc4,zc4}}=rzcorn;\n   For [i=1,i<=Ner+1,i++, \
For [j=1,j<=Nez+1,j++,\n      \[Xi]=N[2*(i-1)/Ner-1]; \
\[Eta]=N[2*(j-1)/Nez-1];\n      N1=(1-\[Xi])*(1-\[Eta])/4; N2=(1+\[Xi])*(1-\
\[Eta])/4;\n      N3=(1+\[Xi])*(1+\[Eta])/4; N4=(1-\[Xi])*(1+\[Eta])/4;\n     \
 r=rc1*N1+rc2*N2+rc3*N3+rc4*N4;\n      z=zc1*N1+zc2*N2+zc3*N3+zc4*N4;\n      \
nodrzc[[++k]]={r,z};\n      ]];\n   Return[nodrzc]];\n       \n\
GenQuad8NodeCoordinates[rzcorn_,Ner_,Nez_]:= Module[\n   \
{k=0,i,j,dr,dz,numnod,N1,N2,N3,N4,\[Xi],\[Eta],\n   \
rc1,rc2,rc3,rc4,yc1,yc2,yc3,yc4,r,z,nodrzc}, \n   \
numnod=(2*Ner+1)*(2*Nez+1)-Ner*Nez; \n   nodrzc=Table[{0,0},{numnod}];\n   \
{{rc1,zc1},{rc2,zc2},{rc3,zc3},{rc4,zc4}}=rzcorn;\n   For \
[i=1,i<=2*Ner+1,i++, For [j=1,j<=2*Nez+1,j++,\n      If \
[Mod[i,2]==0&&Mod[j,2]==0,Continue[]];\n      \[Xi]=N[(i-Ner-1)/Ner]; \
\[Eta]=N[(j-Nez-1)/Nez];\n      N1=(1-\[Xi])*(1-\[Eta])/4; N2=(1+\[Xi])*(1-\
\[Eta])/4;\n      N3=(1+\[Xi])*(1+\[Eta])/4; N4=(1-\[Xi])*(1+\[Eta])/4;\n     \
 r=rc1*N1+rc2*N2+rc3*N3+rc4*N4;\n      z=zc1*N1+zc2*N2+zc3*N3+zc4*N4; \n      \
If [Mod[i,2]==0&&Mod[j,2]==0,Continue[]];\n      nodrzc[[++k]]={r,z};\n      \
]];\n   Return[nodrzc]]; \n   \n\
GenGradedQuad4NodeCoordinates[rzcorn_,Ner_,Nez_,sr_,sz_]:= \n\
Module[{k=0,i,j,dr,dz,numnod,N1,N2,N3,N4,\[Xi],\[Eta],\n   \
rc1,rc2,rc3,rc4,yc1,yc2,yc3,yc4,r,z,elenod}, \n   numnod=(Ner+1)*(Nez+1); \
rzcoor=Table[0,{numnod}];\n   \
{{rc1,zc1},{rc2,zc2},{rc3,zc3},{rc4,zc4}}=rzcorn;\n   For [i=1,i<=Ner+1,i++,\n\
       If [i<=Length[sr],\[Xi]=sr[[i]],\[Xi]=2*(i-1)/Ner-1];\n       For \
[j=1,j<=Nez+1,j++,\n           If \
[j<=Length[sz],\[Eta]=sz[[j]],\[Eta]=2*(j-1)/Nez-1];\n      N1=(1-\[Xi])*(1-\
\[Eta])/4; N2=(1+\[Xi])*(1-\[Eta])/4;\n      N3=(1+\[Xi])*(1+\[Eta])/4; \
N4=(1-\[Xi])*(1+\[Eta])/4;\n      r=rc1*N1+rc2*N2+rc3*N3+rc4*N4;\n      \
z=zc1*N1+zc2*N2+zc3*N3+zc4*N4;\n      rzcoor[[++k]]={r,z};\n      ]];\n   \
Return[rzcoor]];\n\nGenQuad4ElemNodes[Ner_,Nez_]:= Module[\n   \
{k,i,j,c1,c2,numele,elenod},numele=Ner*Nez; \n   \
elenod=Table[{0,0,0,0},{numele}]; k=0;\n   For [i=1,i<=Ner,i++, For \
[j=1,j<=Nez,j++,\n       c1=(Nez+1)*(i-1)+j; c2=c1+Nez+1; \n       \
elenod[[++k]]={c1,c2,c2+1,c1+1}\n       ]]; Return[elenod]];\n       \n\
GenQuad8ElemNodes[Ner_,Nez_]:= Module[\n  {k,i,j,c1,c2,c3,numele,elenod}, \
numele=Ner*Nez; \n   elenod=Table[{0,0,0,0,0,0,0,0},{numele}]; k=0;\n   For \
[i=1,i<=Ner,i++, For [j=1,j<=Nez,j++,\n       c1=(3*Nez+2)*(i-1)+2*j-1;  \n   \
    c3=c1+2*Nez-j+2; c2=c3+Nez+j;\n       \
elenod[[++k]]={c1,c2,c2+2,c1+2,c3,c2+1,c3+1,c1+1}\n       ]]; \
Return[elenod]]; \n       ",
  StyleBox["        \n(* Ner=2;Nez=2;  rzcorn={{0,0},{1,0},{1,1},{0,1}};\n\
ElemNodes=GenQuad4ElemNodes[Ner,Nez];\n\
NodeCoordinates=GenQuad4NodeCoordinates[rzcorn,Ner,Nez];\n\
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,1,\n    \"Quad4 \
mesh\",True,True];\nElemNodes=GenQuad8ElemNodes[Ner,Nez];\n\
NodeCoordinates=GenQuad8NodeCoordinates[rzcorn,Ner,Nez];\n\
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,1,\n    \"Quad8 \
mesh\",True,True]; *)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["Cell 9C.  Stress contour plot driver", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ContourPlotStresses[nodrzc_,elenod_,nodsig_,whichones_,\n  \
showrange_,plotopt_,legend_,aspect_]:=Module[\n  {sigeps=10.^(-3),nbands=10,\
\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz,n,\n  \
numnod=Length[nodrzc],\n  options={True,True,True,False,True,True},\n  \
\[Sigma]rrmax,\[Sigma]zzmax,\[Sigma]\[Theta]\[Theta]max,\[Sigma]rzmax,\[Sigma]\
rrmin,\[Sigma]zzmin,\[Sigma]\[Theta]\[Theta]min,\[Sigma]rzmin,\n  \
\[Sigma]rrinc,\[Sigma]zzinc,\[Sigma]\[Theta]\[Theta]inc,\[Sigma]rzinc,plotrr,\
plotzz,plot\[Theta]\[Theta],plotrz},\n  \
\[Sigma]rr=Table[nodsig[[n,1]],{n,numnod}]; \n  \
\[Sigma]zz=Table[nodsig[[n,2]],{n,numnod}];\n  \
\[Sigma]\[Theta]\[Theta]=Table[nodsig[[n,3]],{n,numnod}];\n  \
\[Sigma]rz=Table[nodsig[[n,4]],{n,numnod}];\n  {\[Sigma]rrmax,\[Sigma]zzmax,\
\[Sigma]\[Theta]\[Theta]max,\[Sigma]rzmax}=Abs[{Max[\[Sigma]rr],\n       Max[\
\[Sigma]zz],Max[\[Sigma]\[Theta]\[Theta]],Max[\[Sigma]rz]}]+sigeps; \n  {\
\[Sigma]rrmin,\[Sigma]zzmin,\[Sigma]\[Theta]\[Theta]min,\[Sigma]rzmin}=Abs[{\
Min[\[Sigma]rr],\n       Min[\[Sigma]zz],Min[\[Sigma]\[Theta]\[Theta]],Min[\
\[Sigma]rz]}]+sigeps;\n   \[Sigma]rrmax=Max[\[Sigma]rrmax,\[Sigma]rrmin];  \
\[Sigma]rrmin=-\[Sigma]rrmax;\n   \
\[Sigma]zzmax=Max[\[Sigma]zzmax,\[Sigma]zzmin];  \
\[Sigma]zzmin=-\[Sigma]zzmax;\n   \[Sigma]\[Theta]\[Theta]max=Max[\[Sigma]\
\[Theta]\[Theta]max,\[Sigma]\[Theta]\[Theta]min];  \
\[Sigma]\[Theta]\[Theta]min=-\[Sigma]\[Theta]\[Theta]max;\n   \
\[Sigma]rzmax=Max[\[Sigma]rzmax,\[Sigma]rzmin];  \
\[Sigma]rzmin=-\[Sigma]rzmax;\n  \
{\[Sigma]rrinc,\[Sigma]zzinc,\[Sigma]\[Theta]\[Theta]inc,\[Sigma]rzinc}={\
\[Sigma]rrmax-\[Sigma]rrmin,\[Sigma]zzmax-\[Sigma]zzmin,\n      \[Sigma]\
\[Theta]\[Theta]max-\[Sigma]\[Theta]\[Theta]min,\[Sigma]rzmax-\[Sigma]rzmin}/\
nbands;\n  {plotrr,plotzz,plot\[Theta]\[Theta],plotrz}=whichones;\n  If \
[Length[plotopt]>0,options=plotopt];\n  If [plotrr, If [showrange, Print[\"\
\[Sigma]rr min,max,inc=\",\n                  {\[Sigma]rrmin,\[Sigma]rrmax,\
\[Sigma]rrinc}]];\n     ContourBandPlotNodeFuncOver2DMesh[nodrzc,\n       \
elenod,\[Sigma]rr,{\[Sigma]rrmin,\[Sigma]rrmax,\[Sigma]rrinc},options,\n      \
 legend,aspect,\"Radial stress \[Sigma]rr\"]];\n  If [plotzz, If [showrange, \
Print[\"\[Sigma]zz min,max,inc=\",\n                  \
{\[Sigma]zzmin,\[Sigma]zzmax,\[Sigma]zzinc}]];\n     \
ContourBandPlotNodeFuncOver2DMesh[nodrzc,\n       \
elenod,\[Sigma]zz,{\[Sigma]zzmin,\[Sigma]zzmax,\[Sigma]zzinc},options,\n      \
 legend,aspect,\"Axial stress \[Sigma]zz\"]];\n  If [plot\[Theta]\[Theta], If \
[showrange, Print[\"\[Sigma]\[Theta]\[Theta] min,max,inc=\",\n                \
  {\[Sigma]\[Theta]\[Theta]min,\[Sigma]\[Theta]\[Theta]max,\[Sigma]\[Theta]\
\[Theta]inc}]];\n     ContourBandPlotNodeFuncOver2DMesh[nodrzc,\n       \
elenod,\[Sigma]\[Theta]\[Theta],{\[Sigma]\[Theta]\[Theta]min,\[Sigma]\[Theta]\
\[Theta]max,\[Sigma]\[Theta]\[Theta]inc},options,\n       \
legend,aspect,\"Hoop stress \[Sigma]\[Theta]\[Theta]\"]];\n  If [plotrz, If \
[showrange, Print[\"\[Sigma]rz min,max,inc=\",\n                  \
{\[Sigma]rzmin,\[Sigma]rzmax,\[Sigma]rzinc}]];\n     \
ContourBandPlotNodeFuncOver2DMesh[nodrzc,\n       \
elenod,\[Sigma]rz,{\[Sigma]rzmin,\[Sigma]rzmax,\[Sigma]rzinc},options,\n      \
 legend,aspect,\"Shear stress \[Sigma]rz\"]]; \n  \
ClearAll[\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz]; \
Return[]];",
  StyleBox[" ",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 9D: Analytical (exact) solutions for some benchmark problems, \
also
modules to extract tables and do FEM vs exact plots.  Source: Timoshenko's
Theory of Elasticity Ch 2 and Timoshenko's Theory of Plates and Shells Ch \
4.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PressThickCylinder[{a_,b_},{Em_,\[Nu]_},p_,r_,numer_]:=\n   Module[{ur,uz,\
\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz,err,ezz,e\[Theta]\
\[Theta],grz,ueg}, \n     \[Sigma]rr=-p*a^2(b^2/r^2-1)/(b^2-a^2);\n     \
\[Sigma]\[Theta]\[Theta]= p*a^2(b^2/r^2+1)/(b^2-a^2);\n     \
\[Sigma]zz=2*\[Nu]*a^2*p/(b^2-a^2); \[Sigma]rz=0; ezz=grz=0;\n     err= \
a^2*p*(1+\[Nu])*(b^2-r^2*(1-2*\[Nu]))/((a^2-b^2)*Em*r^2); \n     e\[Theta]\
\[Theta]=-a^2*p*(1+\[Nu])*(b^2+r^2*(1-2*\[Nu]))/((a^2-b^2)*Em*r^2); \n     \
ur=((1+\[Nu])/Em)*a^2*p/(b^2-a^2)((1-2*\[Nu])*r +b^2/r); uz=0;\n     \
ueg={ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz};\n     \
If [!numer, ueg=Simplify[ueg]]; If [numer, ueg=N[ueg]]; \n     \
Return[ueg]];",
  StyleBox[" \n     \n",
    FontColor->RGBColor[0, 0, 1]],
  "RotatingThinDisk[{a_,b_,h_},{Em_,\[Nu]_,\[Rho]_},\[Omega]_,{r_,z_},numer_]:\
=\n   Module[{ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz,\
err,ezz,e\[Theta]\[Theta],grz,ueg}, \n     \
ur=((a^2*(3+\[Nu])*(-(r^2*(-1+\[Nu]))+b^2*(1+\[Nu]))-r^2*(-1+\[Nu])*\n        \
 (-(r^2*(1+\[Nu]))+b^2*(3+\[Nu])))*\[Rho]*\[Omega]^2)/(8*Em*r);\n     uz=\
\[Rho]*\[Omega]^2*z*\[Nu]*(1+\[Nu])*(2*\[Nu]-1)*(-2*r^2*(1+\[Nu])+a^2*(3+\[Nu]\
)+\n        b^2*(3+\[Nu]))/(4*Em*(1+3*\[Nu]));\n     \[Sigma]rr=\[Rho]*\
\[Omega]^2*(3+\[Nu])/8*(b^2+a^2-a^2*b^2/r^2-r^2);\n     \[Sigma]\[Theta]\
\[Theta]=\[Rho]*\[Omega]^2*(3+\[Nu])/8*(b^2+a^2+a^2*b^2/r^2-(1+3*\[Nu])*r^2/(\
3+\[Nu]));\n     \[Sigma]zz=\[Sigma]rz=0;\n     \
ueg={ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz};\n     \
If [!numer, ueg=Simplify[ueg]]; If [numer, ueg=N[ueg]]; \n     \
Return[ueg]];",
  StyleBox[" \n\n",
    FontColor->RGBColor[0, 0, 1]],
  "PointLoadCircPlate[{R_,h_},{Em_,\[Nu]_},P_,{rr_,z_},numer_]:=\n   \
Module[{Dm,Mrr,M\[Theta]\[Theta],c,ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\
\[Theta]\[Theta],\[Sigma]rz,duzdr,r,usig},",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "    r=Max[R/1000,rr];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "If [Head[rr]===Symbol,r=rr];\n    Dm=Em*h^3/(12*(1-\[Nu]^2)); \
c=(3+\[Nu])/(1+\[Nu]);\n    uz=P/(16*Pi*Dm)*(c*(R^2-r^2)+2*r^2*Log[r/R]);\n   \
 ur=\[InvisibleSpace](P*r*z*(c-1-2*Log[r/R]))/(8*Dm*Pi);\n    If [rr==0, \
uz=c*P*R^2/(16*Pi*Dm); ur=0];\n    Mrr=-(P/(4*Pi))*(1+\[Nu])*Log[R/r];\n    M\
\[Theta]\[Theta]=-(P/(4*Pi))*((1+\[Nu])*Log[R/r]+1-\[Nu]);\n    \
\[Sigma]rr=-12*Mrr*z/h^3; \
\[Sigma]\[Theta]\[Theta]=-12*M\[Theta]\[Theta]*z/h^3; \
\[Sigma]zz=\[Sigma]rz=0;\n    usig=Simplify[{ur,uz,\[Sigma]rr,\[Sigma]zz,\
\[Sigma]\[Theta]\[Theta],\[Sigma]rz}];\n    If [!numer, usig=Simplify[usig]]; \
If [numer,usig=N[usig]];\n    Return[usig]];",
  StyleBox[" ",
    FontColor->RGBColor[0, 0, 1]],
  "\n",
  StyleBox["\n",
    FontColor->RGBColor[0, 0, 1]],
  "ExactSolution[nodrzc_,dimen_,{Em_,\[Nu]_,\[Rho]_},load_,\n  \
solution_,numer_]:=Module[{n,r,z,p,P,R,\n  ur,uz,\[Sigma]rr,\[Sigma]zz,\
\[Sigma]\[Theta]\[Theta],\[Sigma]rz,numnod=Length[nodrzc],\n  \
exactnoddis,exactnodsig,modname=\"ExactSolution: \"},\n  \
exactnoddis=Table[{0,0},{numnod}];\n  exactnodsig=Table[{0,0,0,0},{numnod}]; \
\n  If [solution==\"PressThickCylinder\", p=load; {a,b}=dimen;\n     For \
[n=1,n<=numnod,n++, {r,z}=nodrzc[[n]];\n         \
{ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}=\n          \
PressThickCylinder[{a,b},{Em,\[Nu]},p,r,numer];\n          \
exactnoddis[[n]]={ur,uz};\n          exactnodsig[[n]]={\[Sigma]rr,\[Sigma]zz,\
\[Sigma]\[Theta]\[Theta],\[Sigma]rz}];\n     \
Return[{exactnoddis,exactnodsig}]];\n  If [solution==\"RotatingThinDisk\", \
\[Omega]=load; {a,b,h}=dimen;\n     For [n=1,n<=numnod,n++, \
{r,z}=nodrzc[[n]];\n         {ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\
\[Theta],\[Sigma]rz}=\n          RotatingThinDisk[{a,b,h},{Em,\[Nu],\[Rho]},\
\[Omega],{r,z},numer];\n          exactnoddis[[n]]={ur,uz};\n          \
exactnodsig[[n]]={\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}];\
\n     Return[{exactnoddis,exactnodsig}]];\n  If \
[solution==\"PointLoadCircPlate\", P=load; {R,h}=dimen;\n     For \
[n=1,n<=numnod,n++, {r,z}=nodrzc[[n]];\n         \
{ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}=\n          \
PointLoadCircPlate[{R,h},{Em,\[Nu]},-P,{r,z},numer];\n          \
exactnoddis[[n]]={ur,uz};\n          exactnodsig[[n]]={\[Sigma]rr,\[Sigma]zz,\
\[Sigma]\[Theta]\[Theta],\[Sigma]rz}];\n     \
Return[{exactnoddis,exactnodsig}]]; \n  Print [modname,solution,\" not \
available\"];\n  Return[{exactnoddis,exactnodsig}]]; "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 9E. Driver to get FEM result vs exact solution r-plots for \
benchmark problems of Chapter 13.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "FillUpFEMResultPlotArray[etype_,nodrzc_,noddis_,nodsig_,\n   \
Ner_,Nez_,iwhat_,ilayer_]:=Module[{n,nn,n1,n2,n3,n4,\n   \
er,ez,k,ir,nr,cn,ncn,ntab,rtab,vtab,vFEM,m=ilayer,\n   \
numnod=Length[nodrzc],num9nod},\n",
  StyleBox["   (*Print[\"enter FillUpFEMResultPlotArray\"]; \
Print[\"numnod=\",numnod];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   If [etype==\"Quad4\",\n       nr=Ner+1; \
ntab=Table[(Nez+1)*n-Nez+m,{n,Ner+1}]];\n   If [etype==\"Quad8\", \
num9nod=(2*Ner+1)*(2*Nez+1);\n       nr=2*Ner+1; \
ntab=Table[(2*Nez+1)*(ir-1)+1+m,{ir,1,nr}];\n       \
cn=Table[2*Nez+3+2*(ez-1)+(4*Nez+2)*(er-1),{er,Ner},{ez,Nez}];\n     ",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  " cn=Flatten[cn];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "ncn=Table[0,{num9nod}]; k=0;\n       For [n=1,n<=num9nod,n++, If \
[MemberQ[cn,n],k++];\n            ncn[[n]]=k];",
  StyleBox["\n       (*Print[\"cn=\",cn];  Print[\"ncn=\",ncn];",
    FontColor->RGBColor[1, 0, 0]],
  "\n       ",
  StyleBox["Print[\"original ntab=\",ntab];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n       For [n=1,n<=nr,n++, nn=ntab[[n]]; k=ncn[[nn]];\n            If \
[MemberQ[cn,nn],ntab[[n]]=-nn; Continue[]];\n            \
ntab[[n]]=ntab[[n]]-k]]; ",
  StyleBox[" \n       (*Print[\"updated ntab=\",ntab];*)\n",
    FontColor->RGBColor[1, 0, 0]],
  "       If [Mod[m,2]==1, ",
  StyleBox["\n       ",
    FontColor->RGBColor[1, 0, 0]],
  "For [n=1,n<=nr,n++, nn=ntab[[n]]; If [nn>0, Continue[]]; \n            \
nn=Abs[nn]; ntab[[n]]={ntab[[n-1]],ntab[[n+1]],\n              \
nn-1-ncn[[nn-1]],nn+1-ncn[[nn+1]]}; ",
  StyleBox["\n           ",
    FontColor->RGBColor[1, 0, 0]],
  " ]]; \n   vFEM=Table[0,{nr}]; \n   rtab=Table[nodrzc[[n,1]],{n,numnod}];\n \
  If [iwhat<=2, vtab=Table[noddis[[n,iwhat]],  {n,numnod}]];\n   If \
[iwhat>=3, vtab=Table[nodsig[[n,iwhat-2]],{n,numnod}]];\n",
  StyleBox["   (*Print[\"rtab=\",rtab]; Print[\"vtab=\",vtab];*)\n",
    FontColor->RGBColor[1, 0, 0]],
  "   For [ir=1,ir<=nr,ir++, n=ntab[[ir]];\n        If [Length[n]==2, \
{n1,n2}=n; \n            vFEM[[ir]]={rtab[[n1]]+rtab[[n2]],\n             \
vtab[[n1]]+vtab[[n2]]}/2; Continue[]];\n        If [Length[n]==4, \
{n1,n2,n3,n4}=n; \n            vFEM[[ir]]={rtab[[n1]]+rtab[[n2]]+\n           \
 rtab[[n3]]+rtab[[n4]],vtab[[n1]]+vtab[[n2]]+\n            \
vtab[[n3]]+vtab[[n4]]}/4; Continue[]];    \n        \
vFEM[[ir]]={rtab[[n]],vtab[[n]]}; \n        ]; \
ClearAll[cn,ncn,ntab,rtab,vtab];\n",
  StyleBox["   (*Print[\"vFEM=\",vFEM//MatrixForm];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   Return[vFEM]];\n",
  StyleBox[" ",
    FontColor->RGBColor[0, 0, 1]],
  "   \nRadialPlotFEMvsExact[etype_,nodrzc_,noddis_,nodsig_,\n  dimen_,{Em_,\
\[Nu]_,\[Rho]_},load_,Nerz_,prob_,what_,ilayer_,\n  \
numer_]:=Module[{n,i,nr,plotwhat,rbeg,rend,iwhat,r,z,\n   Ner,Nez=1,ur,uz,\
\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz,Q4layer,Q8layer,\n  \
 a,b,h,hh,R,p,P,\[Omega],vFEM,vexact,probOK,whatOK,name,\n   \
zlayer,modname=\"RadialPlotFEMvsExact: \"}, \n   ClearAll[r,z];\n   \
probOK={\"PressThickCylinder\",\"PointLoadCircPlate\",\n           \
\"RotatingThinDisk\"};\n   \
whatOK={\"ur\",\"uz\",\"\[Sigma]rr\",\"\[Sigma]zz\",\"\[Sigma]\[Theta]\[Theta]\
\",\"\[Sigma]rz\"};\n   Q4zlayer=(hh/2)*{{-1,1},{-1,0,1},{-1,-1/3,1/3,1},\n   \
      {-1,-1/2,0,1/2,1},{-1,-3/5,-1/5,1/5,3/5,1},\n         \
{-1,-2/3,-1/3,0,1/3,2/3,1}};\n   Q8zlayer=(hh/2)*{{-1,0,1},{-1,-1/2,0,1/2,1},\
\n         {-1,-2/3,-1/3,0,1/3,2/3,1}};  \n   Ner=Nerz; If \
[Length[Nerz]==2,{Ner,Nez}=Nerz];\n   If [!MemberQ[probOK,prob], \
Print[modname,\n       \"exact sol of problem \",prob,\n       \" not \
implemented\"]; Return[]];\n   If [!MemberQ[whatOK,what], Print[modname,\n    \
   \"plot of component \",what,\n       \" not implemented\"]; Return[]]; \n  \
 If [prob==\"PressThickCylinder\", p=load; \n       {a,b}=Take[dimen,2]; \
{rbeg,rend}={a,b};\n       {ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\
\[Theta],\[Sigma]rz}=PressThickCylinder[\n                           \
{a,b},{Em,\[Nu]},p,r,numer];\n        title=\"Pressurized thick cyl: \"];\n   \
If [prob==\"RotatingThinDisk\", \[Omega]=load; \n       {a,b,h}=dimen; \
{rbeg,rend}={a,b};\n       {ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\
\[Theta],\[Sigma]rz}= RotatingThinDisk[\n                  {a,b,h},{Em,\[Nu],\
\[Rho]},\[Omega],{r,z},numer];\n        title=\"Rotating Thin Disk: \"];\n   \
If [prob==\"PointLoadCircPlate\", P=load; \n        {R,h}=Take[dimen,2]; \
{rbeg,rend}={R/1000,R};\n        \
{ur,uz,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}=\
PointLoadCircPlate[\n                        \
{R,h},{Em,\[Nu]},-P,{r,z},numer];\n        If [etype==\"Quad4\",\n           \
zlayer=Q4zlayer[[Nez,ilayer+1]]/.hh->h];\n        If [etype==\"Quad8\",\n     \
      zlayer=Q8zlayer[[Nez,ilayer+1]]/.hh->h];\n",
  StyleBox["        (*Print[\"ilayer=\",ilayer]; Print[\"zlayer=\",zlayer];*) \
  ",
    FontColor->RGBColor[1, 0, 0]],
  "\n        \
{ur,\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}={ur,\[Sigma]rr,\
\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}/.z->zlayer;\n        \
title=\"Point loaded circ plate: \"];\n   {{iwhat}}=Position[whatOK,what];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "name=whatOK[[iwhat]];\n   If [iwhat<=2,  plotwhat=\"disp \"<>name;\n       \
vexact={ur,uz}[[iwhat]] ]; \n   If [iwhat>=3,  plotwhat=\"stress \"<>name;\n  \
     vexact={\[Sigma]rr,\[Sigma]zz,\[Sigma]\[Theta]\[Theta],\[Sigma]rz}[[\
iwhat-2]] ];  \n   vFEM=FillUpFEMResultPlotArray[etype,nodrzc,noddis,nodsig,\n\
            Ner,Nez,iwhat,ilayer];   \n   If [numer, \
{vFEM,rbeg,rend}=N[{vFEM,rbeg,rend}]];\n",
  StyleBox["   (*Print[\"vFEM=\",Transpose[vFEM]//MatrixForm];\n     \
Print[\"name=\",name,\"  exact @ rbeg=\",N[vexact/.r->rbeg]];*)\n",
    FontColor->RGBColor[1, 0, 0]],
  "   If   [$VersionNumber<9.0,\n     \
pFEM=ListPlot[vFEM,PlotJoined->True,PlotStyle->{RGBColor[1,0,0]},\n        \
DisplayFunction->Identity],",
  "\n     pFEM=ListPlot[vFEM,Joined->True,PlotStyle->{RGBColor[1,0,0]},\n     \
   DisplayFunction->Identity]\n        ];\n   \
pexact=Plot[vexact,{r,rbeg,rend},PlotStyle->{RGBColor[0,0,0]},\n        \
DisplayFunction->Identity];\n   Show[Graphics[AbsoluteThickness[2]],pexact,  \
\n        pFEM,TextStyle->{FontFamily->\"Times\",FontSize->9},\n        \
GridLines->Automatic,Axes->True,PlotRange->All,\n        \
PlotLabel->title<>plotwhat<>\" (black=exact,red=FEM)\",\n        \
AspectRatio->1/GoldenRatio,\n        Frame->True, \
DisplayFunction->DisplayChannel]; \n   ClearAll[vFEM]; Return[]]; \n   \n",
  StyleBox["(*etype=\"Quad8\"; Ner=2; Nez=2; iwhat=5; zlayer=3;\n\
nodrzc=Table[{r[i],z[i]},{i,1,25}];\nnoddis=Table[{ur[i],uz[i]},{i,1,25}];\n\
nodsig=Table[{\[Sigma]rr[i],\[Sigma]zz[i],\[Sigma]\[Theta]\[Theta][i],\[Sigma]\
rz[i]},{i,1,25}];\nvFEM=FillUpFEMResultPlotArray[etype,nodrzc,noddis,nodsig,\n\
        Ner,Nez,zlayer,iwhat];\n\
Print[\"vFEM=\",Transpose[vFEM]//MatrixForm];*)",
    FontColor->RGBColor[0, 0, 1]],
  "  "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 10:  Driver for static analysis of axisymmetric solid problem.\
\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "RingAnalysisDriver[nodrzc_,eletyp_,elenod_,elemat_,elebfor_,\n  \
eletfor_,doftag_,dofval_,defopt_]:= Module[{n,ns,j,\n  \
numnod=Length[nodrzc],numele=Length[elenod],eps=10.^(-7),\n  \
f,fb,ft,K,Kmod,u,supdof,supval,noddis,nodfor,nodsig,\n  \
modname=\"RingAnalysisDriver: \"}, \n  noddis=nodfor=nodsig={}; \
fb=ft=Table[0,{2*numnod}];\n  \
K=RingMasterStiffness[nodrzc,eletyp,elenod,elemat,defopt];\n",
  StyleBox["  (*Print[\"Master Stiff K:\\n\",K//MatrixForm];\n  Print[\"eigs \
of K=\",Chop[Eigenvalues[N[K]]]];\n  Print[\"doftag=\",doftag,\" \
dofval=\",dofval]; *)\n  ",
    FontColor->RGBColor[1, 0, 0]],
  "ns=0; Do [Do [If[doftag[[n,j]]>0,ns++],{j,1,2}],{n,1,numnod}];\n  \
supdof=supval=Table[0,{ns}];\n  k=0; Do [Do \
[If[doftag[[n,j]]>0,k++;supdof[[k]]=2*(n-1)+j;\n         \
supval[[k]]=dofval[[n,j]]],{j,1,2}],{n,numnod}];\n",
  StyleBox["  (*Print[\"supdof=\",supdof];*)\n  ",
    FontColor->RGBColor[1, 0, 0]],
  "If [Length[elebfor]>0, ",
  StyleBox["\n      ",
    FontColor->RGBColor[1, 0, 0]],
  "fb=RingBodyForces[nodrzc,eletyp,elenod,elebfor,defopt]];\n  If \
[Length[eletfor]>0,  Print[modname,",
  StyleBox["\n      ",
    FontColor->RGBColor[1, 0, 0]],
  "\" surface traction force lumping unavailable ",
  "- ",
  "sorry\"]];\n  f=ModifiedNodeForces[supdof,supval,K,Flatten[dofval],fb+ft]; \
    \n  Kmod=ModifiedMasterStiffness[supdof,K];\n",
  StyleBox["  (*Print[\"Force vector f=\",f];*)\n  (* Print[\"eigs of \
Kmod=\",Eigenvalues[N[Kmod]]];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n  u=LinearSolve[Kmod,f]; u=Chop[u,eps]; noddis=Partition[u,2];\n  \
f=Simplify[K.u]; f=Chop[f,eps]; nodfor=Partition[f,2];\n  \
nodsig=RingNodalStresses[nodrzc,eletyp,elenod,elemat,\n       noddis,defopt]; \
nodsig=Chop[nodsig,eps]; \n  ClearAll[K,Kmod,u,f,fb,ft]; \n  \
Return[{noddis,nodfor,nodsig}]; \n  ];   "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
 Cell 11:   Pressurized thick cylinder with 2 Quad4 elements in the
radial direction. Cell used for the driver script example of Chapter \
13.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[Em,\[Nu],th,a,b,d,p,numer];  
Em=1000.; \[Nu]=0; etype=\"Quad4\"; numer=True;
Kfac=1; a=4; b=10; d=2; p=10; aspect=d/(b-a);

(*  Define FEM model *)

NodeCoordinates=N[{{a,0},{a,d},{(a+b)/2,0},
    {(a+b)/2,d},{b,0},{b,d}}];          
ElemNodes={{1,3,4,2},{3,5,6,4}};
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes]; 
ElemType=Table[etype,{numele}]; 
ElemMaterial=Table[N[{Em,\[Nu]}],{numele}]; 
FreedomValues=Table[{0,0},{numnod}];
FreedomTags=Table[{0,1},{numnod}];  pfor=N[Kfac*p*a*d];
FreedomValues[[1]]=FreedomValues[[2]]={pfor/2,0}; 
ElemBodyForces={}; ElemTractionForces={};  
Defaultoptions={numer};
PrintRingNodeCoordinates[NodeCoordinates,\"Node Coordinates\",{2,4}];
PrintRingElementTypeNodes[ElemType,ElemNodes,\"Element type & nodes\",{}];
PrintRingElementMaterials[ElemMaterial,\"Material data\",{}]; 
PrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{}];
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"Pressurized thick cylinder mesh\",True,True];

(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[
    NodeCoordinates,ElemType,ElemNodes,
    ElemMaterial,ElemBodyForces,ElemTractionForces,
    FreedomTags,FreedomValues,Defaultoptions];
PrintRingAnalysisSolution[NodeDisplacements,NodeForces,
     NodeStresses,\"Computed solution\",{6,6}];  
{ExactNodeDisplacements,ExactNodeStresses}=
     ExactSolution[NodeCoordinates,{a,b},{Em,\[Nu],\[Rho]},p,
     \"PressThickCylinder\",numer];
PrintRingNodeDispStresses[ExactNodeDisplacements,
    ExactNodeStresses,\"Exact (Lame) solution\",{}];  
 
(*  Plot Node Stress Distributions *)

legend={(a+2*b)/3,0.8*d};
ContourPlotStresses[NodeCoordinates,ElemNodes,NodeStresses,
  {True,False,True,False},True,{},legend,aspect];
    
(* Through-wall plots comparing FEM vs exact solutions *)

pwhat={\"ur\",\"\[Sigma]rr\",\"\[Sigma]zz\",\"\[Sigma]\[Theta]\[Theta]\"}; If \
[\[Nu]==0,pwhat={\"ur\",\"\[Sigma]rr\",\"\[Sigma]\[Theta]\[Theta]\"}];
For [ip=1,ip<=Length[pwhat],ip++, what=pwhat[[ip]];
     RadialPlotFEMvsExact[etype,NodeCoordinates,NodeDisplacements,
       NodeStresses, {a,b},{Em,\[Nu],\[Rho]},p,{2,1},
       \"PressThickCylinder\",what,0,numer]  ];\
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 12:   Pressurized thick cylinder with Quad4 elems: Ner >=1  \
elements in 
radial direction, Nez=1 in the axial direction. a,b inner & outer radii, d \
slice thickness\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[Em,\[Nu],th,a,b,d,p,Ner,Nez,numer];  
Em=1000.; \[Nu]=0; etype=\"Quad4\"; numer=True;
Kfac=1; a=4; b=10; d=2; p=10; Ner=8; Nez=1; aspect=d/(b-a);

(*  Define FEM model using mesh generation *)

CornerNodeCoordinates=N[{{a,0},{b,0},{b,d},{a,d}}]; 
NodeCoordinates=
  GenQuad4NodeCoordinates[CornerNodeCoordinates,Ner,Nez];
ElemNodes=GenQuad4ElemNodes[Ner,Nez];
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes]; 
ElemType=Table[etype,{numele}]; 
ElemMaterial=Table[N[{Em,\[Nu]}],{numele}]; 
FreedomValues=Table[{0,0},{numnod}];
FreedomTags=Table[{0,1},{numnod}];  pfor=N[Kfac*p*a*d];
FreedomValues[[1]]=FreedomValues[[2]]={pfor/2,0}; 
ElemBodyForces={}; ElemTractionForces={};  
Defaultoptions={numer};
PrintRingNodeCoordinates[NodeCoordinates,\"Node Coordinates\",{2,4}];
PrintRingElementTypeNodes[ElemType,ElemNodes,\"Element type & nodes\",{}];
PrintRingElementMaterials[ElemMaterial,\"Material data\",{}]; 
PrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{}];
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];

Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"Pressurized thick cylinder Quad4 mesh\",True,True];

(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[
    NodeCoordinates,ElemType,ElemNodes,
    ElemMaterial,ElemBodyForces,ElemTractionForces,
    FreedomTags,FreedomValues,Defaultoptions];
PrintRingAnalysisSolution[NodeDisplacements,NodeForces,
     NodeStresses,\"Computed solution\",{6,6}];  
{ExactNodeDisplacements,ExactNodeStresses}=
     ExactSolution[NodeCoordinates,{a,b},{Em,\[Nu],\[Rho]},p,
     \"PressThickCylinder\",numer];
PrintRingNodeDispStresses[ExactNodeDisplacements,
    ExactNodeStresses,\"Exact (Lame) solution\",{}];  
 
(*  Contourplot Node Stress Distributions *)

legend={(a+b)/2,0.75*d}; whichones={True,True,True,False};
If [\[Nu]==0, whichones={True,False,True,False}]; 
ContourPlotStresses[NodeCoordinates,ElemNodes,NodeStresses,
  whichones,True,{},legend,aspect];
  
(* Through-wall plots comparing FEM vs exact solutions *)

pwhat={\"ur\",\"\[Sigma]rr\",\"\[Sigma]zz\",\"\[Sigma]\[Theta]\[Theta]\"}; If \
[\[Nu]==0,pwhat={\"ur\",\"\[Sigma]rr\",\"\[Sigma]\[Theta]\[Theta]\"}];
For [ip=1,ip<=Length[pwhat],ip++, what=pwhat[[ip]];
     RadialPlotFEMvsExact[etype,NodeCoordinates,NodeDisplacements,
       NodeStresses,{a,b},{Em,\[Nu],\[Rho]},p,{Ner,Nez},
      \"PressThickCylinder\",what,0,numer]  ];\
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 13:   Pressurized thick cylinder accepting both Quad4 or Quad8 \
elems; Ner >=1 elements in 
radial direction, Nez (assumed 1) in the axial direction. a,b inner & outer \
radii, d slice thickness\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ClearAll[Em,\[Nu],th,a,b,d,p,Ner,Nez];  \nEm=1000.; \[Nu]=0.0; \
etype=\"Quad4\"; numer=True; Ner=4; Nez=1;\nKfac=1; a=4; b=10; d=2; \
aspect=d/(b-a); p=10; \n\n(*  Define FEM model *)\n\n\
MeshCorners=N[{{a,0},{b,0},{b,d},{a,d}}]; \nIf [etype==\"Quad4\",\n    \
NodeCoordinates=GenQuad4NodeCoordinates[MeshCorners,Ner,Nez];          \n    \
ElemNodes=GenQuad4ElemNodes[Ner,Nez]];\nIf [etype==\"Quad8\",\n    \
NodeCoordinates=GenQuad8NodeCoordinates[MeshCorners,Ner,Nez];          \n    \
ElemNodes=GenQuad8ElemNodes[Ner,Nez]];\nnumnod=Length[NodeCoordinates]; \
numele=Length[ElemNodes]; \nElemType=Table[etype,{numele}]; \n\
ElemMaterial=Table[N[{Em,\[Nu]}],{numele}]; \n\
FreedomTags=Table[{0,1},{numnod}];\nFreedomValues=Table[{0,0},{numnod}]; \
pfor=N[Kfac*p*a*d];\nIf [etype==\"Quad4\",\n    \
FreedomValues[[1]]=FreedomValues[[2]]={pfor/2,0}]; \nIf [etype==\"Quad8\",\n  \
  FreedomValues[[1]]=FreedomValues[[3]]={pfor/6,0};\n    FreedomValues[[2]] \
={2*pfor/3,0}];  \nElemBodyForces= ElemTractionForces={}; \
Defaultoptions={numer};\n(* Problem data print statements removed *)\n\
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,\n      \"Pressurized \
thick cylinder\",True,True]; \n      \n(*  Solve problem and print results *)\
\n\n{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[\n    \
NodeCoordinates,ElemType,ElemNodes,\n    \
ElemMaterial,ElemBodyForces,ElemTractionForces,\n    \
FreedomTags,FreedomValues,",
  "Defaultoptions",
  "];\nPrintRingAnalysisSolution[NodeDisplacements,NodeForces,\n     \
NodeStresses,\"Computed solution\",{}]; \n\
{ExactNodeDisplacements,ExactNodeStresses}=\n     \
ExactSolution[NodeCoordinates,{a,b},{Em,\[Nu],\[Rho]},p,\n     \
\"PressThickCylinder\",numer];\n\
PrintRingNodeDispStresses[ExactNodeDisplacements,\n    \
ExactNodeStresses,\"Exact (Lame) solution\",{}];   \n\n(*  Plot Node Stress \
Distributions *)\n\nlegend={(a+b)/2,0.75*d}; \
whichones={True,True,True,False}; \nIf [\[Nu]==0, \
whichones={True,False,True,False}]; \n\
ContourPlotStresses[NodeCoordinates,ElemNodes,NodeStresses,\n  \
whichones,True,{},legend,aspect];\n  \n(* Through-wall plots comparing FEM vs \
exact solutions *)\n\npwhat={\"ur\",\"\[Sigma]rr\",\"\[Sigma]zz\",\"\[Sigma]\
\[Theta]\[Theta]\"};",
  StyleBox["  ",
    FontColor->RGBColor[1, 0, 0]],
  "\nFor [ip=1,ip<=Length[pwhat],ip++, what=pwhat[[ip]];\n     \
RadialPlotFEMvsExact[etype,NodeCoordinates,NodeDisplacements,\n       \
NodeStresses,{a,b},{Em,\[Nu],\[Rho]},p,{Ner,Nez},\n      \"PressThickCylinder\
\",what,1,numer]  ];   "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 14.  Benchmark #2:  Rotating thin disk.  Inner radius a, \
exterior radius b, z-thickness h, mass
density \[Rho]. Disk rotates about z at \[Omega] rads/sec. Accepts both Quad4 \
and Quad8, just change etype.
Body force computed at element corner nodes. Ner arbitrary but Nez=1;\
\>", \
"Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ClearAll[Em,\[Nu],a,b,h,Kfac,\[Rho],\[Omega],Ner,Nez,numer];   \nEm=1000.; \
\[Nu]=N[1/3];  Ner=2; Nez=1; etype=\"Quad8\";\nKfac=1; a=4; b=10; h=1; \
aspect=h/(b-a); \[Rho]=3.0; \[Omega]=0.5;\nnumer=True;\n\n(*  Define FEM \
model *)\n\nMeshCorners=N[{{a,0},{b,0},{b,h},{a,h}}];\nIf [etype==\"Quad4\",\n\
    NodeCoordinates=GenQuad4NodeCoordinates[MeshCorners,Ner,Nez];          \n \
   ElemNodes=  GenQuad4ElemNodes[Ner,Nez]];\nIf [etype==\"Quad8\",\n    \
NodeCoordinates=GenQuad8NodeCoordinates[MeshCorners,Ner,Nez];          \n    \
ElemNodes=  GenQuad8ElemNodes[Ner,Nez]];\nnumnod=Length[NodeCoordinates]; \
numele=Length[ElemNodes];\nElemType=      Table[etype,{numele}]; \n\
ElemMaterial=  Table[{Em,\[Nu]},{numele}]; \n\
ElemBodyForces=Table[{0,0},{numele}];\nFor [e=1,e<=numele,e++, \
enl=ElemNodes[[e]];\n     ncoor=Table[NodeCoordinates[[enl[[i]]]],{i,4}];\n   \
  {{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;\n     ElemBodyForces[[e]]=\[Rho]*\
\[Omega]^2*{{r1,0},{r2,0},{r3,0},{r4,0}}]; ",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "FreedomTags=FreedomValues=Table[{0,0},{numnod}];\nIf [etype==\"Quad4\",  \n\
    For [n=1,n<=numnod-Nez,n=n+Nez+1, FreedomTags[[n]]={0,1}]];\nIf \
[etype==\"Quad8\",  \n    For [n=1,n<=numnod-2*Nez,n=n+3*Nez+2, \
FreedomTags[[n+1]]={0,1}]];\nElemTractionForces={}; Defaultoptions={True};\n\
PrintRingNodeCoordinates[NodeCoordinates,\"Node Coordinates\",{4,2}];\n\
PrintRingElementTypeNodes[ElemType,ElemNodes,\"Element type & nodes\",{}];\n\
PrintRingElementMaterials[ElemMaterial,\"Material data\",{}]; \n\
PrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{5,3}];\n\
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];\n\
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];\n\
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,\n    \"Rotating disk \
mesh\",True,True]; \n  \n(*  Solve problem and print results *)\n\n\
{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[\n    \
NodeCoordinates,ElemType,ElemNodes,\n    \
ElemMaterial,ElemBodyForces,ElemTractionForces,\n    \
FreedomTags,FreedomValues,",
  "Defaultoptions",
  "];\nPrintRingAnalysisSolution[NodeDisplacements,NodeForces,\n    \
NodeStresses,\"Computed solution\",{}]; \n\
{ExactNodeDisplacements,ExactNodeStresses}=\n    \
ExactSolution[NodeCoordinates,{a,b,h},{Em,\[Nu],\[Rho]},\[Omega],\n     \
\"RotatingThinDisk\",numer];\n\
PrintRingNodeDispStresses[ExactNodeDisplacements,\n    \
ExactNodeStresses,\"Exact solution\",{}]; \n    \n(*  Plot Node Stress \
Distributions *)\n\nlegend={(a+b)/2,0.75*h}; \
whichones={True,False,True,False}; \n\
ContourPlotStresses[NodeCoordinates,ElemNodes,NodeStresses,\n  \
whichones,True,{},legend,aspect];\n  \n(* Through-wall plots comparing FEM vs \
exact solutions *)\n\npwhat={\"ur\",\"\[Sigma]rr\",\"\[Sigma]zz\",\"\[Sigma]\
\[Theta]\[Theta]\"};\nFor [ip=1,ip<=Length[pwhat],ip++, what=pwhat[[ip]];\n   \
  RadialPlotFEMvsExact[etype,NodeCoordinates,NodeDisplacements,\n       \
NodeStresses,{a,b,h},{Em,\[Nu],\[Rho]},\[Omega],{Ner,Nez},\n       \
\"RotatingThinDisk\",what,0,numer]  ]; "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 15: Benchmark #3. Thin circular plate under central point \
load.  Exact results taken from
 Timoshenko's Theory of Plates and Shells, Chapter 4.  Both Quad4 and Quad8 \
done by same
 script; just change etypr. Ner elements along radius, Nez=2 (Quad4) or 1 \
(Quad8) \
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ClearAll[Em,\[Nu],a,b,R,h,Kfac,P,Ner,Nez,numer];   \nEm=1000.; \
\[Nu]=N[1/3];  Ner=4; etype=\"Quad4\";\nKfac=1; R=10; h=1; aspect=h/R; P=10.;\
\nNez=2; If [etype==\"Quad8\", Nez=1]; numer=True;\n\n(*  Define FEM model *)\
\n\nMeshCorners=N[{{0,-h/2},{R,-h/2},{R,h/2},{0,h/2}}];\nIf \
[etype==\"Quad4\",\n    \
NodeCoordinates=GenQuad4NodeCoordinates[MeshCorners,Ner,Nez];          \n    \
ElemNodes=  GenQuad4ElemNodes[Ner,Nez]];\nIf [etype==\"Quad8\",\n    \
NodeCoordinates=GenQuad8NodeCoordinates[MeshCorners,Ner,Nez];          \n    \
ElemNodes=  GenQuad8ElemNodes[Ner,Nez]];\nnumnod=Length[NodeCoordinates]; \
numele=Length[ElemNodes];\nElemType=      Table[etype,{numele}]; \n\
ElemMaterial=  Table[{Em,\[Nu]},{numele}]; \n\
ElemBodyForces=ElemTractionForces={}; Defaultoptions={True};",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "FreedomTags=FreedomValues=Table[{0,0},{numnod}];\n\
FreedomTags[[1]]=FreedomTags[[2]]=FreedomTags[[3]]={1,0};\n\
FreedomTags[[numnod-1]]={0,1}; Ps=N[P*Kfac/(2*Pi)];\nIf [etype==\"Quad4\",  \n\
   FreedomValues[[1]]=FreedomValues[[3]]={0,-Ps/4}; \n   \
FreedomValues[[2]]={0,-Ps/2}]; \nIf [etype==\"Quad8\",  \n   \
FreedomValues[[1]]=FreedomValues[[3]]={0,-Ps/6}; \n   \
FreedomValues[[2]]={0,-2*Ps/3}]; \n\
PrintRingNodeCoordinates[NodeCoordinates,\"Node Coordinates\",{4,2}];\n\
PrintRingElementTypeNodes[ElemType,ElemNodes,\"Element type & nodes\",{}];\n\
PrintRingElementMaterials[ElemMaterial,\"Material data\",{}]; \n\
PrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{5,3}];\n\
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];\n\
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];\n\
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,\n    \"Point loaded \
circular plate mesh\",True,True]; \n    \n(*  Solve problem and print results \
*)\n\n{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[\n    \
NodeCoordinates,ElemType,ElemNodes,\n    \
ElemMaterial,ElemBodyForces,ElemTractionForces,\n    \
FreedomTags,FreedomValues,",
  "Defaultoptions",
  "];\nPrintRingAnalysisSolution[NodeDisplacements,NodeForces,\n    \
NodeStresses,\"Computed solution\",{}]; \n\
{ExactNodeDisplacements,ExactNodeStresses}=\n    \
ExactSolution[NodeCoordinates,{R,h},{Em,\[Nu],\[Rho]},P,\n     \
\"PointLoadCircPlate\",numer];\n\
PrintRingNodeDispStresses[ExactNodeDisplacements,\n    \
ExactNodeStresses,\"Exact solution\",{}]; \n   \n(*  Plot Node Stress \
Distributions *)\n\nwhichones={True,False,True,False}; \n\
ContourPlotStresses[NodeCoordinates,ElemNodes,NodeStresses,\n  \
whichones,True,{},{},aspect];\n  \n(* Through-wall plots comparing FEM vs \
exact solutions *)\n\npwhat={\"uz\",\"ur\",\"\[Sigma]rr\",\"\[Sigma]zz\",\"\
\[Sigma]\[Theta]\[Theta]\"};\nFor [ip=1,ip<=Length[pwhat],ip++, \
what=pwhat[[ip]];\n     RadialPlotFEMvsExact[etype,NodeCoordinates,\n       \
NodeDisplacements,NodeStresses, {R,h},{Em,\[Nu],\[Rho]},\n       P,{Ner,Nez},\
\"PointLoadCircPlate\",what,0,numer]  ];\n "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 16: reserved for benchmark problem\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\

\
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 17.  Concrete pile embebded in elastic soil - very coarse \
mesh.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[Epile,\[Nu]pile,Esoil,\[Nu]soil,d,H,P,L,Kfac];   
Epile=300000.; Esoil=Epile/50; \[Nu]pile=0.1; \[Nu]soil=0.4; 
P=500000; d=10; L=250; H=1.2*L; Ner=3; Nez=3; R=8*d; 
Kfac=1; etype=\"Quad4\"; aspect=H/(2*R);  

(*  Define FEM model *)

MeshCorners=N[{{0,-H},{R,-H},{R,0},{0,0}}];
NodeCoordinates=GenGradedQuad4NodeCoordinates[MeshCorners,
     Ner,Nez,{-1,d/R-1,-0.4,1},{-1,1-2*L/H,0.1,1}];        
ElemNodes=GenQuad4ElemNodes[Ner,Nez];
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes];
ElemType=      Table[etype,{numele}]; 
ElemMaterial=   Table[{Esoil,\[Nu]soil},{numele}]; 
ElemMaterial[[2]]=ElemMaterial[[3]]={Epile,\[Nu]pile}; 
ElemBodyForces=ElemTractionForces={};
FreedomTags=        Table[{0,0},{numnod}];  
For[n=2,n<=4,n++,   FreedomTags[[n]]={1,0}];
For[n=14,n<=16,n++, FreedomTags[[n]]={1,1}];
For[n=1,n<=13,n=n+4,FreedomTags[[n]]={1,1}];
FreedomValues=      Table[{0,0},{numnod}];
FreedomValues[[4]]={0,N[-P*Kfac/(2*Pi)]};  
Defaultoptions={True};
PrintRingNodeCoordinates[NodeCoordinates,\"Node Coordinates\",{4,2}];
PrintRingElementTypeNodes[ElemType,ElemNodes,\"Element type & nodes\",{}];
PrintRingElementMaterials[ElemMaterial,\"Material data\",{}]; 
PrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{}];
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"Pile in soil - very coarse mesh\",True,True]; 
    
(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[
    NodeCoordinates,ElemType,ElemNodes,
    ElemMaterial,ElemBodyForces,ElemTractionForces,
    FreedomTags,FreedomValues,Defaultoptions];
PrintRingAnalysisSolution[NodeDisplacements,NodeForces,
     NodeStresses,\"Computed solution\",{6,3}];  
 
(*  Contour plot stress distributions *)
    
ContourPlotStresses[NodeCoordinates,ElemNodes,NodeStresses,
  {True,True,True,True},True,{},{.4*R,-H/4},aspect];\
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 18.  A rotating sphere with the rough dimensions and \
approximate properties of Mother Earth. 
Kg-force/m/sec unit system. The Earth spins with angular velocity  \[Omega]=2\
\[Pi] rad/(24hr)=
2\[Pi]/86400 (rad/sec) about z.  The planet radius is R= 6370Km=6370000m.  
For E take 1/3 of the rigidity of steel =(21000 Kgf/mm^2)/3=(7 10^9 Kgf/m^2, \
as recommended
in Love's Theory of Elasticity, Poisson's ratio \[Nu]=0.3 (my own guess), and \
mass density \[Rho]=5.52 
times the water density.  Watch out for units: the centrifugal body force \
\[Rho]\[Omega]^2r  should 
come up in Kg/m^3.  Very coarse mesh: 4 elements, constructed by eyeballing. \
\
\>", "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[Em,\[Nu],a,b,h,Kfac,\[Rho],\[Gamma],g,\[Omega]];   
Em=N[100]; \[Nu]=0; R=1.; \[Gamma]=5.52*1000; g=9.81; P=10000.
\[Rho]=\[Gamma]/g; \[Omega]=0*N[2*Pi/(3600*24)]; Kfac=1;  aspect=1; \
etype=\"Quad4\";
Print[\"Actual polar flattening of Earth =\",0.001*R/303,\" Km\"];

(*  Define FEM model *)

NodeCoordinates={{0,0},{0.5,0},{1,0},{0,0.5},{0.4,0.4},
                 {0.916,0.4},{0.71,0.71},{0.4,0.916},{0,1}}*R;        
ElemNodes={{1,2,5,4},{2,3,6,5},{4,5,8,9},{5,6,7,8}};
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes];
ElemType=       Table[etype,{numele}]; 
ElemMaterial=   Table[{Em,\[Nu]},{numele}]; 
ElemBodyForces= Table[{0,0},{numele}];
FreedomTags=    Table[{0,0},{numnod}];  
FreedomValues=  Table[{0,0},{numnod}]; 
FreedomTags[[1]]={1,1};FreedomTags[[2]]=FreedomTags[[3]]={0,1};
FreedomTags[[4]]=FreedomTags[[9]]={1,0};
FreedomValues[[3]]={-P,0};
FreedomValues[[9]]={0,-P}/(2*N[Pi]);
For [e=1,e<=numele,e++, enl=ElemNodes[[e]];
     ncoor=Table[NodeCoordinates[[enl[[i]]]],{i,4}];
     {{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;
      ElemBodyForces[[e]]=\[Rho]*\[Omega]^2*{{r1,0},{r2,0},{r3,0},{r4,0}};
    ]; 
Defaultoptions={True};
PrintRingNodeCoordinates[NodeCoordinates,\"Node Coordinates\",{2,4}];
PrintRingElementTypeNodes[ElemType,ElemNodes,\"Element type & nodes\",{}];
PrintRingElementMaterials[ElemMaterial,\"Material data\",{}]; 
PrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{4,2}];
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"Spinning Mother Earth - very coarse mesh\",True,True];

(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[
    NodeCoordinates,ElemType,ElemNodes,
    ElemMaterial,ElemBodyForces,ElemTractionForces,
    FreedomTags,FreedomValues,Defaultoptions];
PrintRingAnalysisSolution[0.001*NodeDisplacements,
     10.^(-6)*NodeForces, 10.^(-6)*NodeStresses,
    \"Solution: disp in Km, forces in MN, stresses in MPa\",{}]; 
    
(*  Contour plot stress distributions *)
  
ContourPlotStresses[NodeCoordinates,ElemNodes,
  10.^(-6)*NodeStresses,
  {True,True,True,True},True,{},{R/3,R/2},aspect]; 
  \
\>", "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[StyleBox["18a. More refined mesh, from a student HW",
  FontSize->14]], "Text",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ClearAll[Em,\[Nu],a,b,h,Kfac,\[Rho],\[Gamma],g,\[Omega]];   \nEm=10^7; \
\[Nu]=0.3; R=1.; P=100000;\n aspect=1; etype=\"Quad4\";\n \n(*  Define FEM \
model *)\n \nNodeCoordinates={{0,0},{0,1}/5,{Cos[3*Pi/8],Sin[3*Pi/8]}/5,\n\
{Cos[Pi/4],Sin[Pi/4]}/5,{Cos[Pi/8],Sin[Pi/8]}/5,{1,0}/5,{0,1}*2/5,\n\
{Cos[3*Pi/8],Sin[3*Pi/8]}*2/5,{Cos[Pi/4],Sin[Pi/4]}*2/5,{Cos[Pi/8],Sin[Pi/8]}*\
2/5,\n{1,0}*2/5,{0,1}*3/5,{Cos[3*Pi/8],Sin[3*Pi/8]}*3/5,{Cos[Pi/4],Sin[Pi/4]}*\
3/5,\n{Cos[Pi/8],Sin[Pi/8]}*3/5,{1,0}*3/5,{0,1}*4/5,{Cos[3*Pi/8],Sin[3*Pi/8]}*\
4/5,\n{Cos[Pi/4],Sin[Pi/4]}*4/5,{Cos[Pi/8],Sin[Pi/8]}*4/5,{1,0}*4/5,{0,1},\n\
{Cos[3*Pi/8],Sin[3*Pi/8]},{Cos[Pi/4],Sin[Pi/4]},{Cos[Pi/8],Sin[Pi/8]},{1,0}}*\
R;        \nElemNodes=Table[{0,0,0,0},{n,1,18}];\nElemNodes[[1]]={1,4,3,2};\n\
ElemNodes[[2]]={1,6,5,4};\nFor[n=3,n<=6,n++,ElemNodes[[n]]={n-1,n,n+4,n+5}];\n\
For[n=7,n<=10,n++,ElemNodes[[n]]={n,n+5,n+6,n+1}];\n\
For[n=11,n<=14,n++,ElemNodes[[n]]={n+1,n+6,n+7,n+2}];\n\
For[n=15,n<=18,n++,ElemNodes[[n]]={n+2,n+7,n+8,n+3}];\n\
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes];\nElemType=       \
Table[etype,{numele}]; \nElemMaterial=   Table[{Em,\[Nu]},{numele}]; \n\
ElemBodyForces= Table[{0,0},{numele}];\nElemTractionForces={};\nFreedomTags=  \
  Table[{0,0},{numnod}];  \nFreedomValues=  Table[{0,0},{numnod}]; \n\
FreedomValues[[22]]={0,N[-P/(2*Pi)]}; \n\
FreedomTags[[1]]={1,1};FreedomTags[[2]]=FreedomTags[[7]]=FreedomTags[[12]]=\
FreedomTags[[17]]=FreedomTags[[22]]={1,0};\n\
FreedomTags[[6]]=FreedomTags[[11]]=FreedomTags[[16]]=FreedomTags[[21]]=\
FreedomTags[[26]]={0,1};\nFor [e=1,e<=numele,e++, enl=ElemNodes[[e]];\n     \
ncoor=Table[NodeCoordinates[[enl[[i]]]],{i,4}];\n     \
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;\n      \
ElemBodyForces[[e]]={{0,0},{0,0},{0,0},{0,0}};\n    ]; \n\
Defaultoptions={True};\nPrintRingNodeCoordinates[NodeCoordinates,\"Node \
Coordinates\",{2,4}];\nPrintRingElementTypeNodes[ElemType,ElemNodes,\"Element \
type & nodes\",{}];\nPrintRingElementMaterials[ElemMaterial,\"Material \
data\",{}]; \nPrintRingElemBodyForces[ElemBodyForces,\"Body forces\",{4,2}];\n\
PrintRingElemTracForces[ElemTractionForces,\"Traction forces\",{}];\n\
PrintRingFreedomActivity[FreedomTags,FreedomValues,\"BC data\",{}];\n\
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,\n    \"Spinning \
Mother Earth - very coarse mesh\",True,True];\n \n(*  Solve problem and print \
results *)\n \n\
{NodeDisplacements,NodeForces,NodeStresses}=RingAnalysisDriver[\n    \
NodeCoordinates,ElemType,ElemNodes,\n    \
ElemMaterial,ElemBodyForces,ElemTractionForces,\n    \
FreedomTags,FreedomValues,Defaultoptions];\n\
PrintRingAnalysisSolution[NodeDisplacements,\n     NodeForces, NodeStresses,\n\
    \"Solution: \",{}]; \n    \n(*  Contour plot stress distributions *)\n  \n\
ContourPlotStresses[NodeCoordinates,ElemNodes,\n  10.^(-6)*NodeStresses,\n  \
{True,True,True,True},True,{},{R/3,R/2},aspect]; ",
  StyleBox["\n ",
    FontFamily->"Verdana",
    FontWeight->"Plain"],
  "\n  ",
  "\n  "
}], "Input",
  CellFrame->True,
  CellMargins->{{15, 80}, {Inherited, Inherited}},
  CellLabelMargins->{{7, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]]
},
FrontEndVersion->"4.2 for Macintosh",
ScreenRectangle->{{0, 1920}, {0, 1180}},
AutoGeneratedPackage->None,
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{1296, 1152},
WindowMargins->{{125, Automatic}, {4, Automatic}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
Magnification->1.5,
MacintoshSystemPageSetup->"\<\
00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001
0000I00000400`<300000BL?00400@00000000000000060801T1T00000000000
00000000000000000000000000000000\>"
]

(*******************************************************************
Cached data follows.  If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of  the file.  The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 1091, 21, 333, "Text"],
Cell[2848, 74, 362, 10, 71, "Text"],
Cell[3213, 86, 373, 10, 102, "Input",
  InitializationCell->True],
Cell[3589, 98, 353, 10, 62, "Input",
  InitializationCell->True],
Cell[3945, 110, 795, 21, 261, "Text",
  InitializationCell->True],
Cell[4743, 133, 2636, 62, 1142, "Input",
  InitializationCell->True],
Cell[7382, 197, 365, 10, 93, "Text"],
Cell[7750, 209, 7913, 141, 2962, "Input",
  InitializationCell->True],
Cell[15666, 352, 368, 10, 93, "Text"],
Cell[16037, 364, 10259, 167, 3882, "Input",
  InitializationCell->True],
Cell[26299, 533, 445, 11, 117, "Text"],
Cell[26747, 546, 3171, 58, 1262, "Input",
  InitializationCell->True],
Cell[29921, 606, 405, 10, 93, "Text"],
Cell[30329, 618, 1391, 24, 602, "Input",
  InitializationCell->True],
Cell[31723, 644, 483, 13, 117, "Text",
  InitializationCell->True],
Cell[32209, 659, 2409, 40, 1002, "Input",
  InitializationCell->True],
Cell[34621, 701, 420, 10, 93, "Text"],
Cell[35044, 713, 3333, 61, 1202, "Input",
  InitializationCell->True],
Cell[38380, 776, 341, 8, 93, "Text"],
Cell[38724, 786, 3685, 64, 1542, "Input",
  InitializationCell->True],
Cell[42412, 852, 358, 9, 93, "Text"],
Cell[42773, 863, 16677, 281, 6122, "Input",
  InitializationCell->True],
Cell[59453, 1146, 245, 5, 69, "Text"],
Cell[59701, 1153, 8677, 196, 3562, "Input",
  InitializationCell->True],
Cell[68381, 1351, 263, 7, 69, "Text"],
Cell[68647, 1360, 3615, 60, 1562, "Input",
  InitializationCell->True],
Cell[72265, 1422, 279, 6, 69, "Text",
  InitializationCell->True],
Cell[72547, 1430, 3700, 60, 922, "Input",
  InitializationCell->True],
Cell[76250, 1492, 450, 11, 117, "Text"],
Cell[76703, 1505, 4828, 84, 1402, "Input",
  InitializationCell->True],
Cell[81534, 1591, 322, 8, 69, "Text"],
Cell[81859, 1601, 7318, 132, 2402, "Input",
  InitializationCell->True],
Cell[89180, 1735, 292, 8, 69, "Text"],
Cell[89475, 1745, 2091, 43, 642, "Input",
  InitializationCell->True],
Cell[91569, 1790, 393, 10, 93, "Text",
  InitializationCell->True],
Cell[91965, 1802, 2572, 60, 1102, "Input"],
Cell[94540, 1864, 395, 10, 93, "Text"],
Cell[94938, 1876, 2748, 63, 1162, "Input"],
Cell[97689, 1941, 423, 10, 93, "Text"],
Cell[98115, 1953, 2781, 49, 1162, "Input"],
Cell[100899, 2004, 507, 12, 117, "Text"],
Cell[101409, 2018, 3341, 55, 1322, "Input"],
Cell[104753, 2075, 495, 12, 117, "Text"],
Cell[105251, 2089, 3234, 55, 1322, "Input"],
Cell[108488, 2146, 262, 7, 69, "Text"],
Cell[108753, 2155, 225, 8, 82, "Input"],
Cell[108981, 2165, 293, 8, 69, "Text"],
Cell[109277, 2175, 2230, 51, 942, "Input"],
Cell[111510, 2228, 879, 19, 213, "Text"],
Cell[112392, 2249, 2499, 58, 1062, "Input"],
Cell[114894, 2309, 292, 6, 72, "Text"],
Cell[115189, 2317, 3248, 54, 1326, "Input"]
}
]
*)



(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)