(***********************************************************************

                    Mathematica-Compatible Notebook

This notebook can be used on any computer system with Mathematica 4.0,
MathReader 4.0, or any compatible application. 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[    136276,       2466]*)
(*NotebookOutlinePosition[    137361,       2501]*)
(*  CellTagsIndexPosition[    137317,       2497]*)
(*WindowFrame->Normal*)



Notebook[{

Cell[CellGroupData[{
Cell[TextData[{
  StyleBox["(*\nC=DECK SM3SHLOCALSYTEM\nC=PURPOSE Define local system of \
3-node triangle in 3D space\nC=AUTHOR C. A. Felippa\nC=VERSION September 1996\
\nC=EQUIPMENT Machine independent\nC=KEYWORDS finite element shell triangle \
local system direction cosines\nC=BLOCK ABSTRACT\nC\nC     SM3SHLOCALSYS \
computes the local corner coordinates\nC     of a flat triangular element in \
3D space and the direction\nC     cosines of the local system\nC\nC=END \
ABSTRACT\nC=BLOCK USAGE\nC\nC     The calling sequence is\nC\nC       call  \
SM3SHLOCALSYS (xg,yg,zg, xl,yl,zl, dcm, area,status)\nC\nC     where the \
input arguments are:\nC\nC       XG,YG,ZG  Corner coordinates of triangle in \
global system\nC\nC     The outputs are:\nC\nC       XL,YL,ZL   Corner \
coordinates of triangle in local system \nC       DCM        Matrix of \
direction cosines of local system\nC       AREA       Signed area of triangle\
\nC       STATUS     Blank if no error detected, else error message\nC\nC     \
The local system is defined as follows:\nC        x' is directed parallel to \
the 2-1 side\nC        z' is the external normal (counterclockwise).\nC       \
 y' computed as z' x x'\nC\nC=END USAGE\nC=BLOCK FORTRAN\n      subroutine  \
SM3SHLOCALSYS \n     &            (xg,yg,zg, xl,yl,zl, dcm, area, status)\nC\n\
C                   A R G U M E N T S\nC\n      character*( * )     status\n  \
    double precision  xg(3), yg(3), zg(3)\n      double precision  xl(3), \
yl(3), zl(3)\n      double precision  dcm(3,3), area\nC\nC                   \
L O C A L   V A R I A B L E S",
    FontColor->RGBColor[0, 0, 1]],
  "\n",
  StyleBox["C\n      double precision  x21, y21, z21, x32, y32, z32\n      \
double precision  xlr, ylr, zlr, x0, y0, z0\n      double precision  dx(3), \
dy(3), dz(3), area\n      integer           i\nC\nC                   L O G I \
C\nC\n      status = ' '\n      x21 =   xg(2) - xg(1)\n      y21 =   yg(2) - \
yg(1)\n      z21 =   zg(2) - zg(1)\n      x32 =   xg(3) - xg(2)\n      y32 =  \
 yg(3) - yg(2)\n      z32 =   zg(3) - zg(2)\n      xlr =   sqrt( \
x21**2+y21**2+z21**2 )\n      if (xlr .eq. 0.0)               then\n         \
status =  'SM3SHLOCALSYTEM: nodes 1-2 coincide'\n         return\n      end \
if\n      dx(1) =  x21/xlr\n      dx(2) =  y21/xlr\n      dx(3) =  z21/xlr\n  \
    dz(1) =  y21*z32 - z21*y32\n      dz(2) =  z21*x32 - x21*z32\n      dz(3) \
=  x21*y32 - y21*x32\n      zlr  =   sqrt( dz(1)**2 + dz(2)**2 + dz(3)**2 )\n \
     if (zlr .eq. 0.0)               then\n         status =  \
'SM3SHLOCALSYTEM: nodes 1-2-3 are colinear'\n         return\n      end if\n  \
    dz(1) =  dz(1)/zlr\n      dz(2) =  dz(2)/zlr\n      dz(3) =  dz(3)/zlr\n  \
    dy(1) =  dz(2) * dx(3) - dz(3) * dx(2)\n      dy(2) =  dz(3) * dx(1) - \
dz(1) * dx(3)\n      dy(3) =  dz(1) * dx(2) - dz(2) * dx(1)\n      ylr  =   \
sqrt( dy(1)**2 + dy(2)**2 + dy(3)**2 )\n      dy(1) =  dy(1)/ylr\n      dy(2) \
=  dy(2)/ylr\n      dy(3) =  dy(3)/ylr\n      x0 =     (xg(1) + xg(2) + \
xg(3))/3.0d0\n      y0 =     (yg(1) + yg(2) + yg(3))/3.0d0\n      z0 =     \
(zg(1) + zg(2) + zg(3))/3.0d0 \n      do 2000  i = 1,3\n        xl(i) =   \
dx(1)*(xg(i)-x0) + dx(2)*(yg(i)-y0) + dx(3)*(zg(i)-z0)\n        yl(i) =   \
dy(1)*(xg(i)-x0) + dy(2)*(yg(i)-y0) + dy(3)*(zg(i)-z0)\n        zl(i) =   \
dz(1)*(xg(i)-x0) + dz(2)*(yg(i)-y0) + dz(3)*(zg(i)-z0)\n        dcm(1,i) =  \
dx(i)\n        dcm(2,i) =  dy(i)\n        dcm(3,i) =  dz(i)\n 2000   continue\
\n      area =  0.5d0*( (xl(2)-xl(1))*(yl(3)-yl(1)) \n     &              - \
(xl(3)-xl(1))*(yl(1)-yl(2)) )\n      return\n      end\nC=END FORTRAN",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)",
  StyleBox["\n    ",
    FontColor->RGBColor[0, 0, 1]],
  "  SM3SHLOCALSYS[xg_,yg_,zg_]:=Module[{\n      \
x21,y21,z21,x32,y32,z32,xlr,ylr,zlr,x0,y0,z0, area,\n      \
dx=dy=dz=xl=yl=zl={0,0,0}, dcm=Table[0,{3},{3}]},\n      x21=xg[[2]]-xg[[1]]; \
y21=yg[[2]]-yg[[1]];\n      z21=zg[[2]]-zg[[1]]; x32=xg[[3]]-xg[[2]];\n      \
y32=yg[[3]]-yg[[2]]; z32=zg[[3]]-zg[[2]];\n      xlr=Sqrt[ x21^2+y21^2+z21^2 \
];\n      If [xlr<=0, Print[\"Nodes 1-2 coincide\"]; Return[Null]];\n      \
dx[[1]]=x21/xlr; dx[[2]]=y21/xlr; dx[[3]]=z21/xlr;\n      \
dz[[1]]=y21*z32-z21*y32; dz[[2]]=z21*x32-x21*z32;\n      \
dz[[3]]=x21*y32-y21*x32;\n      zlr = Sqrt[ dz[[1]]^2+dz[[2]]^2+dz[[3]]^2 ];\n\
      If [zlr<=0, Print[\"Nodes 1-2-3 are colinear\"]; Return[Null]];\n      \
dz[[1]]=dz[[1]]/zlr; dz[[2]]=dz[[2]]/zlr; dz[[3]]=dz[[3]]/zlr;\n      \
dy[[1]]=dz[[2]]*dx[[3]]-dz[[3]]*dx[[2]];\n      \
dy[[2]]=dz[[3]]*dx[[1]]-dz[[1]]*dx[[3]];\n      \
dy[[3]]=dz[[1]]*dx[[2]]-dz[[2]]*dx[[1]];\n      \
ylr=Sqrt[dy[[1]]^2+dy[[2]]^2+dy[[3]]^2];\n      dy[[1]]=dy[[1]]/ylr; \
dy[[2]]=dy[[2]]/ylr; dy[[3]]=dy[[3]]/ylr;\n      \
x0=(xg[[1]]+xg[[2]]+xg[[3]])/3;\n      y0=(yg[[1]]+yg[[2]]+yg[[3]])/3;\n      \
z0=(zg[[1]]+zg[[2]]+zg[[3]])/3;\n      For [i=1,i<=3,i++,\n        \
xl[[i]]=dx[[1]]*(xg[[i]]-x0)+dx[[2]]*(yg[[i]]-y0)+\n                \
dx[[3]]*(zg[[i]]-z0);\n        \
yl[[i]]=dy[[1]]*(xg[[i]]-x0)+dy[[2]]*(yg[[i]]-y0)+\n                \
dy[[3]]*(zg[[i]]-z0);\n        \
zl[[i]]=dz[[1]]*(xg[[i]]-x0)+dz[[2]]*(yg[[i]]-y0)+\n                \
dz[[3]]*(zg[[i]]-z0);\n        dcm[[1,i]]=dx[[i]]; dcm[[2,i]]=dy[[i]]; \
dcm[[3,i]]=dz[[i]];\n      ];\n      \
area=(1/2)*((xl[[2]]-xl[[1]])*(yl[[3]]-yl[[1]]) -\n                  \
(xl[[3]]-xl[[1]])*(yl[[1]]-yl[[2]]) );\n      Return[{xl,yl,zl,dcm,area}]];\n \
",
  StyleBox["     \n \
{xl,yl,zl,dcm,area}=SM3SHLOCALSYS[{0,1,1},{0,0,1},{0,0,1}];\n Print[dcm];  \
Print[N[Transpose[dcm].dcm]];\n Print[N[dcm.Transpose[dcm]]];\n \
Print[Eigenvalues[N[dcm]]]; ",
    FontColor->RGBColor[1, 0, 0]],
  "  \n      "
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({{1, 0, 0}, {0, 1\/\@2, 1\/\@2}, {0, \(-\(1\/\@2\)\), 
        1\/\@2}}\)], "Print"],

Cell[BoxData[
    \({{1.`, 0.`, 0.`}, {0.`, 1.`, 0.`}, {0.`, 0.`, 1.`}}\)], "Print"],

Cell[BoxData[
    \({{1.`, 0.`, 0.`}, {0.`, 1.`, 0.`}, {0.`, 0.`, 1.`}}\)], "Print"],

Cell[BoxData[
    \({\(\(0.7071067811865475`\)\(\[InvisibleSpace]\)\) + 
        0.7071067811865475`\ \[ImaginaryI], \(\(0.7071067811865475`\)\(\
\[InvisibleSpace]\)\) - 0.7071067811865475`\ \[ImaginaryI], 1.`}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  "(*\n",
  StyleBox["C=DECK SM3SHISODB SM3SHISODB FORTRAN\nC=PURPOSE Form isotropic \
bending constitutive matrix DB\nC=AUTHOR C. A. Felippa\nC=BLOCK FORTRAN\n     \
 subroutine  SM3SHISODB (em, nu, h, db)\n      double precision       em, nu, \
h, db(3,3)\n      double precision       c\n      c =          \
em*h**3/(12.d0*(1.-nu**2))\n      db(1,1) =    c\n      db(1,2) =    nu*c\n   \
   db(2,1) =    nu*c\n      db(2,2) =    c\n      db(3,3) =    \
0.5d0*(1.-nu)*c\n      db(1,3) =    0.0\n      db(2,3) =    0.0\n      \
db(3,1) =    0.0\n      db(3,2) =    0.0\n      return\n      end\nC=END \
FORTRAN",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n\nSM3SHISODB[em_,nu_,h_]:=Module[{c,db},\n      \
c=em*h^3/(1-nu^2)/12;\n      db={{c,nu*c,0},{nu*c,c,0},{0,0,c*(1-nu)/2}};\n   \
   Return[db]];\n      \n",
  StyleBox["Print[SM3SHISODB[100,1/4,2]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({{640\/9, 160\/9, 0}, {160\/9, 640\/9, 0}, {0, 0, 80\/3}}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  "(*",
  StyleBox["\nC=DECK SM3SHISODM SM3SHISODM FORTRAN\nC=PURPOSE Form isotropic \
membrane constitutive matrix DM\nC=AUTHOR C. A. Felippa\nC=BLOCK FORTRAN\n    \
  subroutine  SM3SHISODM (em, nu, h, dm)\n      double precision       em, \
nu, h, dm(3,3)\n      double precision       c\n      c =          \
em*h/(1.-nu**2)\n      dm(1,1) =    c\n      dm(1,2) =    nu*c\n      dm(2,1) \
=    nu*c\n      dm(2,2) =    c\n      dm(3,3) =    0.5d0*(1.-nu)*c\n      \
dm(1,3) =    0.0\n      dm(2,3) =    0.0\n      dm(3,1) =    0.0\n      \
dm(3,2) =    0.0\n      return\n      end\nC=END FORTRAN",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\nSM3SHISODM[em_,nu_,h_]:=Module[{c,dm},\n      c=em*h/(1-nu^2);\n     \
 dm={{c,nu*c,0},{nu*c,c,0},{0,0,c*(1-nu)/2}};\n      Return[dm]];\n      \n",
  StyleBox["Print[SM3SHISODM[100,1/4,2]]; ",
    FontColor->RGBColor[1, 0, 0]],
  "     "
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({{640\/3, 160\/3, 0}, {160\/3, 640\/3, 0}, {0, 0, 80}}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  "(*",
  StyleBox["\nC=DECK SM3SHTRANSFORM\nC=PURPOSE Transform stiffness of 3-node \
shell element to global axes\nC=AUTHOR C. A. Felippa\nC=VERSION September \
1996\nC=EQUIPMENT Machine independent\nC=KEYWORDS finite element shell \
triangle local global transform\nC=BLOCK ABSTRACT\nC\nC    SM3SHTRANSFORM \
transforms the local element stiffness to\nC    global coordinates. \nC\nC    \
The stiffness transformation is assumed to have the block \nC    diagonal \
form\nC\nC    [T1' 0  0  0  0  0 ] [S11 S12 S13 S14 S15 S16] [T1  0  0  0  0  \
0]\nC    [ 0 T2' 0  0  0  0 ] [    S22 S23 S24 S25 S26] [ 0 T2  0  0  0  0]\n\
C    [ 0  0 T3' 0  0  0 ] [        S33 S34 S35 S36] [ 0  0 T3  0  0  0]\nC    \
[ 0  0  0 T4' 0  0 ] [            S44 S45 S46] [ 0  0  0 T4  0  0] \nC    [ 0 \
 0  0  0 T5' 0 ] [                S55 S56] [ 0  0  0  0 T5  0]\nC    [ 0  0  \
0  0  0 T6'] [   symm             S66] [ 0  0  0  0  0 T6]\nC\nC    where Ti \
are direction cosines of local directions wrt global.\nC    Each block is 3 x \
3. Indices 1, 3 and 5 pertains to the three \nC    translations at corners 1, \
2 and 3, respectively. Indices 2, 4\nC    and 6 pertain to the three \
rotations at corners 1, 2 and 3.\nC    \nC    The subroutine accounts for the \
possibility that T1, T2, ... T6 \nC    may be different, although they will \
often be the same.\nC\nC    The transformed block Sij is evidently \
Ti'.Sij.Tj, which can be\nC    computed in 54 multiply-add operations.  For \
all blocks this\nC    amounts to about 1000 operations. The implementation \
below uses\nC    50% more operations: 1539, to streamline loops.  A brute \
force   \nC    version ignoring the block-diagonal form would require 2 \
x18^3/2 \nC    ~ 6000  operations.  Tests on the Sun show that the present\nC \
   implementation is 8 times faster than the brute force, because\nC    of \
the careful use of local variables to reduce indexing overhead. \nC     \n\
C=END ABSTRACT\nC=BLOCK USAGE\nC\nC     The calling sequence is\nC\nC       \
call  SM3SHTRANSFORM (sm, T1,T2,T3,T4,T5,T6)\nC\nC     where the input \
arguments are:\nC\nC      SM               Incoming 18 x 18 matrix in local \
system\nC      T1 through T6    Direction cosine matrices - see abstract\nC   \
  \nC\nC     The outputs are:\nC\nC      SM               Output 18 x 18 \
matrix in global system\nC     \nC\nC=END USAGE\nC=BLOCK FORTRAN\n      \
subroutine  SM3SHTRANSFORM \n     &            (sm, t1, t2, t3, t4, t5, t6)\n\
C\nC                   A R G U M E N T S\nC\n      double precision  t1(3,3), \
t2(3,3), t3(3,3)\n      double precision  t4(3,3), t5(3,3), t6(3,3)\n      \
double precision  sm(18,18)\nC\nC                   L O C A L   V A R I A B L \
E S\nC\n      double precision  t(3,3,6), st1(18), st2(18), st3(18)\n      \
double precision  t11, t12, t13, t21, t22, t23, t31, t32, t33\n      double \
precision  tst11, tst12, tst13, tst21, tst22, tst23  \n      double precision \
 tst31, tst32, tst33 \n      integer           i, j, ib, jb\nC",
    FontColor->RGBColor[0, 0, 1]],
  "\n",
  StyleBox["C                   L O G I C\nC\n      do 1400  j = 1,3\n        \
do 1200  i = 1,3\n          t(i,j,1) = t1(i,j)\n          t(i,j,2) = t2(i,j)\n\
          t(i,j,3) = t3(i,j)\n          t(i,j,4) = t4(i,j)\n          \
t(i,j,5) = t5(i,j)\n          t(i,j,6) = t6(i,j)\n 1200   continue\n 1400 \
continue\n      do 3000  jb = 1,6\n        j =     3*(jb-1)\n        t11 =   \
t(1,1,jb)\n        t21 =   t(2,1,jb)\n        t31 =   t(3,1,jb)\n        t12 \
=   t(1,2,jb)\n        t22 =   t(2,2,jb)\n        t32 =   t(3,2,jb)\n        \
t13 =   t(1,3,jb)\n        t23 =   t(2,3,jb)\n        t33 =   t(3,3,jb)\n     \
   do 2200  i = 1,18\n          st1(i) =  sm(i,j+1)*t11 + sm(i,j+2)*t21 + \
sm(i,j+3)*t31\n          st2(i) =  sm(i,j+1)*t12 + sm(i,j+2)*t22 + \
sm(i,j+3)*t32  \n          st3(i) =  sm(i,j+1)*t13 + sm(i,j+2)*t23 + \
sm(i,j+3)*t33\n 2200   continue\n        do 2800  ib = 1,jb",
    FontColor->RGBColor[0, 0, 1]],
  "\n ",
  StyleBox["         i =     3*(ib-1)\n          t11 =   t(1,1,ib)\n          \
t21 =   t(2,1,ib)\n          t31 =   t(3,1,ib)\n          t12 =   t(1,2,ib)\n \
         t22 =   t(2,2,ib)\n          t32 =   t(3,2,ib)\n          t13 =   \
t(1,3,ib)\n          t23 =   t(2,3,ib)\n          t33 =   t(3,3,ib)\n         \
 tst11 = t11*st1(i+1) + t21*st1(i+2) + t31*st1(i+3)\n          tst21 = \
t12*st1(i+1) + t22*st1(i+2) + t32*st1(i+3)  \n          tst31 = t13*st1(i+1) \
+ t23*st1(i+2) + t33*st1(i+3)\n          tst12 = t11*st2(i+1) + t21*st2(i+2) \
+ t31*st2(i+3)\n          tst22 = t12*st2(i+1) + t22*st2(i+2) + t32*st2(i+3)  \
\n          tst32 = t13*st2(i+1) + t23*st2(i+2) + t33*st2(i+3)\n          \
tst13 = t11*st3(i+1) + t21*st3(i+2) + t31*st3(i+3)\n          tst23 = \
t12*st3(i+1) + t22*st3(i+2) + t32*st3(i+3)  \n          tst33 = t13*st3(i+1) \
+ t23*st3(i+2) + t33*st3(i+3)\n          sm(i+1,j+1) =   tst11\n          \
sm(j+1,i+1) =   tst11\n          sm(i+1,j+2) =   tst12\n          sm(j+2,i+1) \
=   tst12\n          sm(i+1,j+3) =   tst13\n          sm(j+3,i+1) =   tst13\n \
         sm(i+2,j+1) =   tst21\n          sm(j+1,i+2) =   tst21\n          \
sm(i+2,j+2) =   tst22\n          sm(j+2,i+2) =   tst22\n          sm(i+2,j+3) \
=   tst23\n          sm(j+3,i+2) =   tst23\n          sm(i+3,j+1) =   tst31\n \
         sm(j+1,i+3) =   tst31\n          sm(i+3,j+2) =   tst32\n          \
sm(j+2,i+3) =   tst32\n          sm(i+3,j+3) =   tst33\n          sm(j+3,i+3) \
=   tst33\n 2800   continue\n 3000 continue\n      return\n      end\nC=END \
FORTRAN\n",
    FontColor->RGBColor[0, 0, 1]],
  "*)\n      SM3SHTRANSFORM[esm_,t1_,t2_,t3_,t4_,t5_,t6_]:=Module[{\n      \
t=Table[0,{3},{3},{6}], st1=st2=st3=Table[0,{18}],\n      t11, t12, t13, t21, \
t22, t23, t31, t32, t33,\n      tst11, tst12, tst13, tst21, tst22, tst23,\n   \
   tst31, tst32, tst33,  i,j,ib,jb, sm=esm},\n\n      For [i=1,i<=3,i++, \
For[j=1,j<=3,j++,\n          t[[i,j,1]]=t1[[i,j]]; t[[i,j,2]]=t2[[i,j]]; \n   \
       t[[i,j,3]]=t3[[i,j]]; t[[i,j,4]]=t4[[i,j]]; \n          \
t[[i,j,5]]=t5[[i,j]]; t[[i,j,6]]=t6[[i,j]] ]];\n\n      For [jb=1,jb<=6,jb++, \
j=3*(jb-1);\n        t11=t[[1,1,jb]]; t21=t[[2,1,jb]]; t31=t[[3,1,jb]];\n     \
   t12=t[[1,2,jb]]; t22=t[[2,2,jb]]; t32=t[[3,2,jb]];\n        \
t13=t[[1,3,jb]]; t23=t[[2,3,jb]]; t33=t[[3,3,jb]];\n        For \
[i=1,i<=18,i++,\n          \
st1[[i]]=sm[[i,j+1]]*t11+sm[[i,j+2]]*t21+sm[[i,j+3]]*t31;\n          \
st2[[i]]=sm[[i,j+1]]*t12+sm[[i,j+2]]*t22+sm[[i,j+3]]*t32; \n          \
st3[[i]]=sm[[i,j+1]]*t13+sm[[i,j+2]]*t23+sm[[i,j+3]]*t33];\n        For \
[ib=1,ib<=jb,ib++, i=3*(ib-1);\n          t11=t[[1,1,ib]]; t21=t[[2,1,ib]]; \
t31=t[[3,1,ib]];\n          t12=t[[1,2,ib]]; t22=t[[2,2,ib]]; \
t32=t[[3,2,ib]];\n          t13=t[[1,3,ib]]; t23=t[[2,3,ib]]; \
t33=t[[3,3,ib]];\n          \
tst11=t11*st1[[i+1]]+t21*st1[[i+2]]+t31*st1[[i+3]];\n          \
tst21=t12*st1[[i+1]]+t22*st1[[i+2]]+t32*st1[[i+3]];  \n          \
tst31=t13*st1[[i+1]]+t23*st1[[i+2]]+t33*st1[[i+3]];\n          \
tst12=t11*st2[[i+1]]+t21*st2[[i+2]]+t31*st2[[i+3]];\n          \
tst22=t12*st2[[i+1]]+t22*st2[[i+2]]+t32*st2[[i+3]];  \n          \
tst32=t13*st2[[i+1]]+t23*st2[[i+2]]+t33*st2[[i+3]];\n          \
tst13=t11*st3[[i+1]]+t21*st3[[i+2]]+t31*st3[[i+3]];\n          \
tst23=t12*st3[[i+1]]+t22*st3[[i+2]]+t32*st3[[i+3]];  \n          \
tst33=t13*st3[[i+1]]+t23*st3[[i+2]]+t33*st3[[i+3]];\n          \
sm[[i+1,j+1]]=tst11; sm[[j+1,i+1]]=tst11;\n          sm[[i+1,j+2]]=tst12; \
sm[[j+2,i+1]]=tst12;\n          sm[[i+1,j+3]]=tst13; sm[[j+3,i+1]]=tst13;\n   \
       sm[[i+2,j+1]]=tst21; sm[[j+1,i+2]]=tst21;\n          \
sm[[i+2,j+2]]=tst22; sm[[j+2,i+2]]=tst22;\n          sm[[i+2,j+3]]=tst23; \
sm[[j+3,i+2]]=tst23;\n          sm[[i+3,j+1]]=tst31; sm[[j+1,i+3]]=tst31;\n   \
       sm[[i+3,j+2]]=tst32; sm[[j+2,i+3]]=tst32;\n          \
sm[[i+3,j+3]]=tst33; sm[[j+3,i+3]]=tst33];\n          ];\n      Return[sm];\n \
     ];\n      \n",
  StyleBox["      \n  esm=N[Table[1,{18},{18}]];\n  \
Print[Chop[Eigenvalues[N[esm]]]];\n  \
t1=t2=t3=t4=t5=t6={{3/5,4/5,0},{-4/5,3/5,0},{0,0,1}};   \n  sm=  \
SM3SHTRANSFORM[esm,t1,t2,t3,t4,t5,t6];\n  Print[sm];\n  \
Print[Chop[Eigenvalues[N[sm]]]];\n  ",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({18.000000000000004`, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0}\)], "Print"],

Cell[BoxData[
    \({{0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\)}, \
{\(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`}, {\(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`, \(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`}, {0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\)}, \
{\(-0.28`\), 1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`}, {\(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`, \(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 
        1.`}, {0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\)}, \
{\(-0.28`\), 1.9599999999999997`, 1.4`, \(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`}, {\(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`, \(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`}, {0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\)}, \
{\(-0.28`\), 1.9599999999999997`, 1.4`, \(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`}, {\(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`, \(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 
        1.`}, {0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.28`\), \(-0.20000000000000007`\)}, \
{\(-0.28`\), 1.9599999999999997`, 1.4`, \(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.28`\), 1.9599999999999997`, 1.4`, \(-0.28`\), 
        1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`, \(-0.2800000000000001`\), 
        1.9599999999999997`, 1.4`}, {\(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`, \(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`}, {0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\), 
        0.04000000000000002`, \(-0.2800000000000001`\), \
\(-0.20000000000000007`\)}, {\(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.28`\), 1.9599999999999997`, 1.4`, \(-0.28`\), 
        1.9599999999999997`, 1.4`, \(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.28`\), 1.9599999999999997`, 
        1.4`, \(-0.2800000000000001`\), 1.9599999999999997`, 
        1.4`}, {\(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`, \(-0.20000000000000007`\), 1.4`, 
        1.`, \(-0.20000000000000007`\), 1.4`, 1.`, \(-0.20000000000000007`\), 
        1.4`, 1.`}}\)], "Print"],

Cell[BoxData[
    \({17.999999999999986`, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  "(*",
  StyleBox["\nC=DECK SM3SHBENDB SM3SHBENDB FORTRAN\nC=PURPOSE Form SM3SHBENDB \
bending stiffness of 3-node, 9 dof Kirchhoff triangle\nC=AUTHOR C. A. \
Felippa, May 1984\nC=VERSION September 1989\nC=EQUIPMENT Machine independent\n\
C=KEYWORDS Kirchhoff thin plate bending\nC=KEYWORDS finite element triangle \
basic stiffness matrix\nC=BLOCK ABSTRACT\nC\nC     SM3SHBENDB forms the basic \
material element stiffness matrix\nC     of a Kirchhoff thin plate bending \
triangle constructed with the \nC     generalized FF/ANS formulation\nC\n\
C=END ABSTRACT\nC=BLOCK USAGE\nC\nC     The calling sequence is\nC\nC       \
call      SM3SHBENDB (x,y,db, f,clr,cqr, ls,sm,m,status)\nC\nC     where the \
input arguments are:\nC\nC       X         (3 x 1) array of x coordinates of \
triangle nodes\nC       Y         (3 x 1) array of y coordinates of triangle \
nodes\nC       DB        (3 x 3) moment-curvature constitutive matrix\nC      \
 F         Factor by which stiffness entries will be multiplied.\nC       \
CLR,CQR   Use CLR*LLR+CQR*LQR for lumping matrix L\nC                 thus \
CLR+CQR  must add up to 1.\nC                 LLR=linear rotation, \
LLR=quadratic rotation (Kirchhoff)\nC       LS        (9 x 1) array of \
stiffness location pointers\nC                 (see output SM).\nC       SM   \
     Incoming material stiffness array.\nC       M         First dimension of \
SM in calling program.\nC\nC     The outputs are:\nC\nC       SM        \
Output stiffness array with basic stiffness\nC                 coefficients \
added in.  the (I,J)-th entry of the\nC                 (9 x 9) element \
bending stiffness is added to\nC                 SM(K,L), where K=LS(I) and \
L=LS(J).\nC       STATUS    Status character variable.  blank if no error \
detected.\nC\nC=END USAGE\nC=BLOCK FORTRAN\n      subroutine    SM3SHBENDB\n  \
   &             ( x, y, db, f, clr, cqr, ls, sm, m, status)\nC\nC            \
       A R G U M E N T S\nC\n      double precision   x(3), y(3), db(3,3), f, \
sm(m,m)\n      double precision   clr, cqr\n      integer            m, ls(9)\
\n      character*( * )      status\nC\nC                   L O C A L   V A R \
I A B L E S\nC\n      double precision   llr(9,3), lqr(9,3), l(9,3)\n      \
double precision   db11, db12, db13, db22, db23, db33\n      double precision \
  x0, y0, cab, a1, a2, a3, b1, b2, b3\n      double precision   x21, x32, \
x13, y21, y32, y13\n      double precision   x12, x23, x31, y12, y23, y31\n   \
   double precision   xl12, xl23, xl31, c12, c23, c31, s12, s23, s31\n      \
double precision   cc12, cc23, cc31, ss12, ss23, ss31\n      double precision \
  cs12, cs23, cs31, s1, s2, s3\n      double precision   area2, c\n      \
integer            i, j, ii, jj\n",
    FontColor->RGBColor[0, 0, 1]],
  " *)     \n      \n      \
SM3SHBENDB[x_,y_,db_,f_,clr_,cqr_,ls_,esm_]:=Module[{\n      x21,x32,x13, \
y21,y32,y13, x12,x23,x31, y12,y23,y31, x0,y0,\n      cab, a1,a2,a3, b1,b2,b3, \
xl12,xl23,xl31,\n      c12,c23,c31,s12,s23,s31, cc12,cc23,cc31, \
ss12,ss23,ss31,\n      cs12,cs23,cs31, s1,s2,s3, area2,c,\n      \
db11,db12,db13,db21,db22,db23,db31,db32,db33,\n      i, j, ii, jj, \n      L, \
Llr=Lqr=Table[0,{9},{3}], sm=esm },\n(*\n",
  StyleBox["C\nC                   L O G I C\nC\n      status =   ' '\n      \
x21 =      x(2) - x(1)\n      x12 =     -x21\n      x32 =      x(3) - x(2)\n  \
    x23 =     -x32\n      x13 =      x(1) - x(3)\n      x31 =     -x13\n      \
y21 =      y(2) - y(1)\n      y12 =     -y21\n      y32 =      y(3) - y(2)\n  \
    y23 =     -y32\n      y13 =      y(1) - y(3)\n      y31 =     -y13\n      \
area2 =    y21*x13 - x21*y13\n      if (area2 .le. 0.0)      then\n        \
status = 'SM3SHBENDB: Negative area'\n        if (area2 .eq. 0.0)   status = \
'SM3SHBENDB: Zero area'\n        return\n      end if",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n      {x21,x32,x13}={x[[2]]-x[[1]],x[[3]]-x[[2]],x[[1]]-x[[3]]};\n   \
   {x12,x23,x31}=-{x21,x32,x13};\n      \
{y21,y32,y13}={y[[2]]-y[[1]],y[[3]]-y[[2]],y[[1]]-y[[3]]};\n      \
{y12,y23,y31}=-{y21,y32,y13};\n      area2=y21*x13-x21*y13;\n      If \
[area2<=0, Print[\"Negative or zero area\"]; Return[Null]];\n(*\n",
  StyleBox["      x0 =      (x(1)+x(2)+x(3))/3.0d0\n      y0 =      \
(y(1)+y(2)+y(3))/3.0d0\n      cab =      3.0d0/area2\n      a1 =      \
-cab*(y(3)-y0)\n      a2 =      -cab*(y(1)-y0)\n      a3 =      \
-cab*(y(2)-y0)\n      b1 =       cab*(x(3)-x0)\n      b2 =       \
cab*(x(1)-x0)\n      b3 =       cab*(x(2)-x0)\n      xl12 =     \
sqrt(x12**2+y12**2)\n      xl23 =     sqrt(x23**2+y23**2)\n      xl31 =     \
sqrt(x31**2+y31**2)\n      do 1200  j = 1,3\n        do 1200  i = 1,9\n       \
   llr(i,j) =  0.0\n          lqr(i,j) =  0.0\n 1200   continue\n      if \
(clr .ne. 0.0)         then\n        llr(3,1) =   y32*0.5d0\n        llr(6,1) \
=   y13*0.5d0\n        llr(9,1) =   y21*0.5d0\n        llr(2,2) =   x32*0.5d0\
\n        llr(5,2) =   x13*0.5d0\n        llr(8,2) =   x21*0.5d0\n        \
llr(2,3) =  -y32*0.5d0\n        llr(3,3) =  -x32*0.5d0\n        llr(5,3) =  \
-y13*0.5d0\n        llr(6,3) =  -x13*0.5d0\n        llr(8,3) =  -y21*0.5d0\n  \
      llr(9,3) =  -x21*0.5d0\n      end if\n",
    FontColor->RGBColor[0, 0, 1]],
  "*)\n   \n",
  StyleBox["    ",
    FontColor->RGBColor[0, 0, 1]],
  "  x0=(x[[1]]+x[[2]]+x[[3]])/3; y0=(y[[1]]+y[[2]]+y[[3]])/3;\n      \
cab=3/area2;\n      a1=-cab*(y[[3]]-y0); a2=-cab*(y[[1]]-y0); \
a3=-cab*(y[[2]]-y0);\n      b1= cab*(x[[3]]-x0); b2= cab*(x[[1]]-x0); b3= \
cab*(x[[2]]-x0);\n      xl12=Sqrt[x12^2+y12^2];\n      \
xl23=Sqrt[x23^2+y23^2];\n      xl31=Sqrt[x31^2+y31^2];",
  StyleBox["\n     ",
    FontColor->RGBColor[0, 0, 1]],
  " Llr[[3,1]]=y32/2;  Llr[[6,1]]=y13/2;  Llr[[9,1]]=y21/2;\n      \
Llr[[2,2]]=x32/2;  Llr[[5,2]]=x13/2;  Llr[[8,2]]=x21/2;\n      \
Llr[[2,3]]=-y32/2; Llr[[3,3]]=-x32/2; Llr[[5,3]]=-y13/2; \n      \
Llr[[6,3]]=-x13/2; Llr[[8,3]]=-y21/2; Llr[[9,3]]=-x21/2;\n(*\n",
  StyleBox["      if (cqr .ne. 0.0)         then\n        c12 =       \
y21/xl12\n        s12 =       x12/xl12\n        c23 =       y32/xl23\n        \
s23 =       x23/xl23\n        c31 =       y13/xl31\n        s31 =       \
x31/xl31\n        cc12 =      c12*c12\n        cc23 =      c23*c23\n        \
cc31 =      c31*c31\n        ss12 =      s12*s12\n        ss23 =      s23*s23\
\n        ss31 =      s31*s31\n        cs12 =      c12*s12\n        cs23 =    \
  c23*s23\n        cs31 =      c31*s31\n        lqr(1,1) =   cs12 - cs31\n    \
    lqr(1,2) =  -lqr(1,1)\n        lqr(1,3) =  (cc31-ss31) - (cc12-ss12)\n    \
    lqr(2,1) =  (cc12*x12 + cc31*x31)*0.5d0\n        lqr(2,2) =  (ss12*x12 + \
ss31*x31)*0.5d0\n        lqr(2,3) =   ss12*y21 + ss31*y13\n        lqr(3,1) = \
-(cc12*y21 + cc31*y13)*0.5d0\n        lqr(3,2) = -0.5d0*lqr(2,3)\n        \
lqr(3,3) =  -2.d0*lqr(2,1)\n        lqr(4,1) =  cs23 - cs12 \n        \
lqr(4,2) =  -lqr(4,1)\n        lqr(4,3) =  (cc12-ss12) - (cc23-ss23)\n        \
lqr(5,1) =  (cc12*x12 + cc23*x23)*0.5d0\n        lqr(5,2) =  (ss12*x12 + \
ss23*x23)*0.5d0\n        lqr(5,3) =   ss12*y21 + ss23*y32\n        lqr(6,1) = \
-(cc12*y21 + cc23*y32)*0.5d0\n        lqr(6,2) = -0.5d0*lqr(5,3)\n        \
lqr(6,3) =  -2.d0*lqr(5,1)\n        lqr(7,1) =  cs31 - cs23\n        lqr(7,2) \
=  -lqr(7,1)\n        lqr(7,3) =  (cc23-ss23) - (cc31-ss31)\n        lqr(8,1) \
=  (cc23*x23 + cc31*x31)*0.5d0\n        lqr(8,2) =  (ss23*x23 + \
ss31*x31)*0.5d0\n        lqr(8,3) =   ss23*y32 + ss31*y13\n        lqr(9,1) = \
-(cc23*y32 + cc31*y13)*0.5d0\n        lqr(9,2) = -0.5d0*lqr(8,3)\n        \
lqr(9,3) =  -2.d0*lqr(8,1)\n      end if\n*     print '(''llr=''/(9F8.4))', \
llr\n*     print '(''lqr=''/(9F8.4))', lqr",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n    ",
  StyleBox["    ",
    FontColor->RGBColor[0, 0, 1]],
  "c12 =y21/xl12; s12 =x12/xl12; c23 =y32/xl23;\n        s23 =x23/xl23; c31 \
=y13/xl31; s31 =x31/xl31;\n        cc12 =c12*c12; cc23 =c23*c23; cc31 \
=c31*c31;\n        ss12 =s12*s12; ss23 =s23*s23; ss31 =s31*s31;\n        cs12 \
=c12*s12; cs23 =c23*s23; cs31 =c31*s31;\n        Lqr[[1,1]]=cs12-cs31;\n      \
  Lqr[[1,2]]=-Lqr[[1,1]];\n        Lqr[[1,3]]=(cc31-ss31)-(cc12-ss12);\n      \
  Lqr[[2,1]]=(cc12*x12+cc31*x31)/2;\n        \
Lqr[[2,2]]=(ss12*x12+ss31*x31)/2;\n        Lqr[[2,3]]=ss12*y21+ss31*y13;\n    \
    Lqr[[3,1]]=-(cc12*y21+cc31*y13)/2;\n        Lqr[[3,2]]=-(1/2)*Lqr[[2,3]];\
\n        Lqr[[3,3]]=-2*Lqr[[2,1]];\n        Lqr[[4,1]]=cs23-cs12; \n        \
Lqr[[4,2]]=-Lqr[[4,1]];\n        Lqr[[4,3]]=(cc12-ss12)-(cc23-ss23);\n        \
Lqr[[5,1]]=(cc12*x12+cc23*x23)/2;\n        Lqr[[5,2]]=(ss12*x12+ss23*x23)/2;\n\
        Lqr[[5,3]]= ss12*y21+ss23*y32;\n        \
Lqr[[6,1]]=-(cc12*y21+cc23*y32)/2;\n        Lqr[[6,2]]=-(1/2)*Lqr[[5,3]];\n   \
     Lqr[[6,3]]=-2*Lqr[[5,1]];\n        Lqr[[7,1]]= cs31-cs23;\n        \
Lqr[[7,2]]=-Lqr[[7,1]];\n        Lqr[[7,3]]=(cc23-ss23)-(cc31-ss31);\n        \
Lqr[[8,1]]=(cc23*x23+cc31*x31)/2;\n        Lqr[[8,2]]=(ss23*x23+ss31*x31)/2;\n\
        Lqr[[8,3]]=ss23*y32+ss31*y13;\n        \
Lqr[[9,1]]=-(cc23*y32+cc31*y13)/2;\n        Lqr[[9,2]]=-(1/2)*Lqr[[8,3]];\n   \
     Lqr[[9,3]]=-2*Lqr[[8,1]];\n(*\n",
  StyleBox["      do 1600  j = 1,9\n        l(j,1) =   clr*llr(j,1) + \
cqr*lqr(j,1)\n        l(j,2) =   clr*llr(j,2) + cqr*lqr(j,2)\n        l(j,3) \
=   clr*llr(j,3) + cqr*lqr(j,3)\n 1600   continue\nc\n      c =        \
2.d0*f/area2\n      db11 =     c*db(1,1)\n      db22 =     c*db(2,2)\n      \
db33 =     c*db(3,3)\n      db12 =     c*db(1,2)\n      db13 =     c*db(1,3)\n\
      db23 =     c*db(2,3)\n      do 4000  j = 1,9\n        jj =     ls(j)\n  \
      s1 =     db11*l(j,1) + db12*l(j,2) + db13*l(j,3)\n        s2 =     \
db12*l(j,1) + db22*l(j,2) + db23*l(j,3)\n        s3 =     db13*l(j,1) + \
db23*l(j,2) + db33*l(j,3)\n        do 3500  i = 1,j\n          ii =      \
ls(i)\n          sm(jj,ii) =  sm(jj,ii) + (s1*l(i,1)+s2*l(i,2)+s3*l(i,3))\n   \
       sm(ii,jj) =  sm(jj,ii)\n 3500     continue\n 4000   continue\n      \
return\n      end\nC=END FORTRAN",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n      L=clr*Llr+cqr*Lqr;\n      \
{{db11,db12,db13},{db21,db22,db23},{db31,db32,db33}}=\n          \
db*2*f/area2;\n      For [j=1,j<=9,j++, jj=ls[[j]];\n          \
s1=db11*L[[j,1]]+db12*L[[j,2]]+db13*L[[j,3]];\n          \
s2=db12*L[[j,1]]+db22*L[[j,2]]+db23*L[[j,3]];\n          \
s3=db13*L[[j,1]]+db23*L[[j,2]]+db33*L[[j,3]];\n          For [i=1,i<=j,i++, \
ii=ls[[i]];\n              sm[[jj,ii]]+=s1*L[[i,1]]+s2*L[[i,2]]+s3*L[[i,3]];\n\
              sm[[ii,jj]]=sm[[jj,ii]] ];\n          ];\n      Return[sm];\n   \
   ];\n      \n \n",
  StyleBox[" \nClearAll[Em]; Em=1;\n\
sm=SM3SHBENDB[{0,1,1/2},{0,0,1},{{Em,0,0},\n  \
{0,Em,0},{0,0,Em/2}},1,.5,0.5,{1,2,3,4,5,6,7,8,9},\n  Table[0,{9},{9}]];\n\
Print[sm];\nPrint[Chop[Eigenvalues[N[sm]]]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({{0.8000000000000002`, \(-0.09999999999999998`\), 
        0.10000000000000002`, \(-0.4800000000000001`\), 
        0.45999999999999996`, \(-0.22000000000000003`\), \
\(-0.32000000000000006`\), \(-0.03999999999999998`\), \(-0.52`\)}, \
{\(-0.09999999999999998`\), 0.38749999999999996`, 0.08125000000000002`, 
        0.45999999999999996`, 
        0.14249999999999996`, \(-0.02875000000000001`\), \(-0.36`\), \
\(-0.16999999999999998`\), 0.22749999999999998`}, {0.10000000000000002`, 
        0.08125000000000002`, 0.415625`, 0.22000000000000003`, 
        0.02875000000000001`, \(-0.40437500000000004`\), \
\(-0.32000000000000006`\), 0.21000000000000002`, 
        0.04875000000000001`}, {\(-0.4800000000000001`\), 
        0.45999999999999996`, 0.22000000000000003`, 
        0.8000000000000002`, \(-0.09999999999999998`\), \
\(-0.10000000000000002`\), \(-0.32000000000000006`\), \(-0.03999999999999998`\
\), 0.52`}, {0.45999999999999996`, 0.14249999999999996`, 
        0.02875000000000001`, \(-0.09999999999999998`\), 
        0.38749999999999996`, \(-0.08125000000000002`\), \(-0.36`\), \
\(-0.16999999999999998`\), \(-0.22749999999999998`\)}, \
{\(-0.22000000000000003`\), \(-0.02875000000000001`\), \
\(-0.40437500000000004`\), \(-0.10000000000000002`\), \(-0.08125000000000002`\
\), 0.415625`, 0.32000000000000006`, \(-0.21000000000000002`\), 
        0.04875000000000001`}, {\(-0.32000000000000006`\), \(-0.36`\), \
\(-0.32000000000000006`\), \(-0.32000000000000006`\), \(-0.36`\), 
        0.32000000000000006`, 0.6400000000000001`, 0.07999999999999996`, 
        0}, {\(-0.03999999999999998`\), \(-0.16999999999999998`\), 
        0.21000000000000002`, \(-0.03999999999999998`\), \
\(-0.16999999999999998`\), \(-0.21000000000000002`\), 0.07999999999999996`, 
        0.26`, 0}, {\(-0.52`\), 0.22749999999999998`, 0.04875000000000001`, 
        0.52`, \(-0.22749999999999998`\), 0.04875000000000001`, 0, 0, 
        0.42250000000000004`}}\)], "Print"],

Cell[BoxData[
    \({1.9587500000000027`, 1.7136315434029563`, 0.856368456597044`, 0, 0, 0, 
      0, 0, 0}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  StyleBox["(*\nC=DECK SM3SHMEMBB\nC=PURPOSE Form basic membrane stiffness of \
9-dof triangle\nC=AUTHOR C. A. Felippa, June 1984\nC=VERSION June 1984\n\
C=EQUIPMENT Machine independent\nC=KEYWORDS finite element membrane plane \
stress\nC=KEYWORDS basic material stiffness matrix\nC=BLOCK ABSTRACT\nC\nC    \
 SM3SHMEMBB forms the material element stiffness matrix associated with\nC    \
 the basic displacement modes (rigid modes + constant strain\nC     modes) of \
a 9-dof plane-stress triangle based on the\nC     free formulation of Bergan \
and Nygard.\nC\nC=END ABSTRACT\nC=BLOCK USAGE\nC\nC     The calling sequence \
is\nC\nC       CALL      SM3SHMEMBB (X, Y, DM, ALPHA, F, LS, SM, M, STATUS)\n\
C\nC     where the input arguments are\nC\nC       X         (3 x 1) array of \
x coordinates of triangle nodes.\nC       Y         (3 x 1) array of y \
coordinates of triangle nodes.\nC       DM        (3 x 3) matrix relating \
in-plane forces to strains.\nC       ALPHA     Rotational lumping factor; if \
zero form CST.\nC       F         Factor by which stiffness entries will be \
multiplied.\nC       LS        (9 x 1) array of stiffness location pointers, \
see output SM\nC       SM        Incoming material stiffness array.\nC       \
M         First dimension of SM in calling program.\nC\nC     The outputs \
are:\nC\nC       SM        Output stiffness array with basic stiffness\nC     \
            coefficients added in.  The (i,j)-th entry of the\nC              \
   basic element stiffness is added to SM(K,L),\nC                 where \
K=LS(I) and L=LS(J).\nC\nC       STATUS    Status character variable.  Blank \
if no error\nC                 detected.\nC\nC     Note: the internal \
ordering of degrees of freedom is\nC\nC          ux1,uy1, ux2,uy2, ux3,uy3, \
thetaz1,thetaz2,thetaz3\nC\nC     This simplifies the stiffness formation \
loop if alpha=0.\nC     \nC\nC=END USAGE\nC=BLOCK FORTRAN\n      subroutine   \
 SM3SHMEMBB\n     &             (x, y, dm, alpha, f, ls, sm, m, status)\nC\nC \
                  A R G U M E N T S\nC\n      character*( * )      status\n   \
   integer            m, ls(9)\n      double precision   x(3), y(3), dm(3,3), \
alpha, f, sm(m,m)\nC\nC                   L O C A L   V A R I A B L E S\nC\n  \
    double precision   area2, c, p(9,3)\n      double precision   dm11, dm12, \
dm13, dm22, dm23, dm33\n      double precision   x21, x32, x13, y21, y32, y13\
\n      double precision   x12, x23, x31, y12, y23, y31\n      double \
precision   s1, s2, s3\n      integer            i, j, k, l, n",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n      SM3SHMEMBB[x_,y_,dm_,alpha_,f_,ls_,esm_]:=Module[{\n      \
x21,x32,x13, y21,y32,y13, x12,x23,x31, y12,y23,y31,\n      area2, c, p, \
s1,s2,s3,\n      dm11,dm12,dm13,dm21,dm22,dm23,dm31,dm32,dm33,\n      i, j, \
k, l, n, sm=esm },\n(*\n",
  StyleBox["C\nC                   L O G I C\nC\n      status =   ' '\n      \
x21 =      x(2) - x(1)\n      x12 =     -x21\n      x32 =      x(3) - x(2)\n  \
    x23 =     -x32\n      x13 =      x(1) - x(3)\n      x31 =     -x13\n      \
y21 =      y(2) - y(1)\n      y12 =     -y21\n      y32 =      y(3) - y(2)\n  \
    y23 =     -y32\n      y13 =      y(1) - y(3)\n      y31 =     -y13\n      \
area2 =    y21*x13 - x21*y13\n      if (area2 .le. 0.0)      then\n        \
status = 'SM3SHMEMBB: Negative area'\n        if (area2 .eq. 0.0)   status = \
'SM3SHMEMBB: Zero area'\n        return\n      end if",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)      \n      \
{x21,x32,x13}={x[[2]]-x[[1]],x[[3]]-x[[2]],x[[1]]-x[[3]]};\n      \
{x12,x23,x31}=-{x21,x32,x13};\n      \
{y21,y32,y13}={y[[2]]-y[[1]],y[[3]]-y[[2]],y[[1]]-y[[3]]};\n      \
{y12,y23,y31}=-{y21,y32,y13};\n      area2=y21*x13-x21*y13;\n      If \
[area2<=0, Print[\"Negative or zero area\"]; Return[Null]];\n(*",
  StyleBox["\n      p(1,1) =   y23\n      p(2,1) =   0.0\n      p(3,1) =   \
y31\n      p(4,1) =   0.0\n      p(5,1) =   y12\n      p(6,1) =   0.0\n      \
p(1,2) =   0.0\n      p(2,2) =   x32\n      p(3,2) =   0.0\n      p(4,2) =   \
x13\n      p(5,2) =   0.0\n      p(6,2) =   x21\n      p(1,3) =   x32\n      \
p(2,3) =   y23\n      p(3,3) =   x13\n      p(4,3) =   y31\n      p(5,3) =   \
x21\n      p(6,3) =   y12\n      n =        6\n      if (alpha .ne. 0.0)      \
   then\n        c =      alpha/6.0d0\n        p(7,1) =   y23*(y13-y21)    *c\
\n        p(7,2) =   x32*(x31-x12)    *c\n        p(7,3) =  (x31*y13-x12*y21) \
*c*2.0d0\n        p(8,1) =   y31*(y21-y32)    *c\n        p(8,2) =   \
x13*(x12-x23)    *c\n        p(8,3) =  (x12*y21-x23*y32) *c*2.0d0\n        \
p(9,1) =   y12*(y32-y13)    *c\n        p(9,2) =   x21*(x23-x31)    *c\n      \
  p(9,3) =  (x23*y32-x31*y13) *c*2.0d0\n        n =        9\n      end if",
    FontColor->RGBColor[0, 0, 1]],
  "\n*) \n      p={{y23,0,x32},{0,x32,y23},{y31,0,x13},{0,x13,y31},\n         \
{y12,0,x21},{0,x21,y12},{y23*(y13-y21)*alpha/6,\n          \
x32*(x31-x12)*alpha/6,(x31*y13-x12*y21)*alpha/3},\n         \
{y31*(y21-y32)*alpha/6,x13*(x12-x23)*alpha/6,\n         \
(x12*y21-x23*y32)*alpha/3},{y12*(y32-y13)*alpha/6,\n          \
x21*(x23-x31)*alpha/6,(x23*y32-x31*y13)*alpha/3}};",
  StyleBox["\n(*       \n      c =        0.5d0*f/area2\n      dm11 =     c * \
dm(1,1)\n      dm22 =     c * dm(2,2)\n      dm33 =     c * dm(3,3)\n      \
dm12 =     c * dm(1,2)\n      dm13 =     c * dm(1,3)\n      dm23 =     c * \
dm(2,3)\n      do 3000  j = 1,n\n        l =      ls(j)\n        s1 =     \
dm11*p(j,1) + dm12*p(j,2) + dm13*p(j,3)\n        s2 =     dm12*p(j,1) + \
dm22*p(j,2) + dm23*p(j,3)\n        s3 =     dm13*p(j,1) + dm23*p(j,2) + \
dm33*p(j,3)\n        do 2500  i = 1,j\n          k =      ls(i)\n          \
sm(k,l) =  sm(k,l) + (s1*p(i,1) + s2*p(i,2) + s3*p(i,3))\n          sm(l,k) = \
 sm(k,l)\n 2500     continue\n 3000   continue",
    FontColor->RGBColor[0, 0, 1]],
  "\n",
  StyleBox["      return\n      end\nC=END FORTRAN\n",
    FontColor->RGBColor[0, 0, 1]],
  "*)\n\n     {{dm11,dm12,dm13},{dm21,dm22,dm23},{dm31,dm32,dm33}}=\n         \
 dm*(f/(2*area2)); \n      For [j=1,j<=9,j++, l=ls[[j]]; \n          \
s1=dm11*p[[j,1]]+dm12*p[[j,2]]+dm13*p[[j,3]];\n          \
s2=dm12*p[[j,1]]+dm22*p[[j,2]]+dm23*p[[j,3]];\n          \
s3=dm13*p[[j,1]]+dm23*p[[j,2]]+dm33*p[[j,3]];\n          For [i=1,i<=j,i++, \
k=ls[[i]];\n            sm[[k,l]]+=s1*p[[i,1]]+s2*p[[i,2]]+s3*p[[i,3]];\n     \
       sm[[l,k]]=sm[[k,l]] ];\n          ];\n      Return[sm];\n     ];\n     \
\n",
  StyleBox["     \nClearAll[Em,nu,h0];  Em=4000; nu=1/3; h0=0.2; \n\
x=N[{0,2,0}]; y=N[{0,0,2}];\nc=Em*h0/(1-nu^2);\n\
Dm={{c,nu*c,0},{nu*c,c,0},{0,0,c*(1-nu)/2}};Print[\"Dm=\",Dm];\n\
sm=SM3SHMEMBB[x,y,Dm,3/2,1,{1,2, 4,5, 7,8, 3,6,9},Table[0,{9},{9}]];\n\
Print[sm];\nPrint[ssm-sm];\nPrint[Chop[Eigenvalues[N[sm]]]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    InterpretationBox[\("Dm="\[InvisibleSpace]{{900.`, 300.`, 0}, {300.`, 
            900.`, 0}, {0, 0, 300.`}}\),
      SequenceForm[ 
      "Dm=", {{900.0, 300.0, 0}, {300.0, 900.0, 0}, {0, 0, 300.0}}],
      Editable->False]], "Print"],

Cell[BoxData[
    \({{600.`, 300.`, \(-150.`\), \(-450.`\), \(-150.`\), 
        375.`, \(-150.`\), \(-150.`\), \(-225.`\)}, {300.`, 600.`, 
        150.`, \(-150.`\), \(-150.`\), 
        225.`, \(-150.`\), \(-450.`\), \(-375.`\)}, {\(-150.`\), 150.`, 
        150.`, 150.`, 0.`, \(-75.`\), 
        0.`, \(-150.`\), \(-75.`\)}, {\(-450.`\), \(-150.`\), 150.`, 450.`, 
        0.`, \(-225.`\), 0.`, 150.`, 75.`}, {\(-150.`\), \(-150.`\), 0.`, 
        0.`, 150.`, \(-150.`\), 150.`, 0.`, 150.`}, {375.`, 
        225.`, \(-75.`\), \(-225.`\), \(-150.`\), 
        262.5`, \(-150.`\), \(-75.`\), \(-187.5`\)}, {\(-150.`\), \(-150.`\), 
        0.`, 0.`, 150.`, \(-150.`\), 150.`, 0.`, 
        150.`}, {\(-150.`\), \(-450.`\), \(-150.`\), 150.`, 0.`, \(-75.`\), 
        0.`, 450.`, 225.`}, {\(-225.`\), \(-375.`\), \(-75.`\), 75.`, 
        150.`, \(-187.5`\), 150.`, 225.`, 262.5`}}\)], "Print"],

Cell[BoxData[
    \({{\(-600.`\) + ssm, \(-300.`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(450.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-375.`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(225.`\)\(\[InvisibleSpace]\)\) + ssm}, {\(-300.`\) + 
          ssm, \(-600.`\) + ssm, \(-150.`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-225.`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(450.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(375.`\)\(\[InvisibleSpace]\)\) + 
          ssm}, {\(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(-150.`\) + ssm, \(-150.`\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(75.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(75.`\)\(\[InvisibleSpace]\)\) + 
          ssm}, {\(\(450.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(-450.`\) + ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(225.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(-75.`\) + ssm}, {\(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm}, {\(-375.`\) + ssm, \(-225.`\) + 
          ssm, \(\(75.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(225.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-262.5`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(75.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(187.5`\)\(\[InvisibleSpace]\)\) + 
          ssm}, {\(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm}, {\(\(150.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(450.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(150.`\)\(\[InvisibleSpace]\)\) + ssm, \(-150.`\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(75.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(0.`\)\(\[InvisibleSpace]\)\) + ssm, \(-450.`\) + 
          ssm, \(-225.`\) + ssm}, {\(\(225.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(375.`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(\(75.`\)\(\[InvisibleSpace]\)\) + ssm, \(-75.`\) + 
          ssm, \(-150.`\) + ssm, \(\(187.5`\)\(\[InvisibleSpace]\)\) + 
          ssm, \(-150.`\) + ssm, \(-225.`\) + ssm, \(-262.5`\) + 
          ssm}}\)], "Print"],

Cell[BoxData[
    \({1799.9999999999993`, 824.9999999999995`, 450.00000000000017`, 0, 0, 0, 
      0, 0, 0}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  "(*",
  StyleBox["\nC=DECK SM3SHMEMBH\nC=PURPOSE Form H.O. material stiffness of \
9-dof ANDES membrane triangle\nC=AUTHOR C. A. Felippa, June 1991\nC=VERSION \
July 1991\nC=EQUIPMENT Machine independent\nC=KEYWORDS finite element\n\
C=KEYWORDS material stiffness matrix high-order\nC=KEYWORDS triangle membrane \
assumed natural deviatoric strain\nC=BLOCK ABSTRACT\nC\nC     SM3MANDES forms \
the higher order element stiffness matrix\nC     of a 9-dof membrane triangle \
based on the ANDES formulation.\nC     Implementation moderately optimized \
for speed.\nC\nC=END ABSTRACT\nC=BLOCK USAGE\nC\nC     The calling sequence \
is\nC\nC       call SM3SHMEMBH (X, Y, DM, F, LS, SM, M, STATUS)\nC\nC     The \
inputs are:\nC\nC     X         (3 x 1) array of x coordinates of triangle \
nodes\nC     Y         (3 x 1) array of y coordinates of triangle nodes\nC    \
 DM        (3 x 3) membrane constitutive matrix already\nC                 \
integrated through the thickness\nC     F         Factor by which all \
stiffness entries will be multiplied.\nC     SM        Incoming material \
stiffness array.\nC     LS        (9 x 1) array of stiffness location \
pointers\nC               (see examples in SM3SHMEMBB)\nC     M         First \
dimension of SM in calling program.\nC\nC     The outputs are:\nC\nC     SM   \
     Output stiffness array with higher order stiffness\nC               \
coefficients added in.\nC               The (i,j)-th entry of the basic \
element stiffness is added\nC               to SM(K,L), where K=LS(I) and \
L=LS(J).\nC     STATUS    Status character variable.  Blank if no error\nC    \
           detected.\nC\nC=END USAGE\nC=BLOCK FORTRAN\n      subroutine    \
SM3SHMEMBH\n     &              (x, y, dm, f, ls, sm, m, status)\nC\nC        \
           A R G U M E N T S\nC\n      integer           ls(9), m\n      \
double precision  x(3),y(3), dm(3,3), f, sm(m,m)\n      character         \
status* ( * )\nC\nC                   L O C A L   V A R I A B L E S\nC\n      \
double precision  x12, x21, x23, x32, x31, x13\n      double precision  y12, \
y21, y23, y32, y31, y13\n      double precision  l21,l32,l13\n      double \
precision  chi213,chi321,chi132\n      double precision  area, area2, area43\n\
      double precision  c(3,3), e(3,3), et(3), d(3), qm(3,3,3)\n      double \
precision  t(3,3), tfac, kth(3,3)\n      double precision  s(3), xyij(6), \
sum, w(3), wfac\n      integer           i, j, k, l",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n\n   SM3SHMEMBH[x_,y_,dm_,f_,ls_,esm_]:=Module[{\n      x21, x32, \
x13, y21, y32, y13, x12, x23, x31, y12, y23, y31,\n      ll21,ll32,ll13, \
chi213,chi321,chi132, area,area2,area43,\n      mu11, mu12, mu13, mu22, mu23, \
mu33,\n      c11,c12,c13,c21,c22,c23, c31,c32,c33, sum,wfac,tfac, \n      \
et1,et2,et3, i,j,k,l,\n      d=s=w={0,0,0}, e=c=t=kth=Table[0,{3},{3}], xyij,\
\n      qm=Table[0,{3},{3},{3}], sm=esm },\n(*\n",
  StyleBox["C\nC                   L O G I C\nC\n      status =   ' '\n      \
if (f .eq. 0.0)         return\n      x12 =      x(1) - x(2)\n      x21 =     \
-x12\n      x23 =      x(2) - x(3)\n      x32 =     -x23\n      x31 =      \
x(3) - x(1)\n      x13 =     -x31\n      y12 =      y(1) - y(2)\n      y21 =  \
   -y12\n      y23 =      y(2) - y(3)\n      y32 =     -y23\n      y31 =      \
y(3) - y(1)\n      y13 =     -y31\n      area2 =    x21*y31-x31*y21\n      if \
(area2 .le. 0.0)      then\n        status = 'SM3SHMEMBH: Negative area'\n    \
    if (area2 .eq. 0.0)   status = 'SM3SHMEMBH: Zero area'\n        return\n  \
    end if\n      area  =     0.5d0*area2",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)         \n      \
{x21,x32,x13}={x[[2]]-x[[1]],x[[3]]-x[[2]],x[[1]]-x[[3]]};\n      \
{x12,x23,x31}=-{x21,x32,x13};\n      \
{y21,y32,y13}={y[[2]]-y[[1]],y[[3]]-y[[2]],y[[1]]-y[[3]]};\n      \
{y12,y23,y31}=-{y21,y32,y13};\n      area2=y21*x13-x21*y13;\n      If \
[area2<=0, Print[\"Negative or zero area\"]; Return[Null]];\n      \
area=area2/2;\n",
  StyleBox["(*\n      l21 =       sqrt(x21**2+y21**2)\n      l32 =       \
sqrt(x32**2+y32**2)\n      l13 =       sqrt(x13**2+y13**2)\n      tfac =      \
0.25d0/area**2\n      t(1,1) =    tfac*y23*y13*l21**2\n      t(1,2) =    \
tfac*y31*y21*l32**2\n      t(1,3) =    tfac*y12*y32*l13**2\n      t(2,1) =    \
tfac*x23*x13*l21**2\n      t(2,2) =    tfac*x31*x21*l32**2\n      t(2,3) =    \
tfac*x12*x32*l13**2\n      t(3,1) =    tfac*(y23*x31+x32*y13)*l21**2\n      \
t(3,2) =    tfac*(y31*x12+x13*y21)*l32**2\n      t(3,3) =    \
tfac*(y12*x23+x21*y32)*l13**2\n      wfac =      0.75d0*f*area\n      e(1,1) \
=    wfac*dm(1,1)\n      e(1,2) =    wfac*dm(1,2)\n      e(1,3) =    \
wfac*dm(1,3)\n      e(2,1) =    wfac*dm(2,1)\n      e(2,2) =    wfac*dm(2,2)\n\
      e(2,3) =    wfac*dm(2,3)\n      e(3,1) =    wfac*dm(3,1)\n      e(3,2) \
=    wfac*dm(3,2)\n      e(3,3) =    wfac*dm(3,3)\n      do 1600  j = 1,3\n   \
     do 1400  i = 1,3\n          et(i) =  \
e(i,1)*t(1,j)+e(i,2)*t(2,j)+e(i,3)*t(3,j)\n 1400     continue\n        do \
1500  i = 1,3\n          c(i,j) = t(1,i)*et(1)+t(2,i)*et(2)+t(3,i)*et(3)\n \
1500     continue\n 1600   continue",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n      ll21=x21^2+y21^2; ll32=x32^2+y32^2; ll13=x13^2+y13^2;\n      \
tfac=1/(4*area^2);\n      t={{y23*y13*ll21,y31*y21*ll32,y12*y32*ll13},\n      \
   {x23*x13*ll21,x31*x21*ll32,x12*x32*ll13},\n         \
{(y23*x31+x32*y13)*ll21,(y31*x12+x13*y21)*ll32,\n          \
(y12*x23+x21*y32)*ll13}}*tfac; \n      wfac=(3/4)*f*area; e=dm*wfac; \n      \
For [j=1,j<=3,j++,\n          \
et1=e[[1,1]]*t[[1,j]]+e[[1,2]]*t[[2,j]]+e[[1,3]]*t[[3,j]];\n          \
et2=e[[2,1]]*t[[1,j]]+e[[2,2]]*t[[2,j]]+e[[2,3]]*t[[3,j]];\n          \
et3=e[[3,1]]*t[[1,j]]+e[[3,2]]*t[[2,j]]+e[[3,3]]*t[[3,j]];\n          For \
[i=1,i<=3,i++,\n              \
c[[i,j]]=t[[1,i]]*et1+t[[2,i]]*et2+t[[3,i]]*et3];\n          ];   \n(*\n",
  StyleBox["      area43 =     (2.d0/3.d0)*area2\n      chi213 =     \
area43/l21**2\n      chi321 =     area43/l32**2\n      chi132 =     \
area43/l13**2\n      qm(1,1,1) =  -0.25d0*chi213\n      qm(1,2,1) =  \
-qm(1,1,1)\n      qm(1,3,1) =   0.0\n      qm(2,1,1) =   0.25d0*chi321\n      \
qm(2,2,1) =   0.50d0*chi321\n      qm(2,3,1) =   qm(2,1,1)\n      qm(3,1,1) = \
 -0.50d0*chi132\n      qm(3,2,1) =  -0.25d0*chi132\n      qm(3,3,1) =   \
qm(3,2,1)\n      qm(1,1,2) =  -0.25d0*chi213\n      qm(1,2,2) =  \
-0.50d0*chi213\n      qm(1,3,2) =   qm(1,1,2)\n      qm(2,1,2) =   0.0\n      \
qm(2,2,2) =  -0.25d0*chi321\n      qm(2,3,2) =  -qm(2,2,2)\n      qm(3,1,2) = \
  0.25d0*chi132\n      qm(3,2,2) =   qm(3,1,2)\n      qm(3,3,2) =   \
0.50d0*chi132\n      qm(1,1,3) =   0.50d0*chi213\n      qm(1,2,3) =   \
0.25d0*chi213\n      qm(1,3,3) =   qm(1,2,3)\n      qm(2,1,3) =  \
-0.25d0*chi321\n      qm(2,2,3) =   qm(2,1,3)\n      qm(2,3,3) =  \
-0.50d0*chi321\n      qm(3,1,3) =   0.25d0*chi132\n      qm(3,2,3) =   0.0\n  \
    qm(3,3,3) =  -qm(3,1,3)\n      kth(1,1) =   0.0\n      kth(1,2) =   0.0\n \
     kth(1,3) =   0.0\n      kth(2,2) =   0.0\n      kth(2,3) =   0.0\n      \
kth(3,3) =   0.0\n      do 2800  k = 1,3\n        do 2600  j = 1,3\n          \
d(1) = c(1,1)*qm(1,j,k)+c(1,2)*qm(2,j,k)+c(1,3)*qm(3,j,k)\n          d(2) = \
c(2,1)*qm(1,j,k)+c(2,2)*qm(2,j,k)+c(2,3)*qm(3,j,k)\n          d(3) = \
c(3,1)*qm(1,j,k)+c(3,2)*qm(2,j,k)+c(3,3)*qm(3,j,k)\n          do 2500  i = \
1,j\n            kth(i,j) = kth(i,j) +\n     $                 \
qm(1,i,k)*d(1)+qm(2,i,k)*d(2)+qm(3,i,k)*d(3)\n            kth(j,i) = kth(i,j)\
\n 2500       continue\n 2600     continue\n 2800   continue",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n      area43=(2/3)*area2;    chi213=area43/ll21;\n      \
chi321=area43/ll32;    chi132=area43/ll13;\n      qm[[1,1,1]]=-(1/4)*chi213; \
qm[[1,2,1]]=-qm[[1,1,1]];\n      qm[[2,1,1]]= (1/4)*chi321; qm[[2,2,1]]= \
(1/2)*chi321;\n      qm[[2,3,1]]= qm[[2,1,1]];  qm[[3,1,1]]=-(1/2)*chi132;\n  \
    qm[[3,2,1]]=-(1/4)*chi132; qm[[3,3,1]]= qm[[3,2,1]];\n      \
qm[[1,1,2]]=-(1/4)*chi213; qm[[1,2,2]]=-(1/2)*chi213;\n      qm[[1,3,2]]= \
qm[[1,1,2]];  qm[[2,2,2]]=-(1/4)*chi321;\n      qm[[2,3,2]]=-qm[[2,2,2]];  \
qm[[3,1,2]]= (1/4)*chi132;\n      qm[[3,2,2]]= qm[[3,1,2]];  qm[[3,3,2]]= \
(1/2)*chi132;\n      qm[[1,1,3]]= (1/2)*chi213; qm[[1,2,3]]= (1/4)*chi213;\n  \
    qm[[1,3,3]]= qm[[1,2,3]];  qm[[2,1,3]]=-(1/4)*chi321;\n      qm[[2,2,3]]= \
qm[[2,1,3]];  qm[[2,3,3]]=-(1/2)*chi321;\n      qm[[3,1,3]]= (1/4)*chi132; \
qm[[3,3,3]]=-qm[[3,1,3]];\n      \
{{c11,c12,c13},{c21,c22,c23},{c31,c32,c33}}=c;\n      kth=Table[0,{3},{3}];\n \
     For [k=1,k<=3,k++,\n          For [j=1,j<=3,j++,\n             \
d={c11*qm[[1,j,k]]+c12*qm[[2,j,k]]+c13*qm[[3,j,k]],\n                \
c21*qm[[1,j,k]]+c22*qm[[2,j,k]]+c23*qm[[3,j,k]],\n                \
c31*qm[[1,j,k]]+c32*qm[[2,j,k]]+c33*qm[[3,j,k]]};\n             For \
[i=1,i<=j,i++,\n              \
kth[[i,j]]+=qm[[1,i,k]]*d[[1]]+qm[[2,i,k]]*d[[2]]+\n                          \
qm[[3,i,k]]*d[[3]]; kth[[j,i]]=kth[[i,j]]\n          ]  ]  ];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "    \n         \n(*\n     ",
  StyleBox[" s(1) =   kth(1,1) + kth(1,2) + kth(1,3)\n      s(2) =   kth(2,1) \
+ kth(2,2) + kth(2,3)\n      s(3) =   kth(3,1) + kth(3,2) + kth(3,3)\n      \
xyij(1) =  0.25d0*x32/area\n      xyij(2) =  0.25d0*y32/area\n      xyij(3) = \
 0.25d0*x13/area\n      xyij(4) =  0.25d0*y13/area\n      xyij(5) =  \
0.25d0*x21/area\n      xyij(6) =  0.25d0*y21/area\n      do 4000  j = 1,9\n   \
     l =    ls(j)\n        do 3600  i = 1,3\n          if (j .le. 6)      \
then\n            w(i) =  s(i)*xyij(j)\n          else\n            w(i) =  \
kth(i,j-6)\n          end if\n 3600     continue\n        sum =    w(1) + \
w(2) + w(3)\n        do 3700  i = 1,j\n          k =      ls(i)\n          if \
(i .le. 6)      then\n             sm(k,l) =  sm(k,l) + sum*xyij(i)\n         \
 else\n             sm(k,l) =  sm(k,l) + w(i-6)\n          end if\n          \
sm(l,k) =  sm(k,l)\n 3700     continue\n 4000   continue\n      return\n      \
end\nC=END FORTRAN",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)  \n    s={kth[[1,1]]+kth[[1,2]]+kth[[1,3]],\n       \
kth[[2,1]]+kth[[2,2]]+kth[[2,3]],\n       kth[[3,1]]+kth[[3,2]]+kth[[3,3]]};\n\
    xyij={x32,y32,x13,y13,x21,y21}/(4*area);\n    For [j=1,j<=9,j++, \
l=ls[[j]];\n        For [i=1,i<=3,i++,\n            If \
[j<=6,w[[i]]=s[[i]]*xyij[[j]],w[[i]]=kth[[i,j-6]]]];\n        \
sum=w[[1]]+w[[2]]+w[[3]];\n        For [i=1,i<=j,i++, k=ls[[i]];\n        If \
[i <=6, sm[[k,l]]+=sum*xyij[[i]],sm[[k,l]]+=w[[i-6]]];\n        \
sm[[l,k]]=sm[[k,l]] ];\n        ];\n    Return[sm]];\n",
  StyleBox["     \nClearAll[Em,nu,h0];  Em=4000; nu=1/3; h0=0.2; \n\
x=N[{0,2,0}]; y=N[{0,0,2}];\nc=Em*h0/(1-nu^2);\n\
Dm={{c,nu*c,0},{nu*c,c,0},{0,0,c*(1-nu)/2}};Print[\"Dm=\",Dm];\n\
betab=(1-4*nu^2)/2; Print[\"betab=\",betab];\n\
sm=SM3SHMEMBH[x,y,Dm,betab,{1,2, 4,5, 7,8, 3,6,9},Table[0,{9},{9}]];\n\
Print[sm];\nPrint[Chop[Eigenvalues[N[sm]]]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    InterpretationBox[\("Dm="\[InvisibleSpace]{{900.`, 300.`, 0}, {300.`, 
            900.`, 0}, {0, 0, 300.`}}\),
      SequenceForm[ 
      "Dm=", {{900.0, 300.0, 0}, {300.0, 900.0, 0}, {0, 0, 300.0}}],
      Editable->False]], "Print"],

Cell[BoxData[
    InterpretationBox[\("betab="\[InvisibleSpace]5\/18\),
      SequenceForm[ "betab=", 
        Rational[ 5, 18]],
      Editable->False]], "Print"],

Cell[BoxData[
    \({{62.5`, \(-62.5`\), \(-118.05555555555554`\), 0.`, 
        62.5`, \(-65.97222222222221`\), \(-62.5`\), 
        0.`, \(-65.97222222222221`\)}, {\(-62.5`\), 62.5`, 
        118.05555555555554`, 0.`, \(-62.5`\), 65.97222222222221`, 62.5`, 0.`, 
        65.97222222222221`}, {\(-118.05555555555554`\), 118.05555555555554`, 
        256.94444444444446`, 0.`, \(-118.05555555555554`\), 
        107.63888888888887`, 118.05555555555554`, 0.`, 
        107.63888888888889`}, {0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 
        0.`}, {62.5`, \(-62.5`\), \(-118.05555555555554`\), 0.`, 
        62.5`, \(-65.97222222222221`\), \(-62.5`\), 
        0.`, \(-65.97222222222221`\)}, {\(-65.97222222222221`\), 
        65.97222222222221`, 107.63888888888887`, 
        0.`, \(-65.97222222222221`\), 90.27777777777777`, 65.97222222222221`, 
        0.`, 65.97222222222221`}, {\(-62.5`\), 62.5`, 118.05555555555554`, 
        0.`, \(-62.5`\), 65.97222222222221`, 62.5`, 0.`, 
        65.97222222222221`}, {0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 
        0.`}, {\(-65.97222222222221`\), 65.97222222222221`, 
        107.63888888888889`, 0.`, \(-65.97222222222221`\), 
        65.97222222222221`, 65.97222222222221`, 0.`, 
        90.27777777777777`}}\)], "Print"],

Cell[BoxData[
    \({614.8814040291797`, 48.31304041526528`, 24.305555555555557`, 0, 0, 0, 
      0, 0, 0}\)], "Print"]
}, Open  ]],

Cell["n", "Input",
  ImageRegion->{{0, 1}, {0, 1}}],

Cell["\<\
SM3SHBENDH forms the higher order material stiffness matrix of a
9-dof thin-plate-bending triangle obtained by using linear
curvatures over the sides.  Implementation optimized for speed\
\>", "Text",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}}],

Cell[CellGroupData[{

Cell[TextData[{
  "(*",
  StyleBox["\nC=DECK SM3SHBENDH\nC=PURPOSE Form Kh for 3-node Kirchhoff ANDES \
plate bending elem\nC=AUTHOR C. Militello and C. Felippa, May 1989\nC=VERSION \
July 1989\nC=EQUIPMENT Machine independent\nC=KEYWORDS thin plate bending\n\
C=KEYWORDS finite element triangle higher order stiffness matrix\nC=BLOCK \
ABSTRACT\nC\nC     SM3SHBENDH forms the higher order material stiffness \
matrix of a\nC     9-dof thin-plate-bending triangle obtained by using linear\
\nC     curvatures over the sides.  Implementation optimized for speed\nC\n\
C=END ABSTRACT\nC=BLOCK USAGE\nC\nC     The calling sequence is\nC\nC       \
CALL      SM3SHBENDH (X, Y, DB, F, LS, SM, M, STATUS)\nC\nC     where the \
input arguments are\nC\nC       X         (3 x 1) array of x coordinates of \
triangle nodes\nC       Y         (3 x 1) array of y coordinates of triangle \
nodes\nC       DB        (3 x 3) moment-curvature matrix.\nC       F         \
Factor by which stiffness entries will be multiplied.\nC       LS        (9 x \
1) array of stiffness location pointers\nC                 (see Output SM).\n\
C       SM        Incoming material stiffness array.\nC       M         First \
dimension of SM in calling program.\nC\nC     The outputs are:\nC\nC       SM \
       Output stiffness array with higher order stiffness\nC                 \
coefficients added in.  The (i,j)-th entry of the\nC                 (9 by 9) \
element bending stiffness is added to\nC                 SM(K,L), where \
K=LS(I) and L=LS(J).\nC       STATUS    Status character variable.  Blank if \
no error\nC                 detected.\nC\nC=END USAGE\nC=BLOCK FORTRAN\n      \
subroutine  SM3SHBENDH (x, y, db, f, ls, sm, m, status)\nC\nC                 \
  A R G U M E N T S\nC\n      character* ( * )     status\n      integer      \
     ls(9), m\n      double precision  x(3), y(3), db(3,3), f\n      double \
precision  sm(m,m)\nC\nC                   L O C A L   V A R I A B L E S\nC\n \
     double precision  x21, x32, x13, y21, y32, y13\n      double precision  \
x12, x23, x31, y12, y23, y31\n      double precision  area, area2, tfac, cfac\
\n      double precision  lam12, lam23, lam31\n      double precision  dbt11, \
dbt12, dbt13, dbt21, dbt22, dbt23\n      double precision  dbt31, dbt32, \
dbt33\n      double precision  t(3,3), sh(9,9)\n      double precision  mu11, \
mu12, mu13, mu22, mu23, mu33\n      double precision  c11, c12, c13, c22, \
c23, c33\n      integer           i, j, k, l\n",
    FontColor->RGBColor[0, 0, 1]],
  "*)\n      SM3SHBENDH[x_,y_,db_,f_,ls_,esm_]:=Module[{\n      x21, x32, \
x13, y21, y32, y13, x12, x23, x31, y12, y23, y31,\n      area, area2, tfac, \
cfac,lam12, lam23, lam31,\n      dbt11, dbt12, dbt13, dbt21, dbt22, dbt23, \
dbt31, dbt32, dbt33,\n      mu11, mu12, mu13, mu22, mu23, mu33,\n      c11, \
c12, c13, c22, c23, c33, i, j, k, l,  \n      t=Table[0,{3},{3}],  \
sh=Table[0,{9},{9}], sm=esm },\n      ",
  StyleBox["\n",
    FontWeight->"Plain",
    FontColor->RGBColor[1, 0, 0]],
  "(*",
  StyleBox["\nC                   L O G I C\nC\n      status =   ' '\n      \
x21 =  x(2) - x(1)\n      x12 =  -x21\n      x32 =  x(3) - x(2)\n      x23 =  \
-x32\n      x13 =  x(1) - x(3)\n      x31 =  -x13\n      y21 =  y(2) - y(1)\n \
     y12 =  -y21\n      y32 =  y(3) - y(2)\n      y23 =  -y32\n      y13 =  \
y(1) - y(3)\n      y31 =  -y13\n      area2 =    y21*x13 - x21*y13\n      if \
(area2 .le. 0.0)      then\n        status = 'SM3SHBENDH: Negative area'\n    \
    if (area2 .eq. 0.0)   status = 'SM3SHBENDH: Zero area'\n        return\n  \
    end if",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)    \n      {x21,x32,x13}={x[[2]]-x[[1]],x[[3]]-x[[2]],x[[1]]-x[[3]]};\
\n      {x12,x23,x31}=-{x21,x32,x13};\n      \
{y21,y32,y13}={y[[2]]-y[[1]],y[[3]]-y[[2]],y[[1]]-y[[3]]};\n      \
{y12,y23,y31}=-{y21,y32,y13};\n      area2=y21*x13-x21*y13;\n      If \
[area2<=0, Print[\"Negative or zero area\"]; Return[Null]];\n(*      \n    ",
  StyleBox["  lam12 =  (x12*x13+y12*y13)/(x21*x21+y21*y21)\n      lam23 =  \
(x23*x21+y23*y21)/(x32*x32+y32*y32)\n      lam31 =  \
(x31*x32+y31*y32)/(x13*x13+y13*y13)\n      mu11 =   2.d0*(lam12**2+1-lam12)\n \
     mu22 =   2.d0*(lam23**2+1-lam23)\n      mu33 =   2.d0*(lam31**2+1-lam31)\
\n      mu12 =   2.d0*lam12-1-(1+lam12)*lam23\n      mu23 =   \
2.d0*lam23-1-(1+lam23)*lam31\n      mu13 =   2.d0*lam31-1-(1+lam31)*lam12",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)               \n      lam12=(x12*x13+y12*y13)/(x21*x21+y21*y21);\n    \
  lam23=(x23*x21+y23*y21)/(x32*x32+y32*y32);\n      \
lam31=(x31*x32+y31*y32)/(x13*x13+y13*y13);\n      mu11=2*(lam12^2+1-lam12);\n \
     mu22=2*(lam23^2+1-lam23);\n      mu33=2*(lam31^2+1-lam31);\n      \
mu12=2*lam12-1-(1+lam12)*lam23;\n      mu23=2*lam23-1-(1+lam23)*lam31;\n      \
mu13=2*lam31-1-(1+lam31)*lam12; ",
  StyleBox["\n",
    FontWeight->"Plain"],
  "(*      \n",
  StyleBox["      tfac =     1.d0/area2**2\n      t(1,1) =   tfac* y23*y13 \n \
     t(1,2) =   tfac* y31*y21 \n      t(1,3) =   tfac* y12*y32\n      t(2,1) \
=   tfac* x23*x13 \n      t(2,2) =   tfac* x31*x21 \n      t(2,3) =   tfac* \
x12*x32 \n      t(3,1) =   tfac* (x31*y23+x32*y13)\n      t(3,2) =   tfac* \
(x12*y31+x13*y21) \n      t(3,3) =   tfac* (x23*y12+x21*y32)\n      dbt11 =   \
db(1,1) *t(1,1)+db(1,2) *t(2,1)+db(1,3) *t(3,1)\n      dbt12 =   db(1,1) \
*t(1,2)+db(1,2) *t(2,2)+db(1,3) *t(3,2)\n      dbt13 =   db(1,1) \
*t(1,3)+db(1,2) *t(2,3)+db(1,3) *t(3,3)\n      dbt21 =   db(2,1) \
*t(1,1)+db(2,2) *t(2,1)+db(2,3) *t(3,1)\n      dbt22 =   db(2,1) \
*t(1,2)+db(2,2) *t(2,2)+db(2,3) *t(3,2)\n      dbt23 =   db(2,1) \
*t(1,3)+db(2,2) *t(2,3)+db(2,3) *t(3,3)\n      dbt31 =   db(3,1) \
*t(1,1)+db(3,2) *t(2,1)+db(3,3) *t(3,1)\n      dbt32 =   db(3,1) \
*t(1,2)+db(3,2) *t(2,2)+db(3,3) *t(3,2)\n      dbt33 =   db(3,1) \
*t(1,3)+db(3,2) *t(2,3)+db(3,3) *t(3,3)\n",
    FontColor->RGBColor[0, 0, 1]],
  "*) \n      t={{y23*y13,y31*y21,y12*y32},{x23*x13,x31*x21,x12*x32},\n       \
 {x31*y23+x32*y13,x12*y31+x13*y21,x23*y12+x21*y32}}/area2^2;  \n      \
dbt11=db[[1,1]]*t[[1,1]]+db[[1,2]]*t[[2,1]]+db[[1,3]]*t[[3,1]];\n      \
dbt12=db[[1,1]]*t[[1,2]]+db[[1,2]]*t[[2,2]]+db[[1,3]]*t[[3,2]];\n      \
dbt13=db[[1,1]]*t[[1,3]]+db[[1,2]]*t[[2,3]]+db[[1,3]]*t[[3,3]];\n      \
dbt21=db[[2,1]]*t[[1,1]]+db[[2,2]]*t[[2,1]]+db[[2,3]]*t[[3,1]];\n      \
dbt22=db[[2,1]]*t[[1,2]]+db[[2,2]]*t[[2,2]]+db[[2,3]]*t[[3,2]];\n      \
dbt23=db[[2,1]]*t[[1,3]]+db[[2,2]]*t[[2,3]]+db[[2,3]]*t[[3,3]];\n      \
dbt31=db[[3,1]]*t[[1,1]]+db[[3,2]]*t[[2,1]]+db[[3,3]]*t[[3,1]];\n      \
dbt32=db[[3,1]]*t[[1,2]]+db[[3,2]]*t[[2,2]]+db[[3,3]]*t[[3,2]];\n      \
dbt33=db[[3,1]]*t[[1,3]]+db[[3,2]]*t[[2,3]]+db[[3,3]]*t[[3,3]];\n(*\n    ",
  StyleBox["  area =  0.5d0*area2\n      cfac =  area*f\n      c11=    \
cfac*mu11* (t(1,1) *dbt11+t(2,1) *dbt21+t(3,1) *dbt31)\n      c12=    \
cfac*mu12* (t(1,2) *dbt11+t(2,2) *dbt21+t(3,2) *dbt31)\n      c13=    \
cfac*mu13* (t(1,3) *dbt11+t(2,3) *dbt21+t(3,3) *dbt31)\n      c22=    \
cfac*mu22* (t(1,2) *dbt12+t(2,2) *dbt22+t(3,2) *dbt32)\n      c23=    \
cfac*mu23* (t(1,3) *dbt12+t(2,3) *dbt22+t(3,3) *dbt32)\n      c33=    \
cfac*mu33* (t(1,3) *dbt13+t(2,3) *dbt23+t(3,3) *dbt33)",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n      area=area2/2; cfac=area*f;\n      \
c11=cfac*mu11*(t[[1,1]]*dbt11+t[[2,1]]*dbt21+t[[3,1]]*dbt31);\n      \
c12=cfac*mu12*(t[[1,2]]*dbt11+t[[2,2]]*dbt21+t[[3,2]]*dbt31);\n      \
c13=cfac*mu13*(t[[1,3]]*dbt11+t[[2,3]]*dbt21+t[[3,3]]*dbt31);\n      \
c22=cfac*mu22*(t[[1,2]]*dbt12+t[[2,2]]*dbt22+t[[3,2]]*dbt32);\n      \
c23=cfac*mu23*(t[[1,3]]*dbt12+t[[2,3]]*dbt22+t[[3,3]]*dbt32);\n      \
c33=cfac*mu33*(t[[1,3]]*dbt13+t[[2,3]]*dbt23+t[[3,3]]*dbt33);\n      \n(*",
  StyleBox["\n      sh(1,1) =  4.d0*(c33-c13-c13+c11)\n      sh(1,2) =  \
2.d0*((c11-c13)*y21+(c13-c33)*y13)\n      sh(1,3) =  \
2.d0*((c13-c11)*x21+(c33-c13)*x13)\n      sh(1,4) =  4.d0*(-c23+c13+c12-c11)\n\
      sh(1,5) =  2.d0*((c12-c23)*y32+(c11-c13)*y21)\n      sh(1,6) =  \
2.d0*((c23-c12)*x32+(c13-c11)*x21)\n      sh(1,7) =  4.d0*(-c33+c23+c13-c12)\n\
      sh(1,8) =  2.d0*((c12-c23)*y32+(c13-c33)*y13)\n      sh(1,9) =  \
2.d0*((c23-c12)*x32+(c33-c13)*x13)\n      sh(2,2) =  \
c11*y21**2+2.d0*c13*y13*y21+c33*y13**2\n      sh(2,3) =  \
(-c11*x21-c13*x13)*y21+(-c13*x21-c33*x13)*y13\n      sh(2,4) =  \
2.d0*((c12-c11)*y21+(c23-c13)*y13)\n      sh(2,5) =  \
(c12*y21+c23*y13)*y32+c11*y21**2+c13*y13*y21\n      sh(2,6) =  \
(-c12*x32-c11*x21)*y21+(-c23*x32-c13*x21)*y13\n      sh(2,7) =  \
2.d0*((c13-c12)*y21+(c33-c23)*y13)\n      sh(2,8) =  \
(c12*y21+c23*y13)*y32+c13*y13*y21+c33*y13**2\n      sh(2,9) =  \
(-c12*x32-c13*x13)*y21+(-c23*x32-c33*x13)*y13\n      sh(3,3) =  \
c11*x21**2+2.d0*c13*x13*x21+c33*x13**2\n      sh(3,4) =  \
2.d0*((c11-c12)*x21+(c13-c23)*x13)\n      sh(3,5) =  \
(-c12*x21-c23*x13)*y32+(-c11*x21-c13*x13)*y21\n      sh(3,6) =  \
(c12*x21+c23*x13)*x32+c11*x21**2+c13*x13*x21\n      sh(3,7) =  \
2.d0*((c12-c13)*x21+(c23-c33)*x13)\n      sh(3,8) =  \
(-c12*x21-c23*x13)*y32+(-c13*x21-c33*x13)*y13\n      sh(3,9) =  \
(c12*x21+c23*x13)*x32+c13*x13*x21+c33*x13**2\n      sh(4,4) =  \
4.d0*(c22-c12-c12+c11)\n      sh(4,5) =  2.d0*((c22-c12)*y32+(c12-c11)*y21)\n \
     sh(4,6) =  2.d0*((c12-c22)*x32+(c11-c12)*x21)\n      sh(4,7) =  \
4.d0*(c23-c22-c13+c12)\n      sh(4,8) =  2.d0*((c22-c12)*y32+(c23-c13)*y13)\n \
     sh(4,9) =  2.d0*((c12-c22)*x32+(c13-c23)*x13)\n      sh(5,5) =  \
c22*y32**2+2.*c12*y21*y32+c11*y21**2\n      sh(5,6) =  \
(-c22*x32-c12*x21)*y32+(-c12*x32-c11*x21)*y21\n      sh(5,7) =  \
2.d0*((c23-c22)*y32+(c13-c12)*y21)\n      sh(5,8) =  \
c22*y32**2+(c12*y21+c23*y13)*y32+c13*y13*y21\n      sh(5,9) =  \
(-c22*x32-c23*x13)*y32+(-c12*x32-c13*x13)*y21\n      sh(6,6) =  \
c22*x32**2+2.*c12*x21*x32+c11*x21**2\n      sh(6,7) =  \
2.d0*((c22-c23)*x32+(c12-c13)*x21)\n      sh(6,8) =  \
(-c22*x32-c12*x21)*y32+(-c23*x32-c13*x21)*y13\n      sh(6,9) =  \
c22*x32**2+(c12*x21+c23*x13)*x32+c13*x13*x21\n      sh(7,7) =  \
4.d0*(c33-c23-c23+c22)\n      sh(7,8) =  2.d0*((c23-c22)*y32+(c33-c23)*y13)\n \
     sh(7,9) =  2.d0*((c22-c23)*x32+(c23-c33)*x13)\n      sh(8,8) =  \
c22*y32**2+2.*c23*y13*y32+c33*y13**2\n      sh(8,9) =  \
(-c22*x32-c23*x13)*y32+(-c23*x32-c33*x13)*y13\n      sh(9,9) =  \
c22*x32**2+2.*c23*x13*x32+c33*x13**2",
    FontColor->RGBColor[0, 0, 1]],
  "\n*)\n",
  StyleBox["  ",
    FontColor->RGBColor[0, 0, 1]],
  "    sh[[1,1]]=4*(c33-c13-c13+c11);\n      \
sh[[1,2]]=2*((c11-c13)*y21+(c13-c33)*y13);\n      \
sh[[1,3]]=2*((c13-c11)*x21+(c33-c13)*x13);\n      \
sh[[1,4]]=4*(-c23+c13+c12-c11);\n      \
sh[[1,5]]=2*((c12-c23)*y32+(c11-c13)*y21);\n      \
sh[[1,6]]=2*((c23-c12)*x32+(c13-c11)*x21);\n      \
sh[[1,7]]=4*(-c33+c23+c13-c12);\n      \
sh[[1,8]]=2*((c12-c23)*y32+(c13-c33)*y13);\n      \
sh[[1,9]]=2*((c23-c12)*x32+(c33-c13)*x13);\n      \
sh[[2,2]]=c11*y21^2+2*c13*y13*y21+c33*y13^2;\n      \
sh[[2,3]]=(-c11*x21-c13*x13)*y21+(-c13*x21-c33*x13)*y13;\n      \
sh[[2,4]]=2*((c12-c11)*y21+(c23-c13)*y13);\n      \
sh[[2,5]]=(c12*y21+c23*y13)*y32+c11*y21^2+c13*y13*y21;\n      \
sh[[2,6]]=(-c12*x32-c11*x21)*y21+(-c23*x32-c13*x21)*y13;\n      \
sh[[2,7]]=2*((c13-c12)*y21+(c33-c23)*y13);\n      \
sh[[2,8]]=(c12*y21+c23*y13)*y32+c13*y13*y21+c33*y13^2;\n      \
sh[[2,9]]=(-c12*x32-c13*x13)*y21+(-c23*x32-c33*x13)*y13;\n      \
sh[[3,3]]=c11*x21^2+2*c13*x13*x21+c33*x13^2;\n      \
sh[[3,4]]=2*((c11-c12)*x21+(c13-c23)*x13);\n      \
sh[[3,5]]=(-c12*x21-c23*x13)*y32+(-c11*x21-c13*x13)*y21;\n      \
sh[[3,6]]=(c12*x21+c23*x13)*x32+c11*x21^2+c13*x13*x21;\n      \
sh[[3,7]]=2*((c12-c13)*x21+(c23-c33)*x13);\n      \
sh[[3,8]]=(-c12*x21-c23*x13)*y32+(-c13*x21-c33*x13)*y13;\n      \
sh[[3,9]]=(c12*x21+c23*x13)*x32+c13*x13*x21+c33*x13^2;\n      \
sh[[4,4]]=4*(c22-c12-c12+c11);\n      \
sh[[4,5]]=2*((c22-c12)*y32+(c12-c11)*y21);\n      \
sh[[4,6]]=2*((c12-c22)*x32+(c11-c12)*x21);\n      \
sh[[4,7]]=4*(c23-c22-c13+c12);\n      \
sh[[4,8]]=2*((c22-c12)*y32+(c23-c13)*y13);\n      \
sh[[4,9]]=2*((c12-c22)*x32+(c13-c23)*x13);\n      \
sh[[5,5]]=c22*y32^2+2*c12*y21*y32+c11*y21^2;\n      \
sh[[5,6]]=(-c22*x32-c12*x21)*y32+(-c12*x32-c11*x21)*y21;\n      \
sh[[5,7]]=2*((c23-c22)*y32+(c13-c12)*y21);\n      \
sh[[5,8]]=c22*y32^2+(c12*y21+c23*y13)*y32+c13*y13*y21;\n      \
sh[[5,9]]=(-c22*x32-c23*x13)*y32+(-c12*x32-c13*x13)*y21;\n      \
sh[[6,6]]=c22*x32^2+2*c12*x21*x32+c11*x21^2;\n      \
sh[[6,7]]=2*((c22-c23)*x32+(c12-c13)*x21);\n      \
sh[[6,8]]=(-c22*x32-c12*x21)*y32+(-c23*x32-c13*x21)*y13;\n      \
sh[[6,9]]=c22*x32^2+(c12*x21+c23*x13)*x32+c13*x13*x21;\n      \
sh[[7,7]]=4*(c33-c23-c23+c22);\n      \
sh[[7,8]]=2*((c23-c22)*y32+(c33-c23)*y13);\n      \
sh[[7,9]]=2*((c22-c23)*x32+(c23-c33)*x13);\n      \
sh[[8,8]]=c22*y32^2+2*c23*y13*y32+c33*y13^2;\n      \
sh[[8,9]]=(-c22*x32-c23*x13)*y32+(-c23*x32-c33*x13)*y13;\n      \
sh[[9,9]]=c22*x32^2+2*c23*x13*x32+c33*x13^2;\n(*",
  StyleBox["\n      do 3500  i = 1,9\n        k =   ls(i)\n        do 3400  j \
= i,9\n          l =   ls(j)\n            sm(k,l) =   sm(k,l) + sh(i,j)\n     \
       sm(l,k) =   sm(k,l)\n 3400      continue\n 3500    continue\n      \
return\n      end\nC=END FORTRAN\n",
    FontColor->RGBColor[0, 0, 1]],
  "*) \n      For [i=1,i<=9,i++, k=ls[[i]]; For [j=i,j<=9,j++,\n           \
l=ls[[j]]; sm[[l,k]]=sm[[k,l]]+=sh[[i,j]] ]];\n      Return[sm];\n ];\n \n",
  StyleBox[" \nClearAll[Em]; Em=1;\n\
sm=SM3SHBENDH[{0,1,1/2},{0,0,1},{{Em,0,0},\n  \
{0,Em,0},{0,0,Em/2}},1,{1,2,3,4,5,6,7,8,9},\n  Table[0,{9},{9}]];\nPrint[sm];\
\nPrint[Chop[Eigenvalues[N[sm]]]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({{2067\/400, 
        213\/200, \(-\(1641\/800\)\), \(-\(1363\/400\)\), \(-\(37\/200\)\), \
\(-\(1289\/800\)\), \(-\(44\/25\)\), 22\/25, \(-\(5\/8\)\)}, {213\/200, 
        57\/100, \(-\(99\/400\)\), \(-\(37\/200\)\), \(-\(13\/100\)\), \
\(-\(11\/400\)\), \(-\(22\/25\)\), 
        11\/25, \(-\(7\/20\)\)}, {\(-\(1641\/800\)\), \(-\(99\/400\)\), 
        1443\/1600, 1289\/800, 11\/400, 1267\/1600, 11\/25, \(-\(11\/50\)\), 
        11\/80}, {\(-\(1363\/400\)\), \(-\(37\/200\)\), 1289\/800, 2067\/400, 
        213\/200, 1641\/800, \(-\(44\/25\)\), 22\/25, 
        5\/8}, {\(-\(37\/200\)\), \(-\(13\/100\)\), 11\/400, 213\/200, 
        57\/100, 99\/400, \(-\(22\/25\)\), 11\/25, 
        7\/20}, {\(-\(1289\/800\)\), \(-\(11\/400\)\), 1267\/1600, 1641\/800, 
        99\/400, 1443\/1600, \(-\(11\/25\)\), 11\/50, 
        11\/80}, {\(-\(44\/25\)\), \(-\(22\/25\)\), 
        11\/25, \(-\(44\/25\)\), \(-\(22\/25\)\), \(-\(11\/25\)\), 
        88\/25, \(-\(44\/25\)\), 0}, {22\/25, 11\/25, \(-\(11\/50\)\), 
        22\/25, 11\/25, 11\/50, \(-\(44\/25\)\), 22\/25, 
        0}, {\(-\(5\/8\)\), \(-\(7\/20\)\), 11\/80, 5\/8, 7\/20, 11\/80, 0, 
        0, 7\/20}}\)], "Print"],

Cell[BoxData[
    \({10.412639708663404`, 6.709999999999994`, 0.9061102913365954`, 0, 0, 0, 
      0, 0, 0}\)], "Print"]
}, Open  ]],

Cell[CellGroupData[{

Cell[TextData[{
  "\n",
  StyleBox["(*\nC=DECK SM3SHELL SM3SHELL FORTRAN\nC=PURPOSE Form  stiffness \
of 3-node, 18-dof flat triangle shell\nC=AUTHOR C. A. Felippa, October 1996\n\
C=VERSION October 1996\nC=EQUIPMENT Machine independent\nC=KEYWORDS Kirchhoff \
thin shell isotropic homogeneous\nC=KEYWORDS finite element triangle \
stiffness matrix\nC=BLOCK ABSTRACT\nC\nC     SM3SHBENDB forms the material \
element stiffness matrix\nC     of a 3-node, 18-dof triangular shell \
constructed with the \nC     generalized FF/ANS formulation.  Isotropic \
material,\nC     homogeneous wall construction.\nC\nC=END ABSTRACT\nC=BLOCK \
USAGE\nC\nC     The calling sequence is\nC\nC     call SM3SHELL (x, y, z, em, \
nu, h, rho, sm, mm, m, status)\nC\nC     The input arguments are:\nC\nC       \
X         (3 x 1) array of x coordinates of triangle nodes\nC       Y         \
(3 x 1) array of y coordinates of triangle nodes\nC       Z         (3 x 1) \
array of z coordinates of triangle nodes\nC       EM        Elastic modulus\n\
C       NU        Poisson's ratio\nC       H         Corner thicknesses\nC    \
   RHO       Dummy\nC       M         First dimension of SM in calling \
program.\nC\nC     The outputs are:\nC\nC       SM        Output stiffness \
array with basic stiffness\nC                 coefficients added in.  the \
(I,J)-th entry of the\nC                 (9 x 9) element bending stiffness is \
added to\nC                 SM(K,L), where K=LS(I) and L=LS(J).\nC       MM   \
     Mass matrix - a dummy for now\nC       STATUS    Status character \
variable. Blank if no error detected.\nC\nC=END USAGE\n \n      subroutine \
SM3SHELL\n     &      (xg, yg, zg, em, nu, h, rho, esm, mm, m, status)\n\nC\n\
C                   A R G U M E N T S\nC\n      character*( * )     status\n  \
    integer           m\n      double precision  xg(3), yg(3), zg(3) \n      \
double precision  em, nu, rho, h(3)\n      double precision  esm(m,m), \
mm(m,m)\nC\nC                   L O C A L   V A R I A B L E S\nC\n      \
double precision  sm(18,18), db(3,3), dm(3,3)\n      double precision  \
xlp(3), ylp(3), zlp(3), dcm(3,3)\n      double precision  area, h0, clr, cqr, \
betab\n      integer           lb(9), le(9)\n      integer           i, j, nd\
\nC\nC                   D A T A\nC\n      data       lb    /3,4,5, 9,10,11, \
15,16,17/\n      data       le    /1,2, 7,8, 13,14, 6,12,18/",
    FontColor->RGBColor[0, 0, 1]],
  "\n",
  StyleBox["\nC\nC                   L O G I C\nC\n      status =  ' '\nC\nC  \
             Establish local system\nC\n      call SM3SHLOCALSYS (xg,yg,zg, \
xlp,ylp,zlp, dcm,area, status)\n      if (status(1:1) .ne. ' ')       return\n\
      if (area .le. 0.0)      then\n        status = 'SM3SHELLST: Negative \
area'\n        if (area .eq. 0.0)   status = 'SM3SHELLST: Zero area'\n        \
return\n      end if\nC\nC               Form constitutive matrices\nC\n      \
h0 =    (h(1) + h(2) + h(3))/3.d0\n      call   SM3SHISODM (em, nu, h0, dm)\n \
     call   SM3SHISODB (em, nu, h0, db)\nC\nC               Form local \
stiffnesses\nC\n      do 1800  i = 1,18\n        do 1800  j = 1,18\n          \
sm(i,j)  = 0.0\n 1800   continue\n      betab =  max(0.01,0.5*(1.-4.*nu**2))\n\
      clr =    0.0d0\n      cqr =    1.0d0\n      nd =     18\n      call   \
SM3SHMEMBB (xlp, ylp, dm, 1.5d0,1.0d0, le, sm, nd, status)\n      call   \
SM3SHMEMBH (xlp, ylp, dm, betab, le, sm, nd, status)\n      call   SM3SHBENDB \
(xlp, ylp, db,1.0d0,clr,cqr, lb, sm, nd, status)\n      call   SM3SHBENDH \
(xlp, ylp, db,1.0d+0, lb, sm, nd, status)\nC\nC              Transform to \
global coordinates\nC\n      call   SM3SHTRANSFORM (sm, \
dcm,dcm,dcm,dcm,dcm,dcm)\n      do 3000  i = 1,18\n        do 2500  j = 1,18\n\
          esm(i,j)  = sm(i,j)\n 2500   continue\n 3000 continue\n      return\
\n      end\nC=END FORTRAN\n*)\n\n",
    FontColor->RGBColor[0, 0, 1]],
  "SM3SHELL[xg_,yg_,zg_,em_,nu_,h_,rho_]:=Module[{\n  db,dm,xlp,ylp,zlp,dcm, \
area,h0,betab,lb,le,i,j, esm},\n      lb={3,4,5, 9,10,11, 15,16,17};\n      \
le={1,2, 7,8, 13,14, 6,12,18};\n      \
{xlp,ylp,zlp,dcm,area}=SM3SHLOCALSYS[xg,yg,zg];\n      If \
[Length[h]>0,h0=(h[[1]]+h[[2]]+h[[3]])/3,h0=h];\n      \
dm=SM3SHISODM[em,nu,h0]; \n      db=SM3SHISODB[em,nu,h0];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      esm=Table[0,{18},{18}];\n      betab=Max[1/100,(1-4*nu^2)/2]; \
alphab=3/2;\n      esm=SM3SHMEMBB[xlp,ylp, dm, alphab, 1, le, esm];\n      \
esm=SM3SHMEMBH[xlp,ylp, dm, betab,  le, esm];\n      esm=SM3SHBENDB[xlp,ylp, \
db, 1, 0, 1, lb,esm];\n      esm=SM3SHBENDH[xlp,ylp, db, 1, lb, esm];  \n     \
 esm=SM3SHTRANSFORM[esm, dcm,dcm,dcm,dcm,dcm,dcm];  \n      Return[esm];\n    \
  ];\n",
  StyleBox["\n\nClearAll[Em,nu]; Em=1; nu=0;\n\
sm=SM3SHELL[{0,1,1/2},{0,0,1},{0,0,0},Em,nu,\n  {1,1,1},rho];\nPrint[sm];\n\
Print[\"eigs=\",Chop[Eigenvalues[N[sm]]]];",
    FontColor->RGBColor[1, 0, 0]],
  "\n"
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[BoxData[
    \({{913\/1536, 47\/768, 0, 0, 
        0, \(-\(141\/1024\)\), \(-\(623\/1536\)\), \(-\(47\/768\)\), 0, 0, 0, 
        115\/1024, \(-\(145\/768\)\), 0, 0, 0, 
        0, \(-\(157\/1536\)\)}, {47\/768, 193\/384, 0, 0, 0, 101\/512, 
        47\/768, \(-\(97\/384\)\), 0, 0, 0, 
        53\/512, \(-\(47\/384\)\), \(-\(1\/4\)\), 0, 0, 
        0, \(-\(35\/768\)\)}, {0, 0, 3347\/4800, 
        253\/2400, \(-\(1961\/9600\)\), 0, 0, 0, \(-\(2131\/4800\)\), 
        131\/2400, \(-\(1993\/9600\)\), 0, 0, 0, \(-\(19\/75\)\), 
        7\/75, \(-\(127\/800\)\), 0}, {0, 0, 253\/2400, 
        73\/800, \(-\(13\/1600\)\), 0, 0, 0, 131\/2400, 
        21\/800, \(-\(7\/4800\)\), 0, 0, 0, \(-\(4\/25\)\), 
        17\/400, \(-\(19\/1200\)\), 0}, {0, 
        0, \(-\(1961\/9600\)\), \(-\(13\/1600\)\), 2243\/19200, 0, 0, 0, 
        1993\/9600, 7\/4800, 979\/19200, 0, 0, 0, \(-\(1\/300\)\), 1\/100, 
        61\/1600, 0}, {\(-\(141\/1024\)\), 101\/512, 0, 0, 0, 35\/256, 
        115\/1024, \(-\(53\/512\)\), 0, 0, 0, 61\/3072, 
        13\/512, \(-\(3\/32\)\), 0, 0, 0, 59\/3072}, {\(-\(623\/1536\)\), 
        47\/768, 0, 0, 0, 115\/1024, 913\/1536, \(-\(47\/768\)\), 0, 0, 
        0, \(-\(141\/1024\)\), \(-\(145\/768\)\), 0, 0, 0, 
        0, \(-\(157\/1536\)\)}, {\(-\(47\/768\)\), \(-\(97\/384\)\), 0, 0, 
        0, \(-\(53\/512\)\), \(-\(47\/768\)\), 193\/384, 0, 0, 
        0, \(-\(101\/512\)\), 47\/384, \(-\(1\/4\)\), 0, 0, 0, 35\/768}, {0, 
        0, \(-\(2131\/4800\)\), 131\/2400, 1993\/9600, 0, 0, 0, 3347\/4800, 
        253\/2400, 1961\/9600, 0, 0, 0, \(-\(19\/75\)\), 7\/75, 127\/800, 
        0}, {0, 0, 131\/2400, 21\/800, 7\/4800, 0, 0, 0, 253\/2400, 73\/800, 
        13\/1600, 0, 0, 0, \(-\(4\/25\)\), 17\/400, 19\/1200, 0}, {0, 
        0, \(-\(1993\/9600\)\), \(-\(7\/4800\)\), 979\/19200, 0, 0, 0, 
        1961\/9600, 13\/1600, 2243\/19200, 0, 0, 0, 1\/300, \(-\(1\/100\)\), 
        61\/1600, 0}, {115\/1024, 53\/512, 0, 0, 0, 
        61\/3072, \(-\(141\/1024\)\), \(-\(101\/512\)\), 0, 0, 0, 35\/256, 
        13\/512, 3\/32, 0, 0, 0, 
        59\/3072}, {\(-\(145\/768\)\), \(-\(47\/384\)\), 0, 0, 0, 
        13\/512, \(-\(145\/768\)\), 47\/384, 0, 0, 0, 13\/512, 145\/384, 0, 
        0, 0, 0, 157\/768}, {0, \(-\(1\/4\)\), 0, 0, 0, \(-\(3\/32\)\), 
        0, \(-\(1\/4\)\), 0, 0, 0, 3\/32, 0, 1\/2, 0, 0, 0, 0}, {0, 
        0, \(-\(19\/75\)\), \(-\(4\/25\)\), \(-\(1\/300\)\), 0, 0, 
        0, \(-\(19\/75\)\), \(-\(4\/25\)\), 1\/300, 0, 0, 0, 
        38\/75, \(-\(14\/75\)\), 0, 0}, {0, 0, 7\/75, 17\/400, 1\/100, 0, 0, 
        0, 7\/75, 17\/400, \(-\(1\/100\)\), 0, 0, 0, \(-\(14\/75\)\), 
        61\/600, 0, 0}, {0, 0, \(-\(127\/800\)\), \(-\(19\/1200\)\), 
        61\/1600, 0, 0, 0, 127\/800, 19\/1200, 61\/1600, 0, 0, 0, 0, 0, 
        33\/400, 0}, {\(-\(157\/1536\)\), \(-\(35\/768\)\), 0, 0, 0, 
        59\/3072, \(-\(157\/1536\)\), 35\/768, 0, 0, 0, 59\/3072, 157\/768, 
        0, 0, 0, 0, 185\/1536}}\)], "Print"],

Cell[BoxData[
    InterpretationBox[\("eigs="\[InvisibleSpace]{1.3364927241045919`, 
          1.0716114940520929`, 0.9940655723461995`, 0.9299320278721559`, 
          0.7779726074441162`, 0.5930639792253276`, 0.07551920531784548`, 
          0.07480219689406195`, 0.04455057057756234`, 0.04026577523378098`, 
          0.01727787767045752`, 0.011503260928470871`, 0, 0, 0, 0, 0, 0}\),
      SequenceForm[ 
      "eigs=", {1.3364927241045919, 1.0716114940520929, 0.99406557234619952, 
        0.9299320278721559, 0.77797260744411623, 0.59306397922532761, 
        0.075519205317845478, 0.074802196894061951, 0.044550570577562343, 
        0.040265775233780977, 0.01727787767045752, 0.011503260928470871, 0, 0,
         0, 0, 0, 0}],
      Editable->False]], "Print"]
}, Open  ]],

Cell["\<\

This is a companion test for the MEMBRANE component of
the shell triangle.   Please check that your codes agree
with these.  If not, please contact me.   Next I will
be checking the complete shell program.

The test triangle is the same as for bending stiffness:

                3(x3,y3)
                  o  
                  + +   
                  +   +
                  +     +  
                  +       + 
                  +         +  
                  +           +
                  o+++++++++++++o
              1(x1,y1)         2(x2,y2)
 
 x =   0, 2, 0          y =   0, 0, 2
 Em =  4000, nu = 1/3,  t = 0.2

The degrees of freedom are arranged u1,v1,thetaz1, u2, ... thetaz3

 Constitutive matrix:
          1       2       3
    1  900.00  300.00    0.00
    2  300.00  900.00    0.00
    3    0.00    0.00  300.00
 Matrix P:
          1       2       3
    1   -2.00    0.00   -2.00
    2    0.00   -2.00   -2.00
    3    2.00    0.00    0.00
    4    0.00    0.00    2.00
    5    0.00    0.00    2.00
    6    0.00    2.00    0.00
    7    1.00   -1.00    0.00
    8   -1.00    0.00   -2.00
    9    0.00    1.00    2.00

 Basic membrane stiffness:
          1       2       3       4       5       6       7       8       9
    1  600.00  300.00 -150.00 -450.00 -150.00  375.00 -150.00 -150.00 -225.00
    2  300.00  600.00  150.00 -150.00 -150.00  225.00 -150.00 -450.00 -375.00
    3 -150.00  150.00  150.00  150.00    0.00  -75.00    0.00 -150.00  -75.00
    4 -450.00 -150.00  150.00  450.00    0.00 -225.00    0.00  150.00   75.00
    5 -150.00 -150.00    0.00    0.00  150.00 -150.00  150.00    0.00  150.00
    6  375.00  225.00  -75.00 -225.00 -150.00  262.50 -150.00  -75.00 -187.50
    7 -150.00 -150.00    0.00    0.00  150.00 -150.00  150.00    0.00  150.00
    8 -150.00 -450.00 -150.00  150.00    0.00  -75.00    0.00  450.00  225.00
    9 -225.00 -375.00  -75.00   75.00  150.00 -187.50  150.00  225.00  262.50

 HO membrane stiffness:
          1       2       3       4       5       6       7       8       9
    1   62.50  -62.50 -118.06    0.00   62.50  -65.97  -62.50    0.00  -65.97
    2  -62.50   62.50  118.06    0.00  -62.50   65.97   62.50    0.00   65.97
    3 -118.06  118.06  256.94    0.00 -118.06  107.64  118.06    0.00  107.64
    4    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00
    5   62.50  -62.50 -118.06    0.00   62.50  -65.97  -62.50    0.00  -65.97
    6  -65.97   65.97  107.64    0.00  -65.97   90.28   65.97    0.00   65.97
    7  -62.50   62.50  118.06    0.00  -62.50   65.97   62.50    0.00   65.97
    8    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00
    9  -65.97   65.97  107.64    0.00  -65.97   65.97   65.97    0.00   90.28
    
 Total membrane stiffness:
          1       2       3       4       5       6       7       8       9
    1  662.50  237.50 -268.06 -450.00  -87.50  309.03 -212.50 -150.00 -290.97
    2  237.50  662.50  268.06 -150.00 -212.50  290.97  -87.50 -450.00 -309.03
    3 -268.06  268.06  406.94  150.00 -118.06   32.64  118.06 -150.00   32.64
    4 -450.00 -150.00  150.00  450.00    0.00 -225.00    0.00  150.00   75.00
    5  -87.50 -212.50 -118.06    0.00  212.50 -215.97   87.50    0.00   84.03
    6  309.03  290.97   32.64 -225.00 -215.97  352.78  -84.03  -75.00 -121.53
    7 -212.50  -87.50  118.06    0.00   87.50  -84.03  212.50    0.00  215.97
    8 -150.00 -450.00 -150.00  150.00    0.00  -75.00    0.00  450.00  225.00
    9 -290.97 -309.03   32.64   75.00   84.03 -121.53  215.97  225.00  352.78
    
 Eigenvalues
          1       2       3       4       5       6       7       8       9
    1 1805.46 1025.31  452.79  426.49   36.39   16.06    0.00    0.00    0.00
\
\>", "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  FontSize->10],

Cell[CellGroupData[{

Cell[TextData[{
  "x={0,2,0}; y={0,0,2};\nle={1,2,7, 3,4,8, 5,6,9};\nem=4000; nu=1/3; h0=0.2;\
\nesm=Table[0,{9},{9}];\ndm=SM3SHISODM[em,nu,h0]; ",
  StyleBox["Print[\"dm=\",dm];",
    FontColor->RGBColor[1, 0, 0]],
  "\nbetab=Max[1/100,(1/2)*(1-4*nu^2)];\nesm=SM3SHMEMBB[x,y, dm, 3/2, 1, le, \
esm];\n",
  StyleBox["Print[\" after MEMB esm=\",esm]; \n\
Print[Chop[Eigenvalues[N[esm]]]];",
    FontColor->RGBColor[1, 0, 0]],
  "\nesm=SM3SHMEMBH[x,y, dm, betab,  le, esm];\n",
  StyleBox["Print[\" after MEMH esm=\",esm];  \n\
Print[Chop[Eigenvalues[N[esm]]]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0.720928, 0.889998, 0.213611]],

Cell["\<\
dm={{900., 300., 0}, {300., 900., 0}, {0, 0, 300.}}
p={{-2, 0, -2}, {0, -2, -2}, {2, 0, 0}, {0, 0, 2}, {0, 0, 2}, 
 
   {0, 2, 0}, {1, -1, 0}, {-1, 0, -2}, {0, 1, 2}}
 after MEMB esm={{600., 300., -150., -150., -150., 375., -450., 
 
    -150., -225.}, {300., 600., -150., -150., 150., 225., -150., 
 
    -450., -375.}, {-150., -150., 150., 150., 0, -150., 0, 0, 150.}, 
 
   {-150., -150., 150., 150., 0, -150., 0, 0, 150.}, 
 
   {-150., 150., 0, 0, 150., -75., 150., -150., -75.}, 
 
   {375., 225., -150., -150., -75., 262.5, -225., -75., -187.5}, 
 
   {-450., -150., 0, 0, 150., -225., 450., 150., 75.}, 
 
   {-150., -450., 0, 0, -150., -75., 150., 450., 225.}, 
 
   {-225., -375., 150., 150., -75., -187.5, 75., 225., 262.5}}
{1800., 825., 450., 0, 0, 0, 0, 0, 0}
area=2
 after MEMH esm={{662.5, 237.5, -87.5, -212.5, -268.056, 309.028, 
 
    -450., -150., -290.972}, 
 
   {237.5, 662.5, -212.5, -87.5, 268.056, 290.972, -150., -450., 
 
    -309.028}, {-87.5, -212.5, 212.5, 87.5, -118.056, -215.972, 0, 
 
    0, 84.0278}, {-212.5, -87.5, 87.5, 212.5, 118.056, -84.0278, 0, 
 
    0, 215.972}, {-268.056, 268.056, -118.056, 118.056, 406.944, 
 
    32.6389, 150., -150., 32.6389}, 
 
   {309.028, 290.972, -215.972, -84.0278, 32.6389, 352.778, -225., 
 
    -75., -121.528}, {-450., -150., 0, 0, 150., -225., 450., 150., 
 
    75.}, {-150., -450., 0, 0, -150., -75., 150., 450., 225.}, 
 
   {-290.972, -309.028, 84.0278, 215.972, 32.6389, -121.528, 75., 
 
    225., 352.778}}
{1805.46, 1025.31, 452.789, 426.494, 36.3933, 16.0552, 0, 0, 0}\
\>", "Print",\

  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  FontSize->10,
  Background->RGBColor[1, 1, 0]]
}, Open  ]],

Cell["\<\

             
Re: my previous e-mail.

The test triangle for basic bending AQR stiffness is defined by

                3(x3,y3)
                  o  
                  + +   
                  +   +
                  +     +  
                  +       + 
                  +         +  
                  +           +
                  o+++++++++++++o
              1(x1,y1)         2(x2,y2)
 
 x =   0, 2, 0          y =   0, 0, 2
 Em =  4000, nu = 1/3,  t = 0.2

These are selected to get \"nice\" numbers in the stiffnesses listed below.
Arrangement of DOFs is:  w1,thetax1,thetay1,w2, ... thetay3.

Constitutive matrix:

       3.000       1.000       0.000
       1.000       3.000       0.000
       0.000       0.000       1.000

Incorrect AQR Kb:

  2.0000  0.0000  0.0000 -1.0000  1.0000 -1.0000 -0.5000  1.0000 -1.0000
  0.0000  1.5000 -0.5000  0.5000  0.5000  1.0000  0.2500 -1.0000  0.5000
  0.0000 -0.5000  1.5000  0.5000  0.5000 -1.0000 -0.7500  1.0000  0.5000
 -1.0000  0.5000  0.5000  1.0000  0.0000  0.5000  0.0000 -0.5000  1.0000
  1.0000  0.5000  0.5000  0.0000  1.0000 -0.5000 -0.5000  0.5000  0.0000
 -1.0000  1.0000 -1.0000  0.5000 -0.5000  1.5000  0.7500 -1.5000  0.5000
 -0.5000  0.2500 -0.7500  0.0000 -0.5000  0.7500  0.5000 -0.7500  0.0000
  1.0000 -1.0000  1.0000 -0.5000  0.5000 -1.5000 -0.7500  1.5000 -0.5000
 -1.0000  0.5000  0.5000  1.0000  0.0000  0.5000  0.0000 -0.5000  1.0000

Correct AQR Kb:

  2.0000  0.0000  0.0000 -1.0000  1.0000 -1.0000 -1.0000  1.0000 -1.0000
  0.0000  1.5000 -0.5000  0.5000  0.5000  1.0000 -0.5000 -1.0000  0.5000
  0.0000 -0.5000  1.5000  0.5000  0.5000 -1.0000 -0.5000  1.0000  0.5000
 -1.0000  0.5000  0.5000  1.0000  0.0000  0.5000  0.0000 -0.5000  1.0000
  1.0000  0.5000  0.5000  0.0000  1.0000 -0.5000 -1.0000  0.5000  0.0000
 -1.0000  1.0000 -1.0000  0.5000 -0.5000  1.5000  0.5000 -1.5000  0.5000
 -1.0000 -0.5000 -0.5000  0.0000 -1.0000  0.5000  1.0000 -0.5000  0.0000
  1.0000 -1.0000  1.0000 -0.5000  0.5000 -1.5000 -0.5000  1.5000 -0.5000
 -1.0000  0.5000  0.5000  1.0000  0.0000  0.5000  0.0000 -0.5000  1.0000

Difference correct Kb - incorrect Kb:

  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -0.5000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -0.7500  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2500  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -0.5000  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -0.2500  0.0000  0.0000
 -0.5000 -0.7500  0.2500  0.0000 -0.5000 -0.2500  0.5000  0.2500  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.2500  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Incorrect total AQR stiffness Kb + Kh:

  8.0000  3.0000 -3.0000 -4.0000  1.0000 -4.0000 -3.5000  4.0000 -1.0000
  3.0000  5.5000  0.5000  0.7500 -0.2500  1.2500 -3.0000  2.2500 -0.2500
 -3.0000  0.5000  5.5000  3.7500 -0.2500  2.2500 -1.0000  1.2500 -0.2500
 -4.0000  0.7500  3.7500  4.2500  0.0000  3.7500 -0.2500 -0.2500  1.0000
  1.0000 -0.2500 -0.2500  0.0000  1.7500 -0.5000 -0.5000  0.5000  0.7500
 -4.0000  1.2500  2.2500  3.7500 -0.5000  4.7500  0.5000 -1.2500  0.5000
 -3.5000 -3.0000 -1.0000 -0.2500 -0.5000  0.5000  3.7500 -4.0000  0.0000
  4.0000  2.2500  1.2500 -0.2500  0.5000 -1.2500 -4.0000  4.7500 -0.5000
 -1.0000 -0.2500 -0.2500  1.0000  0.7500  0.5000  0.0000 -0.5000  1.7500
eigs of above:
 18.2642 12.2463  4.5545  2.8897  1.0941  0.9511  0.0000  0.0000  0.0000

Correct total AQR stiffness Kb + Kh:

  8.0000  3.0000 -3.0000 -4.0000  1.0000 -4.0000 -4.0000  4.0000 -1.0000
  3.0000  5.5000  0.5000  0.7500 -0.2500  1.2500 -3.7500  2.2500 -0.2500
 -3.0000  0.5000  5.5000  3.7500 -0.2500  2.2500 -0.7500  1.2500 -0.2500
 -4.0000  0.7500  3.7500  4.2500  0.0000  3.7500 -0.2500 -0.2500  1.0000
  1.0000 -0.2500 -0.2500  0.0000  1.7500 -0.5000 -1.0000  0.5000  0.7500
 -4.0000  1.2500  2.2500  3.7500 -0.5000  4.7500  0.2500 -1.2500  0.5000
 -4.0000 -3.7500 -0.7500 -0.2500 -1.0000  0.2500  4.2500 -3.7500  0.0000
  4.0000  2.2500  1.2500 -0.2500  0.5000 -1.2500 -3.7500  4.7500 -0.5000
 -1.0000 -0.2500 -0.2500  1.0000  0.7500  0.5000  0.0000 -0.5000  1.7500
eigs of above:
 18.5726 12.5217  4.4925  3.0000  0.9783  0.9348  0.0000  0.0000  0.0000
 
The higher order stiffness matrix Kh is that produced by the ANDES
formulation and determined to be optimal.  The scaling coefficient
gamma in  Kb + gamma*Kh is taken as one.

Note that the incorrectness will not show up in the eigenvalue test:
both matrices have the correct rank.
\
\>", "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  FontSize->10],

Cell[CellGroupData[{

Cell[TextData[{
  "x={0,2,0}; y={0,0,2};\nlb={1,2,3, 4,5,6, 7,8,9};\nem=4000; nu=1/3; h0=0.2;\
\nesm=Table[0,{9},{9}];\ndb=SM3SHISODB[em,nu,h0]; ",
  StyleBox["Print[\"db=\",db];",
    FontColor->RGBColor[1, 0, 0]],
  "\n",
  StyleBox["Print[Chop[Eigenvalues[N[esm]]]];",
    FontColor->RGBColor[1, 0, 0]],
  "\nesm=SM3SHBENDB[x,y, db, 1, 0, 1, lb,esm]; ",
  StyleBox["esm=Chop[esm];",
    FontColor->RGBColor[1, 0, 0]],
  "\n",
  StyleBox["Print[\" after BENB esm=\",esm//MatrixForm];  \n\
Print[Chop[Eigenvalues[N[esm]]]];",
    FontColor->RGBColor[1, 0, 0]],
  "\nesm=SM3SHBENDH[x,y, db, 1, lb, esm];  \n",
  StyleBox["esm=Chop[esm]; \nPrint[\" after BENH esm=\",esm//MatrixForm];\n\
Print[Chop[Eigenvalues[N[esm]]]];",
    FontColor->RGBColor[1, 0, 0]]
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0.720928, 0.889998, 0.213611]],

Cell["\<\
db={{3., 1., 0}, {1., 3., 0}, {0, 0, 1.}}
{0, 0, 0, 0, 0, 0, 0, 0, 0}
 after BENB esm=2.     0      0      -1.    1.     -1.    -1.    1.     -1.

                0      1.5    -0.5   0.5    0.5    1.     -0.5   -1.    0.5

                0      -0.5   1.5    0.5    0.5    -1.    -0.5   1.     0.5

                -1.    0.5    0.5    1.     0      0.5    0      -0.5   1.

                1.     0.5    0.5    0      1.     -0.5   -1.    0.5    0

                -1.    1.     -1.    0.5    -0.5   1.5    0.5    -1.5   0.5

                -1.    -0.5   -0.5   0      -1.    0.5    1.     -0.5   0

                1.     -1.    1.     -0.5   0.5    -1.5   -0.5   1.5    -0.5

                -1.    0.5    0.5    1.     0      0.5    0      -0.5   1.
{6., 3., 3., 0, 0, 0, 0, 0, 0}
 after BENH esm=8.      3.      -3.     -4.     1.      -4.     -4.     4.    \
  -1.

                3.      5.5     0.5     0.75    -0.25   1.25    -3.75   2.25  \
  -0.25

                -3.     0.5     5.5     3.75    -0.25   2.25    -0.75   1.25  \
  -0.25

                -4.     0.75    3.75    4.25    0       3.75    -0.25   -0.25 \
  1.

                1.      -0.25   -0.25   0       1.75    -0.5    -1.     0.5   \
  0.75

                -4.     1.25    2.25    3.75    -0.5    4.75    0.25    -1.25 \
  0.5

                -4.     -3.75   -0.75   -0.25   -1.     0.25    4.25    -3.75 \
  0

                4.      2.25    1.25    -0.25   0.5     -1.25   -3.75   4.75  \
  -0.5

                -1.     -0.25   -0.25   1.      0.75    0.5     0       -0.5  \
  1.75
{18.5726, 12.5217, 4.49253, 3., 0.978302, 0.934825, 0, 0, 0}\
\>", "Print",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  FontSize->10,
  Background->RGBColor[1, 1, 0]]
}, Open  ]],

Cell[TextData[{
  "MasterStiffnessShell[nodcoor_,elenod_,elemat_,elefab_,eleopt_]:=\n  \
Module[ {numele=Length[elenod],numnod=Length[nodcoor],\n  \
xg,yg,zb,Em,nu,h,rho=0,num,\n  \
e,eNL,eftab,ni,nj,i,j,ncoor,mprop,fprop,opt,Ke,K},\n  \
K=Table[0,{6*numnod},{6*numnod}]; \n  For [e=1, e<=numele, e++,\n     \
eNL=elenod[[e]]; {ni,nj,nk}=eNL; \n     \
eftab=Flatten[{Table[6*(ni-1)+i,{i,1,6}],\n                    \
Table[6*(nj-1)+i,{i,1,6}],\n                    Table[6*(nk-1)+i,{i,1,6}]}];\n\
     xg={nodcoor[[ni,1]],nodcoor[[nj,1]],nodcoor[[nk,1]]};\n     \
yg={nodcoor[[ni,2]],nodcoor[[nj,2]],nodcoor[[nk,2]]};  \n     \
zg={nodcoor[[ni,3]],nodcoor[[nj,3]],nodcoor[[nk,3]]};          \n     \
Em=elemat[[e,1]]; nu=elemat[[e,2]]; h=elefab[[e]];\n     num=eleopt;\n     If \
[num,{xg,yg,zg,Em,nu,rho}=N[{xg,yg,zg,Em,nu,rho}]];\n    ",
  StyleBox[" (*Print[\"xg=\",xg]; Print[\"yg=\",yg]; Print[\"zg=\",zg];\n     \
Print[\"eftab=\",eftab]; Print[\"Em,nu,h=\",{Em,nu,h}];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n     Ke=SM3SHELL[xg,yg,zg,Em,nu,h,rho];  Ke=Chop[Ke];\n    ",
  StyleBox[" (*Print[\"Eigs of Ke=\",Chop[Eigenvalues[N[Ke]]]];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n     neldof=Length[Ke];\n     For [i=1, i<=neldof, i++, ii=eftab[[i]];\n \
        For [j=i, j<=neldof, j++, jj=eftab[[j]];\n             \
K[[jj,ii]]=K[[ii,jj]]+=Ke[[i,j]] \n         ];\n       ];\n     ];\n  \
Return[Chop[K]]\n  ];\n  \nModifyMasterStiffForDBC[pdof_,K_] := Module[\n  \
{i,j,k,n=Length[K],np=Length[pdof],Kii,Kmod}, Kmod=K; \n  For [k=1,k<=np,k++, \
i=pdof[[k]]; Kii=Kmod[[i,i]]; \n      For [j=1,j<=n,j++, \
Kmod[[i,j]]=Kmod[[j,i]]=0];\n      If [Kii==0, Kii=1]; Kmod[[i,i]]=Kii\n      \
]; \n  Return[Kmod]\n];\n\nModifyNodeForcesForDBC[pdof_,f_] := Module[\n  \
{i,k,np=Length[pdof],fmod}, fmod=f; \n     For [k=1,k<=np,k++, i=pdof[[k]]; \
fmod[[i]]=0];\n     Return[fmod]\n];\n\n\
IntForcesShellElement[nodcoor_,elenod_,elemat_,elefab_,\n  eleopt_,u_]:= \
Module[{numele=Length[elenod],\n  numnod=Length[nodcoor],e,eNL,eftab,ni,nj,i,\
\n  ncoor,mprop,fprop,opt,ue,p},\n  p=Table[0,{numele}]; ue=Table[0,{6}];\n  \
For [e=1, e<=numele, e++, \n      eNL=elenod[[e]]; {ni,nj}=eNL;\n      \
ncoor={nodcoor[[ni]],nodcoor[[nj]]};           \n      mprop=elemat[[e]]; \
fprop=elefab[[e]]; opt=eleopt;\n      \
eftab={3*ni-2,3*ni-1,3*ni,3*nj-2,3*nj-1,3*nj}; \n      For [i=1,i<=6,i++, \
ii=eftab[[i]]; ue[[i]]=u[[ii]]];\n      \
p[[e]]=IntForce2DBeamColumn[ncoor,mprop,fprop,opt,ue]\n      ]; \n  Return[p]\
\n];"
}], "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[CellGroupData[{

Cell["\<\
ClearAll[a,alpha,nex,ney,Lx,Ly,R,Em,nu,h];
R=10; Lx=40; Ly=20; nex=8; ney=8;  Em=100; nu=0; h=1;
S=Sqrt[R^2-(Ly/2)^2];  
If [R==Ly/2, alpha=Pi/2,alpha=ArcSin[(Ly/2)/R],
                        alpha=ArcSin[(Ly/2)/R]];
dalpha=alpha/(nex/2);   a=-alpha-dalpha;
ny=ney+1; nx=ney+1; numnod=nx*ny; numele=2*nex*ney;
NodeCoordinates=Table[0,{numnod}]; n=0; 
For [j=1,j<=nx,j++, a=-alpha-dalpha;   
    For [i=1,i<=ny,i++, a=a+dalpha; n++;
        y=R*Sin[a]; z=R*Cos[a]-S; x=Lx*(nx-j)/nex;
        NodeCoordinates[[n]]=N[{x,y,z}];
        ]
    ];
Print[NodeCoordinates];
numele=4*nex*ney;
ElementNodes=Table[{0,0,0},{numele}]; e=n=0;
For [j=1,j<=nex,j++, n++;  
    For [i=1,i<=ney,i++,
        e++; ElementNodes[[e]]={n,n+1,n+nx+1}; 
        e++; ElementNodes[[e]]={n+nx,n,n+nx+1};
        e++; ElementNodes[[e]]={n,n+1,n+nx}; 
        e++; ElementNodes[[e]]={n+1,n+nx+1,n+nx};  
        n++;
        ]
    ];
Print[ElementNodes];
ElemMaterial=Table[{Em,nu},{numele}];
ElemFabrication=Table[h,{numele}];
K=MasterStiffnessShell[NodeCoordinates,ElementNodes,
         ElemMaterial,ElemFabrication,{True}];

(*Print[\"Eigenvalues of K=\",Chop[Eigenvalues[N[K]]]];*)

ClearAll[P,My];
P=3;  My=0;
NodeDisplacement=FreedomTag=FreedomValue=
       Table[{0,0,0,0,0,0},{numnod}]; 
For[i=1,i<=ny,i++, 
   If [i>1&&i<ny, FreedomValue[[i]]={P,0,0,0,0,0}]; 
   If [i==1,  FreedomValue[[i]]={P/2,0,0,0,My,0}];
   If [i==ny, FreedomValue[[i]]={P/2,0,0,0,My,0}]];  
Print[\"Applied node forces=\"]; Print[FreedomValue];
For[i=1,i<=ny,i++, FreedomTag[[i,4]]=1];
For[i=1,i<=ny,i++, FreedomTag[[(nx-1)*ny+i,1]]=1]; 
For[i=1,i<=nx,i++, FreedomTag[[(i-1)*nx+1,3]]=1]; 
For[i=1,i<=nx,i++, FreedomTag[[(i-1)*nx+nx,3]]=1]; 
For[i=1,i<=nx,i++, FreedomTag[[(i-1)*nx+1,5]]=1]; 
For[i=1,i<=nx,i++, FreedomTag[[(i-1)*nx+nx,5]]=1];
For[i=1,i<=nx,i++, FreedomTag[[(i-1)*nx+nex/2+1,2]]=1];
Print[\"Freedom Tags=\"]; Print[FreedomTag];

f=Flatten[FreedomValue];
pdof={}; 
For[n=1,n<=numnod,n++, For[j=1,j<=6,j++,
    If[FreedomTag[[n,j]]>0, AppendTo[pdof,6*(n-1)+j]]]];
Print[\"Prescribed DOF=\",pdof];
    
Kmod=ModifyMasterStiffForDBC[pdof,K];

(*{ev,vec}=Chop[Eigensystem[N[Kmod]]];

Print[\"Eigenvalues of Kmod=\", ev]; numdof=6*numnod;
Print[\"Null vector of Kmod=\", vec[[numdof]] ];*)

fmod=ModifyNodeForcesForDBC [pdof,f];
u=LinearSolve[Kmod,fmod];  u=Chop[u,.000001];
Print[\"Computed Nodal Displacements:\"]; (*Print[u];*)
For[n=1,n<=numnod,n++, For[i=1,i<=6,i++,
    NodeDisplacement[[n,i]]=u[[6*(n-1)+i]] ]];
Print[NodeDisplacement//MatrixForm];\
\>", "Input",
  CellFrame->True,
  CellMargins->{{13, -14}, {Inherited, Inherited}},
  CellLabelMargins->{{6, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[BoxData[
    \({{40.`, \(-10.`\), 0.`}, {40.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {40.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {40.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {40.`, 0.`, 10.`}, {40.`, 3.826834323650898`, 
        9.238795325112868`}, {40.`, 7.0710678118654755`, 
        7.0710678118654755`}, {40.`, 9.238795325112868`, 
        3.8268343236508984`}, {40.`, 10.`, 0.`}, {35.`, \(-10.`\), 
        0.`}, {35.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {35.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {35.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {35.`, 0.`, 10.`}, {35.`, 3.826834323650898`, 
        9.238795325112868`}, {35.`, 7.0710678118654755`, 
        7.0710678118654755`}, {35.`, 9.238795325112868`, 
        3.8268343236508984`}, {35.`, 10.`, 0.`}, {30.`, \(-10.`\), 
        0.`}, {30.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {30.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {30.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {30.`, 0.`, 10.`}, {30.`, 3.826834323650898`, 
        9.238795325112868`}, {30.`, 7.0710678118654755`, 
        7.0710678118654755`}, {30.`, 9.238795325112868`, 
        3.8268343236508984`}, {30.`, 10.`, 0.`}, {25.`, \(-10.`\), 
        0.`}, {25.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {25.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {25.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {25.`, 0.`, 10.`}, {25.`, 3.826834323650898`, 
        9.238795325112868`}, {25.`, 7.0710678118654755`, 
        7.0710678118654755`}, {25.`, 9.238795325112868`, 
        3.8268343236508984`}, {25.`, 10.`, 0.`}, {20.`, \(-10.`\), 
        0.`}, {20.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {20.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {20.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {20.`, 0.`, 10.`}, {20.`, 3.826834323650898`, 
        9.238795325112868`}, {20.`, 7.0710678118654755`, 
        7.0710678118654755`}, {20.`, 9.238795325112868`, 
        3.8268343236508984`}, {20.`, 10.`, 0.`}, {15.`, \(-10.`\), 
        0.`}, {15.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {15.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {15.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {15.`, 0.`, 10.`}, {15.`, 3.826834323650898`, 
        9.238795325112868`}, {15.`, 7.0710678118654755`, 
        7.0710678118654755`}, {15.`, 9.238795325112868`, 
        3.8268343236508984`}, {15.`, 10.`, 0.`}, {10.`, \(-10.`\), 
        0.`}, {10.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {10.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {10.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {10.`, 0.`, 10.`}, {10.`, 3.826834323650898`, 
        9.238795325112868`}, {10.`, 7.0710678118654755`, 
        7.0710678118654755`}, {10.`, 9.238795325112868`, 
        3.8268343236508984`}, {10.`, 10.`, 0.`}, {5.`, \(-10.`\), 
        0.`}, {5.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {5.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {5.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {5.`, 0.`, 10.`}, {5.`, 3.826834323650898`, 
        9.238795325112868`}, {5.`, 7.0710678118654755`, 
        7.0710678118654755`}, {5.`, 9.238795325112868`, 
        3.8268343236508984`}, {5.`, 10.`, 0.`}, {0.`, \(-10.`\), 
        0.`}, {0.`, \(-9.238795325112868`\), 
        3.8268343236508984`}, {0.`, \(-7.0710678118654755`\), 
        7.0710678118654755`}, {0.`, \(-3.826834323650898`\), 
        9.238795325112868`}, {0.`, 0.`, 10.`}, {0.`, 3.826834323650898`, 
        9.238795325112868`}, {0.`, 7.0710678118654755`, 
        7.0710678118654755`}, {0.`, 9.238795325112868`, 
        3.8268343236508984`}, {0.`, 10.`, 0.`}}\)], "Print"],

Cell[BoxData[
    \({{1, 2, 11}, {10, 1, 11}, {1, 2, 10}, {2, 11, 10}, {2, 3, 12}, {11, 2, 
        12}, {2, 3, 11}, {3, 12, 11}, {3, 4, 13}, {12, 3, 13}, {3, 4, 
        12}, {4, 13, 12}, {4, 5, 14}, {13, 4, 14}, {4, 5, 13}, {5, 14, 
        13}, {5, 6, 15}, {14, 5, 15}, {5, 6, 14}, {6, 15, 14}, {6, 7, 
        16}, {15, 6, 16}, {6, 7, 15}, {7, 16, 15}, {7, 8, 17}, {16, 7, 
        17}, {7, 8, 16}, {8, 17, 16}, {8, 9, 18}, {17, 8, 18}, {8, 9, 
        17}, {9, 18, 17}, {10, 11, 20}, {19, 10, 20}, {10, 11, 19}, {11, 20, 
        19}, {11, 12, 21}, {20, 11, 21}, {11, 12, 20}, {12, 21, 20}, {12, 13, 
        22}, {21, 12, 22}, {12, 13, 21}, {13, 22, 21}, {13, 14, 23}, {22, 13, 
        23}, {13, 14, 22}, {14, 23, 22}, {14, 15, 24}, {23, 14, 24}, {14, 15, 
        23}, {15, 24, 23}, {15, 16, 25}, {24, 15, 25}, {15, 16, 24}, {16, 25, 
        24}, {16, 17, 26}, {25, 16, 26}, {16, 17, 25}, {17, 26, 25}, {17, 18, 
        27}, {26, 17, 27}, {17, 18, 26}, {18, 27, 26}, {19, 20, 29}, {28, 19, 
        29}, {19, 20, 28}, {20, 29, 28}, {20, 21, 30}, {29, 20, 30}, {20, 21, 
        29}, {21, 30, 29}, {21, 22, 31}, {30, 21, 31}, {21, 22, 30}, {22, 31, 
        30}, {22, 23, 32}, {31, 22, 32}, {22, 23, 31}, {23, 32, 31}, {23, 24, 
        33}, {32, 23, 33}, {23, 24, 32}, {24, 33, 32}, {24, 25, 34}, {33, 24, 
        34}, {24, 25, 33}, {25, 34, 33}, {25, 26, 35}, {34, 25, 35}, {25, 26, 
        34}, {26, 35, 34}, {26, 27, 36}, {35, 26, 36}, {26, 27, 35}, {27, 36, 
        35}, {28, 29, 38}, {37, 28, 38}, {28, 29, 37}, {29, 38, 37}, {29, 30, 
        39}, {38, 29, 39}, {29, 30, 38}, {30, 39, 38}, {30, 31, 40}, {39, 30, 
        40}, {30, 31, 39}, {31, 40, 39}, {31, 32, 41}, {40, 31, 41}, {31, 32, 
        40}, {32, 41, 40}, {32, 33, 42}, {41, 32, 42}, {32, 33, 41}, {33, 42, 
        41}, {33, 34, 43}, {42, 33, 43}, {33, 34, 42}, {34, 43, 42}, {34, 35, 
        44}, {43, 34, 44}, {34, 35, 43}, {35, 44, 43}, {35, 36, 45}, {44, 35, 
        45}, {35, 36, 44}, {36, 45, 44}, {37, 38, 47}, {46, 37, 47}, {37, 38, 
        46}, {38, 47, 46}, {38, 39, 48}, {47, 38, 48}, {38, 39, 47}, {39, 48, 
        47}, {39, 40, 49}, {48, 39, 49}, {39, 40, 48}, {40, 49, 48}, {40, 41, 
        50}, {49, 40, 50}, {40, 41, 49}, {41, 50, 49}, {41, 42, 51}, {50, 41, 
        51}, {41, 42, 50}, {42, 51, 50}, {42, 43, 52}, {51, 42, 52}, {42, 43, 
        51}, {43, 52, 51}, {43, 44, 53}, {52, 43, 53}, {43, 44, 52}, {44, 53, 
        52}, {44, 45, 54}, {53, 44, 54}, {44, 45, 53}, {45, 54, 53}, {46, 47, 
        56}, {55, 46, 56}, {46, 47, 55}, {47, 56, 55}, {47, 48, 57}, {56, 47, 
        57}, {47, 48, 56}, {48, 57, 56}, {48, 49, 58}, {57, 48, 58}, {48, 49, 
        57}, {49, 58, 57}, {49, 50, 59}, {58, 49, 59}, {49, 50, 58}, {50, 59, 
        58}, {50, 51, 60}, {59, 50, 60}, {50, 51, 59}, {51, 60, 59}, {51, 52, 
        61}, {60, 51, 61}, {51, 52, 60}, {52, 61, 60}, {52, 53, 62}, {61, 52, 
        62}, {52, 53, 61}, {53, 62, 61}, {53, 54, 63}, {62, 53, 63}, {53, 54, 
        62}, {54, 63, 62}, {55, 56, 65}, {64, 55, 65}, {55, 56, 64}, {56, 65, 
        64}, {56, 57, 66}, {65, 56, 66}, {56, 57, 65}, {57, 66, 65}, {57, 58, 
        67}, {66, 57, 67}, {57, 58, 66}, {58, 67, 66}, {58, 59, 68}, {67, 58, 
        68}, {58, 59, 67}, {59, 68, 67}, {59, 60, 69}, {68, 59, 69}, {59, 60, 
        68}, {60, 69, 68}, {60, 61, 70}, {69, 60, 70}, {60, 61, 69}, {61, 70, 
        69}, {61, 62, 71}, {70, 61, 71}, {61, 62, 70}, {62, 71, 70}, {62, 63, 
        72}, {71, 62, 72}, {62, 63, 71}, {63, 72, 71}, {64, 65, 74}, {73, 64, 
        74}, {64, 65, 73}, {65, 74, 73}, {65, 66, 75}, {74, 65, 75}, {65, 66, 
        74}, {66, 75, 74}, {66, 67, 76}, {75, 66, 76}, {66, 67, 75}, {67, 76, 
        75}, {67, 68, 77}, {76, 67, 77}, {67, 68, 76}, {68, 77, 76}, {68, 69, 
        78}, {77, 68, 78}, {68, 69, 77}, {69, 78, 77}, {69, 70, 79}, {78, 69, 
        79}, {69, 70, 78}, {70, 79, 78}, {70, 71, 80}, {79, 70, 80}, {70, 71, 
        79}, {71, 80, 79}, {71, 72, 81}, {80, 71, 81}, {71, 72, 80}, {72, 81, 
        80}}\)], "Print"],

Cell[BoxData[
    \("Applied node forces="\)], "Print"],

Cell[BoxData[
    \({{3\/2, 0, 0, 0, 0, 0}, {3, 0, 0, 0, 0, 0}, {3, 0, 0, 0, 0, 0}, {3, 0, 
        0, 0, 0, 0}, {3, 0, 0, 0, 0, 0}, {3, 0, 0, 0, 0, 0}, {3, 0, 0, 0, 0, 
        0}, {3, 0, 0, 0, 0, 0}, {3\/2, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}}\)], "Print"],

Cell[BoxData[
    \("Freedom Tags="\)], "Print"],

Cell[BoxData[
    \({{0, 0, 1, 1, 1, 0}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 
        1, 0, 0}, {0, 1, 0, 1, 0, 0}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 1, 0, 
        0}, {0, 0, 0, 1, 0, 0}, {0, 0, 1, 1, 1, 0}, {0, 0, 1, 0, 1, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 1, 0, 1, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 
        0}, {0, 0, 1, 0, 1, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 1, 0, 1, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 
        0, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 
        1, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 1, 0, 
        1, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
        0}, {0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 
        0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0}, {1, 0, 1, 0, 1, 0}, {1, 0, 0, 0, 
        0, 0}, {1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 0, 
        0}, {1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {1, 
        0, 1, 0, 1, 0}}\)], "Print"],

Cell[BoxData[
    InterpretationBox[\("Prescribed DOF="\[InvisibleSpace]{3, 4, 5, 10, 16, 
          22, 26, 28, 34, 40, 46, 51, 52, 53, 57, 59, 80, 105, 107, 111, 113, 
          134, 159, 161, 165, 167, 188, 213, 215, 219, 221, 242, 267, 269, 
          273, 275, 296, 321, 323, 327, 329, 350, 375, 377, 381, 383, 404, 
          429, 431, 433, 435, 437, 439, 445, 451, 457, 458, 463, 469, 475, 
          481, 483, 485}\),
      SequenceForm[ 
      "Prescribed DOF=", {3, 4, 5, 10, 16, 22, 26, 28, 34, 40, 46, 51, 52, 53,
         57, 59, 80, 105, 107, 111, 113, 134, 159, 161, 165, 167, 188, 213, 
        215, 219, 221, 242, 267, 269, 273, 275, 296, 321, 323, 327, 329, 350, 
        375, 377, 381, 383, 404, 429, 431, 433, 435, 437, 439, 445, 451, 457, 
        458, 463, 469, 475, 481, 483, 485}],
      Editable->False]], "Print"],

Cell[BoxData[
    \("Computed Nodal Displacements:"\)], "Print"],

Cell[BoxData[
    TagBox[
      RowBox[{"(", "\[NoBreak]", GridBox[{
            {"0.16048400470903734`", "0.027095109219610093`", "0", "0", "0", 
              "0.017594092966067603`"},
            {"0.160375087841624`", 
              "0.02539003393734254`", \(-0.010778342671650934`\), "0", 
              "0.00661558034632277`", "0.01599414401748093`"},
            {"0.1603493135887319`", 
              "0.019536058167425575`", \(-0.019882729984017522`\), "0", 
              "0.012215306415998632`", "0.012217967215469697`"},
            {"0.1603486968636527`", 
              "0.010597326414386328`", \(-0.025952793799051726`\), "0", 
              "0.015975147467105048`", "0.006617762392124577`"},
            {"0.16034914047021528`", "0", \(-0.028083910892759464`\), "0", 
              "0.017298677180290427`", "0"},
            {
              "0.16034869686365213`", \(-0.010597326414385723`\), \
\(-0.025952793799052905`\), "0", 
              "0.01597514746710498`", \(-0.006617762392124307`\)},
            {
              "0.1603493135887308`", \(-0.0195360581674251`\), \
\(-0.019882729984018525`\), "0", 
              "0.01221530641599849`", \(-0.01221796721546958`\)},
            {
              "0.16037508784162305`", \(-0.025390033937340947`\), \
\(-0.010778342671651299`\), "0", 
              "0.0066155803463227155`", \(-0.0159941440174809`\)},
            {"0.1604840047090365`", \(-0.027095109219606565`\), "0", "0", 
              "0", \(-0.01759409296606761`\)},
            {"0.13779449922650816`", \(-0.004333841320861762`\), 
              "0", \(-0.0006999798944649924`\), 
              "0", \(-0.000497342935703113`\)},
            {"0.13776202168546076`", \(-0.002791358637611063`\), 
              "0.0008293529594021684`", \(-0.0001416780568995586`\), \
\(-0.00021840126405155242`\), \(-0.0005066258147228275`\)},
            {"0.13773759836748448`", \(-0.001998196279619036`\), 
              "0.0016225381380803921`", \(-0.000012246130795166589`\), \
\(-0.0004047347867333915`\), \(-0.00039765468889877337`\)},
            {"0.13772882317372118`", \(-0.0010814586656288287`\), 
              "0.002232490632963065`", 
              "3.470260064669395`*^-6", \(-0.0005209084513355357`\), \
\(-0.00021390887850329408`\)},
            {"0.13772704896415436`", "0", "0.002455956843890499`", 
              "0", \(-0.0005597043901093025`\), "0"},
            {"0.1377288231737203`", "0.0010814586656285746`", 
              "0.002232490632961217`", \(-3.470260064932211`*^-6\), \
\(-0.0005209084513359123`\), "0.00021390887850358486`"},
            {"0.13773759836748337`", "0.0019981962796189676`", 
              "0.0016225381380786984`", 
              "0.000012246130795507266`", \(-0.0004047347867334359`\), 
              "0.0003976546888988189`"},
            {"0.13776202168545978`", "0.002791358637612759`", 
              "0.0008293529594016093`", 
              "0.00014167805690024635`", \(-0.00021840126405153743`\), 
              "0.000506625814722685`"},
            {"0.13779449922650722`", "0.004333841320866263`", "0", 
              "0.0006999798944657742`", "0", "0.0004973429357029727`"},
            {"0.11857855265121806`", \(-0.001641901613340331`\), 
              "0", \(-0.00034746150599325644`\), 
              "0", \(-0.0002852695383264864`\)},
            {
              "0.11861262166397692`", \(-0.0005767490654172114`\), \
\(-0.00009954039693956279`\), \(-0.00017215216154052224`\), \
\(-0.00006824132625295842`\), \(-0.00017104545575165035`\)},
            {
              "0.11859648533055749`", \(-0.0002058594672967346`\), \
\(-0.00021713952710917153`\), \(-0.000038582948864451384`\), \
\(-0.0001489438010119438`\), \(-0.00014129110554950047`\)},
            {
              "0.11857475954299995`", \(-0.0000970442789607837`\), \
\(-0.0001948515862460819`\), 
              "3.2297610609434805`*^-6", \(-0.00020892394134541055`\), \
\(-0.00008193814122549756`\)},
            {"0.11856728264399523`", "0", \(-0.00016216717252666018`\), 
              "0", \(-0.0002282701363443303`\), "0"},
            {"0.11857475954299938`", 
              "0.00009704427896029768`", \(-0.0001948515862482581`\), \
\(-3.2297610612887638`*^-6\), \(-0.00020892394134540515`\), 
              "0.00008193814122556476`"},
            {"0.11859648533055651`", 
              "0.00020585946729619842`", \(-0.0002171395271114405`\), 
              "0.00003858294886476382`", \(-0.0001489438010121888`\), 
              "0.00014129110554976444`"},
            {"0.11861262166397603`", 
              "0.0005767490654188451`", \(-0.00009954039694038621`\), 
              "0.00017215216154152466`", \(-0.0000682413262530634`\), 
              "0.00017104545575186036`"},
            {"0.11857855265121729`", "0.0016419016133463792`", "0", 
              "0.00034746150599440656`", "0", "0.0002852695383263222`"},
            {"0.09939983711387006`", \(-0.001065742988769653`\), 
              "0", \(-0.00019350720078273914`\), 
              "0", \(-0.000021459662360646185`\)},
            {
              "0.09942299030044724`", \(-0.00038573851422467895`\), \
\(-0.0001551556096080986`\), \(-0.00014885593284465697`\), 
              "0.000014393315364745138`", "0.000025425352482395273`"},
            {
              "0.09942193415226697`", \(-0.00004586173690446714`\), \
\(-0.0004025677960573501`\), \(-0.00006593136912375268`\), 
              "0.000019358426173006834`", "0.000022457824155902515`"},
            {"0.09940430063266675`", 
              "0.000016718339485505745`", \(-0.000528877857610043`\), \
\(-0.000014120960974671171`\), "9.89987521770421`*^-6", 
              "7.977789528470244`*^-6"},
            {"0.09939535130052089`", "0", \(-0.0005512093969470917`\), "0", 
              "4.348600755546752`*^-6", "0"},
            {
              "0.09940430063266646`", \(-0.000016718339485887395`\), \
\(-0.0005288778576117523`\), "0.000014120960974319339`", 
              "9.899875217777837`*^-6", \(-7.977789528438233`*^-6\)},
            {"0.09942193415226638`", 
              "0.0000458617369038124`", \(-0.0004025677960594896`\), 
              "0.000065931369123901`", 
              "0.000019358426173175505`", \(-0.000022457824155984874`\)},
            {"0.09942299030044642`", 
              "0.0003857385142257878`", \(-0.0001551556096090628`\), 
              "0.00014885593284560592`", 
              "0.000014393315364752549`", \(-0.000025425352482316445`\)},
            {"0.09939983711386922`", "0.0010657429887754788`", "0", 
              "0.0001935072007840835`", "0", "0.00002145966236075464`"},
            {"0.08018153245598719`", \(-0.000998672971124753`\), 
              "0", \(-0.00015152626075726186`\), 
              "0", \(-0.000011161568238490528`\)},
            {
              "0.08018925806455864`", \(-0.0004520505318572291`\), \
\(-0.00011098932886013218`\), \(-0.0001290002074468328`\), 
              "2.0852583328090755`*^-6", "0"},
            {
              "0.08019198656497696`", \(-0.00011512483975602373`\), \
\(-0.00033704708151145865`\), \(-0.00007548235871568516`\), 
              "3.7097572031159286`*^-6", "3.7769326888358565`*^-6"},
            {
              "0.08018521777952123`", \(-7.413333678016997`*^-6\), \
\(-0.0004991402225555232`\), \(-0.000027521917990850317`\), "0", 
              "1.13958883237509`*^-6"},
            {"0.080180361135488`", "0", \(-0.0005468048839899072`\), 
              "0", \(-3.8017371815413625`*^-6\), "0"},
            {"0.08018521777952096`", 
              "7.413333677728894`*^-6", \(-0.0004991402225568139`\), 
              "0.00002752191799059835`", "0", \(-1.1395888323191762`*^-6\)},
            {"0.0801919865649765`", 
              "0.00011512483975553291`", \(-0.00033704708151306506`\), 
              "0.00007548235871578972`", 
              "3.7097572031827655`*^-6", \(-3.7769326888163502`*^-6\)},
            {"0.08018925806455793`", 
              "0.0004520505318580261`", \(-0.00011098932886091765`\), 
              "0.00012900020744755657`", "2.0852583328717094`*^-6", "0"},
            {"0.08018153245598626`", "0.000998672971129379`", "0", 
              "0.00015152626075844557`", "0", "0.000011161568238722308`"},
            {"0.06096305928303822`", \(-0.0009462498586949745`\), 
              "0", \(-0.0001530809252387379`\), 
              "0", \(-6.2033239685912265`*^-6\)},
            {
              "0.060958479417914734`", \(-0.0003967064746987928`\), \
\(-0.00012903014865529101`\), \(-0.0001249294478803969`\), \
\(-8.857565204936893`*^-6\), \(-0.00002169645230849977`\)},
            {
              "0.06096361584721887`", \(-0.00008630469368882174`\), \
\(-0.0003578874443480921`\), \(-0.00007239024748363946`\), \
\(-0.000010450741465010426`\), \(-0.000012307687555785934`\)},
            {"0.06096368854374677`", 
              "7.4716829643846466`*^-6", \(-0.0005319620157042892`\), \
\(-0.00003145950905451617`\), \(-0.00001227208305619935`\), \
\(-5.363961296149957`*^-6\)},
            {"0.060961229387726974`", "0", \(-0.0005941199463460002`\), 
              "0", \(-0.00001539678588004887`\), "0"},
            {
              "0.06096368854374669`", \(-7.471682964595456`*^-6\), \
\(-0.0005319620157051718`\), 
              "0.00003145950905433513`", \(-0.000012272083056104402`\), 
              "5.363961296152434`*^-6"},
            {"0.06096361584721854`", 
              "0.00008630469368844877`", \(-0.0003578874443491788`\), 
              "0.00007239024748371355`", \(-0.00001045074146486581`\), 
              "0.000012307687555737812`"},
            {"0.060958479417914074`", 
              "0.0003967064746993036`", \(-0.000129030148655801`\), 
              "0.00012492944788092708`", \(-8.85756520488059`*^-6\), 
              "0.000021696452308580865`"},
            {"0.060963059283037246`", "0.00094624985869828`", "0", 
              "0.00015308092523955217`", "0", "6.203323969022683`*^-6"},
            {"0.041785110598164005`", \(-0.0013472575633650173`\), 
              "0", \(-0.00023422661274861623`\), "0", 
              "0.00023820370936774595`"},
            {
              "0.0417786125178582`", \(-0.0006298270751260415`\), \
\(-0.000029466152459459617`\), \(-0.00011603327220161539`\), 
              "0.00007900766209283046`", "0.00018549609291552746`"},
            {
              "0.04179030279377126`", \(-0.0003108188064376271`\), \
\(-0.00011668708043911932`\), \(-0.00006128976563578508`\), 
              "0.00015913247532828173`", "0.00015660391685536716`"},
            {
              "0.04178591285863299`", \(-0.00011847897366866571`\), \
\(-0.00022071335719785787`\), \(-0.00003730663329293441`\), 
              "0.00020094052521837512`", "0.00008433926051641334`"},
            {"0.04178056914298976`", "0", \(-0.0002778985822716775`\), "0", 
              "0.0002111259692041739`", "0"},
            {"0.041785912858632944`", 
              "0.00011847897366855922`", \(-0.00022071335719823438`\), 
              "0.000037306633292837406`", 
              "0.00020094052521848945`", \(-0.00008433926051643193`\)},
            {"0.04179030279377104`", 
              "0.00031081880643743184`", \(-0.0001166870804396053`\), 
              "0.00006128976563582944`", 
              "0.0001591324753283573`", \(-0.0001566039168553662`\)},
            {"0.04177861251785772`", 
              "0.0006298270751262155`", \(-0.00002946615245968837`\), 
              "0.00011603327220183573`", 
              "0.00007900766209289385`", \(-0.00018549609291547116`\)},
            {"0.0417851105981632`", "0.0013472575633663407`", "0", 
              "0.00023422661274892498`", "0", \(-0.0002382037093672892`\)},
            {"0.022581850640945984`", \(-0.003726256942343171`\), 
              "0", \(-0.00041328229162844915`\), "0", 
              "0.0004082744941190682`"},
            {"0.02264416464752724`", \(-0.002989628102579266`\), 
              "0.000987312252177279`", \(-0.00005850696666226051`\), 
              "0.0002379573432509847`", "0.0005509695506482787`"},
            {"0.02263708063130135`", \(-0.0021796293468348614`\), 
              "0.0017604430341790814`", \(-0.00007180631396729338`\), 
              "0.00040300248497019396`", "0.0004068347559387164`"},
            {"0.022622404951431113`", \(-0.0011097843875752969`\), 
              "0.002145865102282747`", \(-0.000056273339510864485`\), 
              "0.0005061571399995519`", "0.00021275273248978709`"},
            {"0.022617541233341656`", "0", "0.002246841712703079`", "0", 
              "0.000543510108419775`", "0"},
            {"0.022622404951431096`", "0.0011097843875752978`", 
              "0.0021458651022828274`", "0.00005627333951086397`", 
              "0.0005061571399995968`", \(-0.00021275273248979882`\)},
            {"0.022637080631301276`", "0.002179629346834828`", 
              "0.0017604430341791246`", "0.00007180631396729176`", 
              "0.0004030024849703141`", \(-0.00040683475593878446`\)},
            {"0.02264416464752696`", "0.002989628102579196`", 
              "0.0009873122521773395`", "0.00005850696666220143`", 
              "0.00023795734325101055`", \(-0.0005509695506482112`\)},
            {"0.022581850640945526`", "0.0037262569423426374`", "0", 
              "0.00041328229162831666`", "0", \(-0.00040827449411882786`\)},
            {"0", "0.027984918720361858`", "0", "0.0009292779942258843`", 
              "0", \(-0.017606170653594447`\)},
            {"0", 
              "0.02490575594629221`", \(-0.010542133812582165`\), \
\(-0.000052710464122062345`\), \(-0.0066332279490009015`\), \
\(-0.015973757310564864`\)},
            {"0", 
              "0.019341370333582197`", \(-0.01977639600504438`\), \
\(-0.00014226598456721146`\), \(-0.012215259484430795`\), \
\(-0.012220191224550543`\)},
            {"0", 
              "0.010590056573301929`", \(-0.026139282149654035`\), \
\(-0.00007624019829848485`\), \(-0.0159952111767461`\), \
\(-0.006625558817937257`\)},
            {"0", "0", \(-0.028393153990399968`\), 
              "0", \(-0.017321011361983586`\), "0"},
            {"0", \(-0.010590056573301816`\), \(-0.026139282149653685`\), 
              "0.00007624019829855939`", \(-0.01599521117674605`\), 
              "0.006625558817937233`"},
            {"0", \(-0.019341370333582`\), \(-0.019776396005043967`\), 
              "0.00014226598456715335`", \(-0.01221525948443071`\), 
              "0.012220191224550458`"},
            {"0", \(-0.02490575594629247`\), \(-0.010542133812582063`\), 
              "0.000052710464121906735`", \(-0.006633227949000904`\), 
              "0.01597375731056486`"},
            {"0", \(-0.027984918720362805`\), 
              "0", \(-0.0009292779942260736`\), "0", "0.01760617065359438`"}
            }], "\[NoBreak]", ")"}],
      (MatrixForm[ #]&)]], "Print"]
}, Open  ]]
},
FrontEndVersion->"4.0 for Macintosh",
ScreenRectangle->{{0, 1152}, {0, 850}},
AutoGeneratedPackage->None,
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{739, 622},
WindowMargins->{{9, Automatic}, {Automatic, 21}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
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[CellGroupData[{
Cell[1739, 51, 5969, 94, 2233, "Input",
  InitializationCell->True],
Cell[7711, 147, 104, 2, 42, "Print"],
Cell[7818, 151, 84, 1, 22, "Print"],
Cell[7905, 154, 84, 1, 22, "Print"],
Cell[7992, 157, 223, 3, 22, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[8252, 165, 1150, 23, 463, "Input",
  InitializationCell->True],
Cell[9405, 190, 90, 1, 38, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[9532, 196, 1141, 22, 448, "Input",
  InitializationCell->True],
Cell[10676, 220, 87, 1, 38, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[10800, 226, 8518, 126, 3043, "Input",
  InitializationCell->True],
Cell[19321, 354, 112, 2, 22, "Print"],
Cell[19436, 358, 6451, 113, 367, "Print"],
Cell[25890, 473, 112, 2, 22, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[26039, 480, 11163, 168, 4273, "Input",
  InitializationCell->True],
Cell[37205, 650, 1981, 31, 142, "Print"],
Cell[39189, 683, 120, 2, 22, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[39346, 690, 7102, 106, 2893, "Input",
  InitializationCell->True],
Cell[46451, 798, 253, 5, 22, "Print"],
Cell[46707, 805, 897, 14, 142, "Print"],
Cell[47607, 821, 3286, 56, 202, "Print"],
Cell[50896, 879, 120, 2, 22, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[51053, 886, 11308, 168, 4348, "Input",
  InitializationCell->True],
Cell[62364, 1056, 253, 5, 22, "Print"],
Cell[62620, 1063, 163, 4, 38, "Print"],
Cell[62786, 1069, 1267, 20, 142, "Print"],
Cell[64056, 1091, 119, 2, 22, "Print"]
}, Open  ]],
Cell[64190, 1096, 51, 1, 27, "Input"],
Cell[64244, 1099, 376, 8, 81, "Text"],

Cell[CellGroupData[{
Cell[64645, 1111, 13855, 223, 4363, "Input",
  InitializationCell->True],
Cell[78503, 1336, 1196, 19, 276, "Print"],
Cell[79702, 1357, 120, 2, 22, "Print"]
}, Open  ]],

Cell[CellGroupData[{
Cell[79859, 1364, 5185, 80, 2188, "Input",
  InitializationCell->True],
Cell[85047, 1446, 2976, 45, 616, "Print"],
Cell[88026, 1493, 767, 12, 37, "Print"]
}, Open  ]],
Cell[88808, 1508, 3944, 87, 1081, "Input"],

Cell[CellGroupData[{
Cell[92777, 1599, 828, 19, 208, "Input"],
Cell[93608, 1620, 1798, 55, 635, "Print"]
}, Open  ]],
Cell[95421, 1678, 4850, 106, 1328, "Input"],

Cell[CellGroupData[{
Cell[100296, 1788, 989, 24, 223, "Input"],
Cell[101288, 1814, 1874, 54, 518, "Print"]
}, Open  ]],
Cell[103177, 1871, 2758, 46, 928, "Input",
  InitializationCell->True],

Cell[CellGroupData[{
Cell[105960, 1921, 2745, 77, 1093, "Input"],
Cell[108708, 2000, 3912, 63, 322, "Print"],
Cell[112623, 2065, 4028, 52, 577, "Print"],
Cell[116654, 2119, 55, 1, 22, "Print"],
Cell[116712, 2122, 1875, 25, 306, "Print"],
Cell[118590, 2149, 48, 1, 22, "Print"],
Cell[118641, 2152, 1869, 25, 262, "Print"],
Cell[120513, 2179, 839, 13, 67, "Print"],
Cell[121355, 2194, 64, 1, 22, "Print"],
Cell[121422, 2197, 14838, 266, 1361, "Print"]
}, Open  ]]
}
]
*)




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