PROGRAMA DE ELEMENTOS FINITOS COMPLETO PARA EL PROBLEMA DE TENSION PLANA - v2018 - m8-a2a

Cell 0.  Mathematica version - dependent settings (to get plots displayed correctly).

DisplayChannel=$DisplayFunction;
If [$VersionNumber>=6.0, DisplayChannel=Print]; 
  (* fix for Mathematica 6 & later *)

Off[General::spell1]; Off[General::spell];

Cell  1.  Modules QuadGaussRuleInfo, TrigGaussRuleInfo and LineGaussRuleInfo provide 
quadrature information for  iso - P quadrilaterals, triangles, and lines, respectively.
    
   For quadrilaterals the rules implemented are 1 x1 through 5 x5.  These are identified by p = 1, ... 5.
See Ch 23 for argument rules. For triangles the rules are p = 1, 3, -3, 6, -6, 7;  see Ch 24.

This code is used by all quadrilateral  stiffness and stress modules, and all triangle modules 
except the 3 - node linear element. Hence it must be initialized before their appearance.

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

TrigGaussRuleInfo[{rule_,numer_},point_]:= Module[
 {zeta,p=rule,i=point,g1,g2,info={{Null,Null,Null},0}  },
  If [p== 1, info={{1/3,1/3,1/3},1}];
  If [p== 3, info={{1,1,1}/6,1/3}; info[[1,i]]=2/3];
  If [p==-3, info={{1,1,1}/2,1/3}; info[[1,i]]=0  ];
  If [p== 6,  g1=(8-Sqrt[10]+Sqrt[38-44*Sqrt[2/5]])/18;
              g2=(8-Sqrt[10]-Sqrt[38-44*Sqrt[2/5]])/18;
     If [i<4, info={{g1,g1,g1},(620+Sqrt[213125-
              53320*Sqrt[10]])/3720}; info[[1,i]]=1-2*g1];
     If [i>3, info={{g2,g2,g2},(620-Sqrt[213125-
              53320*Sqrt[10]])/3720}; info[[1,i-3]]=1-2*g2]];
  If [p== -6, 
     If [i<4, info={{1,1,1}/6,3/10}; info[[1,i]]=2/3];  
     If [i>3, info={{1,1,1}/2,1/30}; info[[1,i-3]]=0]]; 
  If [p== 7,  g1=(6-Sqrt[15])/21; g2=(6+Sqrt[15])/21;
     If [i<4,      info={{g1,g1,g1},(155-Sqrt[15])/1200};
                   info[[1,i]]=  1-2*g1];
     If [i>3&&i<7, info={{g2,g2,g2},(155+Sqrt[15])/1200};
                   info[[1,i-3]]=1-2*g2];
     If [i==7,     info={{1/3,1/3,1/3},9/40} ]];    
  If [numer, Return[N[info]], Return[Simplify[info]]];
];

LineGaussRuleInfo[{rule_,numer_},point_]:= Module[
  {g2={-1,1}/Sqrt[3],w3={5/9,8/9,5/9}, 
   g3={-Sqrt[3/5],0,Sqrt[3/5]}, 
   w4={(1/2)-Sqrt[5/6]/6, (1/2)+Sqrt[5/6]/6,
       (1/2)+Sqrt[5/6]/6, (1/2)-Sqrt[5/6]/6},
   g4={-Sqrt[(3+2*Sqrt[6/5])/7],-Sqrt[(3-2*Sqrt[6/5])/7],
        Sqrt[(3-2*Sqrt[6/5])/7], Sqrt[(3+2*Sqrt[6/5])/7]},
   g5={-Sqrt[5+2*Sqrt[10/7]],-Sqrt[5-2*Sqrt[10/7]],0, 
        Sqrt[5-2*Sqrt[10/7]], Sqrt[5+2*Sqrt[10/7]]}/3,
   w5={322-13*Sqrt[70],322+13*Sqrt[70],512,
       322+13*Sqrt[70],322-13*Sqrt[70]}/900,
   i=point,p=rule,info={Null,0}}, 
  If [p==1, info={0,2}];
  If [p==2, info={g2[[i]],1}];
  If [p==3, info={g3[[i]],w3[[i]]}]; 
  If [p==4, info={g4[[i]],w4[[i]]}];
  If [p==5, info={g5[[i]],w5[[i]]}];
  If [numer, Return[N[info]], Return[Simplify[info]]];
];

Cell 2 A :  Modules for 4 - node bilinear quad in plane stress, 
Uses QuadGaussRuleInfo, which must have been initialized.

Quad4IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= 
  Module[{i,k,p=2,numer=False,h=th,qcoor,c,w,Nf,
    dNx,dNy,Jdet,Be,Ke=Table[0,{8},{8}]},  
  If [Length[options]==1, {numer}=options];
  If [Length[options]==2, {numer,p}=options];
  If [p<1||p>5, Print["p out of range"]; Return[K]];
  For [k=1, k<=p*p, k++, 
       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];
       {Nf,dNx,dNy,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
       If [!numer,{Nf,dNx,dNy,Jdet}=Simplify[{Nf,dNx,dNy,Jdet}]];
       If [Length[th]==4, h=th.Nf]; c=w*Jdet*h;
       Be= {Flatten[Table[{dNx[[i]],       0},{i,4}]],
            Flatten[Table[{0,       dNy[[i]]},{i,4}]],
            Flatten[Table[{dNy[[i]],dNx[[i]]},{i,4}]]}; 
       Ke+=c*Transpose[Be].(Emat.Be);   
      ]; If [!numer,Ke=Simplify[Ke]]; Return[Ke]
   ];

Quad4IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= 
  Module[{i,k,numer=False,qcoor,Nf,dNx,dNy,Jdet,Be,
  qctab,ue=udis,sige=Table[{0,0,0},{4}]},  
  qctab={{-1,-1},{1,-1},{1,1},{-1,1}}; 
  If [Length[options]>0, numer=options[[1]]]; 
  If [Length[udis]==4, ue=Flatten[udis]]; 
  For [k=1, k<=Length[sige], k++, 
       qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; 
       {Nf,dNx,dNy,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
       If [!numer,{Nf,dNx,dNy,Jdet}=Simplify[{Nf,dNx,dNy,Jdet}]];
       Be= {Flatten[Table[{dNx[[i]],       0},{i,4}]],
            Flatten[Table[{0,       dNy[[i]]},{i,4}]],
            Flatten[Table[{dNy[[i]],dNx[[i]]},{i,4}]]}; 
       sige[[k]]=Emat.(Be.ue); 
      ]; If [!numer,sige=Simplify[sige]]; Return[sige]
   ];

Quad4IsoPShapeFunDer[ncoor_,qcoor_]:= Module[
  {Nf,dNx,dNy,dNξ,dNη,i,J11,J12,J21,J22,Jdet,ξ,η,x,y},
  {ξ,η}=qcoor; 
   Nf={(1-ξ)*(1-η),(1+ξ)*(1-η),(1+ξ)*(1+η),(1-ξ)*(1+η)}/4;
   dNξ ={-(1-η), (1-η),(1+η),-(1+η)}/4;
   dNη= {-(1-ξ),-(1+ξ),(1+ξ), (1-ξ)}/4;
   x=Table[ncoor[[i,1]],{i,4}]; y=Table[ncoor[[i,2]],{i,4}];
   J11=dNξ.x; J12=dNξ.y; J21=dNη.x; J22=dNη.y;
   Jdet=Simplify[J11*J22-J12*J21];
   dNx=( J22*dNξ-J12*dNη)/Jdet; dNy=(-J21*dNξ+J11*dNη)/Jdet;
   Return[{Nf,dNx,dNy,Jdet}]
];

Cell 2C:  Modules for 8-node serendipity quad in plane stress,
Uses QuadGaussRuleInfo,which must have been initialized.

Quad8IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= 
  Module[{i,k,p=3,numer=False,h=th,qcoor,c,w,Nf,
    dNx,dNy,Jdet,Be,Ke=Table[0,{16},{16}]},  
  If [Length[options]==1, {numer}=options];
  If [Length[options]==2, {numer,p}=options];
  If [p<1||p>5, Print["p out of range"]; Return[K]];
  For [k=1, k<=p*p, k++,  
       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];
       {Nf,dNx,dNy,Jdet}=Quad8IsoPShapeFunDer[ncoor,qcoor];
        If [Length[th]==8, h=th.Nf]; c=w*Jdet*h;
        Be={Flatten[Table[{dNx[[i]],       0},{i,8}]],
            Flatten[Table[{0,       dNy[[i]]},{i,8}]],
            Flatten[Table[{dNy[[i]],dNx[[i]]},{i,8}]]}; 
        Ke+=c*Transpose[Be].(Emat.Be);   
      ]; If [!numer,Ke=Simplify[Ke]]; Return[Ke]
   ];

