Matrix extension

src/ASM/ASMs2D.C

@@ -578,8 +578,12 @@ bool ASMs2D::assignNodeNumbers (BlockNodes& nodes, int basis)
 #if SP_DEBUG > 2
   if (basis > 0) std::cout <<"\nBasis "<< basis <<":";
   for (int i = inod-n1*n2; i < inod; i++)
-    std::cout <<"\nNode "<< i+1 <<"\t: "<< nodeInd[i].I <<" "<< nodeInd[i].J
-	      <<"\tglobal no. "<< MLGN[i];
+  {
+    std::cout <<"\nNode "<< i+1 <<"\t: ";
+    if (!nodeInd.empty())
+      std::cout << nodeInd[i].I <<" "<< nodeInd[i].J;
+    std::cout <<"\tglobal no. "<< MLGN[i];
+  }
   std::cout << std::endl;
 #endif
   return true;

src/ASM/ASMs3D.C

@@ -604,9 +604,12 @@ bool ASMs3D::assignNodeNumbers (BlockNodes& nodes, int basis)
 #if SP_DEBUG > 2
   if (basis > 0) std::cout <<"\nBasis "<< basis <<":";
   for (int i = inod-n1*n2*n3; i < inod; i++)
-    std::cout <<"\nNode "<< i+1 <<"\t: "
-	      << nodeInd[i].I <<" "<< nodeInd[i].J <<" "<< nodeInd[i].K
-	      <<"\tglobal no. "<< MLGN[i];
+  {
+    std::cout <<"\nNode "<< i+1 <<"\t: ";
+    if (!nodeInd.empty())
+      std::cout << nodeInd[i].I <<" "<< nodeInd[i].J <<" "<< nodeInd[i].K;
+    std::cout <<"\tglobal no. "<< MLGN[i];
+  }
   std::cout << std::endl;
 #endif
   return true;

src/ASM/IntegrandBase.h

@@ -59,7 +59,7 @@ class IntegrandBase : public Integrand
 
   //! \brief Defines the solution mode before the element assembly is started.
   virtual void setMode(SIM::SolutionMode mode) { m_mode = mode; }
-  //! \brief Obtain current solution mode
+  //! \brief Returns current solution mode.
   SIM::SolutionMode getMode() const { return m_mode; }
   //! \brief Initializes an integration parameter for the integrand.
   virtual void setIntegrationPrm(unsigned short int, double) {}
@@ -114,7 +114,7 @@ class IntegrandBase : public Integrand
                            const FiniteElement& fe,
                            const Vec3& X0, size_t nPt, LocalIntegral& elmInt);
   //! \brief Initializes current element for numerical integration (mixed).
-  //! \param[in] MNPC Nodal point correspondance for the bases
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] elem_sizes Size of each basis on the element
   //! \param[in] basis_sizes Size of each basis on the patch
   //! \param elmInt Local integral for element
@@ -129,7 +129,7 @@ class IntegrandBase : public Integrand
   virtual bool initElementBou(const std::vector<int>& MNPC,
                               LocalIntegral& elmInt);
   //! \brief Initializes current element for boundary integration (mixed).
-  //! \param[in] MNPC Nodal point correspondance for the bases
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] elem_sizes Size of each basis on the element
   //! \param[in] basis_sizes Size of each basis on the patch
   //! \param elmInt Local integral for element
@@ -153,7 +153,7 @@ class IntegrandBase : public Integrand
   //! \param[in] X Cartesian coordinates of current point
   //! \param[in] MNPC Nodal point correspondance for the basis function values
   virtual bool evalSol(Vector& s, const FiniteElement& fe, const Vec3& X,
-		       const std::vector<int>& MNPC) const;
+                       const std::vector<int>& MNPC) const;
 
   //! \brief Evaluates the secondary solution at a result point (mixed problem).
   //! \param[out] s The solution field values at current point
