Isoparametric Truss 2d element

class Truss2d : public StructuralElement
{
/*
 This class implements a two-node truss bar element 
 for two-dimensional analysis.
 A truss bar element is characterized by its 'length' 
 and its 'pitch'. The pitch is the angle in radians 
  between the X-axis anf the axis of the element 
  (oriented node1 to node2).
 */

 protected :
   double        length ;
   double        pitch ;

 public :
   Truss2d (int,Domain*) ;    // constructor
   ~Truss2d ()   {}           // destructor

   // mass matrix coputations
   void computeLumpedMassMatrix (FloatMatrix& answer, 
                                 TimeStep* tStep) ;
   // general mass service overloaded 
   void computeMassMatrix (FloatMatrix& answer, TimeStep* tStep) 
        {computeLumpedMassMatrix(answer, tStep);}

       // DOF management
   virtual int            computeNumberOfDofs (EquationID ut) {return 4;}
   virtual void           giveDofManDofIDMask  (int inode, EquationID, IntArray& ) const;

       double computeVolumeAround (GaussPoint*) ;

 // 
 // definition & identification
 //
   char* giveClassName (char* s) const 
         { return strcpy(s,"Truss2d") ;}
   classType giveClassID () const { return Truss2dClass; } 

   IRResultType initializeFrom (InputRecord* ir);

 protected:
   // computes geometrical matrix 
   void computeBmatrixAt (GaussPoint*, FloatMatrix&, 
                          int=1, int=ALL_STRAINS) ;
   // computes interpolation matrix
   void computeNmatrixAt (GaussPoint*, FloatMatrix&) ;
   // initialize element's integration rules
   void computeGaussPoints () ;
   // transformation from global->local c.s.
   int           computeGtoLRotationMatrix (FloatMatrix&);

   double giveLength () ;
   double givePitch () ;
 } ;

This is a minimal-functionality implementation. It does not include, for example, any support for element loading or material nonlinearity. The implementation starts with standard constructor and destructor.

 Truss2d :: Truss2d (int n, Domain* aDomain) : 
            StructuralElement (n,aDomain)
 // Constructor.
 {
    numberOfDofMans       = 2 ;
    rotationMatrix      = NULL ;
    length              = 0. ;
    pitch               = 10. ;      // a dummy value

 }

Next, the services for computing interpolation and geometrical matrices follow. The implementation of the computeNmatrixAt is not necessary for the current purpose, since the mass matrix computation is overloaded and does not require this service, but it is added for completeness. The both methods compute response at a given integration point, which is passed as a parameter. The computeBmatrixAt has two additional parameters which determine the range of strain components for which response is assembled. This has something to do with support for reduced/selective integration and is not important in presented simple case.

Recently, the interpolation classes have been added, that can significantly facilitate the element implementation. They provide shape functions, their derivatives, transformation jacobians out of the box.

 void
 Truss2d :: computeNmatrixAt (GaussPoint* aGaussPoint, 
                              FloatMatrix& answer) 
 // Returns the displacement interpolation matrix {N} 
 // of the receiver, evaluated at aGaussPoint.
 {
    double       ksi,n1,n2 ;

    ksi = aGaussPoint -> giveCoordinate(1) ;
    n1  = (1. - ksi) * 0.5 ;
    n2  = (1. + ksi) * 0.5 ;

    answer.resize (2,4);
    answer.zero();

    answer.at(1,1) = n1 ;
    answer.at(1,3) = n2 ;
    answer.at(2,2) = n1 ;
    answer.at(2,4) = n2 ;

    return  ;
 }


 void
 Truss2d :: computeBmatrixAt (GaussPoint* aGaussPoint, 
                              FloatMatrix& answer, int li, int ui)
 // 
 // Returns linear part of geometrical 
 // equations of the receiver at gp.
 // Returns the linear part of the B matrix
 //
 {
   double coeff,l;

   answer.resize(1,4);
   l = this->giveLength();
   coeff = 1.0/l;

   answer.at(1,1) =-coeff;
   answer.at(1,2) = 0.0;
   answer.at(1,3) = coeff;
   answer.at(1,4) = 0.0;

   return;
 }

The following two functions compute the basic geometric characteristics of a bar element - its length and pitch, defined as the angle between global x-axis and the local element x-axis (oriented from node1 to node2).

 double  Truss2d :: giveLength ()
    // Returns the length of the receiver.
 {
    double dx,dz ;
    Node   *nodeA,*nodeB ;

    if (length == 0.) {
       nodeA = this->giveNode(1) ;
       nodeB = this->giveNode(2) ;
       dx    = nodeB->giveCoordinate(1)-nodeA->giveCoordinate(1);
       dz    = nodeB->giveCoordinate(3)-nodeA->giveCoordinate(3);
       length= sqrt(dx*dx + dz*dz) ;}

    return length ;
 }


 double  Truss2d :: givePitch ()
    // Returns the pitch of the receiver.
 {
    double xA,xB,zA,zB ;
    Node   *nodeA,*nodeB ;

    if (pitch == 10.) {  // 10. : dummy initialization value
       nodeA  = this -> giveNode(1) ;
       nodeB  = this -> giveNode(2) ;
       xA     = nodeA->giveCoordinate(1) ;
       xB     = nodeB->giveCoordinate(1) ;
       zA     = nodeA->giveCoordinate(3) ;
       zB     = nodeB->giveCoordinate(3) ;
       pitch  = atan2(zB-zA,xB-xA) ;}

    return pitch ;
 }