Quad8IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= 
  Module[{i,k,numer=False,qcoor,Nf,dNx,dNy,Jdet,Be,
  qctab,ue=udis,sige=Table[{0,0,0},{8}]},  
  qctab={{-1,-1},{1,-1},{1,1},{-1,1},{0,-1},{1,0},{0,1},{-1,0}}; 
  If [Length[options]>0, numer=options[[1]]]; 
  If [Length[udis]==8, ue=Flatten[udis]];
  For [k=1, k<=Length[sige], k++, 
       qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; 
       {Nf,dNx,dNy,Jdet}=Quad8IsoPShapeFunDer[ncoor,qcoor];
        Be={Flatten[Table[{dNx[[i]],       0},{i,8}]],
            Flatten[Table[{0,       dNy[[i]]},{i,8}]],
            Flatten[Table[{dNy[[i]],dNx[[i]]},{i,8}]]}; 
        sige[[k]]=Emat.(Be.ue);   
      ]; If [!numer,sige=Simplify[sige]]; Return[sige]
   ];

Quad8IsoPShapeFunDer[ncoor_,qcoor_]:= Module[
  {Nf,dNx,dNy,dNξ,dNη,J11,J12,J21,J22,Jdet,ξ,η,x,y},
  {ξ,η}=qcoor;  
   Nf= {-(1-ξ)*(1-η)*(1+ξ+η)/4,-(1+ξ)*(1-η)*(1-ξ+η)/4,
        -(1+ξ)*(1+η)*(1-ξ-η)/4,-(1-ξ)*(1+η)*(1+ξ-η)/4,
         (1+ξ)*(1-ξ)*(1-η)/2, (1+ξ)*(1+η)*(1-η)/2,
         (1+ξ)*(1-ξ)*(1+η)/2, (1-ξ)*(1+η)*(1-η)/2};
   dNξ={ (1-η)*(2*ξ+η)/4,   (1-η)*(2*ξ-η)/4, 
         (1+η)*(2*ξ+η)/4,   (1+η)*(2*ξ-η)/4, 
         ξ*(η-1),  (1+η)*(1-η)/2, -ξ*(1+η), -(1+η)*(1-η)/2};
   dNη={ (1-ξ)*(ξ+2*η)/4,   -(1+ξ)*(ξ-2*η)/4, 
         (1+ξ)*(ξ+2*η)/4,   -(1-ξ)*(ξ-2*η)/4, 
        -(1+ξ)*(1-ξ)/2, -(1+ξ)*η, (1+ξ)*(1-ξ)/2,-(1-ξ)*η};
   x=Table[ncoor[[i,1]],{i,8}]; y=Table[ncoor[[i,2]],{i,8}];
   J11=dNξ.x; J21=dNξ.y; J12=dNη.x; J22=dNη.y;
   Jdet=Simplify[J11*J22-J12*J21];
   dNx=( J22*dNξ-J21*dNη)/Jdet; dNy=(-J12*dNξ+J11*dNη)/Jdet;  
   Return[{Nf,dNx,dNy,Jdet}]
];

Cell 2B:  Modules for 9-node biquadratic quad in plane stress,
Uses QuadGaussRuleInfo,which must have been initialized.

Quad9IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= 
  Module[{i,k,p=3,numer=False,h=th,qcoor,c,w,Nf,
    dNx,dNy,Jdet,Be,Ke=Table[0,{18},{18}]},  
  If [Length[options]>0, numer=options[[1]]]; 
  If [Length[options]==2, {numer,p}=options];
  If [p<1||p>5, Print["p out of range"]; Return[K]];
  For [k=1, k<=p*p, k++,  
       {qcoor,w}= QuadGaussRuleInfo[{p,numer},k];
       {Nf,dNx,dNy,Jdet}=Quad9IsoPShapeFunDer[ncoor,qcoor];
        If [Length[th]==9, h=th.Nf]; c=w*Jdet*h;
        Be={Flatten[Table[{dNx[[i]],       0},{i,9}]],
            Flatten[Table[{0,       dNy[[i]]},{i,9}]],
            Flatten[Table[{dNy[[i]],dNx[[i]]},{i,9}]]}; 
        Ke+=c*Transpose[Be].(Emat.Be);   
      ]; If [!numer,Ke=Simplify[Ke]]; Return[Ke]
   ];

Quad9IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= 
  Module[{i,k,numer=False,qcoor,Nf,dNx,dNy,Jdet,Be,
    qctab,ue=udis,sige=Table[{0,0,0},{9}]},  
  qctab={{-1,-1},{1,-1},{1,1},{-1,1},{0,-1},{1,0},{0,1},{-1,0},{0,0}}; 
  If [Length[options]==1, {numer}=options];
  If [Length[options]==2, {numer,p}=options];
  If [Length[udis]==9, ue=Flatten[udis]];
  For [k=1, k<=Length[sige], k++, 
       qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; 
       {Nf,dNx,dNy,Jdet}=Quad9IsoPShapeFunDer[ncoor,qcoor];
        Be={Flatten[Table[{dNx[[i]],       0},{i,9}]],
            Flatten[Table[{0,       dNy[[i]]},{i,9}]],
            Flatten[Table[{dNy[[i]],dNx[[i]]},{i,9}]]}; 
        sige[[k]]=Emat.(Be.ue);   
      ]; If [!numer,sige=Simplify[sige]]; Return[sige]
   ];

Quad9IsoPShapeFunDer[ncoor_,qcoor_]:= Module[
  {Nf,dNx,dNy,dNξ,dNη,J11,J12,J21,J22,Jdet,ξ,η,x,y},
  {ξ,η}=qcoor;  
   Nf= {(1-ξ)*ξ*(1-η)*η/4,-(1+ξ)*ξ*(1-η)*η/4,
        (1+ξ)*ξ*(1+η)*η/4,-(1-ξ)*ξ*(1+η)*η/4,
       -(1-ξ^2)*(1-η)*η/2, (1+ξ)*ξ*(1-η^2)/2,
        (1-ξ^2)*(1+η)*η/2,-(1-ξ)*ξ*(1-η^2)/2,(1-ξ^2)*(1-η^2)};
   dNξ={(1-2*ξ)*η*(1-η)/4,-(1+2*ξ)*η*(1-η)/4,
        (1+2*ξ)*η*(1+η)/4,-(1-2*ξ)*η*(1+η)/4,
         ξ*η*(1-η), (1/2+ξ)*(1-η^2),
        -ξ*η*(1+η),-(1/2-ξ)*(1-η^2), 2*ξ*(η^2-1)};     
   dNη={ξ*(1-ξ)*(1-2*η)/4, -ξ*(1+ξ)*(1-2*η)/4,
        ξ*(1+ξ)*(1+2*η)/4, -ξ*(1-ξ)*(1+2*η)/4,
        -(1/2-η)*(1-ξ^2),  -ξ*(1+ξ)*η,         
         (1/2+η)*(1-ξ^2),   ξ*(1-ξ)*η,2*(ξ^2-1)*η};
   x=Table[ncoor[[i,1]],{i,9}]; y=Table[ncoor[[i,2]],{i,9}];
   J11=dNξ.x; J21=dNξ.y; J12=dNη.x; J22=dNη.y;
   Jdet=Simplify[J11*J22-J12*J21];
   dNx= ( J22*dNξ-J21*dNη)/Jdet;  dNx=Simplify[dNx];
   dNy= (-J12*dNξ+J11*dNη)/Jdet;  dNy=Simplify[dNy];
   Return[{Nf,dNx,dNy,Jdet}]
];

Cell 2C:  Modules for 3-node linear triangle in plane stress.  Self-contained.

Trig3IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:=
  Module[{x1,x2,x3,y1,y2,y3,x21,x13,x32,y12,y31,y23,
  numer=False,h=th,A,Be,Ke}, {{x1,y1},{x2,y2},{x3,y3}}=ncoor;
  If [Length[th]==3, h=(th[[1]]+th[[2]]+th[[3]])/3];
  If [Length[options]>0, numer=options[[1]]]; 
  {x21,x13,x32,y12,y31,y23}={x2-x1,x1-x3,x3-x2,
       y1-y2,y3-y1,y2-y3};
  A=Simplify[x21*y31-x13*y12]/2; If [!numer,A=Simplify[A]];
  Be={{y23,0,y31,0,y12,0},{0,x32,0,x13,0,x21},
      {x32,y23,x13,y31,x21,y12}}/(2*A);
  Ke=A*h*Transpose[Be].Emat.Be; If [!numer, Ke=Simplify[Ke]];
  Return[Ke]];

