Input file for Plane Stress Analysis of Cantilever

/*
 *  =============================================================================
 *  Plane stress analysis of rectangular cantilever beam subject to an end moment
 *
 *  Written By : Clara Popescu and Mark Austin                          1997-1998
 *  =============================================================================
 */ 

/* Define problem specific parameters */

NDimension         = 2;
NDofPerNode        = 2;
MaxNodesPerElement = 4;

StartMesh();

/* Generate a 3 by 11 grid of nodes */

   beam_length  =  200 cm;
   beam_width   =    5 cm;
   beam_height  =   30 cm;

   node = 0;
   x = 0 cm ; y = 0 cm; 

   while(x <= beam_length) {
      node = node + 1;
      AddNode(node, [x, y]);

      y = y + beam_height/2;
      node = node + 1;
      AddNode(node, [x, y]);

      y = y + beam_height/2;
      node = node + 1;
      AddNode(node, [x, y]);
	
      x = x + beam_length/10;
      y = y - beam_height;
   }

/* Attach elements to grid of nodes in anticlockwise fashion */

   i = 1; elmtno = 1;
   while(elmtno <= 19) {

      AddElmt( elmtno, [   i, i+3, i+4, i+1 ], "cantilever_elmt_attr");
      elmtno = elmtno + 1;
	
      AddElmt( elmtno, [ i+1, i+4, i+5, i+2 ], "cantilever_elmt_attr");
      elmtno = elmtno + 1;

      i = i + 3;
   }

/* Define element section and material properties */

   ElementAttr("cantilever_elmt_attr") { type     = "PLANE_STRESS";
                                         section  = "cantilever_section";
                                         material = "cantilever_material"; }

   SectionAttr("cantilever_section") { depth = beam_height;
                                       width = beam_width; }

   MaterialAttr("cantilever_material") { poisson = 1/3;   
                                         E       = 20000000 kN/m^2; }

/* Setup boundary conditions */ 

   FixNode( 1, [1,1] );
   FixNode( 2, [1,0] );
   FixNode( 3, [1,0] );

/* Add point nodal loads to end of cantilever */

   Fx = 10.0 kN; Fy = 0.0 kN;

   NodeLoad( 31, [  Fx,  Fy]);
   NodeLoad( 33, [ -Fx,  Fy]);

/* Compile and print finite element mesh */

   EndMesh();
   PrintMesh();

/* Compute stiffness and external load matrices; static analysis */

   stiff  = Stiff();
   eload  = ExternalLoad();
   displ  = Solve(stiff, eload);

/* Print the beam displacements */

   PrintDispl(displ);

/* Compute nodal coordinates of structure with scaled displacements */

   print "\n";
   print "Coordinates of Displaced Cantilever \n";
   print "=================================== \n\n";

   scale = 1000.0;
   for (ii = 1; ii <= 33; ii = ii + 1 ) {
        coord     = GetCoord([ii]);
        nodal_dof = GetDof([ii]);

        if( nodal_dof[1][1] > 0 ) {
            coord[1][1] = coord[1][1] + scale*displ[ nodal_dof[1][1] ][1];
        }
        if( nodal_dof[1][2] > 0 ) {
            coord[1][2] = coord[1][2] + scale*displ[ nodal_dof[1][2] ][1];
        }

        print ii, coord[1][1], coord[1][2], "\n";
   }

/* Setup matrix for exptrapolation of stresses */

   M = Zero([4,4]);

   M[1][1] = (sqrt(3) + 1)^2;
   M[1][2] = 2;
   M[1][3] = (sqrt(3) - 1)^2;
   M[1][4] = 2;
   M[2][1] = 2;
   M[2][2] = (sqrt(3) + 1)^2;
   M[2][3] = 2;
   M[2][4] = (sqrt(3) - 1)^2;
   M[3][1] = (sqrt(3) - 1)^2;
   M[3][2] = 2;
   M[3][3] = (sqrt(3) + 1)^2;
   M[3][4] = 2;
   M[4][1] = 2;
   M[4][2] = (sqrt(3) - 1)^2;
   M[4][3] = 2;
   M[4][4] = (sqrt(3) + 1)^2;

   PrintMatrix(M);
   T = Inverse((1/12)*M);
  
/* Systematically retrieve stresses from individual elements */

   print "\n";
   print "Element Stresses (sigma_xx at the nodal points)\n";
   print "===============================================\n\n";

   for( ii = 1; ii <= 20 ; ii = ii + 1 ) {

        /* retrieve stresses and extrapolate out to nodal coordinates    */

        actions = GetStress( [ii], displ );
        extrap  = T*[ actions[1][3];
                      actions[2][3];
                      actions[3][3];
                      actions[4][3] ];

        /* map extrapolated stresses onto nodal coordinate system    */

        extrap = [ 0, 0, 1, 0;
                   0, 0, 0, 1;
                   1, 0, 0, 0;
                   0, 1, 0, 0 ] * extrap;
       
        /* print extrapolated stresses in format for MATLAB plotting */

        print ii;
        print QDimenLess(extrap[1][1]);
        print QDimenLess(extrap[2][1]);
        print QDimenLess(extrap[3][1]);
        print QDimenLess(extrap[4][1]);
        print "\n";
   }
  
   quit;