/*
 *  =====================================================================
 *  Finite element model of Shear-beam using 4 node shell finite element.
 *
 *  Written By : Lanheng Jin                                  April, 1994
 *  =====================================================================
 */

print "*** DEFINE PROBLEM SPECIFIC PARAMETERS \n\n";

NDimension         = 3;
NDofPerNode        = 6;
MaxNodesPerElement = 4;

/*
 *  ===============
 *  Define the Mesh
 *  ===============
 */

StartMesh();

print "*** GENERATE GRID OF GEOMETRY FOR FE MODEL \n\n";

L   =  48 in;   /* Length     */
b   =  12 in;   /* Width      */
t   =   1 in;   /* Thickness  */

print "*** GENERATE GRID OF NODES FOR FE MODEL \n\n";

y  = 0 in;

AddNode(1,[0.0 in,y,6.0 in]);  
AddNode(2,[0.0 in,y,0.0 in]);  
AddNode(3,[0.0 in,y,-6.0 in]);  
AddNode(4,[6.0 in,y,6.0 in]);  
AddNode(5,[7.0 in,y,0.0 in]);  
AddNode(6,[8.0 in,y,-6.0 in]);  
AddNode(7,[12.0 in,y,6.0 in]);  
AddNode(8,[14.0 in,y,0.0 in]);  
AddNode(9,[16.0 in,y,-6.0 in]);  
AddNode(10,[18.0 in,y,6.0 in]);  
AddNode(11,[18.0 in,y,0.0 in]);  
AddNode(12,[18.0 in,y,-6.0 in]);  
AddNode(13,[24.0 in,y,6.0 in]);  
AddNode(14,[22.0 in,y,0.0 in]);  
AddNode(15,[20.0 in,y,-6.0 in]);  
AddNode(16,[30.0 in,y,6.0 in]);  
AddNode(17,[27.0 in,y,0.0 in]);  
AddNode(18,[24.0 in,y,-6.0 in]);  
AddNode(19,[36.0 in,y,6.0 in]);  
AddNode(20,[32.0 in,y,0.0 in]);  
AddNode(21,[28.0 in,y,-6.0 in]);  
AddNode(22,[42.0 in,y,6.0 in]);  
AddNode(23,[40.0 in,y,0.0 in]);  
AddNode(24,[38.0 in,y,-6.0 in]);  
AddNode(25,[48.0 in,y,6.0 in]);  
AddNode(26,[48.0 in,y,0.0 in]);  
AddNode(27,[48.0 in,y,-6.0 in]);  

print "*** ATTACH ELEMENTS TO GRID OF NODES \n\n";

nn  = 9;

for(j=1; j<nn; j=j+1) {
     for(i=1; i<=2; i=i+1) {
          elmtno = i+2*(j-1);
          a = 3*(j-1)+i+1;
          b = 3*j+i+1;
          c = 3*j+i;
          d = 3*(j-1)+i;
          node_connec = [a, b, c, d];
          AddElmt(elmtno, node_connec, "name_of_elmt_attr");
      }
}

/*
 *  ===============================================
 *  Define Element, Section and Material Properties
 *  ===============================================
 */

print "*** DEFINE ELEMENT, SECTION AND MATERIAL PROPERTIES \n\n";

ElementAttr("name_of_elmt_attr") { type     = "SHELL_4NQ";
                                   section  = "mysection";
                                   material = "ELASTIC";
                                  }

MaterialAttr("ELASTIC") { poisson = 0.25;   
	                  E       = 30000 ksi;
                         }

SectionAttr("mysection") { thickness =   1 in; }

/*
 *  =========================
 *  Setup Boundary Conditions
 *  =========================
 */

print "*** SET UP BOUNDARY CONDITIONS \n\n";

/* full fixity condition  */

u_id  = 1; v_id  = 1; w_id  = 1;
rx_id = 1; ry_id = 1; rz_id = 1;
bc_fc = [u_id,v_id,w_id,rx_id,ry_id,rz_id];

/* apply full fixity to corner nodes  */

for (i=1; i<=3; i=i+1) {
     FixNode(i, bc_fc);
}

/*
 *  ==================
 *  Add external loads
 *  ==================
 */

print "*** APPLY EXTERNAL LOADS \n\n";

/* [0] load for end nodes  */

Fx = 0 lbf;   Fy = 0 lbf;   Fz = 10000 lbf;
Mx = 0 lbf*in; My = 0 lbf*in; Mz = 0 lbf*in;
NodeLoad(25,[Fx, Fy, Fz, Mx, My, Mz]);
NodeLoad(27,[Fx, Fy, Fz, Mx, My, Mz]);

Fx = 0 lbf;   Fy = 0 lbf;   Fz = 20000 lbf;
Mx = 0 lbf*in; My = 0 lbf*in; Mz = 0 lbf*in;
NodeLoad(26,[Fx, Fy, Fz, Mx, My, Mz]);

/*
 *  =====================================
 *  Compile and Print Finite Element Mesh
 *  =====================================
 */

EndMesh();
PrintMesh();

/*
 *  ==========================
 *  Compute Stiffness Matrices
 *  ==========================
 */

print "\n*** COMPUTE STIFFNESS MATRICES \n\n";

SetUnitsType("US");
stiff  = Stiff();
eload  = ExternalLoad();
lu     = Decompose(Copy(stiff));
displ  = Substitution(lu, eload);

/*
 *  ================================
 *  Print displacements and stresses
 *  ================================
 */

PrintDispl(displ);
PrintStress(displ);

quit;