Trig3IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:=
  Module[{x1,x2,x3,y1,y2,y3,x21,x13,x32,y12,y31,y23,
  numer=False,h=th,A,Be,ue=udis,sige=Table[{0,0,0},{3}]}, 
  {{x1,y1},{x2,y2},{x3,y3}}=ncoor;
  If [Length[th]==3, h=(th[[1]]+th[[2]]+th[[3]])/3];
  If [Length[options]>0, numer=options[[1]]];
  If [Length[udis]==3, ue=Flatten[udis]];  
  {x21,x13,x32,y12,y31,y23}={x2-x1,x1-x3,x3-x2,
       y1-y2,y3-y1,y2-y3};
  A=Simplify[x21*y31-x13*y12]/2; If [!numer,A=Simplify[A]];
  Be={{y23,0,y31,0,y12,0},{0,x32,0,x13,0,x21},
      {x32,y23,x13,y31,x21,y12}}/(2*A);
  sige[[1]]=sige[[2]]=sige[[3]]=Emat.(Be.ue); 
  If [!numer, sige=Simplify[sige]];
  Return[sige]];

Cell 2E:  Modules for 6-node quadratic triangle in plane stress,
Uses TrigGaussRuleInfo,which must have been initialized.

Trig6IsoPMembraneStiffness[ncoor_,Emat_,th_,options_]:= 
 Module[{i,k,p=3,numer=False,h=th,tcoor,w,c,
  Nf,dNx,dNy,Jdet,Be,Ke=Table[0,{12},{12}]},
  If [Length[options]==1, {numer}  =options];
  If [Length[options]==2, {numer,p}=options];
  If [!MemberQ[{1,-3,3,6,7},p], Print["Illegal p"]; Return[Ke]];
  For [k=1, k<=Abs[p], k++,  
      {tcoor,w}= TrigGaussRuleInfo[{p,numer},k];
      {Nf,dNx,dNy,Jdet}= Trig6IsoPShapeFunDer[ncoor,tcoor];
      If [numer, {Nf,dNx,dNy,Jdet}=N[{Nf,dNx,dNy,Jdet}]];
      If [Length[th]==6, h=th.Nf]; c=w*Jdet*h/2;
      Be={Flatten[Table[{dNx[[i]],0       },{i,6}]],
          Flatten[Table[{0,       dNy[[i]]},{i,6}]],
          Flatten[Table[{dNy[[i]],dNx[[i]]},{i,6}]]};
      Ke+=c*Transpose[Be].(Emat.Be); 
      ]; If[!numer,Ke=Simplify[Ke]]; Return[Ke]
 ];

Trig6IsoPMembraneStresses[ncoor_,Emat_,th_,udis_,options_]:= 
  Module[{i,k,p=3,numer=False,h=th,qcoor,Nf,dNx,dNy,Jdet,
  qctab,ue=udis,Be,sige=Table[{0,0,0},{6}]},  
  qctab={{1,0,0},{0,1,0},{0,0,1},{1/2,1/2,0},{0,1/2,1/2},{1/2,0,1/2}};
  If [Length[options]>0, numer=options[[1]]];
  If [Length[udis]==6, ue=Flatten[udis]];
  For [k=1, k<=Length[sige], k++, 
       qcoor=qctab[[k]]; If [numer, qcoor=N[qcoor]]; 
      {Nf,dNx,dNy,Jdet}=Trig6IsoPShapeFunDer[ncoor,qcoor];
       Be={Flatten[Table[{dNx[[i]], 0      },{i,6}]],
           Flatten[Table[{0,       dNy[[i]]},{i,6}]],
           Flatten[Table[{dNy[[i]],dNx[[i]]},{i,6}]]}; 
       sige[[k]]=Emat.(Be.ue); 
      ]; If [!numer, sige=Simplify[sige]]; Return[sige]
 ];

Trig6IsoPShapeFunDer[ncoor_,tcoor_]:= Module[
  {ζ1,ζ2,ζ3,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,
  dx4,dx5,dx6,dy4,dy5,dy6,Jx21,Jx32,Jx13,Jy12,Jy23,Jy31,
  Nf,dNx,dNy,Jdet}, {ζ1,ζ2,ζ3}=tcoor;  
  {{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6}}=ncoor;
  dx4=x4-(x1+x2)/2; dx5=x5-(x2+x3)/2; dx6=x6-(x3+x1)/2; 
  dy4=y4-(y1+y2)/2; dy5=y5-(y2+y3)/2; dy6=y6-(y3+y1)/2;
  Nf={ζ1*(2*ζ1-1),ζ2*(2*ζ2-1),ζ3*(2*ζ3-1),4*ζ1*ζ2,4*ζ2*ζ3,4*ζ3*ζ1};
  Jx21= x2-x1+4*(dx4*(ζ1-ζ2)+(dx5-dx6)*ζ3);
  Jx32= x3-x2+4*(dx5*(ζ2-ζ3)+(dx6-dx4)*ζ1);
  Jx13= x1-x3+4*(dx6*(ζ3-ζ1)+(dx4-dx5)*ζ2);
  Jy12= y1-y2+4*(dy4*(ζ2-ζ1)+(dy6-dy5)*ζ3);
  Jy23= y2-y3+4*(dy5*(ζ3-ζ2)+(dy4-dy6)*ζ1);
  Jy31= y3-y1+4*(dy6*(ζ1-ζ3)+(dy5-dy4)*ζ2);
  Jdet = Jx21*Jy31-Jy12*Jx13;
  dNx= {(4*ζ1-1)*Jy23,(4*ζ2-1)*Jy31,(4*ζ3-1)*Jy12,4*(ζ2*Jy23+ζ1*Jy31),
         4*(ζ3*Jy31+ζ2*Jy12),4*(ζ1*Jy12+ζ3*Jy23)}/Jdet;
  dNy= {(4*ζ1-1)*Jx32,(4*ζ2-1)*Jx13,(4*ζ3-1)*Jx21,4*(ζ2*Jx32+ζ1*Jx13),
         4*(ζ3*Jx13+ζ2*Jx21),4*(ζ1*Jx21+ζ3*Jx32)}/Jdet;
  Return[{Nf,dNx,dNy,Jdet}]
];

Cell 3:  Plane stress assembler  of MET type. Stiffness modules must have been initialized.

