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

                    Mathematica-Compatible Notebook

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

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

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

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

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

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

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

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

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     34326,        688]*)
(*NotebookOutlinePosition[     35397,        722]*)
(*  CellTagsIndexPosition[     35353,        718]*)
(*WindowFrame->Normal*)



Notebook[{
Cell["Strain-displacement utilities 8-node brick", "Text",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "\nIsoHexa8ShapeFunCarDer[xyznlist_,\[Xi]\[Eta]\[Mu]list_,numer_]:= \n  \
Module[{xn,yn,zn,dN\[Xi],dN\[Eta],dN\[Mu],J,J11,J21,J31,J12,J22,J32,\n   \
J13,J23,J33,Jinv,Jdet,Bx,By,Bz},{\[Xi],\[Eta],\[Mu]}=\[Xi]\[Eta]\[Mu]list;\n  \
 xn=Table[xyznlist[[i,1]],{i,8}];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   yn=Table[xyznlist[[i,2]],{i,8}]; ",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "   zn=Table[xyznlist[[i,3]],{i,8}];     \n   \
NHexa8[\[Xi]_,\[Eta]_,\[Mu]_]:={ \n       (1-\[Xi])*(1-\[Eta])*(1-\[Mu]),(1+\
\[Xi])*(1-\[Eta])*(1-\[Mu]),\n       \
(1+\[Xi])*(1+\[Eta])*(1-\[Mu]),(1-\[Xi])*(1+\[Eta])*(1-\[Mu]),\n       (1-\
\[Xi])*(1-\[Eta])*(1+\[Mu]),(1+\[Xi])*(1-\[Eta])*(1+\[Mu]),\n       \
(1+\[Xi])*(1+\[Eta])*(1+\[Mu]),(1-\[Xi])*(1+\[Eta])*(1+\[Mu])}/8;\n   \
dNHexa8[\[Xi]_,\[Eta]_,\[Mu]_]:={\n       \
{-(1-\[Eta])*(1-\[Mu]),(1-\[Eta])*(1-\[Mu]),(1+\[Eta])*(1-\[Mu]),-(1+\[Eta])*(\
1-\[Mu]),\n        -(1-\[Eta])*(1+\[Mu]),(1-\[Eta])*(1+\[Mu]),(1+\[Eta])*(1+\
\[Mu]),-(1+\[Eta])*(1+\[Mu])},\n       {-(1-\[Xi])*(1-\[Mu]),-(1+\[Xi])*(1-\
\[Mu]),(1+\[Xi])*(1-\[Mu]),(1-\[Xi])*(1-\[Mu]),\n        \
-(1-\[Xi])*(1+\[Mu]),-(1+\[Xi])*(1+\[Mu]),(1+\[Xi])*(1+\[Mu]),(1-\[Xi])*(1+\
\[Mu])},\n       {-(1-\[Xi])*(1-\[Eta]),-(1+\[Xi])*(1-\[Eta]),-(1+\[Xi])*(1+\
\[Eta]),-(1-\[Xi])*(1+\[Eta]),\n         (1-\[Xi])*(1-\[Eta]),(1+\[Xi])*(1-\
\[Eta]),(1+\[Xi])*(1+\[Eta]),(1-\[Xi])*(1+\[Eta])}}/8;\n   \
{dN\[Xi],dN\[Eta],dN\[Mu]}=dNHexa8[\[Xi],\[Eta],\[Mu]];\n   J11=dN\[Xi].xn; \
J21=dN\[Eta].xn; J31=dN\[Mu].xn;\n   J12=dN\[Xi].yn; J22=dN\[Eta].yn; J32=dN\
\[Mu].yn;\n   J13=dN\[Xi].zn; J23=dN\[Eta].zn; J33=dN\[Mu].zn;  \n   \
J={{J11,J12,J13},{J21,J22,J23},{J31,J32,J33}}; ",
  StyleBox["(*Print[\"J=\",J//MatrixForm];\n   \
J0=Simplify[J/.{\[Xi]->0,\[Eta]->0,\[Mu]->0}];  \
Print[\"J0=\",J0//MatrixForm];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   If [!numer,J=Simplify[J]]; \n   \
{{J11,J12,J13},{J21,J22,J23},{J31,J32,J33}}=J;\n   Jdet= \
J11*J22*J33+J21*J32*J13+J31*J12*J23-\n         \
J31*J22*J13-J11*J32*J23-J21*J12*J33;\n   If [!numer,Jdet=Simplify[Jdet]];\n   \
Jinv={{J22*J33-J32*J23,J32*J13-J12*J33,J12*J23-J22*J13},\n         \
{J31*J23-J21*J33,J11*J33-J31*J13,J21*J13-J11*J23},\n         \
{J21*J32-J31*J22,J31*J12-J11*J32,J11*J22-J21*J12}};\n   If \
[!numer,Jinv=Simplify[Jinv]];         \n   \
{Bx,By,Bz}=Jinv.{dN\[Xi],dN\[Eta],dN\[Mu]}/Jdet;\n   If \
[!numer,{Bx,By,Bz}=Simplify[{Bx,By,Bz}]];\n   Return[{Bx,By,Bz,Jdet}]]; \n    \
 \nIsoHexa8StrainDispMatrix[xyznlist_,\[Xi]\[Eta]\[Mu]list_,numer_]:= \n  \
Module[{n,Bx,By,Bz,Bexx,Beyy,Bezz,Bgxy,Bgyz,Bgzx,Jdet},\n  \
{Bx,By,Bz,Jdet}=IsoHexa8ShapeFunCarDer[xyznlist,\[Xi]\[Eta]\[Mu]list,numer];\n\
   Bexx=Flatten[Table[{Bx[[n]],0,0},{n,8}]];\n   \
Beyy=Flatten[Table[{0,By[[n]],0},{n,8}]];\n   \
Bezz=Flatten[Table[{0,0,Bz[[n]]},{n,8}]];\n   \
Bgxy=Flatten[Table[{By[[n]],Bx[[n]],0},{n,8}]];\n   \
Bgyz=Flatten[Table[{0,Bz[[n]],By[[n]]},{n,8}]];\n   \
Bgzx=Flatten[Table[{Bz[[n]],0,Bx[[n]]},{n,8}]];\n   \
Be={Bexx,Beyy,Bezz,Bgxy,Bgyz,Bgzx};\n   Return[Be]];\n   \n\
CubicHexa8TemplateGrcMatrices[{a_,b_,c_}]:=Module[{Gr,Gc}, \n   Gr=\
\[InvisibleSpace]{{2,0,0,-b,0,c},{0,2,0,a,-c,0},{0,0,2,0,b,-a},\n      \
{2,0,0,-b,0,c},{0,2,0,-a,-c,0},{0,0,2,0,b,a},\n      \
{2,0,0,b,0,c},{0,2,0,-a,-c,0},{0,0,2,0,-b,a},\n      \
{2,0,0,b,0,c},{0,2,0,a,-c,0},{0,0,2,0,-b,-a},\n      \
{2,0,0,-b,0,-c},{0,2,0,a,c,0},{0,0,2,0,b,-a},\n      \
{2,0,0,-b,0,-c},{0,2,0,-a,c,0},{0,0,2,0,b,a},\n      \
{2,0,0,b,0,-c},{0,2,0,-a,c,0},{0,0,2,0,-b,a},\n      \
{2,0,0,b,0,-c},{0,2,0,a,c,0},{0,0,2,0,-b,-a}}/2;\n   \
Gc=\[InvisibleSpace]{{-2*a,0,0,-b,0,-c},{0,-2*b,0,-a,-c,0},{0,0,-2*c,0,-b,-a},\
\n      {2*a,0,0,-b,0,-c},{0,-2*b,0,a,-c,0},{0,0,-2*c,0,-b,a},\n      \
{2*a,0,0,b,0,-c},{0,2*b,0,a,-c,0},{0,0,-2*c,0,b,a},\n      \
{-2*a,0,0,b,0,-c},{0,2*b,0,-a,-c,0},{0,0,-2*c,0,b,-a},\n      \
{-2*a,0,0,-b,0,c},{0,-2*b,0,-a,c,0},{0,0,2*c,0,-b,-a},\n      \
{2*a,0,0,-b,0,c},{0,-2*b,0,a,c,0},{0,0,2*c,0,-b,a},\n      \
{2*a,0,0,b,0,c},{0,2*b,0,a,c,0},{0,0,2*c,0,b,a},\n      \
{-2*a,0,0,b,0,c},{0,2*b,0,-a,c,0},{0,0,2*c,0,b,-a}}/4;\n   Return[{Gr,Gc}]];\n\
 \nCubicHexa8TemplateLMatrix[{a_,b_,c_}]:=Module[{LT}, LT={\n   b*c*{-1,0,0, \
1,0,0, 1,0,0, -1,0,0, -1,0,0, 1,0,0, 1,0,0, -1,0,0},\n   c*a*{0,-1,0, 0,-1,0, \
0,1,0, 0,1,0, 0,-1,0, 0,-1,0, 0,1,0, 0,1,0},\n   a*b*{0,0,-1, 0,0,-1, 0,0,-1, \
0,0,-1, 0,0,1, 0,0,1, 0,0,1, 0,0,1},\n   c*a*{-1,0,0, -1,0,0, 1,0,0, 1,0,0, \
-1,0,0, -1,0,0, 1,0,0, 1,0,0}+\n   b*c*{0,-1,0, 0,1,0, 0,1,0, 0,-1,0, 0,-1,0, \
0,1,0, 0,1,0, 0,-1,0},\n   c*a*{0,0,-1, 0,0,-1, 0,0,1, 0,0,1, 0,0,-1, 0,0,-1, \
0,0,1, 0,0,1}+\n   a*b*{0,-1,0, 0,-1,0, 0,-1,0, 0,-1,0, 0,1,0, 0,1,0, 0,1,0, \
0,1,0},\n   b*c*{0,0,-1, 0,0,1, 0,0,1, 0,0,-1, 0,0,-1, 0,0,1, 0,0,1, 0,0,-1}+\
\n   a*b*{-1,0,0, -1,0,0, -1,0,0, -1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0}}/4;\n   \
Return[Transpose[LT]]];\n        \n\
CubicHexa8TemplateHhMatrix[{a_,b_,c_}]:=Module[{HhT},  \n  HhT=\
\[InvisibleSpace]{{-2*a*b,-2*a*c,0,0,0,0,2*b*c,0,0,-a*b*c,0,0},\n      \
{0,0,-2*a*b,-2*b*c,0,0,0,2*a*c,0,0,-a*b*c,0},\n      \
{0,0,0,0,-2*a*c,-2*b*c,0,0,2*a*b,0,0,-a*b*c},\n      \
{2*a*b,2*a*c,0,0,0,0,2*b*c,0,0,a*b*c,0,0},\n      \
{0,0,2*a*b,-2*b*c,0,0,0,-2*a*c,0,0,a*b*c,0},\n      \
{0,0,0,0,2*a*c,-2*b*c,0,0,-2*a*b,0,0,a*b*c},\n      \
{-2*a*b,2*a*c,0,0,0,0,-2*b*c,0,0,-a*b*c,0,0},\n      \
{0,0,-2*a*b,2*b*c,0,0,0,-2*a*c,0,0,-a*b*c,0},\n      \
{0,0,0,0,2*a*c,2*b*c,0,0,2*a*b,0,0,-a*b*c},\n      \
{2*a*b,-2*a*c,0,0,0,0,-2*b*c,0,0,a*b*c,0,0},\n      \
{0,0,2*a*b,2*b*c,0,0,0,2*a*c,0,0,a*b*c,0},\n      \
{0,0,0,0,-2*a*c,2*b*c,0,0,-2*a*b,0,0,a*b*c},\n      \
{-2*a*b,2*a*c,0,0,0,0,-2*b*c,0,0,a*b*c,0,0},\n      \
{0,0,-2*a*b,2*b*c,0,0,0,-2*a*c,0,0,a*b*c,0},\n      \
{0,0,0,0,2*a*c,2*b*c,0,0,2*a*b,0,0,a*b*c},\n      \
{2*a*b,-2*a*c,0,0,0,0,-2*b*c,0,0,-a*b*c,0,0},\n      \
{0,0,2*a*b,2*b*c,0,0,0,2*a*c,0,0,-a*b*c,0},\n      \
{0,0,0,0,-2*a*c,2*b*c,0,0,-2*a*b,0,0,-a*b*c},\n      \
{-2*a*b,-2*a*c,0,0,0,0,2*b*c,0,0,a*b*c,0,0},\n      \
{0,0,-2*a*b,-2*b*c,0,0,0,2*a*c,0,0,a*b*c,0},\n      \
{0,0,0,0,-2*a*c,-2*b*c,0,0,2*a*b,0,0,a*b*c},\n      \
{2*a*b,2*a*c,0,0,0,0,2*b*c,0,0,-a*b*c,0,0},\n      \
{0,0,2*a*b,-2*b*c,0,0,0,-2*a*c,0,0,-a*b*c,0},\n      \
{0,0,0,0,2*a*c,-2*b*c,0,0,-2*a*b,0,0,-a*b*c}}/8;\n      \
Return[Transpose[HhT]]];\n      \n\
CubicHexa8TemplateWMatrix[{a_,b_,c_}]:=Module[{W},  \n  \
W=DiagonalMatrix[{1/(a*b),1/(a*c),1/(a*b),1/(b*c),\n   \
1/(a*c),1/(b*c),1/(b*c),1/(a*c),1/(a*b),\n   1/(a*b*c),1/(a*b*c),1/(a*b*c)}]; \
 (* First cut by eyeballing *)\n  Return[W]];\n      ",
  StyleBox["\n(* \
ClearAll[a,b,c,\[Xi],\[Eta],\[Mu],x,y,z,exx,eyy,ezz,gxy,gyz,gzx,\[Omega]xy,\
\[Omega]zx,\[Omega]yz]; \n\[Xi]\[Eta]\[Mu]list={\[Xi],\[Eta],\[Mu]}; \
a=3;b=2;c=1;  ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["SeedRandom[3141];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\nxyznlist={{-a,-b,-c},{a,-b,-c},{a,b,-c},{-a,b,-c},\n          \
{-a,-b,c},{a,-b,c},{a,b,c},{-a,b,c}}/2;        \nFor \
[n=1,n<=8,n++,For[i=1,i<=3,i++,xyznlist[[n,i]]",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["+=Rationalize[100*Random[],001]/400]];",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\nPrint[\"xyznlist=\",xyznlist//MatrixForm]; \n\
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},\n \
{x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8}}=xyznlist;\n\
xyz0={x0,y0,z0}={x1+x2+x3+x4+x5+x6+x7+x8,y1+y2+y3+y4+y5+y6+y7+y8,\n           \
 z1+z2+z3+z4+z5+z6+z7+z8}/8; \nFor \
[n=1,n<=8,n++,For[i=1,i<=3,i++,xyznlist[[n,i]]-=xyz0[[i]] ]];\n\
Print[\"xyznlist=\",xyznlist//MatrixForm]; \n\
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},\n \
{x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8}}=xyznlist;\n\
xyz0={x0,y0,z0}={x1+x2+x3+x4+x5+x6+x7+x8,y1+y2+y3+y4+y5+y6+y7+y8,\n           \
 z1+z2+z3+z4+z5+z6+z7+z8}/8; Print[\"check xyz0=\",xyz0];\n\
Be=IsoHexa8StrainDispMatrix[xyznlist,\[Xi]\[Eta]\[Mu]list,False]; \n\
Be0=Simplify[Be/.{\[Xi]->0,\[Eta]->0,\[Mu]->0}];\n\
Print[\"Be=\",Be//MatrixForm]; Print[\"Be0=\",Be0//MatrixForm]; *)\n\n(*\n\
Hh=CubicHexa8TemplateHhMatrix[{a,b,c}]; Print[\"Hh=\",Hh//MatrixForm];\n\
W=CubicHexa8TemplateWMatrix[{a,b,c}]; Print[\"W=\",W//MatrixForm];\n\
Print[Simplify[W.Hh]//MatrixForm]; *)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Forming and testing Gr,Gc matrices with random rational \
perturbation of brick corner coordinates.
Uses modules of previous cell for checks\
\>", "Text",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ClearAll[a,b,c,\[Xi],\[Eta],\[Mu],x,y,z,exx,eyy,ezz,gxy,gyz,gzx,\[Omega]xy,\
\[Omega]zx,\[Omega]yz]; \n\[Xi]\[Eta]\[Mu]list={\[Xi],\[Eta],\[Mu]}; a=10; \
b=3; c=1; ",
  StyleBox[" SeedRandom[314];",
    FontColor->RGBColor[1, 0, 0]],
  "\nxyznlist={{-a,-b,-c},{a,-b,-c},{a,b,-c},{-a,b,-c},\n          \
{-a,-b,c},{a,-b,c},{a,b,c},{-a,b,c}}/2;        \nFor \
[n=1,n<=8,n++,For[i=1,i<=3,i++,xyznlist[[n,i]]",
  StyleBox["+=Rationalize[100*Random[],001]/400]",
    FontColor->RGBColor[1, 0, 0]],
  "];\nPrint[\"xyznlist=\",xyznlist//MatrixForm]; \n\
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},\n \
{x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8}}=xyznlist;\n\
xyz0={x0,y0,z0}={x1+x2+x3+x4+x5+x6+x7+x8,y1+y2+y3+y4+y5+y6+y7+y8,\n           \
 z1+z2+z3+z4+z5+z6+z7+z8}/8; \nFor \
[n=1,n<=8,n++,For[i=1,i<=3,i++,xyznlist[[n,i]]-=xyz0[[i]] ]];\n\
Print[\"xyznlist=\",xyznlist//MatrixForm]; \n\
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},\n \
{x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8}}=xyznlist;\n\
xyz0={x0,y0,z0}={x1+x2+x3+x4+x5+x6+x7+x8,y1+y2+y3+y4+y5+y6+y7+y8,\n           \
 z1+z2+z3+z4+z5+z6+z7+z8}/8; Print[\"check xyz0=\",xyz0];\n            \n(* \
Assumed rc displacement modes *)\n\
ux=ux0+exx*(x-x0)+(gxy/2+\[Omega]xy)*(y-y0)+(gzx/2-\[Omega]zx)*(z-z0); \n\
uy=uy0+(gxy/2-\[Omega]xy)*(x-x0)+eyy*(y-y0)+(gyz/2+\[Omega]yz)*(z-z0);\n\
uz=uz0+(gzx/2+\[Omega]zx)*(x-x0)+(gyz/2-\[Omega]yz)*(y-y0)+ezz*(z-z0);\n\n\
Print[\"strains=\",Simplify[{D[ux,x],D[uy,y],D[uz,z],\n              \
D[uy,x]+D[ux,y],D[uz,y]+D[uy,z],D[ux,z]+D[uz,x]}]];\n\
Print[\"rotations=\",Simplify[{-D[uy,x]+D[ux,y],\n             \
-D[uz,y]+D[uy,z],-D[ux,z]+D[uz,x]}/2]];\nurc=Table[{ux/.{x-> \
xyznlist[[n,1]],y-> xyznlist[[n,2]],z-> xyznlist[[n,3]]},\n           \
uy/.{x-> xyznlist[[n,1]],y-> xyznlist[[n,2]],z-> xyznlist[[n,3]]},\n          \
 uz/.{x-> xyznlist[[n,1]],y-> xyznlist[[n,2]],z-> xyznlist[[n,3]]}},\n        \
   {n,8}]; urc=Simplify[Flatten[urc]]; \n\
Be=IsoHexa8StrainDispMatrix[xyznlist,\[Xi]\[Eta]\[Mu]list,False]; \n\
Print[\"chk Be.u=\",Simplify[Be.urc]];\n\n(* Testing Gr *)\nGr=Transpose[\n   \
{{1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0},\n    {0,1,0, \
0,1,0, 0,1,0, 0,1,0, 0,1,0, 0,1,0, 0,1,0, 0,1,0},\n    {0,0,1, 0,0,1, 0,0,1, \
0,0,1, 0,0,1, 0,0,1, 0,0,1, 0,0,1},\n    \
{y1,-x1,0,y2,-x2,0,y3,-x3,0,y4,-x4,0,y5,-x5,0,y6,-x6,0,y7,-x7,0,y8,-x8,0},\n  \
  {0,z1,-y1,0,z2,-y2,0,z3,-y3,0,z4,-y4,0,z5,-y5,0,z6,-y6,0,z7,-y7,0,z8,-y8},\n\
    {-z1,0,x1,-z2,0,x2,-z3,0,x3,-z4,0,x4,-z5,0,x5,-z6,0,x6,-z7,0,x7,-z8,0,x8}}\
];\nGrT=Table[0,{6},{24}]; \n\
ezero={exx->0,eyy->0,ezz->0,gxy->0,gyz->0,gzx->0};\n\
GrT[[1]]=(urc/.{ux0->1,uy0->0,uz0->0,\[Omega]xy->0,\[Omega]yz->0,\[Omega]zx->\
0})/.ezero;\n\
GrT[[2]]=(urc/.{ux0->0,uy0->1,uz0->0,\[Omega]xy->0,\[Omega]yz->0,\[Omega]zx->\
0})/.ezero;\n\
GrT[[3]]=(urc/.{ux0->0,uy0->0,uz0->1,\[Omega]xy->0,\[Omega]yz->0,\[Omega]zx->\
0})/.ezero;\n\
GrT[[4]]=(urc/.{ux0->0,uy0->0,uz0->0,\[Omega]xy->1,\[Omega]yz->0,\[Omega]zx->\
0})/.ezero;\n\
GrT[[5]]=(urc/.{ux0->0,uy0->0,uz0->0,\[Omega]xy->0,\[Omega]yz->1,\[Omega]zx->\
0})/.ezero;\n\
GrT[[6]]=(urc/.{ux0->0,uy0->0,uz0->0,\[Omega]xy->0,\[Omega]yz->0,\[Omega]zx->\
1})/.ezero;\nPrint[\"GrT=\",GrT//MatrixForm];\n\
Print[\"Gr'=\",Transpose[Gr]//MatrixForm];\nPrint[\"chk \
Gr=\",Simplify[Transpose[Gr]-GrT]//MatrixForm];\n\n(* Testing Gc *)\n\
Gc=Transpose[\n   {{x1,0,0,x2,0,0,x3,0,0,x4,0,0,x5,0,0,x6,0,0,x7,0,0,x8,0,0},\
\n    {0,y1,0,0,y2,0,0,y3,0,0,y4,0,0,y5,0,0,y6,0,0,y7,0,0,y8,0},\n    \
{0,0,z1,0,0,z2,0,0,z3,0,0,z4,0,0,z5,0,0,z6,0,0,z7,0,0,z8},\n    \
{y1,x1,0,y2,x2,0,y3,x3,0,y4,x4,0,y5,x5,0,y6,x6,0,y7,x7,0,y8,x8,0}/2,\n    \
{0,z1,y1,0,z2,y2,0,z3,y3,0,z4,y4,0,z5,y5,0,z6,y6,0,z7,y7,0,z8,y8}/2,\n    \
{z1,0,x1,z2,0,x2,z3,0,x3,z4,0,x4,z5,0,x5,z6,0,x6,z7,0,x7,z8,0,x8}/2}];\n\
GcT=Table[0,{6},{24}];\n\
rzero={ux0->0,uy0->0,uz0->0,\[Omega]xy->0,\[Omega]yz->0,\[Omega]zx->0};\n\
GcT[[1]]=(urc/.{exx->1,eyy->0,ezz->0,gxy->0,gyz->0,gzx->0})/.rzero;\n\
GcT[[2]]=(urc/.{exx->0,eyy->1,ezz->0,gxy->0,gyz->0,gzx->0})/.rzero;\n\
GcT[[3]]=(urc/.{exx->0,eyy->0,ezz->1,gxy->0,gyz->0,gzx->0})/.rzero;\n\
GcT[[4]]=(urc/.{exx->0,eyy->0,ezz->0,gxy->1,gyz->0,gzx->0})/.rzero;\n\
GcT[[5]]=(urc/.{exx->0,eyy->0,ezz->0,gxy->0,gyz->1,gzx->0})/.rzero;\n\
GcT[[6]]=(urc/.{exx->0,eyy->0,ezz->0,gxy->0,gyz->0,gzx->1})/.rzero;\n\
Print[\"GcT=\",GcT//MatrixForm];\nPrint[\"Gc'=\",Transpose[Gc]//MatrixForm];\n\
Print[\"chk Gc=\",Simplify[Transpose[Gc]-GcT]//MatrixForm];"
}], "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell[CellGroupData[{

Cell["\<\
Forming and testing Hh matrices. Uses modules of first cell to get \
Gr, Gc\
\>", "Text",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[a,b,c,\[Xi],\[Eta],\[Mu],x,y,z,exx,eyy,ezz,gxy,gyz,gzx,\
\[Omega]xy,\[Omega]zx,\[Omega]yz]; 
ClearAll[\[Kappa]xy,\[Kappa]xz,\[Kappa]yx,\[Kappa]yz,\[Kappa]zx,\[Kappa]zy,\
\[Phi]x,\[Phi]y,\[Phi]z,\[Psi]x,\[Psi]y,\[Psi]z,A1,A2,A3,B1,B2,B3];
(* a=10; b=3; c=1; *)
xyznlist={{-a,-b,-c},{a,-b,-c},{a,b,-c},{-a,b,-c},
          {-a,-b,c},{a,-b,c},{a,b,c},{-a,b,c}}/2;
{Gr,Gc}=CubicHexa8TemplateGrcMatrices[{a,b,c}]; 
(* Assumed h displacement modes; *)
uxb=-\[Kappa]xy*x*y-\[Kappa]xz*x*z  +\[Kappa]yx*(y^2/2-b^2/8)           +\
\[Kappa]zx*(z^2/2-c^2/8);
uyb= \[Kappa]xy*(x^2/2-a^2/8)-\[Kappa]yx*y*x   -\[Kappa]yz*y*z                \
  +\[Kappa]zy*(z^2/2-c^2/8);
uzb=         \[Kappa]xz*(x^2/2-a^2/8)   \
+\[Kappa]yz*(y^2/2-b^2/8)-\[Kappa]zx*z*x-\[Kappa]zy*z*y;
(* ux\[Phi]= \[Phi]y*y*z-\[Phi]z*z*y; uy\[Phi]=-\[Phi]x*x*z+\[Phi]z*z*x; uz\
\[Phi]= \[Phi]x*x*y-\[Phi]y*y*x; torsion is NG *)
ux\[Psi]= \[Psi]x*y*z; uy\[Psi]= \[Psi]y*x*z; uz\[Psi]= \[Psi]z*x*y;
ux\[Phi]= \[Phi]x*x*y*z; uy\[Phi]=\[Phi]y*x*y*z; uz\[Phi]= \[Phi]z*x*y*z;
ux=uxb+ux\[Phi]+ux\[Psi]; uy=uyb+uy\[Phi]+uy\[Psi]; uz=uzb+uz\[Phi]+uz\[Psi]; \
 
Print[\"strains=\",Simplify[{D[ux,x],D[uy,y],D[uz,z],
              D[uy,x]+D[ux,y],D[uz,y]+D[uy,z],D[ux,z]+D[uz,x]}]];
uh=Table[{ux/.{x-> xyznlist[[n,1]],y-> xyznlist[[n,2]],z-> xyznlist[[n,3]]},
          uy/.{x-> xyznlist[[n,1]],y-> xyznlist[[n,2]],z-> xyznlist[[n,3]]},
          uz/.{x-> xyznlist[[n,1]],y-> xyznlist[[n,2]],z-> xyznlist[[n,3]]}},
           {n,8}]; uh=Simplify[Flatten[uh]]; 
Hh= Table[0,{12},{24}];
Hh[[ 1]]=uh/.{\[Kappa]xy->1,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 2]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->1,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 3]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->1,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 4]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->1,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 5]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->1,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 6]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->1,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 7]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->1,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 8]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->1,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[ 9]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->1,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->0};
Hh[[10]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->1,\[Phi]y->0,\
\[Phi]z->0};
Hh[[11]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->1,\
\[Phi]z->0};
Hh[[12]]=uh/.{\[Kappa]xy->0,\[Kappa]xz->0,\[Kappa]yx->0,\[Kappa]yz->0,\[Kappa]\
zx->0,\[Kappa]zy->0,\[Psi]x->0,\[Psi]y->0,\[Psi]z->0,\[Phi]x->0,\[Phi]y->0,\
\[Phi]z->1};

