Added: Account for element damping matrix as the fourth matrix.

Added: Modal damping when either alpha1 or alpha2 is given.
Added: Check that the expanded solution array size is unchanged.

src/ASM/NewmarkMats.C

@@ -47,10 +47,15 @@ const Matrix& NewmarkMats::getNewtonMatrix () const
   N = A[1];
   N.multiply(alpha_m + alpha_f*alpha1*gamma*h);
   N.add(A[2],alpha_f*(alpha2*gamma + beta*h)*h);
-  if (slvDisp) N.multiply(1.0/(beta*h*h));
+  if (A.size() > 3)
+    N.add(A[3],alpha_f*gamma*h);
+  if (slvDisp)
+    N.multiply(1.0/(beta*h*h));
 #if SP_DEBUG > 2
   std::cout <<"\nElement mass matrix"<< A[1];
   std::cout <<"Element stiffness matrix"<< A[2];
+  if (A.size() > 3)
+    std::cout <<"Element damping matrix"<< A[3];
   std::cout <<"Resulting Newton matrix"<< A[0];
 #endif
 
@@ -66,9 +71,16 @@ const Vector& NewmarkMats::getRHSVector () const
   int iv = vec.size() - 2; // index to element velocity vector (v)
 
 #if SP_DEBUG > 2
+  if (!vec.empty()) std::cout <<"\n";
+  int iu = vec.size() - 3; // index to element displacement vector (u)
+  if (iu >= 0) std::cout <<"u:"<< vec[iu];
+  if (iv >= 0) std::cout <<"v:"<< vec[iv];
+  if (ia >= 0) std::cout <<"a:"<< vec[ia];
   std::cout <<"\nf_ext - f_s"<< dF;
   if (A.size() > 1 && ia >= 0)
     std::cout <<"f_i = M*a"<< A[1]*vec[ia];
+  if (A.size() > 3 && iv >= 0)
+    std::cout <<"f_d = C*v"<< A[3]*vec[iv];
   if (alpha1 > 0.0 && A.size() > 1 && iv >= 0)
     std::cout <<"f_d1/alpha1 = M*v (alpha1="<< alpha1 <<")"<< A[1]*vec[iv];
   if (alpha2 > 0.0 && A.size() > 2 && iv >= 0)
@@ -78,6 +90,9 @@ const Vector& NewmarkMats::getRHSVector () const
   if (A.size() > 1 && ia >= 0)
     dF.add(A[1]*vec[ia],-1.0);    // dF = Fext - M*a
 
+  if (A.size() > 3 && iv >= 0)
+    dF.add(A[3]*vec[iv],-1.0);    // dF -= C*v
+
   if (alpha1 > 0.0 && A.size() > 1 && iv >= 0)
     dF.add(A[1]*vec[iv],-alpha1); // dF -= alpha1*M*v
 

src/ASM/NewmarkMats.h

@@ -30,7 +30,7 @@ public:
   //! \param[in] a2 Stiffness-proportional damping coefficient
   //! \param[in] b Time integration parameter
   //! \param[in] c Time integration parameter