PlaneStressMasterStiffness[nodxyz_,eletyp_,elenod_,elemat_,
  elefab_,prcopt_]:=Module[{numele=Length[elenod],
  numnod=Length[nodxyz],ncoor,type,e,enl,neldof,OKtyp,OKenl,
  i,j,n,ii,jj,eft,Emat,th,numer,Ke,K},
  OKtyp={"Bar2","Bar3","Quad4","Quad4.1","Quad8","Quad8.2","Quad9",
         "Quad9.2","Trig3","Trig6","Trig6.-3","Trig10","Trig10.6"};
  OKenl= {2,3,4,4,8,8,9,9,3,6,6,10,10};
  K=Table[0,{2*numnod},{2*numnod}];  numer=prcopt[[1]];
  For [e=1,e<=numele,e++, type=eletyp[[e]];
       If [!MemberQ[OKtyp,type], Print["Illegal type:",type,
       " element e=",e," Assembly aborted"]; Return[Null] ];
       enl=elenod[[e]]; n=Length[enl]; {{j}}=Position[OKtyp,type]; 
       If [OKenl[[j]]!=n, Print ["Wrong node list length, element=",e,
           " Assembly aborted"]; Return[Null] ];
       eft=Flatten[Table[{2*enl[[i]]-1,2*enl[[i]]},{i,1,n}]];
       ncoor=Table[nodxyz[[enl[[i]]]],{i,1,n}]; 
       (*Print["e=",e," type=",type,"  enl=",enl,"  eft=",eft]; *)        
       If [type=="Bar2",     Em=elemat[[e]]; A=elefab[[e]];
           Ke=PlaneBar2Stiffness[ncoor,Em,A,{numer}] ];
       If [type=="Bar3",     Em=elemat[[e]]; A=elefab[[e]];
           Ke=PlaneBar3Stiffness[ncoor,Em,A,{numer}] ];
       If [type=="Quad4",    Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Quad4IsoPMembraneStiffness[ncoor,Emat,th,{numer,2}] ];
       If [type=="Quad4.1",  Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Quad4IsoPMembraneStiffness[ncoor,Emat,th,{numer,1}] ];
       If [type=="Quad8",    Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Quad8IsoPMembraneStiffness[ncoor,Emat,th,{numer,3}] ];
       If [type=="Quad8.2",  Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Quad8IsoPMembraneStiffness[ncoor,Emat,th,{numer,2}] ];
       If [type=="Quad9",    Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Quad9IsoPMembraneStiffness[ncoor,Emat,th,{numer,3}] ];
       If [type=="Quad9.2",  Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Quad9IsoPMembraneStiffness[ncoor,Emat,th,{numer,2}] ];
       If [type=="Trig3",    Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Trig3IsoPMembraneStiffness[ncoor,Emat,th,{numer}] ]; 
       If [type=="Trig6",    Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Trig6IsoPMembraneStiffness[ncoor,Emat,th,{numer,3}] ]; 
       If [type=="Trig6.-3", Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Trig6IsoPMembraneStiffness[ncoor,Emat,th,{numer,-3}] ]; 
       If [type=="Trig10",   Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Trig10IsoPMembraneStiffness[ncoor,Emat,th,{numer,7}] ]; 
       If [type=="Trig10.6", Emat=elemat[[e]]; th=elefab[[e]];
           Ke=Trig10IsoPMembraneStiffness[ncoor,Emat,th,{numer,6}] ];
       neldof=Length[Ke]; (*Print["Ke=",SetPrecision[Chop[Ke],5]//MatrixForm];*)
       For [i=1,i<=neldof,i++, ii=eft[[i]];
           For [j=i,j<=neldof,j++, jj=eft[[j]];
                K[[jj,ii]]=K[[ii,jj]]+=Ke[[i,j]] ];
       ];
     ];  Return[K]; 
];

Cell 4:  Application of displacement BC  on master stiffness equations.  Code takes only single DOF 
constraints, but accepts NZ prescribed displacements. Self contained.

ModifiedMasterStiffness[nodtag_,K_] := Module[
  {i,j,k,n=Length[K],pdof,np,Kmod=K},
  pdof=PrescDisplacementDOFTags[nodtag]; np=Length[pdof];
  For [k=1,k<=np,k++, i=pdof[[k]]; 
      For [j=1,j<=n,j++, Kmod[[i,j]]=Kmod[[j,i]]=0];
      Kmod[[i,i]]=1]; 
  Return[Kmod]];

ModifiedNodeForces[nodtag_,nodval_,K_,f_]:= Module[
  {i,j,k,n=Length[K],pdof,pval,np,d,c,fmod=f},
  pdof=PrescDisplacementDOFTags[nodtag]; np=Length[pdof];
  pval=PrescDisplacementDOFValues[nodtag,nodval]; c=Table[1,{n}]; 
  For [k=1,k<=np,k++, i=pdof[[k]]; c[[i]]=0];
  For [k=1,k<=np,k++, i=pdof[[k]]; d=pval[[k]]; 
      fmod[[i]]=d; If [d==0, Continue[]];
      For [j=1,j<=n,j++, fmod[[j]]-=K[[i,j]]*c[[j]]*d]; 
      ];  
  Return[fmod]];

PrescDisplacementDOFTags[nodtag_]:= Module [
 {j,n,numnod=Length[nodtag],pdof={},k=0,m},
  For [n=1,n<=numnod,n++, m=Length[nodtag[[n]]];
      For [j=1,j<=m,j++, If [nodtag[[n,j]]>0, 
                AppendTo[pdof,k+j]];
         ]; k+=m;
     ]; Return[pdof]];

PrescDisplacementDOFValues[nodtag_,nodval_]:= Module [
 {j,n,numnod=Length[nodtag],pval={},k=0,m},
  For [n=1,n<=numnod,n++, m=Length[nodtag[[n]]];
      For [j=1,j<=m,j++, If [nodtag[[n,j]]>0, 
                AppendTo[pval,nodval[[n,j]]]];
         ]; k+=m;
     ]; Return[pval]];

Cell 5:  Stress recovery in plate elements &  internal force recovery in bars

PlaneStressPlateStresses[nodxyz_,eletyp_,elenod_,elemat_,
  elefab_,noddis_,prcopt_]:=Module[{numele=Length[elenod],
  numnod=Length[nodxyz],ncoor,type,e,enl,i,k,n,Emat,th,
  numer=True,nodpnc,nodsig}, If [Length[prcopt]>0, numer=prcopt[[1]]]; 
  nodpnc=Table[0,{numnod}]; nodsig=Table[{0,0,0},{numnod}];
  If [Length[prcopt]>=1, numer=prcopt[[1]]]; 
  For [e=1,e<=numele,e++, type=eletyp[[e]];
       enl=elenod[[e]]; k=Length[enl]; 
       ncoor=Table[nodxyz[[enl[[i]]]],{i,1,k}];
       ue=Table[{noddis[[enl[[i]],1]],noddis[[enl[[i]],2]]},{i,1,k}];
       ue=Flatten[ue]; 
       (*Print["e=",e," type=",type,"  enl=",enl,"  ue=",ue]; *)        
       If [StringTake[type,3]=="Bar", Continue[]];
       If [type=="Quad4",    Emat=elemat[[e]]; th=elefab[[e]];
           sige=Quad4IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       If [type=="Quad4.1",  Emat=elemat[[e]]; th=elefab[[e]];
           sige=Quad4IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       If [type=="Quad8",    Emat=elemat[[e]]; th=elefab[[e]];
           sige=Quad8IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       If [type=="Quad8.2",  Emat=elemat[[e]]; th=elefab[[e]];
           sige=Quad8IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       If [type=="Quad9",    Emat=elemat[[e]]; th=elefab[[e]];
           sige=Quad9IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       If [type=="Quad9.2",  Emat=elemat[[e]]; th=elefab[[e]];
           sige=Quad9IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       If [type=="Trig3",    Emat=elemat[[e]]; th=elefab[[e]];
           sige=Trig3IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; 
       If [type=="Trig6",    Emat=elemat[[e]]; th=elefab[[e]];
           sige=Trig6IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; 
       If [type=="Trig6.-3", Emat=elemat[[e]]; th=elefab[[e]];
           sige=Trig6IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; 
       If [type=="Trig10",   Emat=elemat[[e]]; th=elefab[[e]];
           sige=Trig10IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ]; 
       If [type=="Trig10.6", Emat=elemat[[e]]; th=elefab[[e]];
           sige=Trig10IsoPMembraneStresses[ncoor,Emat,th,ue,{numer}] ];
       (*Print["sige=",sige//MatrixForm]; *)     
       Do [n=enl[[i]]; nodpnc[[n]]++; nodsig[[n]]+=sige[[i]], {i,1,k}]; 
      ];
  Do [k=nodpnc[[n]]; If [k>1,nodsig[[n]]/=k], {n,1,numnod}];
  If [numer, nodsig=Chop[nodsig]];
  Return[{nodpnc,nodsig}]; 
];

PlaneStressBarForces[nodxyz_,eletyp_,elenod_,elemat_,
  elefab_,noddis_,prcopt_]:=Module[{numele=Length[elenod],
  numnod=Length[nodxyz],ncoor,type,e,enl,
  i,k,n,Em,A,numer=True,pe,barele={},barfor={}},
  If [Length[prcopt]>=1, numer=prcopt[[1]] ];
  For [e=1,e<=numele,e++, type=eletyp[[e]];
       If [StringTake[type,3]!="Bar", Continue[]];
       enl=elenod[[e]]; k=Length[enl]; 
       ncoor=Table[nodxyz[[enl[[i]]]],{i,1,k}];
       ue=Table[{noddis[[enl[[i]],1]],noddis[[enl[[i]],2]]},{i,1,k}];
       ue=Flatten[ue]; 
       (*Print["e=",e," type=",type,"  enl=",enl,"  ue=",ue];*)     
       If [type=="Bar2", bar=True; Em=elemat[[e]]; A=elefab[[e]];
           pe=PlaneBar2IntForce[ncoor,Em,A,ue,{numer}] ];
       If [type=="Bar3", bar=True; Em=elemat[[e]]; A=elefab[[e]];
           pe=PlaneBar3IntForce[ncoor,Em,A,ue,{numer}] ];
       AppendTo[barele,e]; AppendTo[barfor,pe]; 
      ];
   If [numer, nodsig=Chop[barfor]];
   Return[{barele,barfor}]; 
];

Cell 6:  Mesh plot modules.  Written by C. A. Felippa,1997/98.  Mathematica adaptation of ancient  (1966) Fortran code.

Plot2DElements[xycoor_,elenod_,aspect_,label_]:= 
Module[{enl,n,nc,xyc,poly=sides={},numele=Length[elenod]},
 Do [enl=elenod[[ne]]; nc=Length[enl]; 
    If [nc!=3 && nc!=4,Continue[]]; xyc=Table[0,{nc+1}];
    enl=AppendTo[enl,enl[[1]]];    
    Do [n=enl[[i]]; xyc[[i]]=xycoor[[n]],{i,1,nc+1}];
    AppendTo[poly,Graphics[Polygon[Take[xyc,nc]]]];
    AppendTo[sides,Graphics[Line[xyc]]],
 {ne,1,numele}];
 Show[Graphics[RGBColor[1,1,0]],poly,
      Graphics[Thickness[0.008]],Graphics[RGBColor[0,0,0]],
      sides,AspectRatio->aspect,PlotLabel->label,
      DisplayFunction->DisplayChannel ];
 ClearAll[poly,sides];
 ];

Plot2DElementsAndNodes[xycoor_,elenod_,aspect_,label_,
  elabels_,nlabels_]:= Module[{i,j,k,e,enl,elab,n,nlab,
   x,y,delx,dely,xmin,xmax,ymin,ymax,xy0,eps,ex,ey,b,nb,
   numele=Length[elenod],rade,gside,gpoly,style,
   numnod=Length[xycoor],radn,pnod=pnlab={},
   pefill=pesid=pelab=pecirc=pedisk={}},
   x=Table[xycoor[[i,1]],{i,1,numnod}]; 
   y=Table[xycoor[[i,2]],{i,1,numnod}];
   {xmin,xmax,ymin,ymax}=N[{Min[x],Max[x],Min[y],Max[y]}];
   {delx,dely}={xmax-xmin,ymax-ymin}/Sqrt[N[numnod]];
   radn=0.07*Min[delx,dely]; {ex,ey}={-2.3*radn,1.25*radn}; 
   rade=0.16*Min[delx,dely]; 
   (*Print["delx,dely,rade=",{delx,dely,rade}];*) 
   style=TextStyle->{FontFamily->"Times",FontSize->11,
                     FontWeight->"Plain",FontSlant->"Plain"};
   For [e=1,e<=numele,e++, enl=elenod[[e]]; n=Length[enl];  nb=0;
       If [n==3, b={1,2,3};             nb=3];
       If [n==4, b={1,2,3,4};           nb=4];
       If [n==6, b={1,4,2,5,3,6};       nb=6];
       If [n==8, b={1,5,2,6,3,7,4,8};   nb=8];
       If [n==9, b={1,5,2,6,3,7,4,8};   nb=8];
       If [n==10,b={1,4,5,2,6,7,3,8,9}; nb=9];
       If [nb==0, Continue[]];
       gside=gpoly=Table[xycoor[[enl[[b[[i]]]]]],{i,nb}]; 
       AppendTo[gside,gpoly[[1]]];
       xy0={0,0};Do [xy0+=gside[[i]],{i,nb}];xy0=xy0/nb; 
       AppendTo[pefill,Graphics[Polygon[gpoly]]];
       AppendTo[pesid, Graphics[Line[gside]]];
       If [n==9,xy0+={0,-2*rade}];
       If [elabels, AppendTo[pedisk,Graphics[Disk[xy0,rade]]];
           AppendTo[pecirc,Graphics[Circle[xy0,rade]]];
           elab=ToString[e];
           AppendTo[pelab,Graphics[Text[elab,xy0,style]]]];
          ];
   For [n=1,n<=numnod,n++,
       If [nlabels,  nlab=ToString[n]; AppendTo[
             pnlab,Graphics[Text[nlab,xycoor[[n]]+{ex,ey},style]]]];
       AppendTo[pnod, Graphics[Disk[xycoor[[n]],radn]]]
       ];  eps=0.1*(ymax-ymin);
  Show[Graphics[RGBColor[1,1,0]],pefill,
       Graphics[GrayLevel[0]],
       Graphics[AbsoluteThickness[1.5]],pesid,
       Graphics[GrayLevel[0]],pnod,
       Graphics[AbsoluteThickness[0.8]],
       Graphics[GrayLevel[1]],pedisk,
       Graphics[GrayLevel[0]],pecirc,pelab,
       Graphics[GrayLevel[0]],pnlab,
       PlotLabel->label,
       PlotRange->{{xmin-eps,xmax+eps},{ymin-eps,ymax+eps}},
       Axes->False,AspectRatio->aspect,
       DisplayFunction->DisplayChannel];
 ClearAll[pefill,pesid,pnod,pedisk,pnlab];
 ];

Cell 7:  Contour plot modules.  Written by C. A. Felippa, 1997, redone 2008. 
Mathematica adaptation of ancient  (1966) code.

INICIALIZACIONES

DisplayChannel=$DisplayFunction;
If [$VersionNumber>=6.0, DisplayChannel=Print]; 
  (* fix for Mathematica 6 & later *)

SELECCION COLORES BANDAS

ContourBandColor[f_,fmin_,fmax_]:= Module[{fs,h,s,b=1},
   If [fmax==0 ||fmin>=fmax, (* White if wrong inputs *)
       Return[Graphics[RGBColor[1,1,1]] ]];
   If [f>fmax || f<fmin, (* Black if outside range *)
       Return[Graphics[RGBColor[0,0,0]] ]]; 
   fs=(fmax-f)/(fmax-fmin); h=0.70*fs; s=1-1/(1+100*(fs-1/2)^2); 
   Return[Graphics[Hue[h,s,b]]]];

SELECCION BANDAS EN TRIANGULOS

TrigBandSegment[xs_,ys_,fs_,fv_,eps_]:=Module[{x1,x2,x3,
  y1,y2,y3,f1,f2,f3,x21,y21,x32,y32,x13,y13,f21,f32,f31,
  ζ211,ζ212,ζ311,ζ313,ζ322,ζ323,P1,P3,P21,P32,P13}, 
  {x1,x2,x3}=xs; {y1,y2,y3}=ys; {f1,f2,f3}=fs; P1={x1,y1}; P3={x3,y3};
  If [fv<f1, Return[{{P1,P1},-1}]]; If [fv>f3, Return[{{P3,P3},0}]];
  f21=f2-f1; If [f21<eps,f1=f1-eps;f21=eps];
  f32=f3-f2; If [f32<eps,f3=f2+eps;f32=eps]; f31=f3-f1; 
  ζ212=(fv-f1)/f21; ζ212=Max[Min[ζ212,1],0]; ζ211=1-ζ212;
  ζ313=(fv-f1)/f31; ζ313=Max[Min[ζ313,1],0]; ζ311=1-ζ313;
  P21={x2*ζ212+x1*ζ211,y2*ζ212+y1*ζ211};
  P13={x1*ζ311+x3*ζ313,y1*ζ311+y3*ζ313};
  If [fv<=f2,Return[{{P21,P13},1}]];
  ζ323=(fv-f2)/f32; ζ323=Max[Min[ζ323,1],0]; ζ322=1-ζ323; 
  P32={x3*ζ323+x2*ζ322,y3*ζ323+y2*ζ322};
  Return[{{P32,P13},2}]];

VISUALIZACION BANDAS SOBRE TRIANGULOS 3 NODOS

ContourBandPlotFunctionOnTrig3[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
  lines_]:=Module[{nv,iv,xc1,yc1,xc2,yc2,xc3,yc3,fc1,fc2,fc3,
  xs,ys,fs,x1,x2,x3,y1,y2,y3,f1,f2,f3,firsthit,Pc1,Pc2,Pc3,
  ftop,fbot,favg,Pb,Pt,tb,tt,P1,P2,P3,P4,Q1,Q2,Q3,Q4,ptab,
  flast,Plast,tlast,kb,kl,p,pb,pl}, 
  If [fmax<fmin||finc<=0, Return[{{},{}}]];
  nv=Floor[(fmax-fmin+eps)/finc]; nv=Min[nv,1000]; 
  pb=Table[{},{nv}]; pl=Table[{},{nv+1}]; 
 {{xc1,yc1},{xc2,yc2},{xc3,yc3}}=N[xyc]; {fc1,fc2,fc3}=N[fc];
  {{x1,y1,f1},{x2,y2,f2},{x3,y3,f3}}=Sort[{{xc1,yc1,fc1},
          {xc2,yc2,fc2},{xc3,yc3,fc3}},#1[[3]]<=#2[[3]]&]; 
  xs={x1,x2,x3}; ys={y1,y2,y3}; fs={f1,f2,f3}; 
  Pc1={x1,y1}; Pc2={x2,y2}; Pc3={x3,y3};
  ptab={{{},{Pc1,Pc2,Pc3},{Pc1,Q4,Q3},{Pc1,Pc2,Q3,Q4}},
     {{Null},{},{Null},{Null}},
     {{Null},{Pc3,Pc2,Q1,Q2},{Q1,Q2,Q4,Q3},{Pc2,Q1,Q2,Q4,Q3}},
     {{Null},{Pc3,Q1,Q2},{Null},{Q1,Q2,Q4,Q3}}};
  flast=fmin-1000; firsthit=True; kb=kl=0;
  For [iv=1,iv<=nv,iv++, 
       fbot=N[fmin+(iv-1)*finc]; ftop=N[fmin+iv*finc];
       If [ftop<=f1 || fbot>=f3, Continue[]]; favg=(fbot+ftop)/2;
       If [fbot==flast, {fbot,Pb,tb}={flast,Plast,tlast},
          {Pb,tb}=TrigBandSegment[xs,ys,fs,fbot,eps],
          {Pb,tb}=TrigBandSegment[xs,ys,fs,fbot,eps]]; 
       {Pt,tt}=TrigBandSegment[xs,ys,fs,ftop,eps];
       {P1,P2}=Pb; {P3,P4}=Pt; {flast,Plast,tlast}={ftop,Pt,tt};
       If [lines,
          If [firsthit&&(tb>0),  
              pl[[++kl]]=Graphics[Line[Pb]]; firsthit=False];
	      If [tt>0, pl[[++kl]]=Graphics[Line[Pt]]]
	      ];
       p=ptab[[tb+2,tt+2]]/.{Q1->P1,Q2->P2,Q3->P3,Q4->P4};
       pb[[++kb]]={ContourBandColor[favg,fmin,fmax],Graphics[Polygon[p]]};
	   ];  
  If [kb>0,pb=Take[pb,kb],{},{}]; If [kl>0,pl=Take[pl,kl],{},{}];
  ClearAll[ptab]; Return[{pb,pl}];
 ];

VISUALIZACION BANDAS SOBRE CUADRILATEROS 4 NODOS

ContourBandPlotFunctionOnQuad4[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
   lines_]:=Module[{x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,f1,f2,f3,f4,f5,
   e,xytab,ftab,xyt,ft, pbq,plq,pbt,plt },
   pbq=plq=Table[{},{4}];
   {{x1,y1},{x2,y2},{x3,y3},{x4,y4}}=xyc; {f1,f2,f3,f4}=fc;
   x5=(x1+x2+x3+x4)/4; y5=(y1+y2+y3+y4)/4; f5=(f1+f2+f3+f4)/4;
   xytab={{{x1,y1},{x2,y2},{x5,y5}},{{x2,y2},{x3,y3},{x5,y5}},
          {{x3,y3},{x4,y4},{x5,y5}},{{x4,y4},{x1,y1},{x5,y5}}};
   ftab={{f1,f2,f5},{f2,f3,f5},{f3,f4,f5},{f4,f1,f5}};
   For [e=1,e<=4,e++, xyt=xytab[[e]]; ft=ftab[[e]];
       {pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnTrig3[
            xyt,ft,{fmin,fmax,finc,eps},lines];
       ];       
    Return[{pbq,plq}];
	];

VISUALIZACION BANDAS SOBRE CUADRILATEROS 9 NODOS

ContourBandPlotFunctionOnQuad9[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
   lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,y1,y2,y3,y4,
   y5,y6,y7,y8,y9,f1,f2,f3,f4,f5,f6,f7,f8,f9,
   e,xytab,ftab,xyt,ft, pbq,plq,pbt,plt },
   pbq=plq=Table[{},{4}];
   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},
    {x5,y5},{x6,y6},{x7,y7},{x8,y8},{x9,y9}}=xyc; 
    {f1,f2,f3,f4,f5,f6,f7,f8,f9}=fc;
   xytab={{{x1,y1},{x5,y5},{x9,y9},{x8,y8}},
          {{x2,y2},{x6,y6},{x9,y9},{x5,y5}},
          {{x3,y3},{x7,y7},{x9,y9},{x6,y6}},
          {{x4,y4},{x8,y8},{x9,y9},{x7,y7}}};
   ftab={{f1,f5,f9,f8},{f2,f6,f9,f5},{f3,f7,f9,f6},{f4,f8,f9,f7}};
   For [e=1,e<=4,e++, xye=xytab[[e]]; fe=ftab[[e]];
       {pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnQuad4[
             xye,fe,{fmin,fmax,finc,eps},lines];
       ];       
    Return[{pbq,plq}];
	];

VISUALIZACION BANDAS SOBRE CUADRILATEROS 8 NODOS

ContourBandPlotFunctionOnQuad8[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
   lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,y1,y2,y3,y4,
   y5,y6,y7,y8,y9,f1,f2,f3,f4,f5,f6,f7,f8,f9,
   e,xytab,ftab,xyt,ft, pbq,plq,pbt,plt },
   pbq=plq=Table[{},{4}];
   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},
    {x5,y5},{x6,y6},{x7,y7},{x8,y8}}=xyc; 
    {f1,f2,f3,f4,f5,f6,f7,f8}=fc;
   x9=(2*(x5+x6+x7+x8)-(x1+x2+x3+x4))/4;
   y9=(2*(y5+y6+y7+y8)-(y1+y2+y3+y4))/4;
   f9=(2*(f5+f6+f7+f8)-(f1+f2+f3+f4))/4;
   xytab={{{x1,y1},{x5,y5},{x9,y9},{x8,y8}},
          {{x2,y2},{x6,y6},{x9,y9},{x5,y5}},
          {{x3,y3},{x7,y7},{x9,y9},{x6,y6}},
          {{x4,y4},{x8,y8},{x9,y9},{x7,y7}}};
   ftab={{f1,f5,f9,f8},{f2,f6,f9,f5},{f3,f7,f9,f6},{f4,f8,f9,f7}};
   For [e=1,e<=4,e++, xye=xytab[[e]]; fe=ftab[[e]];
       {pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnQuad4[
                   xye,fe,{fmin,fmax,finc,eps},lines];
       ];       
    Return[{pbq,plq}];
	];

VISUALIZACION BANDAS SOBRE CUADRILATEROS 16 NODOS

ContourBandPlotFunctionOnQuad16[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
   lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,
   x11,x12,x13,x14,x15,x16,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,
   y11,y12,y13,y14,y15,y16,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,
   f11,f12,f13,f14,f15,f16,
   e,xytab,ftab,xyt,ft, pbq,plq,pbt,plt },
   pbq=plq=Table[{},{9}];
   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},
    {x5,y5},{x6,y6},{x7,y7},{x8,y8},
    {x9,y9},{x10,y10},{x11,y11},{x12,y12},
    {x13,y13},{x14,y14},{x15,y15},{x16,y16}}=xyc; 
   {f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,
    f13,f14,f15,f16}=fc;
   xytab={{{x1,y1},  {x5,y5},  {x13,y13},{x12,y12}},
          {{x2,y2},  {x7,y7},  {x14,y14},{x6,y6}},
          {{x3,y3},  {x9,y9},  {x15,y15},{x8,y8}},
          {{x4,y4},  {x11,y11},{x16,y16},{x10,y10}},
          {{x5,y5},  {x6,y6},  {x14,y14},{x13,y13}},
          {{x7,y7},  {x8,y8},  {x15,y15},{x14,y14}},
          {{x9,y9},  {x10,y10},{x16,y16},{x15,y15}},
          {{x11,y11},{x12,y12},{x13,y13},{x16,y16}},
          {{x13,y13},{x14,y14},{x15,y15},{x16,y16}}};
   ftab={{f1,f5,f13,f12},{f2,f7,f14,f6},{f3,f9,f15,f8},
        {f4,f11,f16,f10},{f5,f6,f14,f13},{f7,f8,f15,f14},
        {f9,f10,f16,f15},{f11,f12,f13,f16},{f13,f14,f15,f16}};
   For [e=1,e<=9,e++, xye=xytab[[e]]; fe=ftab[[e]];
       {pbq[[e]],plq[[e]]}=ContourBandPlotFunctionOnQuad4[
                xye,fe,{fmin,fmax,finc,eps},lines];
       ];       
    Return[{pbq,plq}];
	];

VISUALIZACION BANDAS SOBRE TRIANGULOS 6 NODOS

ContourBandPlotFunctionOnTrig6[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
   lines_]:=Module[{x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,
   f1,f2,f3,f4,f5,f6,e,xytab,ftab,xyt,ft,pbt,plt,pbe,ple},
   pbt=plt=Table[{},{4}];
   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6}}=xyc; 
   {f1,f2,f3,f4,f5,f6}=fc;
   xytab={{{x1,y1},{x4,y4},{x6,y6}},{{x2,y2},{x5,y5},{x4,y4}},
          {{x3,y3},{x6,y6},{x5,y5}},{{x4,y4},{x5,y5},{x6,y6}}};
   ftab={{f1,f4,f6},{f2,f5,f4},{f3,f6,f5},{f4,f5,f6}};
   For [e=1,e<=4,e++, xyt=xytab[[e]]; ft=ftab[[e]];
       {pbt[[e]],plt[[e]]}=ContourBandPlotFunctionOnTrig3[
             xyt,ft,{fmin,fmax,finc,eps},lines];
       ];       
    Return[{pbt,plt}];
	];

VISUALIZACION BANDAS SOBRE TRIANGULOS 10 NODOS

ContourBandPlotFunctionOnTrig10[xyc_,fc_,{fmin_,fmax_,finc_,eps_},
   lines_]:=Module[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x0,y1,y2,y3,
   y4,y5,y6,y7,y8,y9,y0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f0,
   e,xytab,ftab,xyt,ft,pbt,plt},
   pbt=plt=Table[{},{9}]; 
   {{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6},
    {x7,y7},{x8,y8},{x9,y9},{x0,y0}}=xyc; 
   {f1,f2,f3,f4,f5,f6,f7,f8,f9,f0}=fc;
   xytab={{{x1,y1},{x4,y4},{x9,y9}},{{x2,y2},{x6,y6},{x5,y5}},
          {{x3,y3},{x8,y8},{x7,y7}},{{x9,y9},{x4,y4},{x0,y0}},
          {{x4,y4},{x5,y5},{x0,y0}},{{x5,y5},{x6,y6},{x0,y0}},
          {{x6,y6},{x7,y7},{x0,y0}},{{x7,y7},{x8,y8},{x0,y0}},
          {{x8,y8},{x9,y9},{x0,y0}}};
   ftab={{f1,f4,f9},{f2,f6,f5},{f3,f8,f7},{f9,f4,f0},
         {f4,f5,f0},{f5,f6,f0},{f6,f7,f0},{f7,f8,f0},
         {f8,f9,f0}};
   For [e=1,e<=9,e++, xyt=xytab[[e]]; ft=ftab[[e]];
       {pbt[[e]],plt[[e]]}=ContourBandPlotFunctionOnTrig3[
            xyt,ft,{fmin,fmax,finc,eps},lines];
       ];     
    Return[{pbt,plt}];
	];

VISUALIZACION LEYENDA CONTORNO DE BANDAS

ContourPlotValueLegend[fmin_,fmax_,xyleg_,mv_]:=Module[
  {f,fs,finc,fa=Abs[fmin],fb=Abs[fmax],xylinL,xylinR,xytext,
   black=Graphics[RGBColor[0,0,0]],baselineskip=10,
   white=Graphics[RGBColor[1,1,1]],preleg,plab,iv,
   thick=Graphics[AbsoluteThickness[10]],xf1,xf2,xf3,xf4,
   thin= Graphics[AbsoluteThickness[1]],color,padOK},
   xf1=Offset[{-30,  10},xyleg]; xf2=Offset[{85,10},xyleg];
   xf3=Offset[{ 85,-10-baselineskip*(mv+1)},xyleg];
   xf4=Offset[{-30,-10-baselineskip*(mv+1)},xyleg];
   preleg={thin,white,
       Graphics[Polygon[{xf1,xf2,xf3,xf4}]], black,
       Graphics[Line[{xf1,xf2,xf3,xf4,xf1}]],thick,
       Graphics[Text["LEGEND",Offset[{5,0},xyleg],{-1,0}]]};
   finc=(fmax-fmin)/mv; plab=Table[{},{mv+1}]; 
   padOK=(Max[fa,fb]<=1000)&&(Min[fa,fb]>=1/100);
   For [iv=1,iv<=mv+1,iv++, f=N[fmin+(iv-1)*finc]; 
        xylinL=Offset[{-20,-baselineskip*iv},xyleg];
        xylinR=Offset[{  5,-baselineskip*iv},xyleg]; 
        xytext=Offset[{ 80,-baselineskip*iv},xyleg];
        If [padOK, fs=PaddedForm[f,4,NumberSigns->{"-","+"}],
               fs=ScientificForm[f,4,NumberSigns->{"-","+"}]];
        color=ContourBandColor[f,fmin,fmax]; 
        plab[[iv]]={color,Graphics[Line[{xylinL,xylinR}]],
                    black,Graphics[Text[fs,xytext,{1,0}]]}]; 
   Return[{preleg,plab}]];

VISUALIZACION CONTORNO DE BANDAS SOBRE MALLA DOMINIO PLANO

ContourBandPlotNodeFuncOver2DMesh[nodxyz_,elenod_,fn_,{fmin_,fmax_,finc_},
  {colors_,mesh_,nodes_,bacg_,lines_,bound_},xyleg_,aspect_,label_]:= Module[
  {enl,i,k,n,ne,xyc,fc,c,nn,pb,pl,eps,prebacg,pbacg,preband,pband,preline,
  bnt,bn,nbn,enb,xyb,bseg,ebs,seg,rseg,nbseg,prebound,pbound,
  pline,premesh,pmesh,prenode,preleg,pleg,pnode,numele=Length[elenod]},
  bnt={{},{},{1,2,3,1},{1,2,3,4,1},{},{1,4,2,5,3,6,1},
      {},{1,5,2,6,3,7,4,8,1},{1,5,2,6,3,7,4,8,1},
      {1,4,5,2,6,7,3,8,9,1},{},{},{},{},{},{1,5,6,2,7,8,3,9,10,4,11,12,1}}; 
  eps=10.^(-8)*Max[Abs[fmax],Abs[fmin]]; pband=pmesh=pnode=pline=pbacg={}; 
  If [colors,pband=Table[{},{numele}]]; If [mesh, pmesh=Table[{},{numele}]]; 
  If [nodes, pnode=Table[{},{numele}]]; If [lines,pline=Table[{},{numele}]]; 
  If [bacg,  pbacg=Table[{},{numele}]]; preleg=pleg=bseg=pbound={};
  prebacg= {Graphics[RGBColor[.7,.7,.7]]}; preband={};
  preline= {Graphics[RGBColor[0,0,0]],Graphics[AbsoluteThickness[1]]};
  premesh= {Graphics[RGBColor[1,0,0]],Graphics[AbsoluteThickness[2]],
            Graphics[AbsoluteDashing[{1,10}]]};
  prenode= {Graphics[RGBColor[0,0,0]],Graphics[AbsoluteDashing[{}]],
            Graphics[AbsolutePointSize[5]]};
  prebound={Graphics[RGBColor[0,0,1]],Graphics[AbsoluteDashing[{}]],
            Graphics[AbsoluteThickness[2]]};
  For [e=1,e<=Length[elenod],e++, 
      enl=elenod[[e]]; ne=Length[enl];  
      If [!MemberQ[{3,4,6,8,9,10,16},ne], Print["ContourBandPlot:",
           "Elem ",e," with ",ne," nodes skipped"]; Continue[]];
      fc=Table[fn[[enl[[i]]]],{i,ne}]; xyc=Table[nodxyz[[enl[[i]]]],{i,ne}];
      bn=bnt[[ne]];  nbn=Length[bn]; xyb=Table[xyc[[bn[[i]]]],{i,nbn}];
      If [bound, enb=Table[enl[[bn[[i]]]],{i,nbn}]; 
          ebs=Table[{enb[[i]],enb[[i+1]]},{i,1,nbn-1}]; 
          For [i=1,i<=nbn-1,i++, seg=ebs[[i]]; rseg=Reverse[seg]; 
              If [MemberQ[bseg,rseg], {{k}}=Position[bseg,rseg,1,1]; 
              bseg=Drop[bseg,{k}]; Continue[]]; 
              AppendTo[bseg,seg]];
         ];
      If [ne==3, {pb,pl}=ContourBandPlotFunctionOnTrig3[
                  xyc,fc,{fmin,fmax,finc,eps},lines]];
      If [ne==4, {pb,pl}=ContourBandPlotFunctionOnQuad4[
                  xyc,fc,{fmin,fmax,finc,eps},lines]];
      If [ne==6, {pb,pl}=ContourBandPlotFunctionOnTrig6[
                  xyc,fc,{fmin,fmax,finc,eps},lines]];
      If [ne==8, {pb,pl}=ContourBandPlotFunctionOnQuad8[
                  xyc,fc,{fmin,fmax,finc,eps},lines]];
      If [ne==9, {pb,pl}=ContourBandPlotFunctionOnQuad9[
                  xyc,fc,{fmin,fmax,finc,eps},lines]];
      If [ne==10,{pb,pl}=ContourBandPlotFunctionOnTrig10[
                  xyc,fc,{fmin,fmax,finc,eps},lines]];
      If [ne==16,{pb,pl}=ContourBandPlotFunctionOnQuad16[
                  xyc,fc,{fmin,fmax,finc,eps},lines]]; 
      If [colors, pband[[e]]=pb]; If [lines, pline[[e]]=pl];
      If [mesh,   pmesh[[e]]=Graphics[Line[xyb]]]; 
      If [bacg,   pbacg[[e]]=Graphics[Polygon[xyb]]]; 
      If [nodes,  pnode[[e]]=Table[Graphics[Point[xyb[[i]] ]],{i,ne}]];   
      ];
 If [!colors,pband={}]; If [!mesh, pmesh={}]; If [!nodes, pnode={}]; 
 If [!lines, pline={}]; If [!bacg, pbacg={}]; pleg={};
 nbseg=Length[bseg]; 
 If [nbseg>0, pbound=Table[0,{nbseg}];  (*Print["bseg=",bseg];*)
     For [k=1,k<=nbseg,k++, {i,j}=bseg[[k]];
         pbound[[k]]=Graphics[Line[{nodxyz[[i]],nodxyz[[j]]}]] ]];
 If [Length[xyleg]==0||!colors, Show[prebacg,pbacg,
     preband,pband,preline,pline,premesh,pmesh,prenode,
     pnode,prebound,pbound,TextStyle->{FontFamily->"Courier",
     FontSlant->"Plain", FontSize->10},
     AspectRatio->aspect,PlotLabel->label,      
     DisplayFunction->DisplayChannel ]];
 If [Length[xyleg]>0&&colors,
     {preleg,pleg}=ContourPlotValueLegend[fmin,fmax,xyleg,6];
     Show[prebacg,pbacg,preband,pband,preline,pline,
          premesh,pmesh,prenode,pnode,prebound,pbound,preleg,pleg, 
          TextStyle->{FontFamily->"Courier",
                  FontSlant->"Plain", FontSize->9},
          AspectRatio->aspect,PlotLabel->label,
          DisplayFunction->DisplayChannel ]];
 ClearAll[pband,pline,pmesh,pnode,pbacg,pbound,pleg]; 
 ];

Cell 8: Print utilities

PrintPlaneStressNodeCoordinates[nodxyz_,title_,digits_]:= Module[
  {numnod=Length[nodxyz],n,x,y,z,d=6,f=6,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++,{x,y}=nodxyz[[n]]; 
       tab[[n]]={ToString[n],PaddedForm[x,{d,f}],PaddedForm[y,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->{Right},   
      TableDirections->{Column,Row},TableSpacing->{0,2},
      TableHeadings ->{None,{"node", "x-coor", "y-coor"}}]];
   ClearAll[tab];
   ];

PrintPlaneStressElementTypeNodes[eletyp_,elenod_,title_,digits_]:=
  Module[{e,numele=Length[elenod],d=6,f=3,tab}, 
   tab=Table[0,{numele}]; If [Length[digits]==2,{d,f}=digits];
   For [e=1,e<=numele,e++, 
        tab[[e]]={ToString[e],eletyp[[e]],ToString[elenod[[e]]]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->Right,
         TableDirections->{Column,Row},TableSpacing->{0,2},
         TableHeadings->{None,{"elem", "type", "node-list"}}]];
   ClearAll[tab];
   ];

PrintPlaneStressElementMatFab[elemat_,elefab_,title_,digits_]:= 
  Module[{e,numele=Length[elemat],mat,fab,d=6,f=3,tab}, 
   tab=Table[0,{numele}]; If [Length[digits]==2,{d,f}=digits];
   For [e=1,e<=numele,e++, mat=elemat[[e]]; fab=elefab[[e]];
        tab[[e]]={ToString[e],ToString[mat],ToString[fab]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab, TableAlignments->Right,
         TableDirections->{Column,Row},TableSpacing->{0,2},
         TableHeadings->{None,{"elem", "material", "fabrication"}}]];
   ClearAll[tab];
   ];

PrintPlaneStressFreedomActivity[nodtag_,nodval_,title_,digits_]:= Module[
  {numnod=Length[nodtag],n,t,v,d=6,f=2,tab}, 
   tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++, t=nodtag[[n]]; v=nodval[[n]];
        tab[[n]]={ToString[n],
          PaddedForm[t[[1]],d],PaddedForm[t[[2]],d],
          PaddedForm[v[[1]],{d,f}],PaddedForm[v[[2]],{d,f}] }];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings->{None,{"node", "x-tag", "y-tag",
                                "x-value", "y-value"}} ]];
   ClearAll[tab];
   ];

PrintPlaneStressNodeDisplacements[noddis_,title_,digits_]:= Module[
  {numnod=Length[noddis],n,ux,uy,d=8,f=6,tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits]; 
   For [n=1,n<=numnod,n++, {ux,uy}=noddis[[n]];
       tab[[n]]={ToString[n],PaddedForm[ux,{d,f}],PaddedForm[uy,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings->{None,{"node", "x-displ", "y-displ"}}]];
   ClearAll[tab];
   ];

PrintPlaneStressNodeForces[nodfor_,title_,digits_]:= Module[
  {numnod=Length[nodfor],n,fx,fy,tab,d=8,f=4}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
   For [n=1,n<=numnod,n++, {fx,fy}=nodfor[[n]]; 
       tab[[n]]={ToString[n],PaddedForm[fx,{d,f}],PaddedForm[fy,{d,f}]}];
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1},
         TableHeadings ->{None,{"node", "x-force", "y-force"}}]];         
   ClearAll[tab];
   ];

PrintPlaneStressPlateNodeStresses[nodpnc_,nodsig_,title_,digits_]:= 
  Module[{numnod=Length[nodpnc],n,k,sxx,syy,sxy,d=8,f=4,tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
  For [n=1,n<=numnod,n++, k=nodpnc[[n]]; {sxx,syy,sxy}=nodsig[[n]];
       If [k>0, tab[[n]]={ToString[n], ToString[k],
                 PaddedForm[sxx,{d,f}],PaddedForm[syy,{d,f}],
                 PaddedForm[sxy,{d,f}]} ]; 
       If [k==0, tab[[n]]={ToString[n], "no plates"," attached"}];
       ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{"node", "cnt", "sigma-xx", 
          "sigma-yy", "sigma-xy"}}]];
   ClearAll[tab];
   ];