Print[\"Hh=\",Hh//MatrixForm];  HH=Simplify[Hh.Transpose[Hh]]; \
Print[\"HH=\",HH//MatrixForm];
Print[NullSpace[Transpose[Hh]]];
Print[\"chk Hh.Gr=0: \",Simplify[Hh.Gr]//MatrixForm];
Print[\"chk Hh.Gc=0: \",Simplify[Hh.Gc]//MatrixForm]; \
\>", "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell["Checking that Grc is OK", "Text",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[a,b,c,\[Xi],\[Eta],\[Mu]]; \[Xi]\[Eta]\[Mu]list={\[Xi],\
\[Eta],\[Mu]}; a=3;b=2;c=1;
xyznlist={{-a,-b,-c},{a,-b,-c},{a,b,-c},{-a,b,-c},
          {-a,-b,c},{a,-b,c},{a,b,c},{-a,b,c}}/2;
{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4},
 {x5,y5,z5},{x6,y6,z6},{x7,y7,z7},{x8,y8,z8}}=xyznlist;
{Bx,By,Bz,Jdet}=IsoHexa8ShapeFunCarDer[xyznlist,\[Xi]\[Eta]\[Mu]list,False];
Bexx=Flatten[Table[{Bx[[n]],0,0},{n,8}]];
Beyy=Flatten[Table[{0,By[[n]],0},{n,8}]];
Bezz=Flatten[Table[{0,0,Bz[[n]]},{n,8}]];
Bgxy=Flatten[Table[{By[[n]],Bx[[n]],0},{n,8}]];
Bgyz=Flatten[Table[{0,Bz[[n]],By[[n]]},{n,8}]];
Bgzx=Flatten[Table[{Bz[[n]],0,Bx[[n]]},{n,8}]];
Be={Bexx,Beyy,Bezz,Bgxy,Bgyz,Bgzx}; 
Be0=Simplify[Be/.{\[Xi]->0,\[Eta]->0,\[Mu]->0}];
(*Print[\"SVD of Be0: \
\",Simplify[Eigensystem[Be0.Transpose[Be0]]]//MatrixForm];*)
Print[\"Be=\",Be//MatrixForm]; Print[\"Be0=\",Be0//MatrixForm]; 
Grc=Transpose[
   {{1,0,0,  1,0,0,  1,0,0,  1,0,0,  1,0,0,  1,0,0,  1,0,0,  1,0,0},
    {0,1,0,  0,1,0,  0,1,0,  0,1,0,  0,1,0,  0,1,0,  0,1,0,  0,1,0},
    {0,0,1,  0,0,1,  0,0,1,  0,0,1,  0,0,1,  0,0,1,  0,0,1,  0,0,1},
    {y1,-x1,0,y2,-x2,0,y3,-x3,0,y4,-x4,0,y5,-x5,0,y6,-x6,0,y7,-x7,0,y8,-x8,0},\

    {0,z1,-y1,0,z2,-y2,0,z3,-y3,0,z4,-y4,0,z5,-y5,0,z6,-y6,0,z7,-y7,0,z8,-y8},\

    {-z1,0,x1,-z2,0,x2,-z3,0,x3,-z4,0,x4,-z5,0,x5,-z6,0,x6,-z7,0,x7,-z8,0,x8},\

    {x1,0,0,x2,0,0,x3,0,0,x4,0,0,x5,0,0,x6,0,0,x7,0,0,x8,0,0},
    {0,y1,0,0,y2,0,0,y3,0,0,y4,0,0,y5,0,0,y6,0,0,y7,0,0,y8,0},
    {0,0,z1,0,0,z2,0,0,z3,0,0,z4,0,0,z5,0,0,z6,0,0,z7,0,0,z8},
    {y1,x1,0,y2,x2,0,y3,x3,0,y4,x4,0,y5,x5,0,y6,x6,0,y7,x7,0,y8,x8,0}/2,
    {0,z1,y1,0,z2,y2,0,z3,y3,0,z4,y4,0,z5,y5,0,z6,y6,0,z7,y7,0,z8,y8}/2,
    {z1,0,x1,z2,0,x2,z3,0,x3,z4,0,x4,z5,0,x5,z6,0,x6,z7,0,x7,z8,0,x8}/2}]; 
Print[\"Grc'.Grc: \",Simplify[Transpose[Grc].Grc]//MatrixForm];
Print[\"Be.Grc=\",Simplify[Be.Grc]//MatrixForm];
Print[\"Be0.Grc=\",Simplify[Be0.Grc]//MatrixForm];
  \
\>", "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell[TextData[StyleBox["Finding possible forms of Hh that satisfy \
orthogonality & rank sufficiency", "Text"]], "Text",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ClearAll[B1,B2,B3,A1,A2,A3,C1,C2,C3,a,b,c]; \nd={{0,0,0},\n   \
{a^2,b^2,c^2},{a^2,b^2,0},{a^2,0,c^2},{0,b^2,c^2},{a^2,0,0},{0,b^2,0},{0,0,c^\
2},\n   {b^2,c^2,a^2},{b^2,c^2,0},{b^2,0,a^2},{0,c^2,a^2},{b^2,0,0},{0,c^2,0},\
{0,0,a^2},\n   \
{c^2,a^2,b^2},{c^2,a^2,0},{c^2,0,b^2},{0,a^2,b^2},{c^2,0,0},{0,a^2,0},{0,0,b^\
2},\n   {a*b,b*c,c*a},{a*b,b*c,0},{a*b,0,c*a},{0,b*c,c*a},{a*b,0,0},{0,b*c,0},\
{0,0,c*a},\n   \
{b*c,c*a,a*b},{b*c,c*a,0},{b*c,0,a*b},{0,c*a,a*b},{b*c,0,0},{0,c*a,0},{0,0,a*\
b},\n   {c*a,a*b,b*c},{c*a,a*b,0},{c*a,0,b*c},{0,a*b,b*c},{c*a,0,0},{0,a*b,0},\
{0,0,b*c}};\n\nmm=Length[d];   Print[\"mm=\",mm]; ksol=0; kmax=10; ABCsol={};\
\nFor[i=1,i<=mm,i++,For [j=1,j<=mm,j++,\n   {A1,A2,A3}=d[[i]]-d[[j]]; \
den=Simplify[A1^2*(A3*B2-A2*B3)^2]; \n   If [den==0,Continue[]]; \n   For \
[k=1,k<=mm,k++,\n      {BB1,B2,B3}=d[[k]]; B1=Simplify[-((A2*B2+A3*B3)/A1)];\n\
      dB1=Simplify[B1-BB1]; If [dB1!=0,Continue[]];\n      \
den=Simplify[A1^2*(A3*B2-A2*B3)^2]; \n      \
{den1,den2,den3,den4}=Simplify[{den/.b->a,den/.c->a,\n                        \
    den/.c->b,den/.{b->a,c->a}}];\n      If \
[den==0||den1==0||den2==0||den3==0||den4==0,Continue[]]; ",
  StyleBox["\n      ",
    FontColor->RGBColor[1, 0, 0]],
  "{C1,C2,C3}= Simplify[{-((A2*B2+A3*B3)/A1),\n         \
-(((A2*B2+A3*B3)*(A2*A3*B2+(A1^2+A3^2)*B3))/\n           \
(A1^2*(A3*B2-A2*B3))),\n           ((A2*B2+A3*B3)*((A1^2+A2^2)*B2+A2*A3*B3))/\
\n           (A1^2*(A3*B2-A2*B3))}];\n      If \
[C1==0&&C2==0&&C3==0,Continue[]];\n      \
{dB1,dC1,dC2,dC3}=Simplify[{Denominator[Together[B1]],\n          \
Denominator[Together[C1]],Denominator[Together[C2]],\n          \
Denominator[Together[C3]]}];\n    ",
  StyleBox["  ",
    FontColor->RGBColor[1, 0, 0]],
  "dBC={{dB1,dC1,dC2,dC3}/.a->0,{dB1,dC1,dC2,dC3}/.b->0,{dB1,dC1,dC2,dC3}/.c->\
0,\n           {dB1,dC1,dC2,dC3}/.{a->0,b->0},{dB1,dC1,dC2,dC3}/.{b->0,c->0},\
\n           {dB1,dC1,dC2,dC3}/.{c->0,a->0}};\n      \
zpos=Position[Flatten[dBC],0]; If [Length[zpos]>0,Continue[]];\n      \
ABC={{C1,C2,C3},{B1,B2,B3},{A1,A2,A3}}; ABClis=Flatten[ABC];\n      \
pos=Position[ABCsol,ABClis]; If [pos=={},ipos=0",
  ",",
  "{{ipos}}=pos];\n      If [ipos>0,Continue[]];\n      \
AppendTo[ABCsol,ABClis]; ",
  StyleBox["Print[\"i,j,k=\",{i,j,k},\" ABC=\",ABC//MatrixForm];",
    FontColor->RGBColor[1, 0, 0]],
  "\n      ksol++; ",
  StyleBox["\n ",
    FontColor->RGBColor[1, 0, 0]],
  "   ]]] ;\n   \n   Print[\"ksol=\",ksol,\" solutions found\"];\n   \
Print[\"ABCsol=\",ABCsol];\n "
}], "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell[CellGroupData[{

Cell["\<\
(* non-invar solution list - discarded *)
ABCsol=\[InvisibleSpace]{{b^2,-2*a^2,-b^2,b^2,0,b^2,a^2,
      b^2,-a^2},{c^2,-c^2,-2*a^2,c^2,c^2,0,a^2,-a^2,
      c^2},{c^2,-c^2,2*b^2,c^2,c^2,0,-b^2,b^2,c^2},{b^2,
      2*c^2,-b^2,b^2,0,b^2,-c^2,b^2,c^2},{c^2,2*a^2,-c^2,c^2,0,
      c^2,a^2,-c^2,-a^2},{b^2,-b^2,2*a^2,b^2,b^2,0,
      a^2,-a^2,-b^2},{a^2,-a^2,-2*b^2,a^2,a^2,0,-b^2,
      b^2,-a^2},{a^2,-2*c^2,-a^2,a^2,0,a^2,-c^2,-a^2,
      c^2},{c^2,-2*b^2,-c^2,c^2,0,c^2,b^2,
      c^2,-b^2},{a^2,-a^2,-2*b^2,a^2,a^2,0,b^2,-b^2,a^2},{c^2,
      2*a^2,-c^2,c^2,0,c^2,-a^2,c^2,a^2},{a^2,-a^2,2*c^2,a^2,a^2,
      0,-c^2,c^2,a^2},{c^2,-c^2,2*b^2,c^2,c^2,0,
      b^2,-b^2,-c^2},{a^2,2*b^2,-a^2,a^2,0,a^2,
      b^2,-a^2,-b^2},{b^2,-b^2,-2*c^2,b^2,b^2,0,-c^2,
      c^2,-b^2},{b^2,-2*a^2,-b^2,b^2,0,b^2,-a^2,-b^2,
      a^2},{a^2,-2*c^2,-a^2,a^2,0,a^2,c^2,
      a^2,-c^2},{b^2,-b^2,-2*c^2,b^2,b^2,0,c^2,-c^2,
      b^2},{b^2,-b^2,2*a^2,b^2,b^2,0,-a^2,a^2,b^2},{a^2,
      2*b^2,-a^2,a^2,0,a^2,-b^2,a^2,b^2},{b^2,2*c^2,-b^2,b^2,0,
      b^2,c^2,-b^2,-c^2},{a^2,-a^2,2*c^2,a^2,a^2,0,
      c^2,-c^2,-a^2},{c^2,-c^2,-2*a^2,c^2,c^2,0,-a^2,
      a^2,-c^2},{c^2,-2*b^2,-c^2,c^2,0,c^2,-b^2,-c^2,
      b^2},{c^2,-2*a*c,-c^2,c^2,0,c^2,a*b,
      b*c,-(a*b)},{a*c,-2*a^2,-(a*c),a*c,0,a*c,a*b,
      b*c,-(a*b)},{b*c,-2*a*b,-(b*c),b*c,0,b*c,a*b,
      b*c,-(a*b)},{a*c,-a^2 - c^2,-c^2,a*c,0,a^2,a*b,
      b*c,-(b*c)},{-(a*c),-c^2,-a^2 - c^2,-(a*c),a^2,0,a*b,
      b*c,-(b*c)},{b*c,-c^2,-b^2 - c^2,b*c,b^2,0,a*b,-(a*c),
      a*c},{-(b*c),-b^2 - c^2,-c^2,-(b*c),0,b^2,a*b,-(a*c),
      a*c},{c^2,-c^2,-2*b*c,c^2,c^2,0,a*b,-(a*b),
      a*c},{b*c,-(b*c),-2*b^2,b*c,b*c,0,a*b,-(a*b),
      a*c},{a*c,-(a*c),-2*a*b,a*c,a*c,0,a*b,-(a*b),a*c},{a^2,-a^2,
      2*a*b,a^2,a^2,0,-(b*c),b*c,a*c},{a*c,-(a*c),2*b*c,a*c,a*c,
      0,-(b*c),b*c,a*c},{a*b,-(a*b),2*b^2,a*b,a*b,0,-(b*c),b*c,
      a*c},{b^2,2*a*b,-b^2,b^2,0,b^2,-(a*c),b*c,a*c},{a*b,
      2*a^2,-(a*b),a*b,0,a*b,-(a*c),b*c,a*c},{b*c,2*a*c,-(b*c),
      b*c,0,b*c,-(a*c),b*c,a*c},{c^2,2*b*c,-c^2,c^2,0,c^2,
      a*b,-(a*c),-(a*b)},{a*c,2*a*b,-(a*c),a*c,0,a*c,
      a*b,-(a*c),-(a*b)},{b*c,2*b^2,-(b*c),b*c,0,b*c,
      a*b,-(a*c),-(a*b)},{c^2,-c^2,2*a*c,c^2,c^2,0,
      a*b,-(a*b),-(b*c)},{b*c,-(b*c),2*a*b,b*c,b*c,0,
      a*b,-(a*b),-(b*c)},{a*c,-(a*c),2*a^2,a*c,a*c,0,
      a*b,-(a*b),-(b*c)},{a^2,-a^2,-2*a*c,a^2,a^2,0,-(b*c),
      b*c,-(a*b)},{a*c,-(a*c),-2*c^2,a*c,a*c,0,-(b*c),
      b*c,-(a*b)},{a*b,-(a*b),-2*b*c,a*b,a*b,0,-(b*c),
      b*c,-(a*b)},{-(a*b),-a^2 - b^2,-b^2,-(a*b),0,a^2,-(a*c),
      b*c,-(b*c)},{a*b,-b^2,-a^2 - b^2,a*b,a^2,0,-(a*c),
      b*c,-(b*c)},{-(a*b),-a^2,-a^2 - b^2,-(a*b),b^2,
      0,-(b*c),-(a*c),a*c},{a*b,-a^2 - b^2,-a^2,a*b,0,
      b^2,-(b*c),-(a*c),a*c},{b^2,-2*b*c,-b^2,b^2,0,
      b^2,-(a*c),-(a*b),a*c},{a*b,-2*a*c,-(a*b),a*b,0,
      a*b,-(a*c),-(a*b),a*c},{b*c,-2*c^2,-(b*c),b*c,0,
      b*c,-(a*c),-(a*b),a*c},{-(a*b),-a^2,-a^2 - b^2,-(a*b),b^2,0,
      b*c,a*c,-(a*c)},{a*b,-a^2 - b^2,-a^2,a*b,0,b^2,b*c,
      a*c,-(a*c)},{a^2,-2*a*b,-a^2,a^2,0,a^2,b*c,
      a*c,-(b*c)},{a*c,-2*b*c,-(a*c),a*c,0,a*c,b*c,
      a*c,-(b*c)},{a*b,-2*b^2,-(a*b),a*b,0,a*b,b*c,
      a*c,-(b*c)},{a^2,-a^2,-2*a*c,a^2,a^2,0,b*c,-(b*c),
      a*b},{a*c,-(a*c),-2*c^2,a*c,a*c,0,b*c,-(b*c),
      a*b},{a*b,-(a*b),-2*b*c,a*b,a*b,0,b*c,-(b*c),
      a*b},{-(a*c),-a^2 - c^2,-a^2,-(a*c),0,c^2,b*c,-(a*b),
      a*b},{a*c,-a^2,-a^2 - c^2,a*c,c^2,0,b*c,-(a*b),a*b},{c^2,
      2*b*c,-c^2,c^2,0,c^2,-(a*b),a*c,a*b},{a*c,2*a*b,-(a*c),a*c,
      0,a*c,-(a*b),a*c,a*b},{b*c,2*b^2,-(b*c),b*c,0,b*c,-(a*b),
      a*c,a*b},{b^2,-b^2,2*b*c,b^2,b^2,0,-(a*c),a*c,
      a*b},{b*c,-(b*c),2*c^2,b*c,b*c,0,-(a*c),a*c,
      a*b},{a*b,-(a*b),2*a*c,a*b,a*b,0,-(a*c),a*c,a*b},{a^2,-a^2,
      2*a*b,a^2,a^2,0,b*c,-(b*c),-(a*c)},{a*c,-(a*c),2*b*c,a*c,
      a*c,0,b*c,-(b*c),-(a*c)},{a*b,-(a*b),2*b^2,a*b,a*b,0,
      b*c,-(b*c),-(a*c)},{a^2,2*a*c,-a^2,a^2,0,a^2,
      b*c,-(a*b),-(b*c)},{a*c,2*c^2,-(a*c),a*c,0,a*c,
      b*c,-(a*b),-(b*c)},{a*b,2*b*c,-(a*b),a*b,0,a*b,
      b*c,-(a*b),-(b*c)},{b*c,-c^2,-b^2 - c^2,b*c,b^2,0,-(a*b),
      a*c,-(a*c)},{-(b*c),-b^2 - c^2,-c^2,-(b*c),0,b^2,-(a*b),
      a*c,-(a*c)},{b^2,-b^2,-2*a*b,b^2,b^2,0,-(a*c),
      a*c,-(b*c)},{b*c,-(b*c),-2*a*c,b*c,b*c,0,-(a*c),
      a*c,-(b*c)},{a*b,-(a*b),-2*a^2,a*b,a*b,0,-(a*c),
      a*c,-(b*c)},{c^2,-2*a*c,-c^2,c^2,0,c^2,-(a*b),-(b*c),
      a*b},{a*c,-2*a^2,-(a*c),a*c,0,a*c,-(a*b),-(b*c),
      a*b},{b*c,-2*a*b,-(b*c),b*c,0,b*c,-(a*b),-(b*c),
      a*b},{b*c,-b^2 - c^2,-b^2,b*c,0,c^2,-(a*c),-(a*b),
      a*b},{-(b*c),-b^2,-b^2 - c^2,-(b*c),c^2,0,-(a*c),-(a*b),
      a*b},{b^2,-2*b*c,-b^2,b^2,0,b^2,a*c,
      a*b,-(a*c)},{a*b,-2*a*c,-(a*b),a*b,0,a*b,a*c,
      a*b,-(a*c)},{b*c,-2*c^2,-(b*c),b*c,0,b*c,a*c,
      a*b,-(a*c)},{b*c,-b^2 - c^2,-b^2,b*c,0,c^2,a*c,
      a*b,-(a*b)},{-(b*c),-b^2,-b^2 - c^2,-(b*c),c^2,0,a*c,
      a*b,-(a*b)},{-(a*b),-a^2 - b^2,-b^2,-(a*b),0,a^2,a*c,-(b*c),
      b*c},{a*b,-b^2,-a^2 - b^2,a*b,a^2,0,a*c,-(b*c),
      b*c},{b^2,-b^2,-2*a*b,b^2,b^2,0,a*c,-(a*c),
      b*c},{b*c,-(b*c),-2*a*c,b*c,b*c,0,a*c,-(a*c),
      b*c},{a*b,-(a*b),-2*a^2,a*b,a*b,0,a*c,-(a*c),b*c},{c^2,-c^2,
      2*a*c,c^2,c^2,0,-(a*b),a*b,b*c},{b*c,-(b*c),2*a*b,b*c,b*c,
      0,-(a*b),a*b,b*c},{a*c,-(a*c),2*a^2,a*c,a*c,0,-(a*b),a*b,
      b*c},{a^2,2*a*c,-a^2,a^2,0,a^2,-(b*c),a*b,b*c},{a*c,
      2*c^2,-(a*c),a*c,0,a*c,-(b*c),a*b,b*c},{a*b,2*b*c,-(a*b),
      a*b,0,a*b,-(b*c),a*b,b*c},{b^2,2*a*b,-b^2,b^2,0,b^2,
      a*c,-(b*c),-(a*c)},{a*b,2*a^2,-(a*b),a*b,0,a*b,
      a*c,-(b*c),-(a*c)},{b*c,2*a*c,-(b*c),b*c,0,b*c,
      a*c,-(b*c),-(a*c)},{b^2,-b^2,2*b*c,b^2,b^2,0,
      a*c,-(a*c),-(a*b)},{b*c,-(b*c),2*c^2,b*c,b*c,0,
      a*c,-(a*c),-(a*b)},{a*b,-(a*b),2*a*c,a*b,a*b,0,
      a*c,-(a*c),-(a*b)},{c^2,-c^2,-2*b*c,c^2,c^2,0,-(a*b),
      a*b,-(a*c)},{b*c,-(b*c),-2*b^2,b*c,b*c,0,-(a*b),
      a*b,-(a*c)},{a*c,-(a*c),-2*a*b,a*c,a*c,0,-(a*b),
      a*b,-(a*c)},{-(a*c),-a^2 - c^2,-a^2,-(a*c),0,c^2,-(b*c),
      a*b,-(a*b)},{a*c,-a^2,-a^2 - c^2,a*c,c^2,0,-(b*c),
      a*b,-(a*b)},{a*c,-a^2 - c^2,-c^2,a*c,0,a^2,-(a*b),-(b*c),
      b*c},{-(a*c),-c^2,-a^2 - c^2,-(a*c),a^2,0,-(a*b),-(b*c),
      b*c},{a^2,-2*a*b,-a^2,a^2,0,a^2,-(b*c),-(a*c),
      b*c},{a*c,-2*b*c,-(a*c),a*c,0,a*c,-(b*c),-(a*c),
      b*c},{a*b,-2*b^2,-(a*b),a*b,0,a*b,-(b*c),-(a*c),b*c}}; 
      
(* investigating candidate Hh matrices *)
ClearAll[a,b,c,C1,C2,C3,B1,B2,B3,A1,A2,A3]; m=Length[ABCsol]; 
For [i=1,i<=m,i++,
  {C1,C2,C3,B1,B2,B3,A1,A2,A3}=ABCsol[[i]]; \
(*Print[{C1,C2,C3,B1,B2,B3,A1,A2,A3}];*)
  Hh789=Simplify[{
       {B1,((A2*A3*B2+(A1^2+A3^2)*B3)*B1)/(A1*(A3*B2-A2*B3)),(A1^2*B2*B1+A2^2*\
B2*B1+
        A2*A3*B3*B1)/(-(A1*A3*B2)+A1*A2*B3),-B1,((A2*A3*B2+(A1^2+A3^2)*B3)*
        B1)/(A1*(-(A3*B2)+A2*B3)),(A1^2*B2*B1+A2^2*B2*B1+A2*A3*B3*B1)/(A1*A3*\
B2-
        A1*A2*B3),B1,((A2*A3*B2+(A1^2+A3^2)*B3)*B1)/(A1*(A3*B2-A2*B3)),(A1^2*\
B2*B1+A2^2*B2*B1+
        A2*A3*B3*B1)/(-(A1*A3*B2)+A1*A2*B3),-B1,((A2*A3*B2+(A1^2+A3^2)*B3)*
        B1)/(A1*(-(A3*B2)+A2*B3)),(A1^2*B2*B1+A2^2*B2*B1+
        A2*A3*B3*B1)/(A1*A3*B2-A1*A2*B3),-B1,((A2*A3*B2+(A1^2+A3^2)*B3)*
        B1)/(A1*(-(A3*B2)+A2*B3)),(A1^2*B2*B1+A2^2*B2*B1+
        A2*A3*B3*B1)/(A1*A3*B2-A1*A2*B3),B1,((A2*A3*B2+(A1^2+A3^2)*B3)*
        B1)/(A1*(A3*B2-A2*B3)),(A1^2*B2*B1+A2^2*B2*B1+A2*A3*B3*B1)/(-(A1*A3*\
B2)+
        A1*A2*B3),-B1,((A2*A3*B2+(A1^2+A3^2)*B3)*B1)/(A1*(-(A3*B2)+A2*B3)),(\
A1^2*B2*B1+A2^2*B2*B1+
        A2*A3*B3*B1)/(A1*A3*B2-A1*A2*B3),B1,((A2*A3*B2+(A1^2+A3^2)*B3)*
        B1)/(A1*(A3*B2-A2*B3)),(A1^2*B2*B1+A2^2*B2*B1+A2*A3*B3*B1)/(-(A1*A3*\
B2)+A1*A2*B3)},
       {-((A2*B2+A3*B3)/A1),B2,B3,(A2*B2+A3*B3)/A1,-B2,-B3,-((A2*B2+A3*B3)/A1)\
,B2,B3,(A2*B2+
        A3*B3)/A1,-B2,-B3,(A2*B2+A3*B3)/A1,-B2,-B3,-((A2*B2+A3*B3)/A1),B2,B3,
        (A2*B2+A3*B3)/A1,-B2,-B3,-((A2*B2+A3*B3)/A1),B2,B3},
       {A1,A2,A3,-A1,-A2,-A3,A1,A2,A3,-A1,-A2,-A3,-A1,-A2,-A3,A1,A2,A3,-A1,-\
A2,-A3,A1,A2,A3}}];
       (*Print[\"Hh789=\",Hh789//MatrixForm]; *)
       DD=Simplify[Hh789.Transpose[Hh789]]; 
       Print[\"orthog chk=\",Table[DD[[i,i]],{i,1,3}]];
       ];\
\>", "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Isoparametric template form (by analytical integration over volume) \
- uses modules of top cell\
\>", "Text",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[a,b,c,\[Xi],\[Eta],\[Mu],x,y,z,Em,\[Nu]]; (* Em=144*6;*) \
(* \[Nu]=0;*) (*b=c=a;*)
Emat=Em/((1+\[Nu])*(1-2*\[Nu]))*\[InvisibleSpace]{{1-\[Nu],\[Nu],\[Nu],0,0,0},\

  {\[Nu],1-\[Nu],\[Nu],0,0,0},{\[Nu],\[Nu],1-\[Nu],0,0,0},{0,0,0,1/2-\[Nu],0,\
0},
  {0,0,0,0,1/2-\[Nu],0},{0,0,0,0,0,1/2-\[Nu]}};  
ue={ux1,uy1,uz1,ux2,uy2,uz2,ux3,uy3,uz3,ux4,uy4,uz4,
    ux5,uy5,uz5,ux6,uy6,uz6,ux7,uy7,uz7,ux8,uy8,uz8};
xyznlist={{-a,-b,-c},{a,-b,-c},{a,b,-c},{-a,b,-c},
          {-a,-b,c},{a,-b,c},{a,b,c},{-a,b,c}}/2;
\[Xi]\[Eta]\[Mu]list={\[Xi],\[Eta],\[Mu]};
Be=IsoHexa8StrainDispMatrix[xyznlist,\[Xi]\[Eta]\[Mu]list,False];
BEB=Simplify[Transpose[Be].Emat.Be]; Print[\"BEB done\"]; 
BEBint=Simplify[Integrate[Integrate[Integrate[
       BEB,{\[Xi],-1,1}],{\[Eta],-1,1}],{\[Mu],-1,1}]]; Print[\"BEBint \
done\"];
J=a*b*c/8;  Ue= (J/2)*ue.BEBint.ue; Print[\"Ue done\"]; 
feiso=Table[D[Ue,ue[[i]] ],{i,24}]; 
Keiso=Table[Table[D[feiso[[i]],ue[[j]]],{i,24}],{j,24}];
Keiso=Simplify[Keiso]; Print[\"Total isoP stiffness done\"];
(*Print[\"eigs of total stiffness=\",Chop[Eigenvalues[N[Keiso]]]];*)
L=CubicHexa8TemplateLMatrix[{a,b,c}]; LT=Transpose[L];
Kb=Simplify[(L.Emat.LT)/(a*b*c)]; Print[\"Basic stiffness done\"];
(*Print[\"eigs of basic stiffness=\",Chop[Eigenvalues[N[Kb]]]];*)
Khiso=Keiso-Kb; Print[\"HO stiffness done\"];
(*Print[\"eigs of HO stiffness=\",Chop[Eigenvalues[N[Khiso]]]];*)
Hh=CubicHexa8TemplateHhMatrix[{a,b,c}]; HhT=Transpose[Hh];
W=CubicHexa8TemplateWMatrix[{a,b,c}]; 
Winv=Inverse[W]; Dhinv=Inverse[Simplify[Hh.HhT]]; 
Riso=Simplify[Dhinv.Winv.(Hh.Khiso.HhT).Winv.Dhinv]; (* get R for isoP brick \
*)
Print[\"Riso=\",Riso//MatrixForm];\
\>", "Input",
  CellFrame->True,
  AspectRatioFixed->True,
  Background->RGBColor[0, 1, 0]]
}, Open  ]]
}, Open  ]]
},
FrontEndVersion->"4.2 for Macintosh",
ScreenRectangle->{{0, 1920}, {0, 1180}},
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{1724, 1144},
WindowMargins->{{12, Automatic}, {Automatic, -5}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
Magnification->1.5,
MacintoshSystemPageSetup->"\<\
00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001
0000I00000400`<300000BL?00400@00000000000000060801T1T00000000000
00000000000000000000000000000000\>"
]

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

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 136, 3, 69, "Text"],
Cell[1893, 56, 8317, 144, 2802, "Input"],
Cell[10213, 202, 242, 7, 93, "Text"],
Cell[10458, 211, 4558, 77, 1482, "Input"],

Cell[CellGroupData[{
Cell[15041, 292, 177, 6, 69, "Text"],
Cell[15221, 300, 3919, 74, 802, "Input"],
Cell[19143, 376, 117, 3, 69, "Text"],
Cell[19263, 381, 2025, 42, 702, "Input"],
Cell[21291, 425, 198, 4, 69, "Text"],
Cell[21492, 431, 2631, 50, 902, "Input"],

Cell[CellGroupData[{
Cell[24148, 485, 8185, 152, 2822, "Input"],
Cell[32336, 639, 198, 6, 69, "Text"],
Cell[32537, 647, 1761, 37, 602, "Input"]
}, Open  ]]
}, Open  ]]
}
]
*)



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