Added: New GlobalIntegral sub-class for integration of nodal forces

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2358 e10b68d5-8a6e-419e-a041-bce267b0401d

src/ASM/GlbForceVec.C

@@ -0,0 +1,103 @@
+// $Id$
+//==============================================================================
+//!
+//! \file GlbForceVec.C
+//!
+//! \date Dec 13 2012
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Storage of a global nodal force vector for a FEM problem.
+//!
+//==============================================================================
+
+#include "GlbForceVec.h"
+#include "ElmMats.h"
+#include "Vec3.h"
+#include "SAM.h"
+#if SP_DEBUG > 1
+#include "Vec3Oper.h"
+#endif
+
+
+bool GlbForceVec::initNodeMap (const std::vector<int>& globalNodes, size_t nfc)
+{
+  nodeMap.clear();
+
+  int illegal = 0, nnod = sam.getNoNodes();
+  for (size_t i = 0; i < globalNodes.size(); i++)
+    if (globalNodes[i] > 0 && globalNodes[i] <= nnod)
+      nodeMap[globalNodes[i]] = i+1;
+    else
+      illegal++;
+
+  F.resize(nfc,nodeMap.size());
+  if (illegal == 0) return true;
+
+  std::cerr <<" *** GlbForceVec::initNodeMap: Detected "<< illegal
+            <<" node numbers out of range."<< std::endl;
+  return false;
+}
+
+
+void GlbForceVec::initialize (bool)
+{
+  F.fill(0.0);
+}
+
+
+bool GlbForceVec::assemble (const LocalIntegral* elmObj, int elmId)
+{
+  const ElmMats* elm = dynamic_cast<const ElmMats*>(elmObj);
+  if (!elm || elm->b.size() < 1) return false;
+
+  const Vector& ES = elm->b.front();
+  const size_t nfc = F.rows();
+  std::vector<int> mnpc;
+  if (!sam.getElmNodes(mnpc,elmId))
+    return false;
+  else if (ES.size() < nfc*mnpc.size())
+  {
+    std::cerr <<" *** GlbForceVec::assemble: Invalid element force vector,"
+              <<" size="<< ES.size() <<" should be (at least) "
+              << nfc*mnpc.size() << std::endl;
+    return false;
+  }
+
+  // Assemble the nodal forces into the Matrix F
+  size_t i, j, k, ninod = 0;
+  std::map<int,size_t>::const_iterator nit;
+  for (i = k = 0; i < mnpc.size(); i++, k += nfc)
+    if ((nit = nodeMap.find(mnpc[i])) == nodeMap.end())
+      ninod++;
+    else for (j = 0; j < nfc; j++)
+      F(j+1,nit->second) += ES[k+j];
+
+  if (ninod < mnpc.size())
+    return true;
+
+  std::cerr <<" *** GlbForceVec::assemble: Element "<< elmId
+            <<" has no nodal force contributions on this boundary"<< std::endl;
+  return false;
+}
+
+
+bool GlbForceVec::finalize (bool)
+{
+  // TODO: Add MPI reduction of the global F here
+  return true;
+}
+
+
+Vec3 GlbForceVec::getForce (int node) const
+{
+  Vec3 force;
+  std::map<int,size_t>::const_iterator nit = nodeMap.find(node);
+  if (nit != nodeMap.end())
+    force = Vec3(F.getColumn(nit->second));
+
+#if SP_DEBUG > 1
+  std::cout <<"Force in node "<< node <<": "<< force << std::endl;
+#endif
+  return force;
+}

src/ASM/GlbForceVec.h

