Input file for Rollup of Cantilever Beam ...

/* 
 *  =====================================================================
 *  Analysis of cantilever with 2D geometrically exact rod finite element
 * 
 *  Theoretical considerations indicate that for a cantilever rod having
 *  length L, Young's Modulus, E, and inertia I, the end moment:
 *
 *  M = 2*pi*(EI/L)
 *
 *  will roll the cantilever into a circle.
 * 
 *  Written By: Mark Austin                                   August 2000
 *  =====================================================================
 */ 

/* Setup problem specific parameters */

   NDimension         = 2;
   NDofPerNode        = 3;
   MaxNodesPerElement = 2;

   StartMesh();

/* Generate line of nodes */

   for(node = 1; node <= 11; node = node + 1) {
       AddNode(node, [ (node-1)*5 m,  0 m ] );
   }

/* Attach beam elements to nodes */

   for(elmtno = 1; elmtno <= 10; elmtno = elmtno + 1) {
       AddElmt( elmtno, [ elmtno, elmtno + 1 ], "rodproperties");
   }

/* Define section and material properties */

   ElementAttr("rodproperties") { type     = "GEXACT_2D";
                                  section  = "rodsection";
                                  material = "rodmaterial";
                                }

   SectionAttr("rodsection") { Izz       =    100 m^4;
                               area      = 0.0001 m^2;
                             }

   MaterialAttr("rodmaterial") { poisson = 0.25;
                                 E       = 10000 Pa;
                               }

/* Apply full-fixity at cantilever base */

   nodeno = 1;
   FixNode( nodeno, [ 1, 1, 1 ]);

/* Add external load to node 4 */

   NodeLoad( 11, [ 0.0 N, 0.0 N, 2*PI*20000 N*m ]);

/* Compile and Print Finite Element Mesh */

   EndMesh();
   PrintMesh();

/* Allocate memory for displacment response */

   noIterations =  6;
   noNode       = 11;
   xResponse = Matrix( [ noNode, noIterations ]);
   yResponse = Matrix( [ noNode, noIterations ]);
   for(i = 1; i <= noIterations; i = i + 1) {
       xResponse = ColumnUnits( xResponse, [m], [i] );
       yResponse = ColumnUnits( yResponse, [m], [i] );
   }

   for (ii = 1; ii <= noNode; ii = ii + 1 ) {
        coord = GetCoord([ii]);
        xResponse[ ii ][ 1 ] = coord [1][1];
        yResponse[ ii ][ 1 ] = coord [1][2];
   }

/* Compute and print initial stiffness, external load, and displacement */

   eload = ExternalLoad();
   y     = L2Norm(eload);

   i  = 0;
   print "*** ====================================== \n";
   print "*** Initial Configuration                  \n";
   print "*** L2norm(unbalance force) = ", y,       "\n";
   print "*** Relative error : err = 1.0             \n";
   print "*** ====================================== \n";

   for (ii = 1; ii <= noNode; ii = ii + 1 ) {
           coord     = GetCoord([ii]);
           print "Node : ", ii, " x = " ,coord [1][1], " y = ", coord[1][2], "\n";
   }

/* Compute and print displacements */

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

/* Main loop for nonlinear iteration */

   x = L2Norm(eload);
   tol = 0.0000001;
   err = 1 + tol;
   while( err > tol ) { 

      i = i + 1;

      /* Save displacement to database      */

      ElmtStateDet( displ );
      UpdateResponse();

      /* Compute tangent stiffness          */

      stiff = Stiff();

      /* Compute internal load matrix       */

      Fs = InternalLoad( displ );

      /* Compute vector of unbalance forces */

      Unbalance = eload - Fs;

      /* Compute increment in displcement   */

      d_displ = Solve(stiff, Unbalance) ;

      /* Update displacement vector */

      displ = displ + d_displ;

      /* Compute relative error in unbalance force */

      y = L2Norm(Unbalance);
      err = y/x;

      /* Compute and print coordinates at ith iteration */

      print "\n";
      print "*** =================================== \n";
      print "*** Iteration No = ", i,               "\n";
      print "*** L2norm(unbalance force) = ", y,    "\n";
      print "*** Relative error : err = ", err,     "\n";
      print "*** =================================== \n";

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

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

           print "Node : ", ii, " x = " ,coord [1][1], " y = ", coord[1][2], "\n";

           if ( i <= 5 ) {
              xResponse[ ii ][ i+1 ] = coord [1][1];
              yResponse[ ii ][ i+1 ] = coord [1][2];
           }
      }
   }

   PrintDispl(displ);

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

   no_node = 11;
   for (ii = 1; ii <= no_node; ii = ii + 1 ) {
        coord     = GetCoord([ii]);
        nodal_dof = GetDof([ii]);

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

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

   /* Summary of time-history response */

   PrintMatrix(xResponse);
   PrintMatrix(yResponse);

   /* Print element level forces */

   PrintStress(displ);

   quit;