@@ -163,7 +163,7 @@ class IntegrandBase : public Integrand
   //! \param[in] elem_sizes Size of each basis on the element
   //! \param[in] basis_sizes Size of each basis on the patch
   virtual bool evalSol(Vector& s, const MxFiniteElement& fe, const Vec3& X,
-		       const std::vector<int>& MNPC,
+                       const std::vector<int>& MNPC,
                        const std::vector<size_t>& elem_sizes,
                        const std::vector<size_t>& basis_sizes) const;
 
@@ -172,21 +172,21 @@ class IntegrandBase : public Integrand
   //! \param[in] asol The analytical solution field (tensor field)
   //! \param[in] X Cartesian coordinates of current point
   virtual bool evalSol(Vector& s,
-		       const TensorFunc& asol, const Vec3& X) const;
+                       const TensorFunc& asol, const Vec3& X) const;
 
   //! \brief Evaluates the analytical secondary solution at a result point.
   //! \param[out] s The solution field values at current point
   //! \param[in] asol The analytical solution field (symmetric tensor field)
   //! \param[in] X Cartesian coordinates of current point
   virtual bool evalSol(Vector& s,
-		       const STensorFunc& asol, const Vec3& X) const;
+                       const STensorFunc& asol, const Vec3& X) const;
 
   //! \brief Evaluates the analytical secondary solution at a result point.
   //! \param[out] s The solution field values at current point
   //! \param[in] asol The analytical solution field (vector field)
   //! \param[in] X Cartesian coordinates of current point
   virtual bool evalSol(Vector& s,
-		       const VecFunc& asol, const Vec3& X) const;
+                       const VecFunc& asol, const Vec3& X) const;
 
   //! \brief Returns an evaluated principal direction vector field for plotting.
   virtual bool getPrincipalDir(Matrix&, size_t, size_t) const { return false; }
@@ -234,6 +234,8 @@ class IntegrandBase : public Integrand
 
   //! \brief Accesses the primary solution vector of current patch.
   Vector& getSolution(size_t n = 0) { return primsol[n]; }
+  //! \brief Accesses the primary solution vectors of current patch.
+  Vectors& getSolutions() { return primsol; }
 
   //! \brief Resets the primary solution vectors.
   void resetSolution();
@@ -255,7 +257,6 @@ class IntegrandBase : public Integrand
     return LinAlg::GENERAL_MATRIX;
   }
 
-protected:
   //! \brief Registers a vector to inject a named field into.
   //! \param[in] name Name of field
   //! \param[in] vec Vector to inject field into
@@ -297,6 +298,7 @@ class NormBase : public Integrand
   //! \brief Sets a vector of LocalIntegrals to be used during norm integration.
   void setLocalIntegrals(LintegralVec* elementNorms) { lints = elementNorms; }
 
+  using Integrand::getLocalIntegral;
   //! \brief Returns a local integral container for the element \a iEl.
   virtual LocalIntegral* getLocalIntegral(size_t, size_t iEl, bool) const;
 
@@ -398,6 +400,7 @@ class ForceBase : public Integrand
   //! \brief Initializes the integrand with the number of integration points.
   virtual void initIntegration(size_t, size_t) {}
 
+  using Integrand::getLocalIntegral;
   //! \brief Returns a local integral container for the element \a iEl.
   virtual LocalIntegral* getLocalIntegral(size_t, size_t iEl,
                                           bool = false) const;
@@ -412,7 +415,8 @@ class ForceBase : public Integrand
   { return false; }
 
   //! \brief Dummy implementation (only boundary integration is relevant).