@@ -0,0 +1,63 @@
+// $Id$
+//==============================================================================
+//!
+//! \file GlbForceVec.h
+//!
+//! \date Dec 13 2012
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Storage of a global nodal force vector for a FEM problem.
+//!
+//==============================================================================
+
+#ifndef _GLB_FORCE_VEC_H
+#define _GLB_FORCE_VEC_H
+
+#include "GlobalIntegral.h"
+#include "MatVec.h"
+#include <map>
+
+class LocalIntegral;
+class Vec3;
+class SAM;
+
+
+/*!
+  \brief Class for storage of a global nodal force vector with assembly methods.
+*/
+
+class GlbForceVec : public GlobalIntegral
+{
+public:
+  //! \brief The constructor only sets its reference to the SAM object.
+  GlbForceVec(const SAM& _sam) : sam(_sam) {}
+  //! \brief Empty destructor.
+  virtual ~GlbForceVec() {}
+
+  //! \brief Initializes the global node map and allocates the force vector.
+  //! \param[in] globalNodes The global node numbers that will have force terms
+  //! \param[in] nfc Number of force components
+  bool initNodeMap(const std::vector<int>& globalNodes, size_t nfc);
+
+  //! \brief Initializes the global nodal force vector to zero.
+  virtual void initialize(bool = false);
+  //! \brief Finalizes the global nodal force vector after element assembly.
+  virtual bool finalize(bool = false);
+
+  //! \brief Adds a set of element nodal forces into the global nodal forces.
+  //! \param[in] elmObj Pointer to the element nodal forces to add into \a *this
+  //! \param[in] elmId Global number of the element associated with \a *elmObj
+  virtual bool assemble(const LocalIntegral* elmObj, int elmId);
+
+  //! \brief Returns the global nodal force vector for a specified node.
+  //! \param[in] node 1-based global node number to return the forces for
+  Vec3 getForce(int node) const;
+
+private:
+  const SAM&           sam;     //!< Data for FE assembly management
+  Matrix               F;       //!< Global nodal forces
+  std::map<int,size_t> nodeMap; //!< Maps from global node number to force index
+};
+
+#endif

src/LinAlg/SAM.h

@@ -57,7 +57,7 @@ public:
   virtual int getNoEquations() const { return neq; }
   //! \brief Returns the Matrix of Accumulated DOFs.
   const int* getMADOF() const { return madof; }
-  //! \brief Returns the Vector of global equation numbers
+  //! \brief Returns the Matrix of EQuation Numbers.
   const int* getMEQN() const { return meqn; }
 
   //! \brief Returns max number of DOF couplings in the model.
@@ -163,7 +163,7 @@ public:
 			      const Real* nS, int inod = 0,
 			      Vector* reactionForces = 0) const;
 
-  //! \brief Finds the matrix of equation numbers for an element.
+  //! \brief Finds the matrix of nodal point correspondance for an element.
   //! \param[out] mnpc Matrix of nodal point correspondance
   //! \param[in] iel Identifier for the element to get the node numbers for
   bool getElmNodes(IntVec& mnpc, int iel) const;
@@ -225,9 +225,10 @@ public:
   //! \details This version is typically used to expand eigenvectors.
   bool expandVector(const Vector& solVec, Vector& dofVec) const;
 
-  //! \brief Apply non-homogenous Dirichlet B.C's to vector
-  //! \param[out] dofVec Degrees of freedom vector, length = NDOF
-  //! \details This is typically used with explicit time integration.
+  //! \brief Applies the non-homogenous Dirichlet BCs to the given vector.
+  //! \param dofVec Degrees of freedom vector, length = NDOF
+  //!
+  //! \details This method is typically used with explicit time integration.
   bool applyDirichlet(Vector& dofVec) const;
 
   //! \brief Computes the dot-product of two vectors of length NDOF.
@@ -308,7 +309,7 @@ protected:
   int*  minex;  //!< Matrix of internal to external node numbers
   int*  meqn;   //!< Matrix of equation numbers
 
-  std::vector<char>   nodeType;     //!< Nodal DOF classification
+  std::vector<char> nodeType; //!< Nodal DOF classification
 
   friend class DenseMatrix;
   friend class SPRMatrix;