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

where is the so-called geometrical matrix, containing derivatives of shape functions and is the material stiffness matrix. If 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

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

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
, where
is the geometrical matrix and
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''.

Borek Patzak 2018-01-02