added mass matrix to link elements

Author: githubgrasp

src/sm/Elements/LatticeElements/latticelink3d.C

@@ -143,6 +143,20 @@ namespace oofem {
         return;
     }
 
+    void
+    LatticeLink3d::computeLumpedMassMatrix(FloatMatrix &answer, TimeStep *tStep)
+    // Returns the lumped mass matrix of the receiver. This expression is
+    // valid in both local and global axes.
+    {
+        GaussPoint *gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
+        double density = static_cast < LatticeCrossSection * > ( this->giveCrossSection() )->give('d', gp);
+        double halfMass = density * computeVolumeAround(gp) / 2.;
+        answer.resize(12, 12);
+        answer.zero();
+        answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = halfMass;
+        answer.at(7, 7) = answer.at(8, 8) = answer.at(9, 9) = halfMass;
+
+    }
 
     void
     LatticeLink3d::computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode,
@@ -328,9 +342,14 @@ namespace oofem {
         FloatArray normal(3), s(3), t(3);
 
         //Calculate normal vector
+        if(this->directionVector.computeNorm() < 1e-12) {
+            OOFEM_ERROR("LatticeLink3d: directionVector is zero. Check input or link setup.");
+        }
         normal = this->directionVector;
         normal.normalize();
 
+
+
         //Construct two perpendicular axis so that n is normal to the plane which they create
         //Check, if one of the components of the normal-direction is zero
         if ( normal.at(1) == 0 ) {

src/sm/Elements/LatticeElements/latticelink3d.h

@@ -89,6 +89,9 @@ class LatticeLink3d : public LatticeStructuralElement
     virtual void giveGPCoordinates(FloatArray &coords);
 
     virtual void computeGeometryProperties();
+    void computeLumpedMassMatrix(FloatMatrix &answer, TimeStep *tStep) override;
+    void computeMassMatrix(FloatMatrix &answer, TimeStep *tStep) override
+    { this->computeLumpedMassMatrix(answer, tStep); }
 
     void giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord = 0) override;
 

src/sm/Materials/LatticeMaterials/latticeslip.C

@@ -117,7 +117,7 @@ namespace oofem {
         //evaluate tempKappa (no elastic strain in axial direction)
         double tempKappa = status->giveKappa() + fabs(totalStrain.at(1) - oldStrain.at(1) );
 
-
+        double deltaSlip = totalStrain.at(1) - oldStrain.at(1);
 
         /*First component is the slip one for which the stress should be limited using plasiticity (frictional slip between fibre and matrix). The other components are kept elastic. */
         FloatArrayF < 6 > stress;
@@ -129,7 +129,7 @@ namespace oofem {
 
         if ( f > 0 ) {//plastic response.
             //Reduced stress by increasing plastic strain.
-            stress.at(1) = sgn(stress.at(1) ) * bondStress;
+            stress.at(1) = sgn(stress.at(1)) * bondStress;
             status->setTempCrackFlag(1);
         }