PrintPlaneStressSolution[noddis_,nodfor_,nodpnc_,nodsig_,title_,digits_]:= 
  Module[{numnod=Length[nodpnc],n,k,sxx,syy,sxy,d=6,f=4,tab}, 
  tab=Table[0,{numnod}]; If [Length[digits]==2,{d,f}=digits];
  For [n=1,n<=numnod,n++,{ux,uy}=noddis[[n]]; {fx,fy}=nodfor[[n]]; 
       k=nodpnc[[n]]; {sxx,syy,sxy}=nodsig[[n]];
       If [k>0, tab[[n]]={ToString[n], 
                 PaddedForm[ux,{d,f}],PaddedForm[uy,{d,f}],
                 PaddedForm[fx,{d,f}],PaddedForm[fy,{d,f}],
                 PaddedForm[sxx,{d,f}],PaddedForm[syy,{d,f}],
                 PaddedForm[sxy,{d,f}]} ]; 
       If [k==0, tab[[n]]={ToString[n], 
                 PaddedForm[ux,{d,f}],PaddedForm[uy,{d,f}],
                 PaddedForm[fx,{d,f}],PaddedForm[fy,{d,f}]} ];
       ]; 
   If [StringLength[title]>0, Print[title]];
   Print[TableForm[tab,TableAlignments->{Right},
         TableDirections->{Column,Row},TableSpacing->{0,1}, 
         TableHeadings ->{None,{"node", "x-displ", "y-displ","x-force",
         "y-force",  "sigma-xx", "sigma-yy", "sigma-xy"}}]];
   ClearAll[tab];
   ];

