Structural Element - Example

This class is the base class for all structural elements. It is derived from the Element class. The basic tasks of the structural element is to compute its contributions to global equilibrium equations (mass and stiffness matrices, various load vectors (due to boundary conditions, force loading, thermal loading, etc.) and computing the corresponding strains and stresses from nodal displacements. Therefore the corresponding virtual services for computing these contributions are declared. These standard contributions can be computed by numerical integration of appropriate terms, which typically depend on element interpolation or material model, over the element volume. Therefore, it is possible to provide general implementations of these services, provided that the corresponding methods for computing interpolation dependent terms and material terms are implemented, and corresponding integration rules are initialized. This concept will be demonstrated on service computing stiffness matrix. Since element stiffness matrix contributes to the global equilibrium, the stiffness will be requested using giveCharacteristicMatrix service. The implementation of the services of the structural element takes into account only material nonlinearity, the geometrical nonlinearity is not taken into account (regarding geometrical nonlinearity, see NLStructuralElement class). The implementation of this service is following:

void
StructuralElement ::  giveCharacteristicMatrix (FloatMatrix& answer, 
                      CharType mtrx, TimeStep *tStep) 
// 
// returns characteristic matrix of receiver according to mtrx
//
{
  if (mtrx == TangentStiffnessMatrix) 
    this -> computeStiffnessMatrix(answer, TangentStiffness, tStep);
  else if (mtrx == SecantStiffnessMatrix) 
    this -> computeStiffnessMatrix(answer, SecantStiffness, tStep);
  else if (mtrx == MassMatrix) 
    this -> computeMassMatrix(answer, tStep);
  else if 
  ....
}

The first parameter is the matrix object to be computed, the mtrx parameter determines the type of contribution and the last parameter, time step, represents time. The element stiffness matrix can be evaluated using

$$\mbox{\boldmath $K$} = \int_{V} \mbox{\boldmath $B$}^T\mbox{\boldmath $D$}\mbox{\boldmath $B$}\;dV,$$

where $\mbox{\boldmath$B$}$ is the so-called geometrical matrix, containing derivatives of shape functions and $\mbox{\boldmath$D$}$ is the material stiffness matrix. If $\mbox{\boldmath$D$}$ is symmetric (which is usually the case) then element stiffness is symmetric, too. The numerical integration is used to evaluate this integral. For numerical integration, the IntegrationRule class is used. The integration rules for a specific element are created during element initialization and are stored in integrationRulesArray attribute, inherited from the Element class. In order to implement the stiffness evaluation, the methods for computing geometrical matrix and material stiffness matrix are declared (as virtual), but not implemented. They have to be implemented by specific elements, because they ``know'' their interpolation and the corresponding material mode. The implementation of computeStiffnessMatrix is following

void
StructuralElement :: computeStiffnessMatrix (FloatMatrix& answer, 
                     MatResponseMode rMode, TimeStep* tStep)
// Computes numerically the stiffness matrix of the receiver.
{
  int j;
  double      dV ;
  FloatMatrix d, bj, dbj;
  GaussPoint  *gp ;
  IntegrationRule* iRule;

  // give reference to integration rule
  iRule = integrationRulesArray[giveDefaultIntegrationRule()];
     
  // loop over integration points
  for (j=0 ; j < iRule->getNumberOfIntegrationPoints() ; j++) {
     gp = iRule->getIntegrationPoint(j) ;
     // compute geometrical matrix of particular element 
     this -> computeBmatrixAt(gp, bj) ;  
     //compute material stiffness
     this -> computeConstitutiveMatrixAt(d, rMode, gp, tStep); 
     // compute jacobian
     dV = this -> computeVolumeAround(gp) ; 
     // evaluate stiffness
     dbj.beProductOf (d, bj) ;
     answer.plusProductSymmUpper(bj,dbj,dV) ; 

  }
 
  answer.symmetrized() ;
  // transform into global coordinate system if necessary
  if (this->updateRotationMatrix()) 
     answer.rotatedWith(*this->rotationMatrix) ;
  return  ;
}

Inside the integration loop, only the upper half of the element stiffness is computed in element local coordinate system. Then, the lower part of the stiffness is initialized from the upper part (answer.symmetrized()) and transformation from element local coordinate system to global coordinate system is done. The transformation matrix is cached in the rotationMatrix attribute declared in the StructuraleElement class. The updateRotationMatrix service updates this matrix (if necessary) and returns pointer to this matrix. The zero value indicates that no transformation is necessary. The other element contributions can be computed using similar procedures. In general, different integration rules can be used for evaluation of different element contributions. For example, the support for the reduced integration of some terms of the stiffness matrix can be implemented - see the implementation of computeStiffnessMatrix in ``structuralelement.C''.

The element strain and stress vectors at a given integration point are computed using computeStrainVector and computeStressVector services. The element strain vector can be evaluated using $\mbox{\boldmath$\varepsilon$}=\mbox{\boldmath$B$}\mbox{\boldmath$u$}$, where $\mbox{\boldmath$B$}$ is the geometrical matrix and $\mbox{\boldmath$u$}$ is element local displacement vector. The stress evaluation is rather simple, since the stress evaluation from a given strain increment and actual state (kept within the integration point) is done at the cross section description level (integration over cross section volume) and material model level (stress evaluation).

void
StructuralElement :: computeStrainVector (FloatArray& answer, 
                     GaussPoint* gp, TimeStep* stepN)
// Computes the vector containing the strains 
// at the Gauss point gp of the receiver, 
// at time step stepN. The nature of these strains depends
// on the element's type.
{
  FloatMatrix b;
  FloatArray  u ;
  
  this -> computeBmatrixAt(gp, b) ;
  // compute vector of element's unknowns
  this -> computeVectorOf(DisplacementVector,
                          UnknownMode_Total,stepN,u) ;
  // transform global unknowns into element local c.s.
  if (this->updateRotationMatrix()) 
     u.rotatedWith(this->rotationMatrix,'n') ;
  answer.beProductOf (b, u) ;
  
  return ;
}


void
StructuralElement :: computeStressVector (FloatArray& answer, 
                     GaussPoint* gp, TimeStep* stepN)
// Computes the vector containing the stresses 
// at the Gauss point gp of the receiver, at time step stepN. 
// The nature of these stresses depends on the element's type.
{

  FloatArray Epsilon ;
  StructuralCrossSection* cs = (StructuralCrossSection*) 
                                this->giveCrossSection();
  Material *mat = this->giveMaterial();
  
  this->computeStrainVector (Epsilon, gp,stepN) ;
  // ask cross section model for real stresses 
  // for given strain increment 
  cs -> giveRealStresses (answer, ReducedForm, gp, Epsilon, stepN);

  return  ;
}

The StructuralElement class overloads also services for context storing/restoring, calling corresponding services for all integration rules (and thus on all integration points belonging to an element).

For further reference see OOFEM library reference manualand files ``structuralelement.h'' and ``structuralelement.C''.