-  virtual bool initElement(const std::vector<int>&, const std::vector<size_t>&,
+  virtual bool initElement(const std::vector<int>&,
+                           const std::vector<size_t>&,
                            const std::vector<size_t>&,
                            LocalIntegral&)
   { return false; }
@@ -421,7 +425,7 @@ class ForceBase : public Integrand
   virtual bool initElementBou(const std::vector<int>& MNPC,
                               LocalIntegral& elmInt);
   //! \brief Initializes current element for boundary integration (mixed).
-  virtual bool initElementBou(const std::vector<int>& MNPC1,
+  virtual bool initElementBou(const std::vector<int>& MNPC,
                               const std::vector<size_t>& elem_sizes,
                               const std::vector<size_t>& basis_sizes,
                               LocalIntegral& elmInt);

src/ASM/NewmarkMats.C

@@ -23,7 +23,7 @@ NewmarkMats::NewmarkMats (double a1, double a2, double b, double c,
 
   if (generalizedAlpha)
   {
-    alpha_m = b;
+    alpha_m = fabs(b);
     alpha_f = c;
     double alpha = alpha_f - alpha_m;
     beta = 0.25*(1.0-alpha)*(1.0-alpha);
@@ -32,9 +32,11 @@ NewmarkMats::NewmarkMats (double a1, double a2, double b, double c,
   else
   {
     alpha_m = alpha_f = 1.0;
-    beta  = b;
+    beta  = fabs(b);
     gamma = c;
   }
+
+  slvDisp = b < 0.0; // Displacement increments are used as primary unknowns
 }
 
 
@@ -45,7 +47,7 @@ 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 SP_DEBUG > 2
   std::cout <<"\nElement mass matrix"<< A[1];
   std::cout <<"Element stiffness matrix"<< A[2];

src/ASM/NewmarkMats.h

@@ -60,6 +60,7 @@ class NewmarkMats : public ElmMats
 private:
   double alpha_m; //!< Generalized-alpha parameter
   double alpha_f; //!< Generalized-alpha parameter
+  bool   slvDisp; //!< If \e true, solve for displacement increments
 };
 
 #endif

src/LinAlg/Test/TestMatrix.C

@@ -17,25 +17,44 @@
 
 TEST(TestMatrix, AddBlock)
 {
-    utl::matrix<int> a(3,3);
-    int data_a[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-    a.fill(data_a, 9);
-
-    utl::matrix<int> b(2,2);
-    int data_b[9] = {1, 2, 3, 4};
-    b.fill(data_b, 4);
-
-    a.addBlock(b, 2, 2, 2, false);
-
-    ASSERT_EQ(a(2,2), 7);
-    ASSERT_EQ(a(3,2), 10);
-    ASSERT_EQ(a(2,3), 14);
-    ASSERT_EQ(a(3,3), 17);
+  utl::matrix<int> a(3,3),b(2,2);
+  std::iota(a.begin(),a.end(),1);
+  std::iota(b.begin(),b.end(),1);
+
+  a.addBlock(b, 2, 2, 2, false);
+  ASSERT_EQ(a(2,2), 7);
+  ASSERT_EQ(a(3,2), 10);
+  ASSERT_EQ(a(2,3), 14);
+  ASSERT_EQ(a(3,3), 17);
+
+  a.addBlock(b, 1, 1, 1, true);
+  ASSERT_EQ(a(1,1), 2);
+  ASSERT_EQ(a(2,1), 5);
+  ASSERT_EQ(a(1,2), 6);
+  ASSERT_EQ(a(2,2), 11);
+}
 
-    a.addBlock(b, 1, 1, 1, true);
 
-    ASSERT_EQ(a(1,1), 2);
-    ASSERT_EQ(a(2,1), 5);
-    ASSERT_EQ(a(1,2), 6);
-    ASSERT_EQ(a(2,2), 11);
+TEST(TestMatrix, AddRows)
+{
+  utl::matrix<int> a(3,5);
+  std::iota(a.begin(),a.end(),1);
+  std::cout <<"A:"<< a;
+
+  a.expandRows(1);
+  std::cout <<"B:"<< a;
+  int fasit = 1;
+  for (size_t j = 1; j <= 5; j++)
+  {
+    for (size_t i = 1; i <= 3; i++, fasit++)
+      ASSERT_EQ(a(i,j), fasit);
+    ASSERT_EQ(a(4,j), 0);
+  }
+
+  a.expandRows(-2);
+  std::cout <<"C:"<< a;
+  fasit = 1;
+  for (size_t j = 1; j <= 5; j++, fasit++)
+    for (size_t i = 1; i <= 2; i++, fasit++)
+      ASSERT_EQ(a(i,j), fasit);
 }