-  //! \param[in] generalizedAlpha If \e true, iterpret \a b and \a c as the
+  //! \param[in] generalizedAlpha If \e true, interpret \a b and \a c as the
   //! generalized alpha parameters \a alpha_m and \a alpha_f, respectively,
   //! otherwise as \a beta and \a gamma
   NewmarkMats(double a1 = 0.0, double a2 = 0.0, double b = 0.0, double c = 0.0,

src/SIM/SIMmodal.C

@@ -80,7 +80,16 @@ bool SIMmodal::swapSystem (AlgEqSystem*& sys, SAM*& sam)
 
 bool SIMmodal::expandSolution (const Vectors& mSol, Vectors& pSol) const
 {
-  pSol.resize(mSol.size());
+  if (pSol.empty())
+    pSol.resize(mSol.size());
+  else if (pSol.size() != mSol.size())
+  {
+    std::cerr <<" *** SIMmodal::expandSolution: Invalid modal solution, "
+              <<" size(mSol) = "<< mSol.size() <<" ("<< pSol.size()
+              <<" expected)."<< std::endl;
+    return false;
+  }
+
   for (size_t i = 0; i < pSol.size(); i++)
     if (mSol[i].size() != myModes.size())
     {
@@ -91,7 +100,19 @@ bool SIMmodal::expandSolution (const Vectors& mSol, Vectors& pSol) const
     }
     else
     {
-      pSol[i].resize(myModes.front().eigVec.size(),true);
+      if (pSol[i].empty())
+        pSol[i].resize(myModes.front().eigVec.size());
+      else if (pSol[i].size() == myModes.front().eigVec.size())
+        pSol[i].fill(0.0);
+      else
+      {
+        std::cerr <<" *** SIMmodal::expandSolution: Logic error, "
+                  <<" size(pSol["<< i <<"]) = "<< pSol[i].size()
+                  <<" != size(eigVec) = "<< myModes.front().eigVec.size()
+                  << std::endl;
+        return false;
+      }
+
       for (size_t j = 0; j < myModes.size(); j++)
         pSol[i].add(myModes[j].eigVec,mSol[i][j]);
     }
@@ -102,7 +123,8 @@ bool SIMmodal::expandSolution (const Vectors& mSol, Vectors& pSol) const
 
 bool SIMmodal::assembleModalSystem (const TimeDomain& time,
                                     const Vectors& pSol, const Vector& Rhs,
-                                    double beta, double gamma)
+                                    double beta, double gamma,
+                                    double alpha1, double alpha2)
 {
   if (!modalSam)
     // Create a diagonal SAM object
@@ -121,7 +143,7 @@ bool SIMmodal::assembleModalSystem (const TimeDomain& time,
     // Create a single-dof element matrix object,
     // the "elements" now being the eigenmodes
     myElmMat = new NewmarkMats(0.0,0.0,beta,gamma);
-    myElmMat->resize(3,1);
+    myElmMat->resize(alpha1 > 0.0 || alpha2 > 0.0 ? 4 : 3, 1);
     myElmMat->redim(1);
   }
 
@@ -144,14 +166,20 @@ bool SIMmodal::assembleModalSystem (const TimeDomain& time,
 
   // Assemble the modal system
   bool ok = true;
+  double omega, damping;
   modalSys->initialize(true);
   for (size_t m = 0; m < myModes.size() && ok; m++)
   {
-    double omega = 2.0*M_PI*myModes[m].eigVal; // Angular eigenfrequency
+    omega = 2.0*M_PI*myModes[m].eigVal; // Angular eigenfrequency
+    if (myModes[m].damping > 0.0)
+      damping = myModes[m].damping;
+    else // Pure Rayleigh damping
+      damping = alpha1 + alpha2*omega*omega;
 
 #if SP_DEBUG > 2
     std::cout <<"\nProcessing mode shape "<< myModes[m].eigNo
-              <<": omega = "<< omega << std::endl;
+              <<": omega = "<< omega
+              <<", damping = "<< damping << std::endl;
 #endif
 
     // Extract current solution associated with the m'th mode shape
@@ -160,6 +188,8 @@ bool SIMmodal::assembleModalSystem (const TimeDomain& time,
 
     myElmMat->A[1].fill(1.0);                        // Modal mass
     myElmMat->A[2].fill(omega*omega);                // Modal stiffness
+    if (myElmMat->A.size() > 3)
+      myElmMat->A[3].fill(damping);                  // Modal damping
     myElmMat->b[0].fill(myModes[m].eigVec.dot(Rhs)); // Modal load
     myElmMat->b[0].add({-omega*omega*pSol[iu][m]});  // Modal residual
     ok = modalSys->assemble(myElmMat,myModes[m].eigNo);

src/SIM/SIMmodal.h

@@ -44,9 +44,12 @@ protected:
   //! \param[in] Rhs Current right-hand-side load vector
   //! \param[in] beta Newmark time integration parameter
   //! \param[in] gamma Newmark time integration parameter
+  //! \param[in] alpha1 Mass-proportional damping factor
+  //! \param[in] alpha2 Stiffness-proportional damping factor
   bool assembleModalSystem(const TimeDomain& time,
                            const Vectors& mSol, const Vector& Rhs,
-                           double beta, double gamma);
+                           double beta, double gamma,
+                           double alpha1 = 0.0, double alpha2 = 0.0);
 
   //! \brief Swaps the modal equation system before/after load vector assembly.
   bool swapSystem(AlgEqSystem*& sys, SAM*& sam);