When an element is created (by the Domain class), the default constructor is called. To initialize the element, according to its record in the input database, the initializeFrom is immediately called after element creation. The element implementation should first call the parent implementation to ensure that attributes declared at parent level are initialized properly. Then the element has to initialize attributes declared by itself and also to set up its integration rules. In our example, special method computeGaussPoints is called to initialize integration rules. In this case, only one integration rule is created. It is of type GaussIntegrationRule, indicating that the Gaussian integration is used. Once integration rule is created, its integration points are created to represent line integral, with 1 integration point. Integration points will be associated to element under consideration and will have 1D material mode (which determines the type of material model response).

 IRResultType
 Truss2d :: initializeFrom (InputRecord* ir)
 {
    this->NLStructuralElement :: initializeFrom (ir);
    this -> computeGaussPoints();
    return IRRT_OK;
 }

 void  Truss2d :: computeGaussPoints ()
    // Sets up the array of Gauss Points of the receiver.
 {

  numberOfIntegrationRules = 1 ;
  integrationRulesArray = new IntegrationRule*;
  integrationRulesArray[0] = new GaussIntegrationRule (1,domain, 1, 2);
  integrationRulesArray[0]->
    setUpIntegrationPoints (_Line, 1, this, _1dMat);

 }

 double  Truss2d :: computeVolumeAround (GaussPoint* aGaussPoint)
 // Returns the volume corresponding to given integration point. 
 {
   double weight  = aGaussPoint -> giveWeight() ;
   return 0.5 * this->giveLength() * weight * 
          this->giveCrossSection()->give('A');
 }

 void
 Truss2d :: computeLumpedMassMatrix (FloatMatrix& answer, TimeStep* tStep)
    // Returns the lumped mass matrix of the receiver. This expression is
    // valid in both local and global axes.
 {
    Material* mat ;
    double    halfMass ;

    mat        = this -> giveMaterial() ;
    halfMass   = mat->give('d') * 
                 this->giveCrossSection()->give('A') * 
                 this->giveLength() / 2.;

    answer.resize (4,4) ; answer.zero();
    answer . at(1,1) = halfMass ;
    answer . at(2,2) = halfMass ;
    answer . at(3,3) = halfMass ;
    answer . at(4,4) = halfMass ;

    if (this->updateRotationMatrix()) 
       answer.rotatedWith(*this->rotationMatrix) ;
    return  ;
 }

The following service computes the part of the element transformation matrix, corresponding to transformation between global and element local coordinate systems. This method is called from the general updateRotationMatrix service, implemented at the StructuralElement level, which computes the element transformation matrix, taking into account further transformations (nodal coordinate system, for example). The default implementation of updateRotationMatrix computes the element transformation matrix only once and stores it in rotationMatrix attribute. Since this service is virtual, particular elements can overload this service to take into account continuously changing local coordinate system (corotational elements in nonlinear analysis).

 int
 Truss2d ::  computeGtoLRotationMatrix (FloatMatrix& rotationMatrix) 
 // computes the rotation matrix of the receiver.
 // r(local) = T * r(global)
 {
   double sine,cosine ;

   sine           = sin (this->givePitch()) ;
   cosine         = cos (pitch) ;

   rotationMatrix.resize(4,4);
   rotationMatrix . at(1,1) =  cosine ;
   rotationMatrix . at(1,2) =  sine   ;
   rotationMatrix . at(2,1) = -sine   ;
   rotationMatrix . at(2,2) =  cosine ;
   rotationMatrix . at(3,3) =  cosine ;
   rotationMatrix . at(3,4) =  sine   ;
   rotationMatrix . at(4,3) = -sine   ;
   rotationMatrix . at(4,4) =  cosine ;

   return 1 ;
 }

 void
 Truss2d ::   giveDofManDofIDMask  (int inode, EquationID, IntArray& answer) const {
 // returns DofId mask array for inode element node.
 // DofId mask array determines the dof ordering requsted from node.
 // DofId mask array contains the DofID constants (defined in cltypes.h)
 // describing physical meaning of particular DOFs.
 //IntArray* answer = new IntArray (2);
  answer.resize (2);

  answer.at(1) = D_u;
  answer.at(2) = D_w;

  return ;
}