Cell 9 :  Miscellaneous utilities

FlatNodePartVector[nv_]:=Flatten[nv];
           
NodePartFlatVector[nfc_,v_]:= Module [
  {i,k,m,n,nv={},numnod},
  If [Length[nfc]==0, nv=Partition[v,nfc]];
  If [Length[nfc]>0,  numnod=Length[nfc]; m=0;
      nv=Table[0,{numnod}];   
      For [n=1,n<=numnod,n++, k=nfc[[n]]; 
          nv[[n]]=Table[v[[m+i]],{i,1,k}]; 
          m+=k]];
  Return[nv]];

Cell 10:       Driver for static analysis of plane stress problem.

PlaneStressSolution[nodxyz_,eletyp_,elenod_,elemat_,elefab_,
   nodtag_,nodval_,prcopt_]:= Module[{K,Kmod,f,fmod,u,numer=True,
   noddis,nodfor,nodpnc,nodsig,barele,barfor},
   If [Length[prcopt]>=1, numer=prcopt[[1]]];
   K=PlaneStressMasterStiffness[nodxyz,eletyp,elenod,elemat,elefab,prcopt];
   If [K==Null,Return[Table[Null,{6}]]];
  (* Print["K=",K//MatrixForm];
   Print["eigs of K=",Chop[Eigenvalues[N[K]]]];*)
   Kmod=ModifiedMasterStiffness[nodtag,K];
   f=FlatNodePartVector[nodval];
   fmod=ModifiedNodeForces[nodtag,nodval,K,f];
   (* Print["eigs of Kmod=",Chop[Eigenvalues[N[Kmod]]]];*)
   u=LinearSolve[Kmod,fmod]; If [numer,u=Chop[u]]; 
   f=K.u;                    If [numer,f=Chop[f]];  
   nodfor=NodePartFlatVector[2,f]; noddis=NodePartFlatVector[2,u];
   {nodpnc,nodsig}=PlaneStressPlateStresses[nodxyz,eletyp,elenod,elemat,
                   elefab,noddis,prcopt];
   {barele,barfor}=PlaneStressBarForces[nodxyz,eletyp,elenod,elemat,
                   elefab,noddis,prcopt];
   ClearAll[K,Kmod];
   Return[{noddis,nodfor,nodpnc,nodsig,barele,barfor}]; 
];