(************** 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[    133196,       2596]*)
(*NotebookOutlinePosition[    134293,       2631]*)
(*  CellTagsIndexPosition[    134249,       2627]*)
(*WindowFrame->Normal*)



Notebook[{
Cell[TextData[{
  StyleBox["To use this program here are the instructions:\n1.  Execute Cells \
0 through 10, in that order to initialize all necessary modules. \n     \
Shortcut: select Kernel ->Evaluation->Evaluate Initialization  \n",
    FontSize->14],
  StyleBox["      Note: some cells in the range 1-10 are not set as \
initialization cells, since they are\n      not needed for the assigned \
problems.",
    FontSize->14,
    FontSlant->"Italic"],
  StyleBox["\n2.  Cells 11 and 12 contain scripts for running the  \
one-element test problems of Figure 27.1. \n    To get confidence in the \
dowloaded file, run them  and check that the answers \n     (displacements \
and stresses) are OK. The contour plots should show only one color.\n4.  If \
the plots produced by  running a main program appear too small,click on the \
top \n     one   (only)  with the mouse,grap a corner \"handle\" and enlarge \
it.  Then rerun \n     the  main program. \n5.  Cell 0 contains ",
    FontSize->14],
  StyleBox["Mathematica",
    FontSize->14,
    FontSlant->"Italic"],
  StyleBox[" version dependent settings",
    FontSize->14]
}], "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {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->{{11, 35}, {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->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, 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->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, 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 p=1, ... 5.
See Ch 23 for argument rules. For triangles the rules are p=1,3,-3,6,-6,7;  \
see Ch 24.

This code is used by all quadrilateral  stiffness and stress modules,and all \
triangle modules 
except the 3-node linear element. Hence it must be initialized before their \
appearance.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, 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->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2A:  Modules for 4-node bilinear quad in plane stress, 
Uses QuadGaussRuleInfo,which must have been initialized. \
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad4IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= \n  \
Module[{i,k,p=2,numer=False,h=th,qcoor,c,w,Nf,\n    \
dNx,dNy,Jdet,Be,Ke=Table[0,{8},{8}]},  \n  If [Length[options]==1, \
{numer}=options];\n  If [Length[options]==2, {numer,p}=options];\n  If \
[p<1||p>5, Print[\"p out of range\"]; Return[K]];\n  For [k=1, k<=p*p, k++, \n\
       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];\n       \
{Nf,dNx,dNy,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];\n       If \
[!numer,{Nf,dNx,dNy,Jdet}=Simplify[{Nf,dNx,dNy,Jdet}]];\n       If \
[Length[th]==4, h=th.Nf]; c=w*Jdet*h;\n       Be= {Flatten[Table[{dNx[[i]],   \
    0},{i,4}]],\n            Flatten[Table[{0,       dNy[[i]]},{i,4}]],\n     \
       Flatten[Table[{dNy[[i]],dNx[[i]]},{i,4}]]}; \n       \
Ke+=c*Transpose[Be].(Emat.Be);   \n      ]; If [!numer,Ke=Simplify[Ke]]; \
Return[Ke]\n   ];\n   \n\
Quad4IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= \n  \
Module[{i,k,numer=False,qcoor,Nf,dNx,dNy,Jdet,Be,\n  \
qctab,ue=udis,sige=Table[{0,0,0},{4}]},  \n  \
qctab={{-1,-1},{1,-1},{1,1},{-1,1}}; \n  If [Length[options]>0, \
numer=options[[1]]]; \n  If [Length[udis]==4, ue=Flatten[udis]]; \n  For \
[k=1, k<=Length[sige], k++, \n       qcoor=qctab[[k]]; If [numer, \
qcoor=N[qcoor]]; \n       \
{Nf,dNx,dNy,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];\n       If \
[!numer,{Nf,dNx,dNy,Jdet}=Simplify[{Nf,dNx,dNy,Jdet}]];\n       Be= \
{Flatten[Table[{dNx[[i]],       0},{i,4}]],\n            Flatten[Table[{0,    \
   dNy[[i]]},{i,4}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,4}]]};",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n       sige[[k]]=Emat.(Be.ue); \n      ]; If \
[!numer,sige=Simplify[sige]]; Return[sige]\n   ];\n      \n\
Quad4IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  \
{Nf,dNx,dNy,dN\[Xi],dN\[Eta],i,J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  {\
\[Xi],\[Eta]}=qcoor; \n   Nf={(1-\[Xi])*(1-\[Eta]),(1+\[Xi])*(1-\[Eta]),(1+\
\[Xi])*(1+\[Eta]),(1-\[Xi])*(1+\[Eta])}/4;\n   dN\[Xi] ={-(1-\[Eta]), (1-\
\[Eta]),(1+\[Eta]),-(1+\[Eta])}/4;\n   dN\[Eta]= {-(1-\[Xi]),-(1+\[Xi]),(1+\
\[Xi]), (1-\[Xi])}/4;\n   x=Table[ncoor[[i,1]],{i,4}]; \
y=Table[ncoor[[i,2]],{i,4}];\n   J11=dN\[Xi].x; J12=dN\[Xi].y; \
J21=dN\[Eta].x; J22=dN\[Eta].y;\n   Jdet=Simplify[J11*J22-J12*J21];\n   dNx=( \
J22*dN\[Xi]-J12*dN\[Eta])/Jdet; dNy=(-J21*dN\[Xi]+J11*dN\[Eta])/Jdet;\n   \
Return[{Nf,dNx,dNy,Jdet}]\n];\n",
  StyleBox[" \nClearAll[Em,\[Nu],a,b,e,h,p,numer,exx,eyy,gxy,x,y];  h=1; a=2; \
b=1;\nncoor={{0,0},{a,0},{a,b},{0,b}};\nEm=96; \[Nu]=1/3; \
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\n\
Print[\"Emat=\",Emat//MatrixForm]; numer=False; p=2;\n\
Ke=Quad4IsoPMembraneStiffness[ncoor,Emat,h,{numer,p}];\n\
Ke=Simplify[Chop[Ke]];  Print[\"Ke=\",Ke//MatrixForm]; \nPrint[\"Eigenvalues \
of Ke=\",Chop[Eigenvalues[N[Ke]],.0000001]];\n\
xx=Table[{ncoor[[i,1]],0},{i,4}]; yy=Table[{0,ncoor[[i,2]]},{i,4}];\n\
xy=Table[{0,ncoor[[i,1]]},{i,4}]; yx=Table[{ncoor[[i,2]],0},{i,4}];\n\
ue=Flatten[xx*exx+yy*eyy+(xy+yx)*gxy/2]; Emat=IdentityMatrix[3];\n\
sige=Quad4IsoPMembraneStresses[ncoor,Emat,h,ue,{False}];\nPrint[\"node strain \
chk=\",Simplify[Transpose[sige]]//MatrixForm]; ",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2B:  Modules for 9-node biquadratic quad in plane stress,
Uses QuadGaussRuleInfo,which must have been initialized.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad9IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= \n  \
Module[{i,k,p=3,numer=False,h=th,qcoor,c,w,Nf,\n    \
dNx,dNy,Jdet,Be,Ke=Table[0,{18},{18}]},  \n  If [Length[options]>0, \
numer=options[[1]]]; \n  If [Length[options]==2, {numer,p}=options];\n  If \
[p<1||p>5, Print[\"p out of range\"]; Return[K]];\n  For [k=1, k<=p*p, k++,  \
\n       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];\n       \
{Nf,dNx,dNy,Jdet}=Quad9IsoPShapeFunDer[ncoor,qcoor];\n        If \
[Length[th]==9, h=th.Nf]; c=w*Jdet*h;\n        Be={Flatten[Table[{dNx[[i]],   \
    0},{i,9}]],\n            Flatten[Table[{0,       dNy[[i]]},{i,9}]],\n     \
       Flatten[Table[{dNy[[i]],dNx[[i]]},{i,9}]]}; \n        \
Ke+=c*Transpose[Be].(Emat.Be);   \n      ]; If [!numer,Ke=Simplify[Ke]]; \
Return[Ke]\n   ];   \n    \n\
Quad9IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= \n  \
Module[{i,k,numer=False,qcoor,Nf,dNx,dNy,Jdet,Be,\n    \
qctab,ue=udis,sige=Table[{0,0,0},{9}]},  \n  \
qctab={{-1,-1},{1,-1},{1,1},{-1,1},{0,-1},{1,0},{0,1},{-1,0},{0,0}}; \n  If \
[Length[options]==1, {numer}=options];\n  If [Length[options]==2, \
{numer,p}=options];\n  If [Length[udis]==9, ue=Flatten[udis]];\n  For [k=1, \
k<=Length[sige], k++, \n       qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; \
\n       {Nf,dNx,dNy,Jdet}=Quad9IsoPShapeFunDer[ncoor,qcoor];\n        \
Be={Flatten[Table[{dNx[[i]],       0},{i,9}]],\n            Flatten[Table[{0, \
      dNy[[i]]},{i,9}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,9}]]}; \n        sige[[k]]=Emat.(Be.ue); \
  \n      ]; If [!numer,sige=Simplify[sige]]; Return[sige]\n   ];  \n\n\
Quad9IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  \
{Nf,dNx,dNy,dN\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  \
{\[Xi],\[Eta]}=qcoor;  \n   Nf= \
{(1-\[Xi])*\[Xi]*(1-\[Eta])*\[Eta]/4,-(1+\[Xi])*\[Xi]*(1-\[Eta])*\[Eta]/4,\n  \
      (1+\[Xi])*\[Xi]*(1+\[Eta])*\[Eta]/4,-(1-\[Xi])*\[Xi]*(1+\[Eta])*\[Eta]/\
4,\n       -(1-\[Xi]^2)*(1-\[Eta])*\[Eta]/2, (1+\[Xi])*\[Xi]*(1-\[Eta]^2)/2,\n\
        (1-\[Xi]^2)*(1+\[Eta])*\[Eta]/2,-(1-\[Xi])*\[Xi]*(1-\[Eta]^2)/2,(1-\
\[Xi]^2)*(1-\[Eta]^2)};\n   dN\[Xi]={(1-2*\[Xi])*\[Eta]*(1-\[Eta])/4,-(1+2*\
\[Xi])*\[Eta]*(1-\[Eta])/4,\n        (1+2*\[Xi])*\[Eta]*(1+\[Eta])/4,-(1-2*\
\[Xi])*\[Eta]*(1+\[Eta])/4,\n         \[Xi]*\[Eta]*(1-\[Eta]), \
(1/2+\[Xi])*(1-\[Eta]^2),\n        -\[Xi]*\[Eta]*(1+\[Eta]),-(1/2-\[Xi])*(1-\
\[Eta]^2), 2*\[Xi]*(\[Eta]^2-1)};     \n   dN\[Eta]={\[Xi]*(1-\[Xi])*(1-2*\
\[Eta])/4, -\[Xi]*(1+\[Xi])*(1-2*\[Eta])/4,\n        \[Xi]*(1+\[Xi])*(1+2*\
\[Eta])/4, -\[Xi]*(1-\[Xi])*(1+2*\[Eta])/4,\n        \
-(1/2-\[Eta])*(1-\[Xi]^2),  -\[Xi]*(1+\[Xi])*\[Eta],         \n         (1/2+\
\[Eta])*(1-\[Xi]^2),   \[Xi]*(1-\[Xi])*\[Eta],2*(\[Xi]^2-1)*\[Eta]};\n   \
x=Table[ncoor[[i,1]],{i,9}]; y=Table[ncoor[[i,2]],{i,9}];\n   J11=dN\[Xi].x; \
J21=dN\[Xi].y; J12=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx= ( J22*dN\[Xi]-J21*dN\[Eta])/Jdet;  \
dNx=Simplify[dNx];\n   dNy= (-J12*dN\[Xi]+J11*dN\[Eta])/Jdet;  \
dNy=Simplify[dNy];\n   Return[{Nf,dNx,dNy,Jdet}]\n];\n",
  StyleBox["\nClearAll[Em,\[Nu],a,b,e,h,p,numer,exx,eyy,gxy];  h=1; a=2; b=1;\
\nncoor={{0,0},{a,0},{a,b},{0,b},{a/2,0},{a,b/2},{a/2,b},{0,b/2},{a/2,b/2}};\n\
Em=96*39*11*55*7; \[Nu]=1/3;  \n\
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\n\
Print[\"Emat=\",Emat//MatrixForm];\np=2;  numer=False;  \n\
Ke=Quad9IsoPMembraneStiffness[ncoor,Emat,h,{numer,p}];\n\
Ke=Simplify[Chop[Ke]];  Print[\"Ke=\",Ke//MatrixForm]; \nPrint[\"Eigenvalues \
of Ke=\",Chop[Eigenvalues[N[Ke]],.00001]];\nxx=Table[{ncoor[[i,1]],0},{i,9}]; \
yy=Table[{0,ncoor[[i,2]]},{i,9}];\nxy=Table[{0,ncoor[[i,1]]},{i,9}]; \
yx=Table[{ncoor[[i,2]],0},{i,9}];\nue=Flatten[xx*exx+yy*eyy+(xy+yx)*gxy/2]; \
Emat=IdentityMatrix[3];\n\
sige=Quad9IsoPMembraneStresses[ncoor,Emat,h,ue,{False}];\nPrint[\"strain chk=\
\",Simplify[Transpose[sige]]//MatrixForm];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2C:  Modules for 8-node serendipity quad in plane stress,
Uses QuadGaussRuleInfo,which must have been initialized.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad8IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= \n  \
Module[{i,k,p=3,numer=False,h=th,qcoor,c,w,Nf,\n    \
dNx,dNy,Jdet,Be,Ke=Table[0,{16},{16}]},  \n  If [Length[options]==1, \
{numer}=options];\n  If [Length[options]==2, {numer,p}=options];\n  If \
[p<1||p>5, Print[\"p out of range\"]; Return[K]];\n  For [k=1, k<=p*p, k++,  \
\n       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];\n       \
{Nf,dNx,dNy,Jdet}=Quad8IsoPShapeFunDer[ncoor,qcoor];\n        If \
[Length[th]==8, h=th.Nf]; c=w*Jdet*h;\n        Be={Flatten[Table[{dNx[[i]],   \
    0},{i,8}]],\n            Flatten[Table[{0,       dNy[[i]]},{i,8}]],\n     \
       Flatten[Table[{dNy[[i]],dNx[[i]]},{i,8}]]}; \n        \
Ke+=c*Transpose[Be].(Emat.Be);   \n      ]; If [!numer,Ke=Simplify[Ke]]; \
Return[Ke]\n   ];   \n    \n\
Quad8IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= \n  \
Module[{i,k,numer=False,qcoor,Nf,dNx,dNy,Jdet,Be,\n  \
qctab,ue=udis,sige=Table[{0,0,0},{8}]},  \n  \
qctab={{-1,-1},{1,-1},{1,1},{-1,1},{0,-1},{1,0},{0,1},{-1,0}}; \n  If \
[Length[options]>0, numer=options[[1]]]; \n  If [Length[udis]==8, \
ue=Flatten[udis]];\n  For [k=1, k<=Length[sige], k++, \n       \
qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; \n       \
{Nf,dNx,dNy,Jdet}=Quad8IsoPShapeFunDer[ncoor,qcoor];\n        \
Be={Flatten[Table[{dNx[[i]],       0},{i,8}]],\n            Flatten[Table[{0, \
      dNy[[i]]},{i,8}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,8}]]}; \n        sige[[k]]=Emat.(Be.ue); \
  \n      ]; If [!numer,sige=Simplify[sige]]; Return[sige]\n   ];  \n\n\
Quad8IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  \
{Nf,dNx,dNy,dN\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  \
{\[Xi],\[Eta]}=qcoor;  \n   Nf= \
{-(1-\[Xi])*(1-\[Eta])*(1+\[Xi]+\[Eta])/4,-(1+\[Xi])*(1-\[Eta])*(1-\[Xi]+\
\[Eta])/4,\n        -(1+\[Xi])*(1+\[Eta])*(1-\[Xi]-\[Eta])/4,-(1-\[Xi])*(1+\
\[Eta])*(1+\[Xi]-\[Eta])/4,\n         (1+\[Xi])*(1-\[Xi])*(1-\[Eta])/2, (1+\
\[Xi])*(1+\[Eta])*(1-\[Eta])/2,\n         (1+\[Xi])*(1-\[Xi])*(1+\[Eta])/2, \
(1-\[Xi])*(1+\[Eta])*(1-\[Eta])/2};\n   dN\[Xi]={ \
(1-\[Eta])*(2*\[Xi]+\[Eta])/4,   (1-\[Eta])*(2*\[Xi]-\[Eta])/4, \n         \
(1+\[Eta])*(2*\[Xi]+\[Eta])/4,   (1+\[Eta])*(2*\[Xi]-\[Eta])/4, \n         \
\[Xi]*(\[Eta]-1),  (1+\[Eta])*(1-\[Eta])/2, -\[Xi]*(1+\[Eta]), \
-(1+\[Eta])*(1-\[Eta])/2};\n   dN\[Eta]={ (1-\[Xi])*(\[Xi]+2*\[Eta])/4,   \
-(1+\[Xi])*(\[Xi]-2*\[Eta])/4, \n         (1+\[Xi])*(\[Xi]+2*\[Eta])/4,   \
-(1-\[Xi])*(\[Xi]-2*\[Eta])/4, \n        -(1+\[Xi])*(1-\[Xi])/2, -(1+\[Xi])*\
\[Eta], (1+\[Xi])*(1-\[Xi])/2,-(1-\[Xi])*\[Eta]};\n   \
x=Table[ncoor[[i,1]],{i,8}]; y=Table[ncoor[[i,2]],{i,8}];\n   J11=dN\[Xi].x; \
J21=dN\[Xi].y; J12=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx=( J22*dN\[Xi]-J21*dN\[Eta])/Jdet; \
dNy=(-J12*dN\[Xi]+J11*dN\[Eta])/Jdet;  \n   Return[{Nf,dNx,dNy,Jdet}]\n];\n",
  StyleBox["\nClearAll[Em,\[Nu],a,b,e,h,p,numer,exx,eyy,gxy];  h=1; a=2; b=1;\
\nncoor={{0,0},{a,0},{a,b},{0,b},{a/2,0},{a,b/2},{a/2,b},{0,b/2}};\nEm=96*5; \
\[Nu]=1/3; \nEmat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\
\nPrint[\"Emat=\",Emat//MatrixForm];\np=3;  numer=True;  \n\
Ke=Quad8IsoPMembraneStiffness[ncoor,Emat,h,{numer,p}];\n\
Ke=Simplify[Chop[Ke]];  Print[\"Ke=\",Ke//MatrixForm]; \nPrint[\"Eigenvalues \
of Ke=\",Chop[Eigenvalues[N[Ke]],.000001]];\n\
xx=Table[{ncoor[[i,1]],0},{i,8}]; yy=Table[{0,ncoor[[i,2]]},{i,8}];\n\
xy=Table[{0,ncoor[[i,1]]},{i,8}]; yx=Table[{ncoor[[i,2]],0},{i,8}];\n\
ue=Flatten[xx*exx+yy*eyy+(xy+yx)*gxy/2]; Emat=IdentityMatrix[3];\n\
sige=Quad8IsoPMembraneStresses[ncoor,Emat,h,ue,{False}];\nPrint[\"strain chk=\
\",Simplify[Transpose[sige]]//MatrixForm];\n",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2X.  16-node bicubic quad with nodes at thirdpoints (cell not \
initialized)\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad16IsoPMembraneStiffness[ncoor_,mprop_,fprop_,options_]:= \n  \
Module[{i,j,k,kk,p=3,numer=False,Emat,th=1,h,qcoor,c,w,Nf,\n   \
dNx,dNy,Jdet,B,Ke=Table[0,{32},{32}]},  \n  Emat=mprop[[1]];\n  If \
[Length[options]==2, {numer,p}=options, {numer}=options];\n  If \
[Length[fprop]>0, th=fprop[[1]]]; Print[\"p=\",p]; \n  If [(p<1||p>5), Print[\
\"p out of range\"]; Return[Null]];\n  kk= If [p>0,p^2,Abs[p]];\n  For [k=1, \
k<=kk, k++, \n       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k]; \n       \
{Nf,dNx,dNy,Jdet}=Quad16IsoPShapeFunDer[ncoor,qcoor];\n        If \
[Length[th]==0, h=th, h=th.Nf]; c=w*Jdet*h;\n        B={ \
Flatten[Table[{dNx[[i]],       0},{i,16}]],\n            Flatten[Table[{0,    \
   dNy[[i]]},{i,16}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,16}]]}; \n  ",
  StyleBox["      (*Print[\"B=\",B//MatrixForm];*) ",
    FontColor->RGBColor[1, 0, 0]],
  "\n        Ke+=Simplify[c*Transpose[B].(Emat.B)];   \n      ]; Return[Ke];\n\
   ];",
  StyleBox["\n   \n",
    FontColor->RGBColor[0, 0, 1]],
  "Quad16IsoPMembraneStresses[ncoor_,mprop_,fprop_,options_,udis_]:= \n  \
Module[{i,k,numer=False,Emat,th=1,h,qcoor,Nf,\n    \
dNx,dNy,Jdet,B,qctab,ue=udis,sige=Table[0,{5},{3}]},  \n  \
qctab={{-1,-1},{1,-1},{1,1},{-1,1},{-1/3,-1},{1/3,-1},\n         \
{1,-1/3},{1,1/3},{1/3,1},{-1/3,1},{1/3,-1},{-1/3,-1},\n         \
{-1/3,-1/3},{1/3,-1/3},{1/3,1/3},{-1/3,1/3}}; \n  Emat=mprop[[1]]; \
numer=options[[1]]; \n  If [Length[udis]==16, ue=Flatten[udis]];\n  For [k=1, \
k<=Length[sige], k++, \n       qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; \
\n       {Nf,dNx,dNy,Jdet}=Quad16IsoPShapeFunDer[ncoor,qcoor];\n        B={ \
Flatten[Table[{dNx[[i]],       0},{i,16}]],\n            Flatten[Table[{0,    \
   dNy[[i]]},{i,16}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,16}]]}; \n        sige[[k]]=Emat.(B.ue); \
  \n      ]; Return[sige]\n   ];\n   \n   \n\
Quad16IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  \
{Nf,dNx,dNy,dN\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  \
{\[Xi],\[Eta]}=qcoor; \n   \
Nf={(1-\[Xi])*(9*\[Xi]^2-1)*(1-\[Eta])*(9*\[Eta]^2-1),\n       (1+\[Xi])*(9*\
\[Xi]^2-1)*(1-\[Eta])*(9*\[Eta]^2-1),\n       (1+\[Xi])*(9*\[Xi]^2-1)*(1+\
\[Eta])*(9*\[Eta]^2-1),\n       \
(1-\[Xi])*(9*\[Xi]^2-1)*(1+\[Eta])*(9*\[Eta]^2-1),\n       \
9*(1-\[Xi]^2)*(1-3*\[Xi])*(1-\[Eta])*(9*\[Eta]^2-1),\n       \
9*(1-\[Xi]^2)*(1+3*\[Xi])*(1-\[Eta])*(9*\[Eta]^2-1),\n       9*(1+\[Xi])*(9*\
\[Xi]^2-1)*(1-\[Eta]^2)*(1-3*\[Eta]),\n       9*(1+\[Xi])*(9*\[Xi]^2-1)*(1-\
\[Eta]^2)*(1+3*\[Eta]),\n       \
9*(1-\[Xi]^2)*(1+3*\[Xi])*(1+\[Eta])*(9*\[Eta]^2-1),\n       \
9*(1-\[Xi]^2)*(1-3*\[Xi])*(1+\[Eta])*(9*\[Eta]^2-1),\n       9*(1-\[Xi])*(9*\
\[Xi]^2-1)*(1-\[Eta]^2)*(1+3*\[Eta]),\n       9*(1-\[Xi])*(9*\[Xi]^2-1)*(1-\
\[Eta]^2)*(1-3*\[Eta]),\n       81*(1-\[Xi]^2)*(1-\[Eta]^2)*(1-3*\[Xi])*(1-3*\
\[Eta]),\n       81*(1-\[Xi]^2)*(1-\[Eta]^2)*(1+3*\[Xi])*(1-3*\[Eta]),\n      \
 81*(1-\[Xi]^2)*(1-\[Eta]^2)*(1+3*\[Xi])*(1+3*\[Eta]),\n       \
81*(1-\[Xi]^2)*(1-\[Eta]^2)*(1-3*\[Xi])*(1+3*\[Eta])}/256;\n   \
dN\[Xi]={-(1+9*\[Xi]*(2-3*\[Xi]))*(1-\[Eta])*(1-9*\[Eta]^2),\n       (1-9*\
\[Xi]*(2+3*\[Xi]))*(1-\[Eta])*(1-9*\[Eta]^2),\n       \
(1-9*\[Xi]*(2+3*\[Xi]))*(1+\[Eta])*(1-9*\[Eta]^2),\n      -(1+9*\[Xi]*(2-3*\
\[Xi]))*(1+\[Eta])*(1-9*\[Eta]^2),\n       \
9*(3+\[Xi]*(2-9*\[Xi]))*(1-\[Eta])*(1-9*\[Eta]^2),\n      -9*(3-\[Xi]*(2+9*\
\[Xi]))*(1-\[Eta])*(1-9*\[Eta]^2),\n      \
-9*(1-9*\[Xi]*(2+3*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1-3*\[Eta]),\n      \
-9*(1-9*\[Xi]*(2+3*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1+3*\[Eta]),\n      -9*(3-\
\[Xi]*(2+9*\[Xi]))*(1+\[Eta])*(1-9*\[Eta]^2),\n       \
9*(3+\[Xi]*(2-9*\[Xi]))*(1+\[Eta])*(1-9*\[Eta]^2),\n       9*(1+9*\[Xi]*(2-3*\
\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1+3*\[Eta]),\n       \
9*(1+9*\[Xi]*(2-3*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1-3*\[Eta]),\n      -81*(3+\
\[Xi]*(2-9*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1-3*\[Eta]),\n       \
81*(3-\[Xi]*(2+9*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1-3*\[Eta]),\n       81*(3-\
\[Xi]*(2+9*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1+3*\[Eta]),\n      \
-81*(3+\[Xi]*(2-9*\[Xi]))*(1-\[Eta])*(1+\[Eta])*(1+3*\[Eta])}/256;\n   dN\
\[Eta]={-(1-\[Xi])*(1-3*\[Xi])*(1+3*\[Xi])*(1+9*\[Eta]*(2-3*\[Eta])),\n       \
-(1+\[Xi])*(1+3*\[Xi])*(1-3*\[Xi])*(1+9*\[Eta]*(2-3*\[Eta])),\n        (1+\
\[Xi])*(1+3*\[Xi])*(1-3*\[Xi])*(1-9*\[Eta]*(2+3*\[Eta])),\n        \
(1-\[Xi])*(1-3*\[Xi])*(1+3*\[Xi])*(1-9*\[Eta]*(2+3*\[Eta])),\n        9*(1-\
\[Xi])*(1+\[Xi])*(1-3*\[Xi])*(1+9*\[Eta]*(2-3*\[Eta])),\n        \
9*(1-\[Xi])*(1+\[Xi])*(1+3*\[Xi])*(1+9*\[Eta]*(2-3*\[Eta])),\n        9*(1+\
\[Xi])*(1+3*\[Xi])*(1-3*\[Xi])*(3+\[Eta]*(2-9*\[Eta])),\n       \
-9*(1+\[Xi])*(1+3*\[Xi])*(1-3*\[Xi])*(3-\[Eta]*(2+9*\[Eta])),\n       -9*(1-\
\[Xi])*(1+\[Xi])*(1+3*\[Xi])*(1-9*\[Eta]*(2+3*\[Eta])),\n       \
-9*(1-\[Xi])*(1+\[Xi])*(1-3*\[Xi])*(1-9*\[Eta]*(2+3*\[Eta])),\n       -9*(1-\
\[Xi])*(1-3*\[Xi])*(1+3*\[Xi])*(3-\[Eta]*(2+9*\[Eta])),\n        \
9*(1-\[Xi])*(1-3*\[Xi])*(1+3*\[Xi])*(3+\[Eta]*(2-9*\[Eta])),\n       -81*(1-\
\[Xi])*(1+\[Xi])*(1-3*\[Xi])*(3-\[Eta]*(-2+9*\[Eta])),\n       \
-81*(1-\[Xi])*(1+\[Xi])*(1+3*\[Xi])*(3+\[Eta]*(2-9*\[Eta])),\n        81*(1-\
\[Xi])*(1+\[Xi])*(1+3*\[Xi])*(3-\[Eta]*(2+9*\[Eta])),\n        \
81*(1-\[Xi])*(1+\[Xi])*(1-3*\[Xi])*(3-\[Eta]*(2+9*\[Eta]))}/256;\n  ",
  StyleBox[" (*Print [\"chk: \",Simplify[dN\[Xi]-D[Nf,\[Xi]]]];\n   Print \
[\"chk: \",Simplify[dN\[Eta]-D[Nf,\[Eta]]]];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   x=Table[ncoor[[i,1]],{i,16}]; y=Table[ncoor[[i,2]],{i,16}];\n   \
J11=dN\[Xi].x; J12=dN\[Xi].y; J21=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx= ( J22*dN\[Xi]-J12*dN\[Eta])/Jdet;  \
dNx=Simplify[dNx];\n   dNy= (-J21*dN\[Xi]+J11*dN\[Eta])/Jdet;  \
dNy=Simplify[dNy];\n   Return[Simplify[{Nf,dNx,dNy,Jdet}]] \n]; \n"
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell  2Y.   Gauss-point  connected (GC) variant of the bicubic quad \
 (cell not initialized)\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Quad16GCIsoPMembraneStiffness[ncoor_,mprop_,fprop_,options_]:= \n  \
Module[{i,j,k,kk,p=3,numer=False,Emat,th=1,h,qcoor,c,w,Nf,\n   \
dNx,dNy,Jdet,B,Ke=Table[0,{32},{32}]},  \n  Emat=mprop[[1]];\n  If \
[Length[options]==2, {numer,p}=options, {numer}=options];\n  If \
[Length[fprop]>0, th=fprop[[1]]]; Print[\"p=\",p]; \n  If [(p<1||p>5), Print[\
\"p out of range\"]; Return[Null]];\n  kk= If [p>0,p^2,Abs[p]];\n  For [k=1, \
k<=kk, k++, \n       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k]; \n       \
{Nf,dNx,dNy,Jdet}=Quad16GCIsoPShapeFunDer[ncoor,qcoor];\n        If \
[Length[th]==0, h=th, h=th.Nf]; c=w*Jdet*h;\n        B={ \
Flatten[Table[{dNx[[i]],       0},{i,16}]],\n            Flatten[Table[{0,    \
   dNy[[i]]},{i,16}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,16}]]}; \n  ",
  StyleBox["      (*Print[\"B=\",B//MatrixForm];*) ",
    FontColor->RGBColor[1, 0, 0]],
  "\n        Ke+=Simplify[c*Transpose[B].(Emat.B)];   \n      ]; Return[Ke];\n\
   ];",
  StyleBox["\n   ",
    FontColor->RGBColor[0, 0, 1]],
  "\n  \nQuad16GCIsoPMembraneStresses[ncoor_,mprop_,fprop_,options_,udis_]:= \
\n  Module[{i,k,numer=False,Emat,\[Alpha]=Sqrt[3]/3,th=1,h,qcoor,Nf,\n    \
dNx,dNy,Jdet,B,qctab,ue=udis,sige=Table[0,{5},{3}]},  \n  \
qctab={{-1,-1},{1,-1},{1,1},{-1,1},{-\[Alpha],-1},{\[Alpha],-1},\n         \
{1,-\[Alpha]},{1,\[Alpha]},{\[Alpha],1},{-\[Alpha],1},{\[Alpha],-1},{-\[Alpha]\
,-1},\n         \
{-\[Alpha],-\[Alpha]},{\[Alpha],-\[Alpha]},{\[Alpha],\[Alpha]},{-\[Alpha],\
\[Alpha]}}; \n  Emat=mprop[[1]]; numer=options[[1]]; \n  If \
[Length[udis]==16, ue=Flatten[udis]];\n  For [k=1, k<=Length[sige], k++, \n   \
    qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; \n       \
{Nf,dNx,dNy,Jdet}=Quad16GCIsoPShapeFunDer[ncoor,qcoor];\n        B={ \
Flatten[Table[{dNx[[i]],       0},{i,16}]],\n            Flatten[Table[{0,    \
   dNy[[i]]},{i,16}]],\n            \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,16}]]}; \n        sige[[k]]=Emat.(B.ue); \
  \n      ]; Return[sige]\n   ];",
  StyleBox["\n",
    FontColor->RGBColor[0, 0, 1]],
  "\nQuad16GCIsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  \
{Nf,dNx,dNy,dN\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],cc,cm,ci,x,y},\
\n  {\[Xi],\[Eta]}=qcoor; \[Alpha]=Sqrt[3]/3; cc=1/(4*(\[Alpha]^2-1)^2);\n   \
cm=1/(4*Sqrt[3]*\[Alpha]*(\[Alpha]^2-1)^2); \
ci=1/(4*\[Alpha]^2*(\[Alpha]^2-1)^2);\n   \
Nf={(1-\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1-\[Eta])*(\[Eta]^2-\[Alpha]^2)*cc,\n     \
  (1+\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1-\[Eta])*(\[Eta]^2-\[Alpha]^2)*cc,\n       \
(1+\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*cc,\n       \
(1-\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*cc,\n       \
Sqrt[3]*(1-\[Xi]^2)*(\[Alpha]-\[Xi])*(1-\[Eta])*(\[Eta]^2-\[Alpha]^2)*cm,\n   \
    Sqrt[3]*(1-\[Xi]^2)*(\[Alpha]+\[Xi])*(1-\[Eta])*(\[Eta]^2-\[Alpha]^2)*cm,\
\n       Sqrt[3]*(1+\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1-\[Eta]^2)*(\[Alpha]-\[Eta])\
*cm,\n       Sqrt[3]*(1+\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1-\[Eta]^2)*(\[Alpha]+\
\[Eta])*cm,\n       \
Sqrt[3]*(1-\[Xi]^2)*(\[Alpha]+\[Xi])*(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*cm,\n   \
    Sqrt[3]*(1-\[Xi]^2)*(\[Alpha]-\[Xi])*(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*cm,\
\n       Sqrt[3]*(1-\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1-\[Eta]^2)*(\[Alpha]+\[Eta])\
*cm,\n       Sqrt[3]*(1-\[Xi])*(\[Xi]^2-\[Alpha]^2)*(1-\[Eta]^2)*(\[Alpha]-\
\[Eta])*cm,\n       \
(1+\[Xi])*(1-\[Xi])*(1+\[Eta])*(1-\[Eta])*(\[Alpha]-\[Xi])*(\[Alpha]-\[Eta])*\
ci,\n       (1+\[Xi])*(1-\[Xi])*(1+\[Eta])*(1-\[Eta])*(\[Alpha]+\[Xi])*(\
\[Alpha]-\[Eta])*ci,\n       (1+\[Xi])*(1-\[Xi])*(1+\[Eta])*(1-\[Eta])*(\
\[Alpha]+\[Xi])*(\[Alpha]+\[Eta])*ci,\n       \
(1+\[Xi])*(1-\[Xi])*(1+\[Eta])*(1-\[Eta])*(\[Alpha]-\[Xi])*(\[Alpha]+\[Eta])*\
ci};\n  dN\[Xi]={-(\[Eta]-1)*(\[Eta]^2-\[Alpha]^2)*(\[Alpha]^2+\[Xi]*(2-3*\
\[Xi]))*cc, \n       (\[Eta]-1)*(\[Eta]^2-\[Alpha]^2)*(\[Alpha]^2-\[Xi]*(2+3*\
\[Xi]))*cc, \n      -(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*(\[Alpha]^2-\[Xi]*(2+3*\
\[Xi]))*cc, \n       (1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*(\[Alpha]^2+\[Xi]*(2-3*\
\[Xi]))*cc, \n       Sqrt[3]*(\[Eta]-1)*(\[Eta]^2-\[Alpha]^2)*(1+2*\[Alpha]*\
\[Xi]-3*\[Xi]^2)*cm, \n      -Sqrt[3]*(\[Eta]-1)*(\[Eta]^2-\[Alpha]^2)*(1-2*\
\[Alpha]*\[Xi]-3*\[Xi]^2)*cm, \n       \
Sqrt[3]*(\[Alpha]-\[Eta])*(\[Eta]^2-1)*(\[Alpha]^2-\[Xi]*(2+3*\[Xi]))*cm, \n  \
     Sqrt[3]*(\[Alpha]+\[Eta])*(\[Eta]^2-1)*(\[Alpha]^2-\[Xi]*(2+3*\[Xi]))*cm,\
 \n       Sqrt[3]*(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*(1-2*\[Alpha]*\[Xi]-3*\[Xi]\
^2)*cm, \n      \
-Sqrt[3]*(1+\[Eta])*(\[Eta]^2-\[Alpha]^2)*(1+2*\[Alpha]*\[Xi]-3*\[Xi]^2)*cm, \
\n      -Sqrt[3]*(\[Alpha]+\[Eta])*(\[Eta]^2-1)*(\[Alpha]^2+(2-3*\[Xi])*\[Xi])\
*cm, \n      -Sqrt[3]*(\[Alpha]-\[Eta])*(\[Eta]^2-1)*(\[Alpha]^2+(2-3*\[Xi])*\
\[Xi])*cm, \n       3*(\[Eta]^2-1)*(\[Alpha]-\[Eta])*(1+2*\[Alpha]*\[Xi]-3*\
\[Xi]^2)*ci/3, \n      \
-3*(\[Eta]^2-1)*(\[Alpha]-\[Eta])*(1-2*\[Alpha]*\[Xi]-3*\[Xi]^2)*ci/3, \n     \
 -3*(\[Eta]^2-1)*(\[Alpha]+\[Eta])*(1-2*\[Alpha]*\[Xi]-3*\[Xi]^2)*ci/3, \n    \
   3*(\[Eta]^2-1)*(\[Alpha]+\[Eta])*(1+2*\[Alpha]*\[Xi]-3*\[Xi]^2)*ci/3}; \n  \
dN\[Eta]={ \
(\[Alpha]^2+\[Eta]*(2-3*\[Eta]))*(\[Alpha]^2-\[Xi]^2)*(-1+\[Xi])*cc, \n       \
-(\[Alpha]^2+\[Eta]*(2-3*\[Eta]))*(\[Alpha]^2-\[Xi]^2)*(1+\[Xi])*cc, \n       \
 (\[Alpha]^2-\[Eta]*(2+3*\[Eta]))*(\[Alpha]^2-\[Xi]^2)*(1+\[Xi])*cc, \n       \
-(\[Alpha]^2-\[Eta]*(2+3*\[Eta]))*(\[Alpha]^2-\[Xi]^2)*(-1+\[Xi])*cc, \n      \
  Sqrt[3]*(\[Alpha]^2+(2-3*\[Eta])*\[Eta])*(\[Alpha]-\[Xi])*(1-\[Xi]^2)*cm, \n\
        Sqrt[3]*(\[Alpha]^2+(2-3*\[Eta])*\[Eta])*(\[Alpha]+\[Xi])*(1-\[Xi]^2)*\
cm, \n        \
Sqrt[3]*(1+2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\[Alpha]^2-\[Xi]^2)*(1+\[Xi])*cm, \n\
       -Sqrt[3]*(1-2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\[Alpha]^2-\[Xi]^2)*(1+\[Xi]\
)*cm, \n        Sqrt[3]*(\[Alpha]^2-\[Eta]*(2+3*\[Eta]))*(\[Alpha]+\[Xi])*(\
\[Xi]^2-1)*cm, \n        Sqrt[3]*(\[Alpha]^2-\[Eta]*(2+3*\[Eta]))*(\[Alpha]-\
\[Xi])*(\[Xi]^2-1)*cm, \n        Sqrt[3]*(1-2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\
\[Alpha]^2-\[Xi]^2)*(-1+\[Xi])*cm, \n       -Sqrt[3]*(1+2*\[Alpha]*\[Eta]-3*\
\[Eta]^2)*(\[Alpha]^2-\[Xi]^2)*(-1+\[Xi])*cm, \n        \
3*(1+2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\[Xi]^2-1)*(\[Alpha]-\[Xi])*ci/3, \n      \
  3*(1+2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\[Xi]^2-1)*(\[Alpha]+\[Xi])*ci/3, \n    \
   -3*(1-2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\[Xi]^2-1)*(\[Alpha]+\[Xi])*ci/3, \n  \
     -3*(1-2*\[Alpha]*\[Eta]-3*\[Eta]^2)*(\[Xi]^2-1)*(\[Alpha]-\[Xi])*ci/3}; \
\n",
  StyleBox["   (*Print [\"chk: \",Simplify[dN\[Xi]-D[Nf,\[Xi]]]];\n   Print [\
\"chk: \",Simplify[dN\[Eta]-D[Nf,\[Eta]]]];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   x=Table[ncoor[[i,1]],{i,16}]; y=Table[ncoor[[i,2]],{i,16}];\n   \
J11=dN\[Xi].x; J12=dN\[Xi].y; J21=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx= ( J22*dN\[Xi]-J12*dN\[Eta])/Jdet;  \
dNx=Simplify[dNx];\n   dNy= (-J21*dN\[Xi]+J11*dN\[Eta])/Jdet;  \
dNy=Simplify[dNy];\n   Return[Simplify[{Nf,dNx,dNy,Jdet}]] \n ];",
  StyleBox["\n",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[StyleBox["Cell SF. Enable for a quick of all test shape \
function modules, or selected ones (cell not initialized)", "Text"]], "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["ClearAll[a,b,\[Xi],\[Eta],exx,eyy,gxy];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\nQuad4IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  {Nf,dNx,dNy,dN\
\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  \
{\[Xi],\[Eta]}=qcoor;  \n   Nf={(1-\[Xi])*(1-\[Eta]),(1+\[Xi])*(1-\[Eta]),(1+\
\[Xi])*(1+\[Eta]),(1-\[Xi])*(1+\[Eta])}/4;\n   dN\[Xi] ={-(1-\[Eta]), (1-\
\[Eta]),(1+\[Eta]),-(1+\[Eta])}/4;\n   dN\[Eta]= {-(1-\[Xi]),-(1+\[Xi]),(1+\
\[Xi]), (1-\[Xi])}/4;\n",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["   Print[Simplify[dN\[Xi]-D[Nf,\[Xi]]]]; \
Print[Simplify[dN\[Eta]-D[Nf,\[Eta]]]];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   x=Table[ncoor[[i,1]],{i,4}]; y=Table[ncoor[[i,2]],{i,4}];\n  \
 J11=dN\[Xi].x; J21=dN\[Xi].y; J12=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx=( J22*dN\[Xi]-J21*dN\[Eta])/Jdet; \
dNy=(-J12*dN\[Xi]+J11*dN\[Eta])/Jdet;\n   Return[{Nf,dNx,dNy,Jdet}]\n];\n\n",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["ncoor={{0,0},{2*a,0},{a,b},{0,b}};\n\
{Nf,dNx,dNy,Jdet}=Simplify[Quad4IsoPShapeFunDer[ncoor,{\[Xi],\[Eta]}]];\n\
Print[{Nf,dNx,dNy,Jdet}];\nx=Table[ncoor[[i,1]],{i,4}]; \
y=Table[ncoor[[i,2]],{i,4}];\n\
Print[Simplify[{dNx.(x*exx),dNy.(y*eyy),dNx.(y*gxy)+dNx.(x*gxy)}]];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n\nQuad9IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  \
{Nf,dNx,dNy,dN\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  \
{\[Xi],\[Eta]}=qcoor;  \n   Nf= \
{(1-\[Xi])*\[Xi]*(1-\[Eta])*\[Eta]/4,-(1+\[Xi])*\[Xi]*(1-\[Eta])*\[Eta]/4,\n  \
      (1+\[Xi])*\[Xi]*(1+\[Eta])*\[Eta]/4,-(1-\[Xi])*\[Xi]*(1+\[Eta])*\[Eta]/\
4,\n       -(1-\[Xi]^2)*(1-\[Eta])*\[Eta]/2, (1+\[Xi])*\[Xi]*(1-\[Eta]^2)/2,\n\
        (1-\[Xi]^2)*(1+\[Eta])*\[Eta]/2,-(1-\[Xi])*\[Xi]*(1-\[Eta]^2)/2,(1-\
\[Xi]^2)*(1-\[Eta]^2)};\n   dN\[Xi]={(1-2*\[Xi])*\[Eta]*(1-\[Eta])/4,-(1+2*\
\[Xi])*\[Eta]*(1-\[Eta])/4,\n        (1+2*\[Xi])*\[Eta]*(1+\[Eta])/4,-(1-2*\
\[Xi])*\[Eta]*(1+\[Eta])/4,\n         \[Xi]*\[Eta]*(1-\[Eta]), \
(1/2+\[Xi])*(1-\[Eta]^2),\n        -\[Xi]*\[Eta]*(1+\[Eta]),-(1/2-\[Xi])*(1-\
\[Eta]^2), 2*\[Xi]*(\[Eta]^2-1)};     \n   dN\[Eta]={\[Xi]*(1-\[Xi])*(1-2*\
\[Eta])/4, -\[Xi]*(1+\[Xi])*(1-2*\[Eta])/4,\n        \[Xi]*(1+\[Xi])*(1+2*\
\[Eta])/4, -\[Xi]*(1-\[Xi])*(1+2*\[Eta])/4,\n        \
-(1/2-\[Eta])*(1-\[Xi]^2),  -\[Xi]*(1+\[Xi])*\[Eta],         \n         (1/2+\
\[Eta])*(1-\[Xi]^2),   \[Xi]*(1-\[Xi])*\[Eta], 2*(\[Xi]^2-1)*\[Eta]};\n",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["   Print[Simplify[dN\[Xi]-D[Nf,\[Xi]]]]; \
Print[Simplify[dN\[Eta]-D[Nf,\[Eta]]]];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   x=Table[ncoor[[i,1]],{i,9}]; y=Table[ncoor[[i,2]],{i,9}];\n  \
 J11=dN\[Xi].x; J21=dN\[Xi].y; J12=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx=( J22*dN\[Xi]-J21*dN\[Eta])/Jdet; \
dNy=(-J12*dN\[Xi]+J11*dN\[Eta])/Jdet;  \n   Return[{Nf,dNx,dNy,Jdet}]\n];\n\n\
",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["ncoor={{0,0},{2*a,0},{a,b},{0,b},{a,0},{a,b/2},{a/2,b},{0,b/2},{3*\
a/4,b/2}}; \n{Nf,dNx,dNy,Jdet}=Simplify[Quad9IsoPShapeFunDer[ncoor,{\[Xi],\
\[Eta]}]];\nPrint[{Nf,dNx,dNy,Jdet}];\nx=Table[ncoor[[i,1]],{i,9}]; \
y=Table[ncoor[[i,2]],{i,9}];\n\
Print[Simplify[{dNx.(x*exx),dNy.(y*eyy),dNx.(y*gxy)+dNx.(x*gxy)}]];\n\n",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["Quad8IsoPShapeFunDer[ncoor_,qcoor_]:= Module[\n  {Nf,dNx,dNy,dN\
\[Xi],dN\[Eta],J11,J12,J21,J22,Jdet,\[Xi],\[Eta],x,y},\n  \
{\[Xi],\[Eta]}=qcoor;  \n   Nf= \
{-(1-\[Xi])*(1-\[Eta])*(1+\[Xi]+\[Eta])/4,-(1+\[Xi])*(1-\[Eta])*(1-\[Xi]+\
\[Eta])/4,\n         -(1+\[Xi])*(1+\[Eta])*(1-\[Xi]-\[Eta])/4,-(1-\[Xi])*(1+\
\[Eta])*(1+\[Xi]-\[Eta])/4,\n          (1+\[Xi])*(1-\[Xi])*(1-\[Eta])/2, (1+\
\[Xi])*(1+\[Eta])*(1-\[Eta])/2,\n          (1+\[Xi])*(1-\[Xi])*(1+\[Eta])/2, \
(1-\[Xi])*(1+\[Eta])*(1-\[Eta])/2};\n   \
dN\[Xi]={(1-\[Eta])*(2*\[Xi]+\[Eta])/4,   (1-\[Eta])*(2*\[Xi]-\[Eta])/4, \n   \
     (1+\[Eta])*(2*\[Xi]+\[Eta])/4,   (1+\[Eta])*(2*\[Xi]-\[Eta])/4, \n       \
  \[Xi]*(\[Eta]-1),  (1+\[Eta])*(1-\[Eta])/2, -\[Xi]*(1+\[Eta]), \
-(1+\[Eta])*(1-\[Eta])/2};\n   dN\[Eta]={(1-\[Xi])*(\[Xi]+2*\[Eta])/4,   -(1+\
\[Xi])*(\[Xi]-2*\[Eta])/4, \n        (1+\[Xi])*(\[Xi]+2*\[Eta])/4,   \
-(1-\[Xi])*(\[Xi]-2*\[Eta])/4, \n       -(1+\[Xi])*(1-\[Xi])/2, -(1+\[Xi])*\
\[Eta], (1+\[Xi])*(1-\[Xi])/2,-(1-\[Xi])*\[Eta]};",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["\n",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["   Print[Simplify[dN\[Xi]-D[Nf,\[Xi]]]]; \
Print[Simplify[dN\[Eta]-D[Nf,\[Eta]]]];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   x=Table[ncoor[[i,1]],{i,8}]; y=Table[ncoor[[i,2]],{i,8}];\n  \
 J11=dN\[Xi].x; J21=dN\[Xi].y; J12=dN\[Eta].x; J22=dN\[Eta].y;\n   \
Jdet=Simplify[J11*J22-J12*J21];\n   dNx=( J22*dN\[Xi]-J21*dN\[Eta])/Jdet; \
dNy=(-J12*dN\[Xi]+J11*dN\[Eta])/Jdet;  \n   Return[{Nf,dNx,dNy,Jdet}]\n];\n\n\
",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["ncoor={{0,0},{2*a,0},{a,b},{0,b},{a,0},{a,b/2},{a/2,b},{0,b/2}}; \
\n{Nf,dNx,dNy,Jdet}=Simplify[Quad8IsoPShapeFunDer[ncoor,{\[Xi],\[Eta]}]];\n\
Print[{Nf,dNx,dNy,Jdet}];\nx=Table[ncoor[[i,1]],{i,8}]; \
y=Table[ncoor[[i,2]],{i,8}];\n\
Print[Simplify[{dNx.(x*exx),dNy.(y*eyy),dNx.(y*gxy)+dNx.(x*gxy)}]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2C:  Modules for 3-node linear triangle in plane stress.  \
Self-contained.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Trig3IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:=\n  \
Module[{x1,x2,x3,y1,y2,y3,x21,x13,x32,y12,y31,y23,\n  \
numer=False,h=th,A,Be,Ke}, {{x1,y1},{x2,y2},{x3,y3}}=ncoor;\n  If \
[Length[th]==3, h=(th[[1]]+th[[2]]+th[[3]])/3];\n  If [Length[options]>0, \
numer=options[[1]]]; \n  {x21,x13,x32,y12,y31,y23}={x2-x1,x1-x3,x3-x2,\n      \
 y1-y2,y3-y1,y2-y3};\n  A=Simplify[x21*y31-x13*y12]/2; If \
[!numer,A=Simplify[A]];\n  Be={{y23,0,y31,0,y12,0},{0,x32,0,x13,0,x21},\n     \
 {x32,y23,x13,y31,x21,y12}}/(2*A);\n  Ke=A*h*Transpose[Be].Emat.Be; If \
[!numer, Ke=Simplify[Ke]];\n  Return[Ke]];\n  \n\
Trig3IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:=\n  \
Module[{x1,x2,x3,y1,y2,y3,x21,x13,x32,y12,y31,y23,\n  \
numer=False,h=th,A,Be,ue=udis,sige=Table[{0,0,0},{3}]}, \n  \
{{x1,y1},{x2,y2},{x3,y3}}=ncoor;\n  If [Length[th]==3, \
h=(th[[1]]+th[[2]]+th[[3]])/3];\n  If [Length[options]>0, \
numer=options[[1]]];\n  If [Length[udis]==3, ue=Flatten[udis]];  \n  \
{x21,x13,x32,y12,y31,y23}={x2-x1,x1-x3,x3-x2,\n       y1-y2,y3-y1,y2-y3};\n  \
A=Simplify[x21*y31-x13*y12]/2; If [!numer,A=Simplify[A]];\n  \
Be={{y23,0,y31,0,y12,0},{0,x32,0,x13,0,x21},\n      \
{x32,y23,x13,y31,x21,y12}}/(2*A);\n  \
sige[[1]]=sige[[2]]=sige[[3]]=Emat.(Be.ue); \n  If [!numer, \
sige=Simplify[sige]];\n  Return[sige]];\n\n",
  StyleBox["ClearAll[Em,\[Nu],a,b,e,h,p,numer,exx,eyy,gxy];  h=1;\n\
ncoor={{0,0},{3,1},{2,2}}; Em=32/3; \[Nu]=1/3; \n\
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\n\
Print[\"Emat=\",Emat//MatrixForm];  options={False};\n\
Ke=Trig3IsoPMembraneStiffness[ncoor,Emat,h,options];  \nPrint[\"Ke = \
\",Chop[Ke]//MatrixForm];\nPrint[\"eigs of Ke = \",Chop[Eigenvalues[N[Ke]]]];\
\nxx=Table[{ncoor[[i,1]],0},{i,3}]; yy=Table[{0,ncoor[[i,2]]},{i,3}];\n\
xy=Table[{0,ncoor[[i,1]]},{i,3}]; yx=Table[{ncoor[[i,2]],0},{i,3}];\n\
ue=Flatten[xx*exx+yy*eyy+(xy+yx)*gxy/2]; Emat=IdentityMatrix[3];\n\
sige=Trig3IsoPMembraneStresses[ncoor,Emat,h,ue,{False}];\nPrint[\"strain \
check = \",Chop[sige]//MatrixForm];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2E:  Modules for 6-node quadratic triangle in plane stress,
Uses TrigGaussRuleInfo,which must have been initialized.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Trig6IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= \n \
Module[{i,k,p=3,numer=False,h=th,tcoor,w,c,\n  \
Nf,dNx,dNy,Jdet,Be,Ke=Table[0,{12},{12}]},\n  If [Length[options]==1, {numer} \
 =options];\n  If [Length[options]==2, {numer,p}=options];\n  If \
[!MemberQ[{1,-3,3,6,7},p], Print[\"Illegal p\"]; Return[Ke]];\n  For [k=1, \
k<=Abs[p], k++,  \n      {tcoor,w}= TrigGaussRuleInfo[{p,numer},k];\n      \
{Nf,dNx,dNy,Jdet}= Trig6IsoPShapeFunDer[ncoor,tcoor];\n      If [numer, \
{Nf,dNx,dNy,Jdet}=N[{Nf,dNx,dNy,Jdet}]];\n      If [Length[th]==6, h=th.Nf]; \
c=w*Jdet*h/2;\n      Be={Flatten[Table[{dNx[[i]],0       },{i,6}]],\n         \
 Flatten[Table[{0,       dNy[[i]]},{i,6}]],\n          \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,6}]]};\n      \
Ke+=c*Transpose[Be].(Emat.Be); \n      ]; If[!numer,Ke=Simplify[Ke]]; \
Return[Ke]\n ];\n  \n\
Trig6IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= \n  \
Module[{i,k,p=3,numer=False,h=th,qcoor,Nf,dNx,dNy,Jdet,\n  \
qctab,ue=udis,Be,sige=Table[{0,0,0},{6}]},  \n  \
qctab={{1,0,0},{0,1,0},{0,0,1},{1/2,1/2,0},{0,1/2,1/2},{1/2,0,1/2}};\n  If \
[Length[options]>0, numer=options[[1]]];\n  If [Length[udis]==6, \
ue=Flatten[udis]];\n  For [k=1, k<=Length[sige], k++, \n       \
qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; \n      \
{Nf,dNx,dNy,Jdet}=Trig6IsoPShapeFunDer[ncoor,qcoor];\n       \
Be={Flatten[Table[{dNx[[i]], 0      },{i,6}]],\n           Flatten[Table[{0,  \
     dNy[[i]]},{i,6}]],\n           \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,6}]]}; \n       sige[[k]]=Emat.(Be.ue); \
\n      ]; If [!numer, sige=Simplify[sige]]; Return[sige]\n ];\n \n\
Trig6IsoPShapeFunDer[ncoor_,tcoor_]:= Module[\n  \
{\[Zeta]1,\[Zeta]2,\[Zeta]3,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,\n  \
dx4,dx5,dx6,dy4,dy5,dy6,Jx21,Jx32,Jx13,Jy12,Jy23,Jy31,\n  Nf,dNx,dNy,Jdet}, {\
\[Zeta]1,\[Zeta]2,\[Zeta]3}=tcoor;  \n  \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6}}=ncoor;\n  dx4=x4-(x1+x2)/2; \
dx5=x5-(x2+x3)/2; dx6=x6-(x3+x1)/2; \n  dy4=y4-(y1+y2)/2; dy5=y5-(y2+y3)/2; \
dy6=y6-(y3+y1)/2;\n  Nf={\[Zeta]1*(2*\[Zeta]1-1),\[Zeta]2*(2*\[Zeta]2-1),\
\[Zeta]3*(2*\[Zeta]3-1),4*\[Zeta]1*\[Zeta]2,4*\[Zeta]2*\[Zeta]3,4*\[Zeta]3*\
\[Zeta]1};\n  Jx21= x2-x1+4*(dx4*(\[Zeta]1-\[Zeta]2)+(dx5-dx6)*\[Zeta]3);\n  \
Jx32= x3-x2+4*(dx5*(\[Zeta]2-\[Zeta]3)+(dx6-dx4)*\[Zeta]1);\n  Jx13= \
x1-x3+4*(dx6*(\[Zeta]3-\[Zeta]1)+(dx4-dx5)*\[Zeta]2);\n  Jy12= y1-y2+4*(dy4*(\
\[Zeta]2-\[Zeta]1)+(dy6-dy5)*\[Zeta]3);\n  Jy23= y2-y3+4*(dy5*(\[Zeta]3-\
\[Zeta]2)+(dy4-dy6)*\[Zeta]1);\n  Jy31= \
y3-y1+4*(dy6*(\[Zeta]1-\[Zeta]3)+(dy5-dy4)*\[Zeta]2);\n  Jdet = \
Jx21*Jy31-Jy12*Jx13;\n  dNx= {(4*\[Zeta]1-1)*Jy23,(4*\[Zeta]2-1)*Jy31,(4*\
\[Zeta]3-1)*Jy12,4*(\[Zeta]2*Jy23+\[Zeta]1*Jy31),\n         4*(\[Zeta]3*Jy31+\
\[Zeta]2*Jy12),4*(\[Zeta]1*Jy12+\[Zeta]3*Jy23)}/Jdet;\n  dNy= \
{(4*\[Zeta]1-1)*Jx32,(4*\[Zeta]2-1)*Jx13,(4*\[Zeta]3-1)*Jx21,4*(\[Zeta]2*Jx32+\
\[Zeta]1*Jx13),\n         4*(\[Zeta]3*Jx13+\[Zeta]2*Jx21),4*(\[Zeta]1*Jx21+\
\[Zeta]3*Jx32)}/Jdet;\n  Return[{Nf,dNx,dNy,Jdet}]\n];\n",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["ClearAll[Em,nu,h,exx,eyy,gxy]; h=1; Em=288; nu=1/3;\n\
ncoor={{0,0},{6,2},{4,4},{3,1},{5,3},{2,2}};\n\
Emat=Em/(1-nu^2)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}};\n\
Print[\"Emat=\",Emat//MatrixForm]; options={False,7};\n\
Ke=Trig6IsoPMembraneStiffness[ncoor,Emat,h,options]; \n\
Print[Chop[Ke]//MatrixForm];\nPrint[\"eigs of \
Ke=\",Chop[Eigenvalues[N[Ke]]]];\nxx=Table[{ncoor[[i,1]],0},{i,6}]; \
yy=Table[{0,ncoor[[i,2]]},{i,6}];\nxy=Table[{0,ncoor[[i,1]]},{i,6}]; \
yx=Table[{ncoor[[i,2]],0},{i,6}];\nue=Flatten[xx*exx+yy*eyy+(xy+yx)*gxy/2]; \
Emat=IdentityMatrix[3];\nsige",
    FontColor->RGBColor[0, 0, 1]],
  "=",
  StyleBox["Trig6IsoPMembraneStresses[ncoor,Emat,h,ue,{False}];\n",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["Print[\"strain chk=\",Simplify[Transpose[sige]]//MatrixForm];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2F:  Modules for 10-node cubic triangle in plane stress,
Uses TrigGaussRuleInfo,which must have been initialized. (cell not \
initialized)\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "Trig10IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= \n \
Module[{i,k,l,p=7,numer=False,h=th,tcoor,w,c,\n  \
Nf,dNx,dNy,Jdet,Be,Ke=Table[0,{20},{20}]}, \n  If [Length[options]>=1, \
numer=options[[1]]];\n  If [Length[options]>=2, p=    options[[2]]];\n  If \
[!MemberQ[{1,3,-3,6,7},p], Print[\"Illegal p\"]; Return[Null]];\n  For [k=1, \
k<=Abs[p], k++,  \n      {tcoor,w}= TrigGaussRuleInfo[{p,numer},k];\n      \
{Nf,dNx,dNy,Jdet}= Trig10IsoPShapeFunDer[ncoor,tcoor];\n      If [numer, \
{Nf,dNx,dNy,Jdet}=N[{Nf,dNx,dNy,Jdet}]];\n      If [Length[th]==10, h=th.Nf]; \
c=w*Jdet*h/2;\n      Be= {Flatten[Table[{dNx[[i]],0       },{i,10}]],\n       \
    Flatten[Table[{0,       dNy[[i]]},{i,10}]],\n           \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,10}]]};\n      \
Ke+=c*Transpose[Be].(Emat.Be); \n      ]; If[!numer,Ke=Simplify[Ke]]; \
Return[Ke]\n ];\n  \n\
Trig10IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:=\n  \
Module[{i,k,p=7,numer=False,h=th,qcoor,Nf,\n  \
dNx,dNy,Jdet,Be,qctab,ue=udis,sige=Table[{0,0,0},{10}]},  \n  \
qctab={{1,0,0},{0,1,0},{0,0,1},{2/3,1/3,0},{1/3,2/3,0},\n         \
{0,2/3,1/3},{0,1/3,2/3},{1/3,0,2/3},{2/3,0,1/3},{1/3,1/3,1/3}}; \n  If \
[Length[options]>0, numer=options[[1]]];\n  If [Length[udis]==10,  \
ue=Flatten[udis]];\n  For [k=1, k<=Length[sige], k++, \n      \
qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; \n     \
{Nf,dNx,dNy,Jdet}=Trig10IsoPShapeFunDer[ncoor,qcoor];\n      Be= \
{Flatten[Table[{dNx[[i]],0       },{i,10}]],\n           Flatten[Table[{0,    \
   dNy[[i]]},{i,10}]],\n           \
Flatten[Table[{dNy[[i]],dNx[[i]]},{i,10}]]};\n      sige[[k]]=Emat.(Be.ue); \n\
      ]; If [!numer, sige=Simplify[sige]]; Return[sige]\n ];\n \n\
Trig10IsoPShapeFunDer[ncoor_,tcoor_]:= Module[\n  \
{\[Zeta]1,\[Zeta]2,\[Zeta]3,x1,x2,x3,x4,x5,x6,x7,x8,x9,x0,y1,y2,y3,y4,y5,y6,\
y7,y8,y9,y0,\n   dx4,dx5,dx6,dx7,dx8,dx9,dx0,dy4,dy5,dy6,dy7,dy8,dy9,dy0,\n  \
Jx21,Jx32,Jx13,Jy12,Jy23,Jy31,Nf,dNx,dNy,Jdet}, \n \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6},{x7,y7},\n  \
{x8,y8},{x9,y9},{x0,y0}}=ncoor;  {\[Zeta]1,\[Zeta]2,\[Zeta]3}=tcoor;     \n  \
dx4=x4-(2*x1+x2)/3; dx5=x5-(x1+2*x2)/3; dx6=x6-(2*x2+x3)/3; \n  \
dx7=x7-(x2+2*x3)/3; dx8=x8-(2*x3+x1)/3; dx9=x9-(x3+2*x1)/3; \n  \
dy4=y4-(2*y1+y2)/3; dy5=y5-(y1+2*y2)/3; dy6=y6-(2*y2+y3)/3; \n  \
dy7=y7-(y2+2*y3)/3; dy8=y8-(2*y3+y1)/3; dy9=y9-(y3+2*y1)/3;\n  \
dx0=x0-(x1+x2+x3)/3; dy0=y0-(y1+y2+y3)/3; \n  Nf={\[Zeta]1*(3*\[Zeta]1-1)*(3*\
\[Zeta]1-2),\[Zeta]2*(3*\[Zeta]2-1)*(3*\[Zeta]2-2),\[Zeta]3*(3*\[Zeta]3-1)*(3*\
\[Zeta]3-2), \n      \
9*\[Zeta]1*\[Zeta]2*(3*\[Zeta]1-1),9*\[Zeta]1*\[Zeta]2*(3*\[Zeta]2-1),9*\
\[Zeta]2*\[Zeta]3*(3*\[Zeta]2-1),\n      \
9*\[Zeta]2*\[Zeta]3*(3*\[Zeta]3-1),9*\[Zeta]3*\[Zeta]1*(3*\[Zeta]3-1),9*\
\[Zeta]3*\[Zeta]1*(3*\[Zeta]1-1),54*\[Zeta]1*\[Zeta]2*\[Zeta]3}/2; \n  \
Jx21=x2-x1+(9/2)*(dx4*(\[Zeta]1*(3*\[Zeta]1-6*\[Zeta]2-1)+\[Zeta]2)+\n    \
dx5*(\[Zeta]2*(1-3*\[Zeta]2+6*\[Zeta]1)-\[Zeta]1)+dx6*\[Zeta]3*(6*\[Zeta]2-1)+\
dx7*\[Zeta]3*(3*\[Zeta]3-1)+\n    \
dx8*\[Zeta]3*(1-3*\[Zeta]3)+dx9*\[Zeta]3*(1-6*\[Zeta]1)+6*dx0*\[Zeta]3*(\
\[Zeta]1-\[Zeta]2)); \n  Jx32=x3-x2+(9/2)*(dx4*\[Zeta]1*(1-3*\[Zeta]1)+dx5*\
\[Zeta]1*(1-6*\[Zeta]2)+\n    \
dx6*(\[Zeta]2*(3*\[Zeta]2-6*\[Zeta]3-1)+\[Zeta]3)+dx7*(\[Zeta]3*(1-3*\[Zeta]3+\
6*\[Zeta]2)-\[Zeta]2)+    \n    dx8*\[Zeta]1*(6*\[Zeta]3-1)+dx9*\[Zeta]1*(3*\
\[Zeta]1-1)+6*dx0*\[Zeta]1*(\[Zeta]2-\[Zeta]3)) ;\n  \
Jx13=x1-x3+(9/2)*(dx4*(6*\[Zeta]1-1)*\[Zeta]2+dx5*\[Zeta]2*(3*\[Zeta]2-1)+\n  \
  dx6*\[Zeta]2*(1-3*\[Zeta]2)+dx7*\[Zeta]2*(1-6*\[Zeta]3)+dx8*(\[Zeta]3*(3*\
\[Zeta]3-6*\[Zeta]1-1)+\[Zeta]1)+\n    \
dx9*(\[Zeta]1*(1-3*\[Zeta]1+6*\[Zeta]3)-\[Zeta]3)+6*dx0*\[Zeta]2*(\[Zeta]3-\
\[Zeta]1));\n  Jy12=y1-y2-(9/2)*(dy4*(\[Zeta]1*(3*\[Zeta]1-6*\[Zeta]2-1)+\
\[Zeta]2)+\n    \
dy5*(\[Zeta]2*(1-3*\[Zeta]2+6*\[Zeta]1)-\[Zeta]1)+dy6*\[Zeta]3*(6*\[Zeta]2-1)+\
dy7*\[Zeta]3*(3*\[Zeta]3-1)+\n    \
dy8*\[Zeta]3*(1-3*\[Zeta]3)+dy9*\[Zeta]3*(1-6*\[Zeta]1)+6*dy0*\[Zeta]3*(\
\[Zeta]1-\[Zeta]2)); \n  Jy23=y2-y3-(9/2)*(dy4*\[Zeta]1*(1-3*\[Zeta]1)+dy5*\
\[Zeta]1*(1-6*\[Zeta]2)+\n    \
dy6*(\[Zeta]2*(3*\[Zeta]2-6*\[Zeta]3-1)+\[Zeta]3)+dy7*(\[Zeta]3*(1-3*\[Zeta]3+\
6*\[Zeta]2)-\[Zeta]2)+    \n    dy8*\[Zeta]1*(6*\[Zeta]3-1)+dy9*\[Zeta]1*(3*\
\[Zeta]1-1)+6*dy0*\[Zeta]1*(\[Zeta]2-\[Zeta]3)) ;\n  \
Jy31=y3-y1-(9/2)*(dy4*(6*\[Zeta]1-1)*\[Zeta]2+dy5*\[Zeta]2*(3*\[Zeta]2-1)+\n  \
  dy6*\[Zeta]2*(1-3*\[Zeta]2)+dy7*\[Zeta]2*(1-6*\[Zeta]3)+dy8*(\[Zeta]3*(3*\
\[Zeta]3-6*\[Zeta]1-1)+\[Zeta]1)+\n    \
dy9*(\[Zeta]1*(1-3*\[Zeta]1+6*\[Zeta]3)-\[Zeta]3)+6*dy0*\[Zeta]2*(\[Zeta]3-\
\[Zeta]1));\n  Jdet = Jx21*Jy31-Jy12*Jx13;\n  dNx={Jy23*(2/9-2*\[Zeta]1+3*\
\[Zeta]1^2),Jy31*(2/9-2*\[Zeta]2+3*\[Zeta]2^2),Jy12*(2/9-2*\[Zeta]3+3*\[Zeta]\
3^2),\n       Jy31*\[Zeta]1*(3*\[Zeta]1-1)+Jy23*\[Zeta]2*(6*\[Zeta]1-1),Jy23*\
\[Zeta]2*(3*\[Zeta]2-1)+Jy31*\[Zeta]1*(6*\[Zeta]2-1), \n       \
Jy12*\[Zeta]2*(3*\[Zeta]2-1)+Jy31*\[Zeta]3*(6*\[Zeta]2-1),Jy31*\[Zeta]3*(3*\
\[Zeta]3-1)+Jy12*\[Zeta]2*(6*\[Zeta]3-1), \n       \
Jy23*\[Zeta]3*(3*\[Zeta]3-1)+Jy12*\[Zeta]1*(6*\[Zeta]3-1),Jy12*\[Zeta]1*(3*\
\[Zeta]1-1)+Jy23*\[Zeta]3*(6*\[Zeta]1-1), \n       \
6*(Jy12*\[Zeta]1*\[Zeta]2+Jy31*\[Zeta]1*\[Zeta]3+Jy23*\[Zeta]2*\[Zeta]3)}/(2*\
Jdet/9);\n  dNy={Jx32*(2/9-2*\[Zeta]1+3*\[Zeta]1^2),Jx13*(2/9-2*\[Zeta]2+3*\
\[Zeta]2^2),Jx21*(2/9-2*\[Zeta]3+3*\[Zeta]3^2), \n       Jx13*\[Zeta]1*(3*\
\[Zeta]1-1)+Jx32*\[Zeta]2*(6*\[Zeta]1-1),Jx32*\[Zeta]2*(3*\[Zeta]2-1)+Jx13*\
\[Zeta]1*(6*\[Zeta]2-1), \n       \
Jx21*\[Zeta]2*(3*\[Zeta]2-1)+Jx13*\[Zeta]3*(6*\[Zeta]2-1),Jx13*\[Zeta]3*(3*\
\[Zeta]3-1)+Jx21*\[Zeta]2*(6*\[Zeta]3-1), \n       \
Jx32*\[Zeta]3*(3*\[Zeta]3-1)+Jx21*\[Zeta]1*(6*\[Zeta]3-1),Jx21*\[Zeta]1*(3*\
\[Zeta]1-1)+Jx32*\[Zeta]3*(6*\[Zeta]1-1), \n       \
6*(Jx21*\[Zeta]1*\[Zeta]2+Jx13*\[Zeta]1*\[Zeta]3+Jx32*\[Zeta]2*\[Zeta]3)}/(2*\
Jdet/9);\n  Return[Simplify[{Nf,dNx,dNy,Jdet}]]\n];\n\n",
  StyleBox["ClearAll[Em,\[Nu],a,b,e,h,exx,eyy,gxy]; \n\
{{x1,y1},{x2,y2},{x3,y3}}={{0,0},{6,2},{4,4}};\nx4=(2*x1+x2)/3; \
x5=(x1+2*x2)/3; y4=(2*y1+y2)/3; y5=(y1+2*y2)/3;\nx6=(2*x2+x3)/3; \
x7=(x2+2*x3)/3; y6=(2*y2+y3)/3; y7=(y2+2*y3)/3;\nx8=(2*x3+x1)/3; \
x9=(x3+2*x1)/3; y8=(2*y3+y1)/3; y9=(y3+2*y1)/3;\nx0=(x1+x2+x3)/3; \
y0=(y1+y2+y3)/3;\nncoor= \
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6},{x7,y7},\n        \
{x8,y8},{x9,y9},{x0,y0}};\nEm=1920; \[Nu]=0; h=1;\nEmat=Em/(1-\[Nu]^2)*{{1,\
\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\nPrint[\"Emat=\",Emat//MatrixForm];\n\
Ke=Trig10IsoPMembraneStiffness[ncoor,Emat,h,{False,7}];  \nKe=Simplify[Ke];  \
Print[Ke//MatrixForm]; \nev=Chop[Eigenvalues[N[Ke]]]; Print[\"eigs of \
Ke=\",ev];\nxx=Table[{ncoor[[i,1]],0},{i,10}]; \
yy=Table[{0,ncoor[[i,2]]},{i,10}];\nxy=Table[{0,ncoor[[i,1]]},{i,10}]; \
yx=Table[{ncoor[[i,2]],0},{i,10}];\nue=Flatten[xx*exx+yy*eyy+(xy+yx)*gxy/2]; \
Emat=IdentityMatrix[3];\n\
sige=Trig10IsoPMembraneStresses[ncoor,Emat,h,ue,{False}];\nPrint[\"strain \
chk=\",Simplify[Transpose[sige]]//MatrixForm]",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[";",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2G - Modules for 2-node plane bar\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PlaneBar2Stiffness[ncoor_,Em_,A_,options_]:= Module[\n \
{x1,x2,y1,y2,x21,y21,EA,numer,L,LL,LLL,Ke}, \n  {{x1,y1},{x2,y2}}=ncoor; \
{x21,y21}={x2-x1,y2-y1};\n  EA=Em*A; {numer}=options; LL=x21^2+y21^2; \
L=Sqrt[LL];\n  If [numer,{x21,y21,EA,LL,L}=N[{x21,y21,EA,LL,L}]]; \n  If \
[!numer, L=PowerExpand[L]]; LLL=Simplify[LL*L];\n  Ke=(Em*A/LLL)*{{ x21*x21, \
x21*y21,-x21*x21,-x21*y21},\n                 { y21*x21, \
y21*y21,-y21*x21,-y21*y21},\n                 {-x21*x21,-x21*y21, x21*x21, \
x21*y21},\n                 {-y21*x21,-y21*y21, y21*x21, y21*y21}}; \n  \
Return[Ke]\n];  \n\nPlaneBar2IntForce[ncoor_,Em_,A_,udis_,options_]:= Module[ \
\n {x1,x2,y1,y2,z1,z2,x21,y21,z21,EA,numer,LL,ue=udis,pe},",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n {{x1,y1},{x2,y2}}=ncoor; {x21,y21}={x2-x1,y2-y1};\n  EA=Em*A; \
{numer}=options;  LL=x21^2+y21^2;\n  If [Length[udis]==2, ue=Flatten[udis]];\n\
  If [numer,{x21,y21,EA,LL}=N[{x21,y21,EA,LL}]];\n  \
pe=(EA/LL)*(x21*(ue[[3]]-ue[[1]])+y21*(ue[[4]]-ue[[2]]));\n  Return[pe]]; \n  \
\n",
  StyleBox["ClearAll[A,Em];\nncoor={{0,0},{30,40}}; Em=1000; A=5; \nKe= \
PlaneBar2Stiffness[ncoor,Em,A,{True}];\nPrint[\"Ke=\",Ke//MatrixForm];\n\
Print[\"Eigenvalues of Ke=\",Chop[Eigenvalues[N[Ke]]]];\npe= \
PlaneBar2IntForce[ncoor,Em,A,{{0,0},{3,4}/10},{True}];\nPrint[\"pe=\",pe];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 2H - Modules for 3-node plane bar (possibly curved in plane)  \
(cell not initialized)\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PlaneBar3Stiffness[ncoor_,Em_,A_,options_]:= \n \
Module[{x1,y1,x2,y2,x3,y3,x21,y21,xm,ym,\[Xi],\[Alpha],\[Beta],LLe,Le,\n      \
   numer,Be,J,JJ,Be1,Be2,J1,J2,Kebar,T,Ke},\n  \
{{x1,y1},{x2,y2},{x3,y3}}=ncoor; numer=options[[1]]; \n   \
{x21,y21}={x2-x1,y2-y1}; {xm,ym}={x1+x2,y1+y2}/2; \n   LLe=x21^2+y21^2; \
Le=Sqrt[LLe];\n   \[Alpha]=( (x3-xm)*x21+(y3-ym)*y21)/LLe;\n   \
\[Beta]=(-(x3-xm)*y21+(y3-ym)*x21)/LLe;\n   If [!numer, Le=PowerExpand[Le]; {\
\[Alpha],\[Beta]}=Simplify[{\[Alpha],\[Beta]}]]; \n   \
JJ=LLe*((1/2-2*\[Alpha]*\[Xi])^2+4*\[Beta]^2*\[Xi]^2);\n   J= Le*Sqrt[(1/2-2*\
\[Alpha]*\[Xi])^2+4*\[Beta]^2*\[Xi]^2];\n   \
Be={{(1/2-2*\[Alpha]*\[Xi])*(\[Xi]-1/2),-2*\[Beta]*\[Xi]*(\[Xi]-1/2),(1/2-2*\
\[Alpha]*\[Xi])*(\[Xi]+1/2),\n      -2*\[Beta]*\[Xi]*(\[Xi]+1/2),-(1/2-2*\
\[Alpha]*\[Xi])*(2*\[Xi]),4*\[Beta]*\[Xi]^2}}*Le/JJ;\n   \
Be1=Be/.\[Xi]->-Sqrt[3]/3;  Be2=Be/.\[Xi]->Sqrt[3]/3; \n   J1 =J \
/.\[Xi]->-Sqrt[3]/3;  J2= J /.\[Xi]->Sqrt[3]/3; \n   If[numer, \
{B1,J1,B2,J2}=N[{B1,J1,B2,J2}]];\n   \
Kebar=Em*A*(J1*Transpose[Be1].Be1+J2*Transpose[Be2].Be2);    \n   T={{x21, \
y21,0,0,0,0},{-y21,x21,0,0,0,0},{0,0, x21,y21,0,0},\n     \
{0,0,-y21,x21,0,0},{0,0,0,0,x21,y21}, {0,0,0,0,-y21,x21}}; \n   \
Ke=(Transpose[T].Kebar.T)/LLe; (* avoids taking Sqrt[LLe] *)\n   If [!numer, \
Ke=Simplify[Ke]]; Return[Ke] ];\n   \n\
PlaneBar3IntForce[ncoor_,Em_,A_,udis_,options_]:= \n \
Module[{x1,y1,x2,y2,x3,y3,x21,y21,xm,ym,\[Xi],\[Alpha],\[Beta],LLe,Le,\n      \
   numer,Be,J,JJ,Be1,Be2,Be3,ue=udis,ubar,T,pe},\n  \
{{x1,y1},{x2,y2},{x3,y3}}=ncoor; numer=options[[1]]; \n   \
{x21,y21}={x2-x1,y2-y1}; {xm,ym}={x1+x2,y1+y2}/2;\n   LLe=x21^2+y21^2; \
Le=Sqrt[LLe];\n   \[Alpha]=( (x3-xm)*x21+(y3-ym)*y21)/LLe;\n   \
\[Beta]=(-(x3-xm)*y21+(y3-ym)*x21)/LLe;\n   If [!numer, Le=PowerExpand[Le]; {\
\[Alpha],\[Beta]}=Simplify[{\[Alpha],\[Beta]}]];\n   T={{x21, \
y21,0,0,0,0},{-y21,x21,0,0,0,0},{0,0, x21,y21,0,0},\n      \
{0,0,-y21,x21,0,0},{0,0,0,0,x21,y21}, {0,0,0,0,-y21,x21}}/Le;  \n   \
JJ=LLe*((1/2-2*\[Alpha]*\[Xi])^2+4*\[Beta]^2*\[Xi]^2);\n   J= Le*Sqrt[(1/2-2*\
\[Alpha]*\[Xi])^2+4*\[Beta]^2*\[Xi]^2];\n   \
Be={{(1/2-2*\[Alpha]*\[Xi])*(\[Xi]-1/2),-2*\[Beta]*\[Xi]*(\[Xi]-1/2),(1/2-2*\
\[Alpha]*\[Xi])*(\[Xi]+1/2),\n      -2*\[Beta]*\[Xi]*(\[Xi]+1/2),-(1/2-2*\
\[Alpha]*\[Xi])*(2*\[Xi]),4*\[Beta]*\[Xi]^2}}*Le/JJ;\n   Be1=Be/.\[Xi]->-1; \
Be2=Be/.\[Xi]->1; Be3=Be/.\[Xi]->0; \n   If [numer, \
{Be1,Be2,Be3}=N[{Be1,Be2,Be3}]];\n   If [Length[udis]==3,ue=Flatten[udis]]; \
ubar=T.ue;\n   pe=Em*A*{Be1.ubar,Be2.ubar,Be3.ubar};\n   If [!numer, \
pe=Simplify[pe]]; \n   Return[pe] ];\n  \n",
  StyleBox["ClearAll[A,Em,\[Epsilon]];\n\
ncoor={{0,0},{30,40},{15,20+\[Epsilon]}}; Em=3000; A=5; \n\
ncoor={{0,0},{50,0},{25,\[Epsilon]}}; Em=3000; A=5; \nKe= \
PlaneBar3Stiffness[ncoor,Em,A,{False}];\nPrint[\"Ke=\",Ke//MatrixForm];\n\
Print[Simplify[Normal[Series[Ke,{\[Epsilon],0,2}]]]//MatrixForm];\n\
(*Print[\"Eigenvalues of Ke=\",Chop[Eigenvalues[N[Ke]]]]",
    FontColor->RGBColor[0, 0, 1]],
  ";*)"
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 3:  Plane stress assembler  of MET type. Stiffness modules \
must have been initialized.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PlaneStressMasterStiffness[nodxyz_,eletyp_,elenod_,elemat_,\n  \
elefab_,prcopt_]:=Module[{numele=Length[elenod],\n  \
numnod=Length[nodxyz],ncoor,type,e,enl,neldof,OKtyp,OKenl,\n  \
i,j,n,ii,jj,eft,Emat,th,numer,Ke,K},\n  \
OKtyp={\"Bar2\",\"Bar3\",\"Quad4\",\"Quad4.1\",\"Quad8\",\"Quad8.2\",\"Quad9\"\
,\n         \
\"Quad9.2\",\"Trig3\",\"Trig6\",\"Trig6.-3\",\"Trig10\",\"Trig10.6\"};\n  \
OKenl= {2,3,4,4,8,8,9,9,3,6,6,10,10};\n  K=Table[0,{2*numnod},{2*numnod}];  \
numer=prcopt[[1]];\n  For [e=1,e<=numele,e++, type=eletyp[[e]];\n       If \
[!MemberQ[OKtyp,type], Print[\"Illegal type:\",type,\n       \" element \
e=\",e,\" Assembly aborted\"]; Return[Null] ];\n       enl=elenod[[e]]; \
n=Length[enl]; {{j}}=Position[OKtyp,type];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n       If [OKenl[[j]]!=n, Print [\"Wrong node list length, element=\",e,\
\n           \" Assembly aborted\"]; Return[Null] ];\n       \
eft=Flatten[Table[{2*enl[[i]]-1,2*enl[[i]]},{i,1,n}]];\n       \
ncoor=Table[nodxyz[[enl[[i]]]],{i,1,n}]; \n",
  StyleBox["       (*Print[\"e=\",e,\" type=\",type,\"  enl=\",enl,\"  \
eft=\",eft]; *)  ",
    FontColor->RGBColor[1, 0, 0]],
  "      \n       If [type==\"Bar2\",     Em=elemat[[e]]; A=elefab[[e]];\n    \
       Ke=PlaneBar2Stiffness[ncoor,Em,A,{numer}] ];\n       If \
[type==\"Bar3\",     Em=elemat[[e]]; A=elefab[[e]];\n           \
Ke=PlaneBar3Stiffness[ncoor,Em,A,{numer}] ];\n       If [type==\"Quad4\",    \
Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Quad4IsoPMembraneStiffness[ncoor,Emat,th,{numer,2}] ];\n       If \
[type==\"Quad4.1\",  Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Quad4IsoPMembraneStiffness[ncoor,Emat,th,{numer,1}] ];\n       If \
[type==\"Quad8\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Quad8IsoPMembraneStiffness[ncoor,Emat,th,{numer,3}] ];\n       If \
[type==\"Quad8.2\",  Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Quad8IsoPMembraneStiffness[ncoor,Emat,th,{numer,2}] ];\n       If \
[type==\"Quad9\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Quad9IsoPMembraneStiffness[ncoor,Emat,th,{numer,3}] ];\n       If \
[type==\"Quad9.2\",  Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Quad9IsoPMembraneStiffness[ncoor,Emat,th,{numer,2}] ];\n       If \
[type==\"Trig3\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Trig3IsoPMembraneStiffness[ncoor,Emat,th,{numer}] ]; \n       If \
[type==\"Trig6\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Trig6IsoPMembraneStiffness[ncoor,Emat,th,{numer,3}] ]; \n       If [type==\
\"Trig6.-3\", Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Trig6IsoPMembraneStiffness[ncoor,Emat,th,{numer,-3}] ]; \n       If \
[type==\"Trig10\",   Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Trig10IsoPMembraneStiffness[ncoor,Emat,th,{numer,7}] ]; \n       If \
[type==\"Trig10.6\", Emat=elemat[[e]]; th=elefab[[e]];\n           \
Ke=Trig10IsoPMembraneStiffness[ncoor,Emat,th,{numer,6}] ];\n       \
neldof=Length[Ke];",
  StyleBox[" (*Print[\"Ke=\",SetPrecision[Chop[Ke],5]//MatrixForm];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n       For [i=1,i<=neldof,i++, ii=eft[[i]];\n           For \
[j=i,j<=neldof,j++, jj=eft[[j]];\n                \
K[[jj,ii]]=K[[ii,jj]]+=Ke[[i,j]] ];\n       ];\n     ];  Return[K]; \n];\n\n\
",
  StyleBox["ClearAll[Em,\[Nu],h,exx,eyy,gxy]; h=1; Em=288; \[Nu]=1/3;\n\
ncoor={{0,0},{6,2},{4,4},{3,1},{5,3},{2,2},{8,8}};\nEmat=Em/(1-\[Nu]^2)*{{1,\
\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\nnodxyz= \
{{0,4},{0,0},{4,4},{4,0},{8,4},{8,0},{12,2}};\nelenod=  \
{{1,2,4,3},{3,4,6,5},{5,6,7}};\neletyp=  {\"Quad4.1\",\"Quad4.1\",\"Trig3\"};\
\nelemat=  {Emat,Emat,Emat};\nelefab=  {1,1,1}; eleopt={True};\n\
K=PlaneStressMasterStiffness[nodxyz,eletyp,elenod,\n          \
elemat,elefab,{True}];\nPrint[\"Master Stiffness Matrix (twice rank deficient \
because of Quad4.1):\"]; \nPrint[Chop[K]//MatrixForm]; \n\
Print[Chop[Eigenvalues[N[K]]]];\n\nClearAll[Em,\[Nu],h,A];\nEm=10000; ",
    FontColor->RGBColor[0, 0, 1]],
  "\[Nu]",
  StyleBox["=1/4; h=3/10; A=2;\nEmat=Em/(1-",
    FontColor->RGBColor[0, 0, 1]],
  "\[Nu]",
  StyleBox["^2)*{{1,",
    FontColor->RGBColor[0, 0, 1]],
  "\[Nu]",
  StyleBox[",0},{",
    FontColor->RGBColor[0, 0, 1]],
  "\[Nu]",
  StyleBox[",1,0},{0,0,(1-",
    FontColor->RGBColor[0, 0, 1]],
  "\[Nu]",
  StyleBox[")/2}};\nnodxyz= {{0,1},{1,0},{5,0},{0,6},{5,6}};\nelenod= \
{{3,5,4},{2,3,4,1},{1,4},{1,2}};\neletyp= \
{\"Trig3\",\"Quad4\",\"Bar2\",\"Bar2\"};\nelemat= {Emat,Emat,Em,Em};\nelefab= \
{h,h,A,A}; prcopt={True};\n\
K=PlaneStressMasterStiffness[nodxyz,eletyp,elenod,elemat,elefab,prcopt];\n\
Print[\"Master Stiffness Matrix:\"]; \n\
Print[SetPrecision[Chop[K],5]//MatrixForm]; \nPrint[\"Eigs of \
K:\",Chop[Eigenvalues[N[K]]]];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 4:  Application of displacement BC  on master stiffness \
equations.  Code takes only single DOF 
constraints, but accepts NZ prescribed displacements. Self contained.\
\>", \
"Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ModifiedMasterStiffness[nodtag_,K_] := Module[\n  \
{i,j,k,n=Length[K],pdof,np,Kmod=K},\n  pdof=PrescDisplacementDOFTags[nodtag]; \
np=Length[pdof];\n  For [k=1,k<=np,k++, i=pdof[[k]]; \n      For \
[j=1,j<=n,j++, Kmod[[i,j]]=Kmod[[j,i]]=0];\n      Kmod[[i,i]]=1]; \n  \
Return[Kmod]];  \n  \nModifiedNodeForces[nodtag_,nodval_,K_,f_]:= Module[\n  \
{i,j,k,n=Length[K],pdof,pval,np,d,c,fmod=f},\n  \
pdof=PrescDisplacementDOFTags[nodtag]; np=Length[pdof];\n  \
pval=PrescDisplacementDOFValues[nodtag,nodval]; c=Table[1,{n}]; \n  For \
[k=1,k<=np,k++, i=pdof[[k]]; c[[i]]=0];\n  For [k=1,k<=np,k++, i=pdof[[k]]; \
d=pval[[k]]; \n      fmod[[i]]=d; If [d==0, Continue[]];\n      For \
[j=1,j<=n,j++, fmod[[j]]-=K[[i,j]]*c[[j]]*d]; \n      ];  \n  Return[fmod]];\n\
  \nPrescDisplacementDOFTags[nodtag_]:= Module [\n \
{j,n,numnod=Length[nodtag],pdof={},k=0,m},\n  For [n=1,n<=numnod,n++, \
m=Length[nodtag[[n]]];\n      For [j=1,j<=m,j++, If [nodtag[[n,j]]>0, \n      \
          AppendTo[pdof,k+j]];\n         ]; k+=m;\n     ]; Return[pdof]]; \n  \
         \nPrescDisplacementDOFValues[nodtag_,nodval_]:= Module [\n \
{j,n,numnod=Length[nodtag],pval={},k=0,m},\n  For [n=1,n<=numnod,n++, \
m=Length[nodtag[[n]]];\n      For [j=1,j<=m,j++, If [nodtag[[n,j]]>0, \n      \
          AppendTo[pval,nodval[[n,j]]]];\n         ]; k+=m;\n     ]; \
Return[pval]];\n\n",
  StyleBox["ClearAll[K,f,v1,v2,v4]; Km=Array[K,{6,6}];\nPrint[\"Master \
Stiffness: \",Km//MatrixForm];\nnodtag={{1,1},{0,1},{0,0}}; \
nodval={{v1,v2},{0,v4},{0,0}};\nKmod=ModifiedMasterStiffness[nodtag,Km];\n\
Print[\"Modified Master Stiffness:\",Kmod//MatrixForm];\nfm=Array[f,{6}]; \
Print[\"Master Force Vector:\",fm];\n\
fmod=ModifiedNodeForces[nodtag,nodval,Km,fm];\nPrint[\"Modified Force Vector:\
\",fmod//MatrixForm];\n\n(*nodtag= {{0,1,0},{0,0,0},{1,0,1},{0,1,1},{0,0,1}};\
\nnodval=-{{1,2,3},{4,5,6},{7,8,9},{10,11,12},{13,14,15}};\n\
Print[PrescDisplacementDOFTags[nodtag]];\n\
Print[PrescDisplacementDOFValues[nodtag,nodval]];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 5:  Stress recovery in plate elements &  internal force \
recovery in bars\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PlaneStressPlateStresses[nodxyz_,eletyp_,elenod_,elemat_,\n  \
elefab_,noddis_,prcopt_]:=Module[{numele=Length[elenod],\n  \
numnod=Length[nodxyz],ncoor,type,e,enl,i,k,n,Emat,th,\n  \
numer=True,nodpnc,nodsig}, If [Length[prcopt]>0, numer=prcopt[[1]]]; \n  \
nodpnc=Table[0,{numnod}]; nodsig=Table[{0,0,0},{numnod}];\n  If \
[Length[prcopt]>=1, numer=prcopt[[1]]]; \n  For [e=1,e<=numele,e++, \
type=eletyp[[e]];\n       enl=elenod[[e]]; k=Length[enl]; \n       \
ncoor=Table[nodxyz[[enl[[i]]]],{i,1,k}];\n       \
ue=Table[{noddis[[enl[[i]],1]],noddis[[enl[[i]],2]]},{i,1,k}];\n       \
ue=Flatten[ue]; \n   ",
  StyleBox["    (*Print[\"e=\",e,\" type=\",type,\"  enl=\",enl,\"  \
ue=\",ue]; *) ",
    FontColor->RGBColor[1, 0, 0]],
  "       \n       If [StringTake[type,3]==\"Bar\", Continue[]];\n       If \
[type==\"Quad4\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Quad4IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n       If \
[type==\"Quad4.1\",  Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Quad4IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n       If \
[type==\"Quad8\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Quad8IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n       If \
[type==\"Quad8.2\",  Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Quad8IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n       If \
[type==\"Quad9\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Quad9IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n       If \
[type==\"Quad9.2\",  Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Quad9IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n       If \
[type==\"Trig3\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Trig3IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; \n       If \
[type==\"Trig6\",    Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Trig6IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; \n       If \
[type==\"Trig6.-3\", Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Trig6IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; \n       If \
[type==\"Trig10\",   Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Trig10IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; \n       If \
[type==\"Trig10.6\", Emat=elemat[[e]]; th=elefab[[e]];\n           \
sige=Trig10IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];\n",
  StyleBox["       (*Print[\"sige=\",sige//MatrixForm]; *)     ",
    FontColor->RGBColor[1, 0, 0]],
  "\n       Do [n=enl[[i]]; nodpnc[[n]]++; nodsig[[n]]+=sige[[i]], {i,1,k}]; \
\n      ];\n  Do [k=nodpnc[[n]]; If [k>1,nodsig[[n]]/=k], {n,1,numnod}];\n  \
If [numer, nodsig=Chop[nodsig]];\n  Return[{nodpnc,nodsig}]; \n];\n\n\
PlaneStressBarForces[nodxyz_,eletyp_,elenod_,elemat_,\n  \
elefab_,noddis_,prcopt_]:=Module[{numele=Length[elenod],\n  \
numnod=Length[nodxyz],ncoor,type,e,enl,\n  \
i,k,n,Em,A,numer=True,pe,barele={},barfor={}},\n  If [Length[prcopt]>=1, \
numer=prcopt[[1]] ];\n  For [e=1,e<=numele,e++, type=eletyp[[e]];\n       If \
[StringTake[type,3]!=\"Bar\", Continue[]];\n       enl=elenod[[e]]; \
k=Length[enl]; \n       ncoor=Table[nodxyz[[enl[[i]]]],{i,1,k}];\n       \
ue=Table[{noddis[[enl[[i]],1]],noddis[[enl[[i]],2]]},{i,1,k}];\n       \
ue=Flatten[ue]; \n",
  StyleBox["       (*Print[\"e=\",e,\" type=\",type,\"  enl=\",enl,\"  \
ue=\",ue];*) ",
    FontColor->RGBColor[1, 0, 0]],
  "    \n       If [type==\"Bar2\", bar=True; Em=elemat[[e]]; A=elefab[[e]];\n\
           pe=PlaneBar2IntForce[ncoor,Em,A,ue,{numer}] ];\n       If \
[type==\"Bar3\", bar=True; Em=elemat[[e]]; A=elefab[[e]];\n           \
pe=PlaneBar3IntForce[ncoor,Em,A,ue,{numer}] ];\n       AppendTo[barele,e]; \
AppendTo[barfor,pe]; \n      ];\n   If [numer, nodsig=Chop[barfor]];\n   \
Return[{barele,barfor}]; \n];\n",
  StyleBox["   \nClearAll[Em,\[Nu],h,exx,eyy,gxy,\[Sigma]xx,\[Sigma]yy,\
\[Sigma]xy]; \nEm=300; \[Nu]=0; h=1; A=10; Em=300; \nEmat=Em/(1-\[Nu]^2)*{{1,\
\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};\n\
{exx,eyy,gxy}=Simplify[Inverse[Emat].{\[Sigma]xx,\[Sigma]yy,\[Sigma]xy}];\n\
nodxyz={{0,5},{0,2},{0,1},{0,0},{3,5},{2,2},{5/2,1},\n        \
{2,-1},{4,4},{4,3/2},{4,-2},{5,-1/2},{6,1},{5,5/2}};\n\
elenod={{1,2,6,5},{5,6,9},{2,4,11,9,3,8,10,6,7},\n        \
{9,11,13,10,12,14},{1,2},{2,6}};\n\
eletyp={\"Quad4\",\"Trig3\",\"Quad9\",\"Trig6\",\"Bar2\",\"Bar2\"};\n\
numnod=Length[nodxyz]; numele=Length[elenod]; prcopt={False};\n\
x=Table[nodxyz[[i,1]],{i,numnod}]; y=Table[nodxyz[[i,2]],{i,numnod}];\n\
noddis=Table[{x[[i]]*exx+y[[i]]*gxy/2,x[[i]]*gxy/2+y[[i]]*eyy},{i,numnod}];\n\
elemat={Emat,Emat,Emat,Emat,Em,Em};  elefab={h,h,h,h,A,A};\n\
{nodpnc,nodsig}=PlaneStressPlateStresses[nodxyz,eletyp,elenod,elemat,\n   \
elefab,noddis,prcopt];\nPrint[\"Plate node counts=\",nodpnc];\nPrint[\"Plate \
stresses=\",Transpose[nodsig]//MatrixForm];\n\
(*{barele,barfor}=PlaneStressBarForces[nodxyz,eletyp,elenod,elemat,\n   \
elefab,noddis,prcopt];\nPrint[\"barele=\",barele,\"  barfor=\",barfor];*) ",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 6:  Mesh plot modules.  Written by C. A. Felippa,1997/98.  \
Mathematica adaptation of ancient  (1966) Fortran code.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {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    \
AppendTo[poly,Graphics[Polygon[Take[xyc,nc]]]];\n    \
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      \
DisplayFunction->DisplayChannel ];\n ClearAll[poly,sides];\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,",
  StyleBox["style,",
    FontColor->RGBColor[0, 0, 1]],
  "\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",
  StyleBox["   (*Print[\"delx,dely,rade=\",{delx,dely,rade}];*)",
    FontColor->RGBColor[1, 0, 0]],
  " \n",
  StyleBox["   style=TextStyle->{FontFamily->\"Times\",FontSize->11,\n        \
             FontWeight->\"Plain\",FontSlant->\"Plain\"};",
    FontColor->RGBColor[0, 0, 1]],
  "\n   For [e=1,e<=numele,e++, enl=elenod[[e]]; n=Length[enl];  nb=0;\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       If [n==10,b={1,4,5,2,6,7,3,8,9}; nb=9];\n       If [nb==0, \
Continue[]];\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       \
AppendTo[pefill,Graphics[Polygon[gpoly]]];\n       AppendTo[pesid, \
Graphics[Line[gside]]];\n       If [n==9,xy0+={0,-2*rade}];\n       If \
[elabels, AppendTo[pedisk,Graphics[Disk[xy0,rade]]];\n           \
AppendTo[pecirc,Graphics[Circle[xy0,rade]]];\n",
  StyleBox["           elab=ToString[e];\n           \
AppendTo[pelab,Graphics[Text[elab,xy0,style]]]];",
    FontColor->RGBColor[0, 0, 1]],
  "\n          ];\n   For [n=1,n<=numnod,n++,\n       If [nlabels, ",
  StyleBox[" nlab=ToString[n]; AppendTo[\n             \
pnlab,Graphics[Text[nlab,xycoor[[n]]+{ex,ey},style]]]];",
    FontColor->RGBColor[0, 0, 1]],
  "\n       AppendTo[pnod, Graphics[Disk[xycoor[[n]],radn]]]\n       ];  \
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       DisplayFunction->DisplayChannel];\n \
ClearAll[pefill,pesid,pnod,pedisk,pnlab];\n ];\n",
  "\n",
  StyleBox["\n\
(*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\
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];*)\n",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

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

Cell[TextData[{
  "DisplayChannel=$DisplayFunction;\nIf [$VersionNumber>=6.0, \
DisplayChannel=Print]; \n  (* fix for Mathematica 6 & later *)\n\n\
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=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}];  ",
  StyleBox["(*Print[\"bseg=\",bseg];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\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",
  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\
sxx=Table[8*nodxyz[[n,1]]/5-4,{n,1,9}]; sxxmax=4; sxxmin=-4;\n sxxinc=.5; \
aspect=6/5; Print[\"sxx=\",sxx];\n\
ContourBandPlotNodeFuncOver2DMesh[nodxyz,elenod,sxx,{sxxmin,sxxmax,sxxinc},\n \
 {True,True,True,True,True,True},{2,2.5},aspect,\"Sigma-xx (white: zero)\"];\n\
syy=Table[4*nodxyz[[n,2]]/6,{n,1,9}]; syymax=4; syymin=-4;\nsyyinc=.5; \
aspect=6/5; Print[\"syy=\",syy];\n\
ContourBandPlotNodeFuncOver2DMesh[nodxyz,elenod,syy,{syymin,syymax,syyinc},\n \
 {True,True,True,True,True,True},{2,2.5},aspect,\"Sigma-yy (white: zero)\"];",
    
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["Cell 8: Print utilities ", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
PrintPlaneStressNodeCoordinates[nodxyz_,title_,digits_]:= Module[
  {numnod=Length[nodxyz],n,x,y,z,d=6,f=6,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++,{x,y}=nodxyz[[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\", \"x-coor\", \"y-coor\"}}]];
   ClearAll[tab];
   ];
   
PrintPlaneStressElementTypeNodes[eletyp_,elenod_,title_,digits_]:=
  Module[{e,numele=Length[elenod],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];
   ];
   
PrintPlaneStressElementMatFab[elemat_,elefab_,title_,digits_]:= 
  Module[{e,numele=Length[elemat],mat,fab,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]]; fab=elefab[[e]];
        tab[[e]]={ToString[e],ToString[mat],ToString[fab]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->Right,
         TableDirections->{Column,Row},TableSpacing->{0,2},
         TableHeadings->{None,{\"elem\", \"material\", \"fabrication\"}}]];
   ClearAll[tab];
   ];
   
PrintPlaneStressFreedomActivity[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\", \"x-tag\", \"y-tag\",
                                \"x-value\", \"y-value\"}} ]];
   ClearAll[tab];
   ];

PrintPlaneStressNodeDisplacements[noddis_,title_,digits_]:= Module[
  {numnod=Length[noddis],n,ux,uy,d=8,f=6,tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits]; 
   For [n=1,n<=numnod,n++, {ux,uy}=noddis[[n]];
       tab[[n]]={ToString[n],PaddedForm[ux,{d,f}],PaddedForm[uy,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings->{None,{\"node\", \"x-displ\", \"y-displ\"}}]];
   ClearAll[tab];
   ];
   
PrintPlaneStressNodeForces[nodfor_,title_,digits_]:= Module[
  {numnod=Length[nodfor],n,fx,fy,tab,d=8,f=4}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++, {fx,fy}=nodfor[[n]]; 
       tab[[n]]={ToString[n],PaddedForm[fx,{d,f}],PaddedForm[fy,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings ->{None,{\"node\", \"x-force\", \"y-force\"}}]];       \
  
   ClearAll[tab];
   ];
   
PrintPlaneStressPlateNodeStresses[nodpnc_,nodsig_,title_,digits_]:= 
  Module[{numnod=Length[nodpnc],n,k,sxx,syy,sxy,d=8,f=4,tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
  For [n=1,n<=numnod,n++, k=nodpnc[[n]]; {sxx,syy,sxy}=nodsig[[n]];
       If [k>0, tab[[n]]={ToString[n], ToString[k],
                 PaddedForm[sxx,{d,f}],PaddedForm[syy,{d,f}],
                 PaddedForm[sxy,{d,f}]} ]; 
       If [k==0, tab[[n]]={ToString[n], \"no plates\",\" attached\"}];
       ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{\"node\", \"cnt\", \"sigma-xx\", 
          \"sigma-yy\", \"sigma-xy\"}}]];
   ClearAll[tab];
   ]; 
   
PrintPlaneStressSolution[noddis_,nodfor_,nodpnc_,nodsig_,title_,digits_]:= 
  Module[{numnod=Length[nodpnc],n,k,sxx,syy,sxy,d=6,f=4,tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
  For [n=1,n<=numnod,n++,{ux,uy}=noddis[[n]]; {fx,fy}=nodfor[[n]]; 
       k=nodpnc[[n]]; {sxx,syy,sxy}=nodsig[[n]];
       If [k>0, tab[[n]]={ToString[n], 
                 PaddedForm[ux,{d,f}],PaddedForm[uy,{d,f}],
                 PaddedForm[fx,{d,f}],PaddedForm[fy,{d,f}],
                 PaddedForm[sxx,{d,f}],PaddedForm[syy,{d,f}],
                 PaddedForm[sxy,{d,f}]} ]; 
       If [k==0, tab[[n]]={ToString[n], 
                 PaddedForm[ux,{d,f}],PaddedForm[uy,{d,f}],
                 PaddedForm[fx,{d,f}],PaddedForm[fy,{d,f}]} ];
       ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{\"node\", \"x-displ\", \
\"y-displ\",\"x-force\",
         \"y-force\",  \"sigma-xx\", \"sigma-yy\", \"sigma-xy\"}}]];
   ClearAll[tab];
   ];  \
\>", "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 9:  Miscellaneous utilities\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "FlatNodePartVector[nv_]:=Flatten[nv];\n           \n\
NodePartFlatVector[nfc_,v_]:= Module [\n  {i,k,m,n,nv={},numnod},\n  If \
[Length[nfc]==0, nv=Partition[v,nfc]];\n  If [Length[nfc]>0,  \
numnod=Length[nfc]; m=0;\n      nv=Table[0,{numnod}];   \n      For \
[n=1,n<=numnod,n++, k=nfc[[n]]; \n          nv[[n]]=Table[v[[m+i]],{i,1,k}]; \
\n          m+=k]];\n  Return[nv]];",
  StyleBox["  ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["\n    ",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 10:       Driver for static analysis of plane stress \
problem.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PlaneStressSolution[nodxyz_,eletyp_,elenod_,elemat_,elefab_,\n   \
nodtag_,nodval_,prcopt_]:= Module[{K,Kmod,f,fmod,u,numer=True,\n   \
noddis,nodfor,nodpnc,nodsig,barele,barfor},\n   If [Length[prcopt]>=1, \
numer=prcopt[[1]]];\n   \
K=PlaneStressMasterStiffness[nodxyz,eletyp,elenod,elemat,elefab,prcopt];\n   \
If [K==Null,Return[Table[Null,{6}]]];\n",
  StyleBox["  (* Print[\"K=\",K//MatrixForm];",
    FontColor->RGBColor[1, 0, 0]],
  "\n",
  StyleBox["   Print[\"eigs of K=\",Chop[Eigenvalues[N[K]]]];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   Kmod=ModifiedMasterStiffness[nodtag,K];\n   \
f=FlatNodePartVector[nodval];\n   fmod=ModifiedNodeForces[nodtag,nodval,K,f];\
\n",
  StyleBox["   (* Print[\"eigs of Kmod=\",Chop[Eigenvalues[N[Kmod]]]];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   u=LinearSolve[Kmod,fmod]; If [numer,u=Chop[u]]; \n   f=K.u;           \
         If [numer,f=Chop[f]];  \n   nodfor=NodePartFlatVector[2,f]; \
noddis=NodePartFlatVector[2,u];\n   \
{nodpnc,nodsig}=PlaneStressPlateStresses[nodxyz,eletyp,elenod,elemat,\n       \
            elefab,noddis,prcopt];\n   \
{barele,barfor}=PlaneStressBarForces[nodxyz,eletyp,elenod,elemat,\n           \
        elefab,noddis,prcopt];\n   ClearAll[K,Kmod];\n   \
Return[{noddis,nodfor,nodpnc,nodsig,barele,barfor}]; \n];"
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 11:  Main program for verifying  4-node quad implementation \
with one element\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell["\<\
ClearAll[Em,\[Nu],th];  
Em=10000; \[Nu]=.25;  th=3; aspect=6/5; Nsub=4;
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};

(*  Define FEM model *)

NodeCoordinates=N[{{0,6},{0,0},{5,6},{5,0}}];
PrintPlaneStressNodeCoordinates[NodeCoordinates,\"\",{6,4}];       
ElemNodes=  {{1,2,4,3}};
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes]; 
ElemTypes=       Table[\"Quad4\",{numele}]; 
PrintPlaneStressElementTypeNodes[ElemTypes,ElemNodes,\"\",{}];
ElemMaterials=   Table[Emat,   {numele}]; 
ElemFabrications=Table[th,     {numele}];
PrintPlaneStressElementMatFab[ElemMaterials,ElemFabrications,\"\",{}];
NodeDOFValues=NodeDOFTags=Table[{0,0},{numnod}];
NodeDOFValues[[1]]=NodeDOFValues[[3]]={0,75}; (* nodal loads *)
NodeDOFTags[[1]]={1,0};     (* vroller @ node 1 *)
NodeDOFTags[[2]]={1,1};     (* fixed node 2     *)
NodeDOFTags[[4]]={0,1};     (* hroller @ node 4 *)
PrintPlaneStressFreedomActivity[NodeDOFTags,NodeDOFValues,\"\",{}];
ProcessOptions={True};
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"One element mesh - 4-node quad\",True,True];
    
(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodePlateCounts,NodePlateStresses,
    ElemBarNumbers,ElemBarForces}= PlaneStressSolution[
    NodeCoordinates,ElemTypes,ElemNodes,
    ElemMaterials,ElemFabrications,
    NodeDOFTags,NodeDOFValues,ProcessOptions];

PrintPlaneStressSolution[NodeDisplacements,NodeForces,NodePlateCounts,
      NodePlateStresses,\"Computed Solution:\",{}];

(*  Plot Averaged Nodal Stresses Distribution *)

sigeps=10.^(-3); nbands=10;
sxx=Table[NodePlateStresses[[n,1]],{n,numnod}]; 
syy=Table[NodePlateStresses[[n,2]],{n,numnod}];
sxy=Table[NodePlateStresses[[n,3]],{n,numnod}];
{sxxmax,syymax,sxymax}=Abs[{Max[sxx],Max[syy],Max[sxy]}]+sigeps; 
{sxxmin,syymin,sxymin}=Abs[{Min[sxx],Min[syy],Min[sxy]}]+sigeps;
sxxmax=Max[sxxmax,sxxmin];  sxxmin=-sxxmax;
syymax=Max[syymax,syymin];  syymin=-syymax;
sxymax=Max[sxymax,sxymin];  sxymin=-sxymax;
{sxxinc,syyinc,sxyinc}={sxxmax-sxxmin,syymax-syymin,sxymax-sxymin}/nbands; 
Print[\"sxxmin,sxxmax,sxxinc=\",{sxxmin,sxxmax,sxxinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxx,{sxxmin,\
sxxmax,sxxinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-xx\"];
Print[\"syymin,syymax,syyinc=\",{syymin,syymax,syyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,syy,{syymin,\
syymax,syyinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-yy\"];
Print[\"sxymin,sxymax,sxyinc=\",{sxymin,sxymax,sxyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxy,{sxymin,\
sxymax,sxyinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-xy\"];\
\>", \
"Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 12:  Main program for verifying  9-node quad implementation \
with one element\
\>", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell["\<\
ClearAll[Em,\[Nu],th];  
Em=10000; \[Nu]=.25;  th=3; aspect=6/5; Nsub=4;
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};

(*  Define FEM model *)

NodeCoordinates=N[{{0,6},{0,3},{0,0},{5/2,6},{5/2,3},
     {5/2,0},{5,6},{5,3},{5,0}}];  
PrintPlaneStressNodeCoordinates[NodeCoordinates,\"\",{6,4}];       
ElemNodes=  {{1,3,9,7,2,6,8,4,5}};
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes]; 
ElemTypes=       Table[\"Quad9.2\",{numele}]; 
PrintPlaneStressElementTypeNodes[ElemTypes,ElemNodes,\"\",{}];
ElemMaterials=   Table[Emat,   {numele}]; 
ElemFabrications=Table[th,     {numele}];
PrintPlaneStressElementMatFab[ElemMaterials,ElemFabrications,\"\",{}];
NodeDOFValues=NodeDOFTags=Table[{0,0},{numnod}];
NodeDOFValues[[1]]=NodeDOFValues[[7]]={0,25}; 
NodeDOFValues[[4]]={0,100};              (* nodal loads *)
NodeDOFTags[[1]]=NodeDOFTags[[2]]={1,0}; (* vroller @ nodes 1,2 *)
NodeDOFTags[[3]]={1,1};                  (* fixed node 3 *)
NodeDOFTags[[6]]=NodeDOFTags[[9]]={0,1}; (* hroller @ nodes 6,9 *)
PrintPlaneStressFreedomActivity[NodeDOFTags,NodeDOFValues,\"\",{}];
ProcessOptions={True};
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"One element mesh - 4-node quad\",True,True];
   
(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodePlateCounts,NodePlateStresses,
    ElemBarNumbers,ElemBarForces}= PlaneStressSolution[
    NodeCoordinates,ElemTypes,ElemNodes,
    ElemMaterials,ElemFabrications,
    NodeDOFTags,NodeDOFValues,ProcessOptions];
    
PrintPlaneStressSolution[NodeDisplacements,NodeForces,NodePlateCounts,
      NodePlateStresses,\"Computed Solution:\",{}]; 
   
(*  Plot Displacement Components Distribution - skipped *)

(* ux=Table[NodeDisplacements[[n,1]],{n,numnod}];
uy=Table[NodeDisplacements[[n,2]],{n,numnod}];
{uxmax,uymax}=Abs[{Max[ux],Max[uy]}];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,ux,
       uxmax,Nsub,aspect,\"Displacement component ux\"];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,uy,
       uymax,Nsub,aspect,\"Displacement component uy\"]; *)

(*  Plot Averaged Nodal Stresses Distribution *)

sigeps=10.^(-3); nbands=10;
sxx=Table[NodePlateStresses[[n,1]],{n,numnod}]; 
syy=Table[NodePlateStresses[[n,2]],{n,numnod}];
sxy=Table[NodePlateStresses[[n,3]],{n,numnod}];
{sxxmax,syymax,sxymax}=Abs[{Max[sxx],Max[syy],Max[sxy]}]+sigeps; 
{sxxmin,syymin,sxymin}=Abs[{Min[sxx],Min[syy],Min[sxy]}]+sigeps;
sxxmax=Max[sxxmax,sxxmin];  sxxmin=-sxxmax;
syymax=Max[syymax,syymin];  syymin=-syymax;
sxymax=Max[sxymax,sxymin];  sxymin=-sxymax;
{sxxinc,syyinc,sxyinc}={sxxmax-sxxmin,syymax-syymin,sxymax-sxymin}/nbands; 
Print[\"sxxmin,sxxmax,sxxinc=\",{sxxmin,sxxmax,sxxinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxx,{sxxmin,\
sxxmax,sxxinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-xx\"];
Print[\"syymin,syymax,syyinc=\",{syymin,syymax,syyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,syy,{syymin,\
syymax,syyinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-yy\"];
Print[\"sxymin,sxymax,sxyinc=\",{sxymin,sxymax,sxyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxy,{sxymin,\
sxymax,sxyinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-xy\"];\
\>", \
"Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Cell 13. Plate with circular hole, 4-node bilinear quad model\
\>", \
"Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell["\<\
ClearAll[Em,\[Nu],th,aspect,Nsub];  
Em=10000; \[Nu]=.25;  th=3; aspect=6/5; Nsub=4;
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};

(*  Define FEM model *)

s={1,0.70,0.48,0.30,0.16,0.07,0.0};
xy1={0,6}; xy7={0,1}; xy8={2.5,6}; xy14={Cos[3*Pi/8],Sin[3*Pi/8]}; 
xy8={2.5,6}; xy21={Cos[Pi/4],Sin[Pi/4]}; xy15={5,6}; 
xy22={5,2};  xy28={Cos[Pi/8],Sin[Pi/8]}; xy29={5,0}; xy35={1,0};
NodeCoordinates=Table[{0,0},{35}];
Do[NodeCoordinates[[n]]=N[s[[n]]   *xy1+(1-s[[n]])    *xy7],{n,1,7}]; 
Do[NodeCoordinates[[n]]=N[s[[n-7]] *xy8+(1-s[[n-7]])  *xy14],{n,8,14}];
Do[NodeCoordinates[[n]]=N[s[[n-14]]*xy15+(1-s[[n-14]])*xy21],{n,15,21}];
Do[NodeCoordinates[[n]]=N[s[[n-21]]*xy22+(1-s[[n-21]])*xy28],{n,22,28}];
Do[NodeCoordinates[[n]]=N[s[[n-28]]*xy29+(1-s[[n-28]])*xy35],{n,29,35}];
PrintPlaneStressNodeCoordinates[NodeCoordinates,\"\",{6,4}];          
ElemNodes=Table[{0,0,0,0},{24}];
ElemNodes[[1]]={1,2,9,8};
Do [ElemNodes[[e]]=ElemNodes[[e-1]]+{1,1,1,1},{e,2,6}];
ElemNodes[[7]]=ElemNodes[[6]]+{2,2,2,2};
Do [ElemNodes[[e]]=ElemNodes[[e-1]]+{1,1,1,1},{e,8,12}];
ElemNodes[[13]]=ElemNodes[[12]]+{2,2,2,2};
Do [ElemNodes[[e]]=ElemNodes[[e-1]]+{1,1,1,1},{e,14,18}];
ElemNodes[[19]]=ElemNodes[[18]]+{2,2,2,2};
Do [ElemNodes[[e]]=ElemNodes[[e-1]]+{1,1,1,1},{e,20,24}];
numnod=Length[NodeCoordinates]; numele=Length[ElemNodes]; 

ElemTypes=       Table[\"Quad4\",{numele}]; 
PrintPlaneStressElementTypeNodes[ElemTypes,ElemNodes,\"\",{}];
ElemMaterials=   Table[Emat,   {numele}]; 
ElemFabrications=Table[th,     {numele}];
(*PrintPlaneStressElementMatFab[ElemMaterials,ElemFabrications,\"\",{}];*)
NodeDOFValues=NodeDOFTags=Table[{0,0},{numnod}];
NodeDOFValues[[1]]=NodeDOFValues[[15]]={0,37.5}; 
NodeDOFValues[[8]]={0,75}; (* nodal loads *)
Do[NodeDOFTags[[n]]={1,0},{n,1,7}];   (* vroller @ nodes 1-7 *)
Do[NodeDOFTags[[n]]={0,1},{n,29,35}]; (* hroller @ node 4 *)
PrintPlaneStressFreedomActivity[NodeDOFTags,NodeDOFValues,\"\",{}];
ProcessOptions={True};
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"One element mesh - 4-node quad\",True,True];
(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodePlateCounts,NodePlateStresses,
    ElemBarNumbers,ElemBarForces}= PlaneStressSolution[
    NodeCoordinates,ElemTypes,ElemNodes,
    ElemMaterials,ElemFabrications,
    NodeDOFTags,NodeDOFValues,ProcessOptions];
   
PrintPlaneStressSolution[NodeDisplacements,NodeForces,NodePlateCounts,
      NodePlateStresses,\"Computed Solution:\",{}]; 
   
(*  Plot Displacement Components Distribution - skipped *)

(* ux=Table[NodeDisplacements[[n,1]],{n,numnod}];
uy=Table[NodeDisplacements[[n,2]],{n,numnod}];
{uxmax,uymax}=Abs[{Max[ux],Max[uy]}];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,ux,
       uxmax,Nsub,aspect,\"Displacement component ux\"];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,uy,
       uymax,Nsub,aspect,\"Displacement component uy\"]; *)

(*  Plot Averaged Nodal Stresses Distribution *)

sigeps=10.^(-3); nbands=10;
sxx=Table[NodePlateStresses[[n,1]],{n,numnod}]; 
syy=Table[NodePlateStresses[[n,2]],{n,numnod}];
sxy=Table[NodePlateStresses[[n,3]],{n,numnod}];
{sxxmax,syymax,sxymax}=Abs[{Max[sxx],Max[syy],Max[sxy]}]+sigeps; 
{sxxmin,syymin,sxymin}=Abs[{Min[sxx],Min[syy],Min[sxy]}]+sigeps;
sxxmax=Max[sxxmax,sxxmin];  sxxmin=-sxxmax;
syymax=Max[syymax,syymin];  syymin=-syymax;
sxymax=Max[sxymax,sxymin];  sxymin=-sxymax;
{sxxinc,syyinc,sxyinc}={sxxmax-sxxmin,syymax-syymin,sxymax-sxymin}/nbands; 
Print[\"sxxmin,sxxmax,sxxinc=\",{sxxmin,sxxmax,sxxinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxx,{sxxmin,\
sxxmax,sxxinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-xx\"];
Print[\"syymin,syymax,syyinc=\",{syymin,syymax,syyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,syy,{syymin,\
syymax,syyinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-yy\"];
Print[\"sxymin,sxymax,sxyinc=\",{sxymin,sxymax,sxyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxy,{sxymin,\
sxymax,sxyinc},
  {True,True,True,False,True,True},{2,2},aspect,\"Stress sigma-xy\"];\
\>", \
"Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["Cell 16 - Bridge with traveling load - animation", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  "ClearAll[Em,\[Nu],th,aspect,Nsub];  \nEm=10000; \[Nu]=.25;  th=3; \
aspect=1/2; Nsub=4;\nEmat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\
\[Nu])/2}};\n\n(*  Define FEM model *)\n\ns=N[{1,0.80,0.60,0.40,0.20,0.0}];\n\
xytop={{-60,64},{-50,64},{-40,64},{-30,64},{-20,64},\n       \
{-12,64},{-5,64},{0,64},{5,64},{12,64},{20,64},\n       \
{30,64},{40,64},{50,64},{60,64}}; \nxybot={{-60,0},{-50,0},{-40,18},{-30,32}, \
{-20,42},\n       {-12,47.12},{-5,49.5},{0,50},{5,49.5},{12,47.12},{20,42},\n \
      {30,32},{40,18},{50,0},{60,0}};\nnx=Length[xytop]; ny=6; \n\
numnod=nx*ny; NodeCoordinates=Table[{0,0},{numnod}];\nFor[k=1,k<=ny,k++,\n    \
For[i=1,i<=nx,i++,c=s[[k]];NodeCoordinates[[nx*(k-1)+i]]=\n                  \
c*xytop[[i]] +(1-c)*xybot[[i]]\n       ]; \n   ]; \
NodeCoordinates=N[NodeCoordinates]; \n\
(*PrintPlaneStressNodeCoordinates[NodeCoordinates,\"\",{6,4}];*)\n\
numele=(nx-1)*(ny-1);      \nElemNodes=Table[{0,0,0,0},{numele}];\n\
ElemNodes[[1]]={1,nx+1,nx+2,2};\nFor[e=2,e<=numele,e++,\n    If \
[Mod[e-1,nx-1]==0,\n        ElemNodes[[e]]=ElemNodes[[e-nx+1]]+{1,1,1,1}*nx,\n\
        ElemNodes[[e]]=ElemNodes[[e-1]]+{1,1,1,1}]];\nElemTypes=       Table[\
\"Quad4\",{numele}]; \n\
(*PrintPlaneStressElementTypeNodes[ElemTypes,ElemNodes,\"\",{}];*)\n\
ElemMaterials=   Table[Emat,   {numele}]; \nElemFabrications=Table[th,     \
{numele}];\nNodeDOFValues=NodeDOFTags=Table[{0,0},{numnod}];\n\n\
Do[NodeDOFTags[[n]]={1,0},{n,1,numnod-nx+1,nx}];  \n\
Do[NodeDOFTags[[n]]={1,0},{n,nx,numnod,nx}];\n\
NodeDOFTags[[numnod-nx+1]]=NodeDOFTags[[numnod]]={1,1}; \n\
NodeDOFTags[[numnod-nx+2]]=NodeDOFTags[[numnod-1]]={0,1};\n\
(*PrintPlaneStressFreedomActivity[NodeDOFTags,NodeDOFValues,\"\",{}];*)\n\
ProcessOptions={True};\nPlot2DElementsAndNodes[NodeCoordinates,ElemNodes,1/2,\
\n    \"Bridge Mesh\",True,True];\n\n  (*  Solve problem and print results *)\
\n  \nPrint[\"Frames for animation:\"];  nloads=nx;",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\nFor[load=1,load<=nloads,load++,\n  NodeDOFValues=Table[{0,0},{numnod}]; \
NodeDOFValues[[load]]={0,-1000}; \n \
{NodeDisplacements,NodeForces,NodePlateCounts,NodePlateStresses,\n    \
ElemBarNumbers,ElemBarForces}= PlaneStressSolution[\n    \
NodeCoordinates,ElemTypes,ElemNodes,\n    ElemMaterials,ElemFabrications,\n   \
 NodeDOFTags,NodeDOFValues,ProcessOptions];\n\n  (*  Plot Stresses \
Distribution *)\n\n  sVM=sig1=sig2=tmax=Table[0,{numnod}]; \n  \
For[n=1,n<=numnod,n++,\n   {sigxx,sigyy,sigxy}=NodePlateStresses[[n]];  \
s0=(sigxx+sigyy)/2; \n    R=N[Sqrt[(sigxx-sigyy)^2+sigxy^2]];\n    \
sig1[[n]]=s0+R; sig2[[n]]=s0-R; tmax[[n]]=R;\n    \
sVM[[n]]=N[Sqrt[3*(sigxx^2+sigyy^2+sigxy^2)/2]]];\n  sVMmax=sg1max=sg2max=0;\n\
  For[n=1,n<=numnod,n++,sg1max=Max[Abs[sig1[[n]]],sg1max];\n                  \
      sg2max=Max[Abs[sig2[[n]]],sg2max];\n                        \
sVMmax=Max[Abs[sVM[[n]]],sVMmax]; ];\n \n  (*Print[\"Max sig1,sig2,sigVM  \
=\",{sg1max,sg2max,sVMmax}];*)\n  sg1max=40;  sg1min=-sg1max; \
sg1inc=(sg1max-sg1min)/10;\n  \
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sig1,{sg1min,\
sg1max,sg1inc},\n  {True,False,False,False,True,True},{},aspect,\"Principal \
stress sig1\"];       \n  ];  "
}], "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  AnimationDisplayTime->0.28561,
  AnimationCycleOffset->1,
  AnimationCycleRepetitions->Infinity,
  Background->RGBColor[0, 1, 0]],

Cell["Cell 14. Circular plate with multiple holes", "Text",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell["\<\
nodecoor=
{{1,{-1.35925,-0.85227}},{2,{1.79965,-1.04346}},{3,{-0.65924,1.05955}},{4,{-1.\
50571,-0.49871}},
{5,{-1.85915,-0.35227}},{6,{-2.21288,-0.49880}},{7,{-2.35925,-0.85227}},{8,{-\
2.21288,-1.20575}},
{9,{-1.85915,-1.35227}},{10,{-1.50571,-1.20584}},{11,{1.65634,-0.46230}},{12,{\
1.25973,-0.01473}},
{13,{0.70048,0.19741}},{14,{0.10654,0.12537}},{15,{-0.38620,-0.21479}},{16,{-\
0.66404,-0.74436}},
{17,{-0.66404,-1.34256}},{18,{-0.38620,-1.87213}},{19,{0.10654,-2.21228}},{20,\
{0.70048,-2.28432}},
{21,{1.25973,-2.07219}},{22,{1.65634,-1.62462}},{23,{-0.79823,1.53275}},{24,{-\
1.17075,1.85548}},
{25,{-1.65854,1.92568}},{26,{-2.10735,1.72074}},{27,{-2.37381,1.30602}},{28,{-\
2.37381,0.81308}},
{29,{-2.10735,0.39836}},{30,{-1.65854,0.19343}},{31,{-1.17075,0.26363}},{32,{-\
0.79823,0.58636}},
{33,{2.65109,0.13232}},{34,{2.62740,0.50858}},{35,{2.55684,0.87839}},{36,{2.\
44042,1.23669}},
{37,{2.28001,1.57758}},{38,{2.07814,1.89567}},{39,{1.83788,2.18609}},{40,{1.\
56334,2.44388}},
{41,{1.25857,2.66530}},{42,{0.92843,2.84680}},{43,{0.57821,2.98546}},{44,{0.\
21348,3.07913}},
{45,{-0.16054,3.12640}},{46,{-0.53728,3.12640}},{47,{-0.91105,3.07918}},{48,{-\
1.27596,2.98549}},
{49,{-1.62612,2.84686}},{50,{-1.95642,2.66528}},{51,{-2.26118,2.44386}},{52,{-\
2.53581,2.18596}},
{53,{-2.77596,1.89567}},{54,{-2.97783,1.57758}},{55,{-3.13836,1.23638}},{56,{-\
3.25469,0.87827}},
{57,{-3.32526,0.50830}},{58,{-3.34891,0.13232}},{59,{-3.32526,-0.24367}},{60,{\
-3.25469,-0.61364}},
{61,{-3.13836,-0.97175}},{62,{-2.97783,-1.31295}},{63,{-2.77596,-1.63104}},{\
64,{-2.53581,-1.92133}},
{65,{-2.26118,-2.17922}},{66,{-1.95642,-2.40065}},{67,{-1.62612,-2.58222}},{\
68,{-1.27596,-2.72085}},
{69,{-0.91105,-2.81455}},{70,{-0.53728,-2.86176}},{71,{-0.16054,-2.86176}},{\
72,{0.21348,-2.81450}},
{73,{0.57821,-2.72083}},{74,{0.92843,-2.58216}},{75,{1.25857,-2.40067}},{76,{\
1.56334,-2.17925}},
{77,{1.83788,-1.92145}},{78,{2.07814,-1.63104}},{79,{2.28001,-1.31295}},{80,{\
2.44042,-0.97206}},
{81,{2.55684,-0.61375}},{82,{2.62740,-0.24395}},{83,{-2.08111,2.16378}},{84,{-\
1.91730,-2.05615}},
{85,{-1.88738,-1.72539}},{86,{-0.32616,0.34708}},{87,{-2.88234,0.64326}},{88,{\
-2.56946,0.36214}},
{89,{-2.38638,0.09053}},{90,{-2.11256,-0.06651}},{91,{-2.56112,-1.24946}},{92,\
{-2.17779,-1.89001}},
{93,{-2.28716,-1.56789}},{94,{-2.97534,0.23861}},{95,{-2.45103,-0.23032}},{96,\
{-2.69152,-0.01102}},
{97,{-3.03458,-0.20330}},{98,{-2.57050,-0.57231}},{99,{-2.73232,-0.32378}},{\
100,{-2.75351,-0.92678}},
{101,{-2.91399,-0.60193}},{102,{-1.62513,-1.91359}},{103,{-0.94916,-1.67251}},\
{104,{-0.34435,-2.53664}},
{105,{-1.59425,-2.22690}},{106,{-1.26905,-1.98716}},{107,{-1.17086,-2.42120}},\
{108,{-0.81763,-2.11673}},
{109,{-1.10044,-1.14419}},{110,{-0.70165,-2.49974}},{111,{-0.42866,-2.24750}},\
{112,{-2.43731,1.79165}},
{113,{-2.74979,1.47774}},{114,{-2.79547,1.05913}},{115,{-1.10402,-0.62002}},{\
116,{-1.45658,-1.57157}},
{117,{-0.77934,0.08361}},{118,{-1.17329,-0.15124}},{119,{-0.82138,-0.32936}},{\
120,{2.20794,-0.05916}},
{121,{2.45321,-0.04498}},{122,{2.39429,-0.29238}},{123,{2.12615,-0.50006}},{\
124,{0.63583,2.58692}},
{125,{0.32174,2.29803}},{126,{-0.02877,2.02782}},{127,{-0.39363,1.77903}},{\
128,{-1.61114,2.41332}},
{129,{-1.19311,2.20807}},{130,{-0.79142,1.99259}},{131,{0.24331,2.69927}},{\
132,{-0.41750,2.18915}},
{133,{-0.08581,2.42559}},{134,{-1.15155,2.67024}},{135,{-0.77530,2.39109}},{\
136,{-0.66814,2.76770}},
{137,{-0.43006,2.51880}},{138,{-0.20358,2.82047}},{139,{0.01857,1.60135}},{\
140,{0.41995,1.44126}},
{141,{0.79972,1.25840}},{142,{1.18778,1.04865}},{143,{1.60187,0.86328}},{144,{\
1.98306,0.70495}},
{145,{2.31237,0.59123}},{146,{2.06177,0.30145}},{147,{2.38563,0.23824}},{148,{\
0.45170,1.05916}},
{149,{0.20798,0.64745}},{150,{0.07437,1.15569}},{151,{-0.34430,1.33047}},{152,\
{-0.30760,0.85443}},
{153,{1.65454,0.39067}},{154,{1.19351,0.54637}},{155,{0.75686,0.79291}},{156,{\
1.82775,-0.05735}},
{157,{0.38419,1.87371}},{158,{0.78193,2.17477}},{159,{0.79203,1.69898}},{160,{\
1.20304,2.31191}},
{161,{1.17559,1.93653}},{162,{1.19129,1.50911}},{163,{0.96156,2.51714}},{164,{\
2.26028,0.91570}},
{165,{1.58493,1.32155}},{166,{2.00892,1.16718}},{167,{1.47757,2.13378}},{168,{\
1.60819,1.79046}},
{169,{1.91204,1.55049}}};
elemnodes={{1,{1,115,4}},{2,{10,109,1}},{3,{115,1,109}},{4,{2,123,11}},
{5,{2,22,78}},{6,{78,79,2}},{7,{79,80,2}},{8,{2,80,123}},
{9,{3,151,23}},{10,{3,32,152}},{11,{151,3,152}},{12,{4,30,5}},
{13,{4,118,30}},{14,{115,118,4}},{15,{5,90,6}},{16,{5,30,90}},
{17,{6,98,7}},{18,{90,95,6}},{19,{95,98,6}},{20,{7,91,8}},
{21,{91,7,100}},{22,{98,100,7}},{23,{8,93,9}},{24,{8,91,93}},
{25,{9,116,10}},{26,{9,93,85}},{27,{9,85,116}},{28,{109,10,116}},
{29,{11,156,12}},{30,{156,123,120}},{31,{123,156,11}},{32,{12,154,13}},
{33,{153,154,12}},{34,{153,12,156}},{35,{13,149,14}},{36,{149,155,148}},
{37,{149,13,155}},{38,{154,155,13}},{39,{14,86,15}},{40,{14,149,86}},
{41,{15,119,16}},{42,{86,117,15}},{43,{15,117,119}},{44,{17,16,109}},
{45,{16,119,115}},{46,{18,17,103}},{47,{109,116,103}},{48,{115,109,16}},
{49,{18,111,19}},{50,{18,108,111}},{51,{103,108,18}},{52,{19,73,20}},
{53,{72,73,19}},{54,{72,19,71}},{55,{19,111,104}},{56,{20,75,21}},
{57,{73,74,20}},{58,{74,75,20}},{59,{21,76,22}},{60,{75,76,21}},
{61,{76,77,22}},{62,{77,78,22}},{63,{23,130,24}},{64,{127,130,23}},
{65,{127,23,151}},{66,{24,129,25}},{67,{129,24,130}},{68,{25,83,26}},
{69,{25,128,83}},{70,{128,25,129}},{71,{26,112,27}},{72,{26,83,112}},
{73,{27,114,28}},{74,{112,113,27}},{75,{113,114,27}},{76,{28,88,29}},
{77,{87,88,28}},{78,{87,28,114}},{79,{29,90,30}},{80,{29,88,89}},
{81,{29,89,90}},{82,{30,118,31}},{83,{31,117,32}},{84,{117,31,118}},
{85,{86,32,117}},{86,{86,152,32}},{87,{33,34,147}},{88,{33,121,82}},
{89,{33,147,121}},{90,{34,35,145}},{91,{145,147,34}},{92,{35,36,164}},
{93,{35,164,145}},{94,{36,166,164}},{95,{37,38,169}},{96,{36,37,166}},
{97,{37,169,166}},{98,{38,39,168}},{99,{38,168,169}},{100,{39,40,167}},
{101,{39,167,168}},{102,{40,41,160}},{103,{40,160,167}},{104,{41,42,163}},
{105,{41,163,160}},{106,{42,43,124}},{107,{42,124,163}},{108,{43,44,131}},
{109,{124,43,131}},{110,{44,45,138}},{111,{44,138,131}},{112,{45,46,138}},
{113,{46,47,136}},{114,{136,138,46}},{115,{136,137,138}},{116,{47,48,134}},
{117,{47,134,136}},{118,{49,134,48}},{119,{49,128,134}},{120,{49,50,128}},
{121,{50,83,128}},{122,{51,52,83}},{123,{51,83,50}},{124,{52,53,112}},
{125,{52,112,83}},{126,{53,54,113}},{127,{112,53,113}},{128,{54,55,113}},
{129,{55,56,114}},{130,{113,55,114}},{131,{56,57,87}},{132,{87,114,56}},
{133,{57,58,94}},{134,{87,57,94}},{135,{58,59,97}},{136,{58,97,94}},
{137,{59,60,97}},{138,{60,61,101}},{139,{97,101,99}},{140,{97,60,101}},
{141,{61,62,100}},{142,{100,101,61}},{143,{62,63,91}},{144,{91,100,62}},
{145,{63,64,93}},{146,{63,93,91}},{147,{64,65,92}},{148,{64,92,93}},
{149,{65,66,84}},{150,{65,84,92}},{151,{66,67,105}},{152,{84,66,105}},
{153,{67,107,105}},{154,{68,69,107}},{155,{67,68,107}},{156,{69,70,110}},
{157,{69,110,107}},{158,{70,71,104}},{159,{70,104,110}},{160,{104,71,19}},
{161,{81,123,80}},{162,{81,122,123}},{163,{81,82,122}},{164,{82,121,122}},
{165,{84,85,92}},{166,{84,102,85}},{167,{84,105,102}},{168,{85,93,92}},
{169,{85,102,116}},{170,{86,149,152}},{171,{87,94,88}},{172,{88,96,89}},
{173,{96,88,94}},{174,{90,89,95}},{175,{95,89,96}},{176,{94,97,96}},
{177,{95,96,99}},{178,{95,99,98}},{179,{96,97,99}},{180,{98,99,101}},
{181,{98,101,100}},{182,{106,116,102}},{183,{106,103,116}},{184,{105,106,102}}\
,
{185,{17,109,103}},{186,{104,111,110}},{187,{105,107,106}},{188,{106,107,108}}\
,
{189,{106,108,103}},{190,{108,107,110}},{191,{108,110,111}},{192,{115,119,118}\
},
{193,{117,118,119}},{194,{120,122,121}},{195,{146,120,147}},{196,{120,123,122}\
},
{197,{146,156,120}},{198,{120,121,147}},{199,{124,131,125}},{200,{124,125,158}\
},
{201,{124,158,163}},{202,{133,126,125}},{203,{126,157,125}},{204,{133,125,131}\
},
{205,{125,157,158}},{206,{126,132,127}},{207,{139,126,127}},{208,{126,133,132}\
},
{209,{126,139,157}},{210,{127,132,130}},{211,{139,127,151}},{212,{128,129,134}\
},
{213,{135,129,130}},{214,{129,135,134}},{215,{135,130,132}},{216,{131,138,133}\
},
{217,{132,133,137}},{218,{132,137,135}},{219,{133,138,137}},{220,{135,136,134}\
},
{221,{135,137,136}},{222,{150,140,139}},{223,{157,139,140}},{224,{139,151,150}\
},
{225,{140,148,141}},{226,{159,140,141}},{227,{140,150,148}},{228,{159,157,140}\
},
{229,{141,155,142}},{230,{162,141,142}},{231,{141,148,155}},{232,{162,159,141}\
},
{233,{143,142,154}},{234,{143,165,142}},{235,{142,155,154}},{236,{142,165,162}\
},
{237,{143,153,144}},{238,{143,144,166}},{239,{143,154,153}},{240,{143,166,165}\
},
{241,{144,146,145}},{242,{144,145,164}},{243,{144,153,146}},{244,{144,164,166}\
},
{245,{145,146,147}},{246,{146,153,156}},{247,{148,150,149}},{248,{150,152,149}\
},
{249,{152,150,151}},{250,{157,159,158}},{251,{159,161,158}},{252,{160,158,161}\
},
{253,{158,160,163}},{254,{161,159,162}},{255,{160,161,167}},{256,{161,162,168}\
},
{257,{161,168,167}},{258,{162,165,168}},{259,{165,166,169}},{260,{165,169,168}\
}};
xy=nodecoor; ijk=elemnodes; (*Print[xy]; Print[ijk];*)
numnod=Length[xy]; numele=Length[ijk];
NodeCoordinates=Table[xy[[i,2]],{i,1,numnod}]; 
ElemNodes=Table[ijk[[i,2]],{i,1,numele}];
(*Print[NodeCoordinates]; Print[ElemNodes];*)

ClearAll[Em,\[Nu],th,x0,y0,aspect,Nsub]; 
Em=1000; \[Nu]=0; th=1; aspect=1; Nsub=4;
Emat=Em/(1-\[Nu]^2)*{{1,\[Nu],0},{\[Nu],1,0},{0,0,(1-\[Nu])/2}};
ElemMaterials=   Table[Emat,{numele}]; 
ElemFabrications=Table[th,  {numele}];
ElemTypes=Table[\"Trig3\",{numele}]; 
NodeDOFValues=NodeDOFTags=Table[{0,0},{numnod}];
NodeDOFTags[[58]]={1,1};  NodeDOFTags[[33]]={0,1};
x0=-0.348906;  y0=0.132318;
For [n=33,n<=82,n++, {xn,yn}=NodeCoordinates[[n]];
     LL=(xn-x0)^2+(yn-y0)^2; 
     NodeDOFValues[[n]]=Chop[3*{-(xn-x0)/LL,-(yn-y0)/LL},0.000001]]; 
ProcessOptions={True};
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    \"Circle Mesh\",True,True];
    
(*  Solve problem and print results *)

{NodeDisplacements,NodeForces,NodePlateCounts,NodePlateStresses,
    ElemBarNumbers,ElemBarForces}= PlaneStressSolution[
    NodeCoordinates,ElemTypes,ElemNodes,
    ElemMaterials,ElemFabrications,
    NodeDOFTags,NodeDOFValues,ProcessOptions];
    
PrintPlaneStressSolution[NodeDisplacements,NodeForces,NodePlateCounts,
      NodePlateStresses,\"Computed Solution:\",{}];

(*  Plot Displacement Components Distribution - skipped *)

(* ux=Table[NodeDisplacements[[n,1]],{n,numnod}];
uy=Table[NodeDisplacements[[n,2]],{n,numnod}];
{uxmax,uymax}=Abs[{Max[ux],Max[uy]}];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,ux,
       uxmax,Nsub,aspect,\"Displacement component ux\"];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,uy,
       uymax,Nsub,aspect,\"Displacement component uy\"]; *)

(*  Plot Principal and VM Nodal Stresses Distribution *)

sVM=sig1=sig2=tmax=Table[0,{numnod}]; 
For[n=1,n<=numnod,n++,
   {sigxx,sigyy,sigxy}=NodePlateStresses[[n]]; s0=(sigxx+sigyy)/2; 
   R=N[Sqrt[(sigxx-sigyy)^2+sigxy^2]];
   sig1[[n]]=s0+R; sig2[[n]]=s0-R; tmax[[n]]=R;
   sVM[[n]]=N[Sqrt[3*(sigxx^2+sigyy^2+sigxy^2)/2]]]; 
Print[\"Max sig1,sig2,VM  =\",{sg1max,sg2max,sVMmax}];

aspect=1;
sigeps=10.^(-3); nbands=40;
sxx=Table[NodePlateStresses[[n,1]],{n,numnod}]; 
syy=Table[NodePlateStresses[[n,2]],{n,numnod}];
sxy=Table[NodePlateStresses[[n,3]],{n,numnod}];
{sxxmax,syymax,sxymax}=Abs[{Max[sxx],Max[syy],Max[sxy]}]+sigeps; 
{sxxmin,syymin,sxymin}=Abs[{Min[sxx],Min[syy],Min[sxy]}]+sigeps;
sxxmax=Max[sxxmax,sxxmin];  sxxmin=-sxxmax;
syymax=Max[syymax,syymin];  syymin=-syymax;
sxymax=Max[sxymax,sxymin];  sxymin=-sxymax;
{sxxinc,syyinc,sxyinc}={sxxmax-sxxmin,syymax-syymin,sxymax-sxymin}/nbands; 
Print[\"sxxmin,sxxmax,sxxinc=\",{sxxmin,sxxmax,sxxinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxx,{sxxmin,\
sxxmax,sxxinc},
  {True,True,True,False,True,True},{0.2,-0.5},aspect,\"Stress sigma-xx\"];
Print[\"syymin,syymax,syyinc=\",{syymin,syymax,syyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,syy,{syymin,\
syymax,syyinc},
  {True,True,True,False,True,True},{0.2,-0.5},aspect,\"Stress sigma-yy\"];
Print[\"sxymin,sxymax,sxyinc=\",{sxymin,sxymax,sxyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxy,{sxymin,\
sxymax,sxyinc},
  {True,True,True,False,True,True},{0.2,-0.5},aspect,\"Stress sigma-xy\"];


(*  Plot Principal and VM Nodal Stresses Distribution *)
  
sVM=sig1=sig2=Table[0,{numnod}]; 
For[n=1,n<=numnod,n++,
   {sigxx,sigyy,sigxy}=NodePlateStresses[[n]]; s0=(sigxx+sigyy)/2; 
   R=N[Sqrt[(sigxx-sigyy)^2+sigxy^2]];
   sig1[[n]]=s0+R; sig2[[n]]=s0-R; tmax[[n]]=R;
   sVM[[n]]=N[Sqrt[3*(sigxx^2+sigyy^2+sigxy^2)/2]]];

{sg1max,sg2max,sVMmax}=Abs[{Max[sig1],Max[sig2],Max[sVM]}]+sigeps; 
{sg1min,sg2min,sVMmin}=Abs[{Min[sig1],Min[sig2],Min[sVM]}]+sigeps;
sg1max=Max[sg1max,sg1min];  sg1min=-sg1max;
sg2max=Max[sg2max,sg2min];  sg2min=-sg2max;
sVMmax=Max[sVMmax,sVMmin];  sVMmin=-sVMmax;
{sg1inc,sg2inc,sVMinc}={sg1max-sg1min,sg2max-sg2min,sVMmax-sVMmin}/nbands; 
Print[\"sg1min,sg1max,sg1inc=\",{sg1min,sg1max,sg1inc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sig1,{sg1min,\
sg1max,sg1inc},
  {True,True,True,False,True,True},{0.2,-0.5},aspect,\"Max principal stress \
sig1\"];
Print[\"sg2min,sg2max,sg2inc=\",{sg2min,sg2max,sg2inc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sig2,{sg2min,\
sg2max,sg2inc},
  {True,True,True,False,True,True},{0.2,-0.5},aspect,\"Min principal stress \
sig2\"];
Print[\"sVMmin,sVMmax,sVMinc=\",{sVMmin,sVMmax,sVMinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sVM,{sVMmin,\
sVMmax,sVMinc},
  {True,True,True,False,True,True},{0.2,-0.5},aspect,\"Von Mises stress\"];\
\
\>", "Input",
  CellFrame->True,
  CellMargins->{{11, 35}, {Inherited, Inherited}},
  CellLabelMargins->{{2, 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->{998, 1030},
WindowMargins->{{371, Automatic}, {Automatic, 59}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
Magnification->1,
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, 1291, 27, 258, "Text"],
Cell[3048, 80, 300, 9, 48, "Text"],
Cell[3351, 91, 374, 11, 88, "Input",
  InitializationCell->True],
Cell[3728, 104, 338, 7, 43, "Input",
  InitializationCell->True],
Cell[4069, 113, 778, 20, 158, "Text"],
Cell[4850, 135, 2636, 62, 853, "Input",
  InitializationCell->True],
Cell[7489, 199, 279, 7, 62, "Text"],
Cell[7771, 208, 3475, 59, 928, "Input",
  InitializationCell->True],
Cell[11249, 269, 280, 7, 62, "Text"],
Cell[11532, 278, 4214, 66, 1093, "Input",
  InitializationCell->True],
Cell[15749, 346, 370, 9, 62, "Text",
  InitializationCell->True],
Cell[16122, 357, 3983, 63, 1048, "Input",
  InitializationCell->True],
Cell[20108, 422, 304, 8, 46, "Text"],
Cell[20415, 432, 6038, 95, 1528, "Input"],
Cell[26456, 529, 316, 8, 46, "Text"],
Cell[26775, 539, 7215, 114, 1543, "Input"],
Cell[33993, 655, 376, 7, 46, "Text",
  InitializationCell->True],
Cell[34372, 664, 5501, 94, 1168, "Input"],
Cell[39876, 760, 241, 7, 46, "Text"],
Cell[40120, 769, 2330, 40, 643, "Input",
  InitializationCell->True],
Cell[42453, 811, 282, 7, 62, "Text"],
Cell[42738, 820, 4167, 70, 1048, "Input",
  InitializationCell->True],
Cell[46908, 892, 304, 8, 62, "Text"],
Cell[47215, 902, 7256, 115, 1528, "Input"],
Cell[54474, 1019, 289, 8, 46, "Text",
  InitializationCell->True],
Cell[54766, 1029, 1629, 29, 463, "Input",
  InitializationCell->True],
Cell[56398, 1060, 342, 9, 46, "Text",
  InitializationCell->True],
Cell[56743, 1071, 3246, 54, 793, "Input"],
Cell[59992, 1127, 254, 7, 46, "Text"],
Cell[60249, 1136, 5075, 94, 1183, "Input",
  InitializationCell->True],
Cell[65327, 1232, 335, 9, 62, "Text"],
Cell[65665, 1243, 2298, 38, 748, "Input",
  InitializationCell->True],
Cell[67966, 1283, 240, 7, 46, "Text"],
Cell[68209, 1292, 5276, 85, 1303, "Input",
  InitializationCell->True],
Cell[73488, 1379, 282, 7, 46, "Text"],
Cell[73773, 1388, 4275, 77, 1228, "Input",
  InitializationCell->True],
Cell[78051, 1467, 288, 8, 62, "Text"],
Cell[78342, 1477, 16722, 283, 4633, "Input",
  InitializationCell->True],
Cell[95067, 1762, 239, 5, 46, "Text"],
Cell[95309, 1769, 5576, 122, 1723, "Input",
  InitializationCell->True],
Cell[100888, 1893, 255, 7, 46, "Text"],
Cell[101146, 1902, 744, 17, 208, "Input",
  InitializationCell->True],
Cell[101893, 1921, 229, 7, 46, "Text"],
Cell[102125, 1930, 1562, 31, 358, "Input",
  InitializationCell->True],
Cell[103690, 1963, 244, 7, 46, "Text"],
Cell[103937, 1972, 2944, 67, 883, "Input"],
Cell[106884, 2041, 244, 7, 46, "Text"],
Cell[107131, 2050, 3523, 79, 1063, "Input"],
Cell[110657, 2131, 224, 7, 46, "Text"],
Cell[110884, 2140, 4366, 94, 1288, "Input"],
Cell[115253, 2236, 201, 4, 46, "Text"],
Cell[115457, 2242, 3554, 59, 1108, "Input"],
Cell[119014, 2303, 196, 4, 46, "Text"],
Cell[119213, 2309, 13979, 285, 3148, "Input"]
}
]
*)



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