From ea0ad9a4feaa20953275f234998927a898f15a26 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 3 Nov 2021 11:11:15 +0100 Subject: [PATCH] Changed: Use std::vector as reaction force container --- src/ASM/AlgEqSystem.C | 4 ++-- src/ASM/AlgEqSystem.h | 13 ++++++++++--- src/LinAlg/SAM.C | 43 ++++++++++++++++++++++--------------------- src/LinAlg/SAM.h | 36 ++++++++++++++++++------------------ src/SIM/NewmarkSIM.C | 2 +- src/SIM/SIMbase.C | 20 +++++++++++--------- src/SIM/SIMbase.h | 4 ++-- src/SIM/SIMoutput.C | 2 +- 8 files changed, 67 insertions(+), 57 deletions(-) diff --git a/src/ASM/AlgEqSystem.C b/src/ASM/AlgEqSystem.C index 7d171b48..249bb595 100644 --- a/src/ASM/AlgEqSystem.C +++ b/src/ASM/AlgEqSystem.C @@ -140,7 +140,7 @@ void AlgEqSystem::initialize (bool initLHS) } #endif - R.fill(0.0); + std::fill(R.begin(),R.end(),0.0); } @@ -159,7 +159,7 @@ bool AlgEqSystem::assemble (const LocalIntegral* elmObj, int elmId) // The algebraic system consists of one system matrix and one RHS-vector. // Extract the element-level Newton matrix and associated RHS-vector for // general time-dependent and/or nonlinear problems. - Vector* reac = R.empty() ? nullptr : &R; + std::vector* reac = R.empty() ? nullptr : &R; status = sam.assembleSystem(*b.front(), elMat->getRHSVector(), elmId, reac); #if SP_DEBUG > 2 for (i = 1; i < b.size() && i < elMat->b.size(); i++) diff --git a/src/ASM/AlgEqSystem.h b/src/ASM/AlgEqSystem.h index 9da8da7b..95e75d01 100644 --- a/src/ASM/AlgEqSystem.h +++ b/src/ASM/AlgEqSystem.h @@ -15,7 +15,14 @@ #define _ALG_EQ_SYSTEM_H #include "GlobalIntegral.h" -#include "SystemMatrix.h" +#include "LinAlgenums.h" +#include "MatVec.h" + +class SAM; +class ProcessAdm; +class LinSolParams; +class SystemMatrix; +class SystemVector; /*! @@ -79,7 +86,7 @@ public: double getScalar(size_t i = 0) { return i < c.size() ? c[i] : 0.0; } //! \brief Returns a pointer to the nodal reaction forces, if any. - const Vector* getReactions() const { return R.empty() ? 0 : &R; } + const std::vector* getReactions() const { return R.empty() ? 0 : &R; } //! \brief Performs static condensation of the indicated equation system. //! \param[out] Ared Reduced System matrix @@ -120,7 +127,7 @@ private: std::vector b; //!< The actual right-hand-side vectors std::vector c; //!< Global scalar quantities std::vector* d; //!< Multithreading buffer for the scalar values - Vector R; //!< Nodal reaction forces + std::vector R; //!< Nodal reaction forces const SAM& sam; //!< Data for FE assembly management const ProcessAdm* adm; //!< Parallel process administrator diff --git a/src/LinAlg/SAM.C b/src/LinAlg/SAM.C index 4dfeb13a..bef4e315 100644 --- a/src/LinAlg/SAM.C +++ b/src/LinAlg/SAM.C @@ -417,7 +417,7 @@ bool SAM::getDofCouplings (std::vector& dofc) const bool SAM::initForAssembly (SystemMatrix& sysK, SystemVector& sysRHS, - Vector* reactionForces, bool dontLockSP) const + RealArray* reactionForces, bool dontLockSP) const { sysK.initAssembly(*this,dontLockSP); return this->initForAssembly(sysRHS,reactionForces); @@ -431,7 +431,8 @@ bool SAM::initForAssembly (SystemMatrix& sysM) const } -bool SAM::initForAssembly (SystemVector& sysRHS, Vector* reactionForces) const +bool SAM::initForAssembly (SystemVector& sysRHS, + RealArray* reactionForces) const { sysRHS.redim(neq); sysRHS.init(); @@ -443,8 +444,8 @@ bool SAM::initForAssembly (SystemVector& sysRHS, Vector* reactionForces) const bool SAM::assembleSystem (SystemMatrix& sysK, SystemVector& sysRHS, - const Matrix& eK, int iel, - Vector* reactionForces) const + const Matrix& eK, int iel, + RealArray* reactionForces) const { if (reactionForces) { @@ -487,8 +488,8 @@ bool SAM::assembleSystem (SystemMatrix& sysM, bool SAM::assembleSystem (SystemVector& sysRHS, - const Matrix& eK, int iel, - Vector* reactionForces) const + const Matrix& eK, int iel, + RealArray* reactionForces) const { IntVec meen; if (!this->getElmEqns(meen,iel,eK.rows())) @@ -527,8 +528,8 @@ bool SAM::assembleSystem (SystemVector& sysRHS, bool SAM::assembleSystem (SystemVector& sysRHS, - const RealArray& eS, int iel, - Vector* reactionForces) const + const RealArray& eS, int iel, + RealArray* reactionForces) const { Real* sysrhsPtr = sysRHS.getPtr(); int ierr = 0; @@ -554,7 +555,7 @@ bool SAM::assembleSystem (SystemVector& sysRHS, bool SAM::assembleSystem (SystemVector& sysRHS, const Real* nS, int inod, - Vector* reactionForces) const + RealArray* reactionForces) const { IntVec mnen; if (!this->getNodeEqns(mnen,inod)) @@ -572,7 +573,7 @@ bool SAM::assembleSystem (SystemVector& sysRHS, const Real* nS, int inod, { int ipR = -msc[j-1]; if (ipR > 0 && (size_t)ipR <= reactionForces->size()) - (*reactionForces)(ipR) += nS[k]; + (*reactionForces)[ipR-1] += nS[k]; } } @@ -635,7 +636,7 @@ void SAM::assembleRHS (Real* RHS, Real value, int ieq) const } -void SAM::assembleReactions (Vector& reac, const RealArray& eS, int iel) const +void SAM::assembleReactions (RealArray& rf, const RealArray& eS, int iel) const { size_t k = 0; int ip = mpmnpc[iel-1]; @@ -649,8 +650,8 @@ void SAM::assembleReactions (Vector& reac, const RealArray& eS, int iel) const for (int j = madof[node-1]; j < madof[node] && k < eS.size(); j++, k++) { int ipR = -msc[j-1]; - if (ipR > 0 && (size_t)ipR <= reac.size()) - reac(ipR) += eS[k]; + if (ipR > 0 && (size_t)ipR <= rf.size()) + rf[ipR-1] += eS[k]; } } } @@ -931,7 +932,7 @@ Real SAM::normInf (const Vector& x, size_t& comp, char dofType) const } -Real SAM::normReact (const Vector& u, const Vector& rf) const +Real SAM::normReact (const RealArray& u, const RealArray& rf) const { Real retVal = Real(0); @@ -939,11 +940,11 @@ Real SAM::normReact (const Vector& u, const Vector& rf) const if (meqn[i] < 0 && msc[i] < 0 && -msc[i] <= (int)rf.size()) if (mpmceq[-meqn[i]] - mpmceq[-meqn[i]-1] == 1) // Only prescribed DOFs { - retVal += u[i]*rf(-msc[i]); + double RF = rf[-msc[i]-1]; + retVal += u[i]*RF; #if SP_DEBUG > 1 std::cout <<"SAM::normReact: idof="<< i+1 <<" SC="<< msc[i] - <<" u="<< u[i] <<" RF="<< rf(-msc[i]) <<" --> " - << retVal << std::endl; + <<" u="<< u[i] <<" RF="<< RF <<" --> "<< retVal << std::endl; #endif } @@ -966,7 +967,7 @@ bool SAM::haveReaction (int dir, const IntVec* nodes) const } -Real SAM::getReaction (int dir, const Vector& rf, const IntVec* nodes) const +Real SAM::getReaction (int dir, const RealArray& rf, const IntVec* nodes) const { Real retVal = Real(0); @@ -977,14 +978,14 @@ Real SAM::getReaction (int dir, const Vector& rf, const IntVec* nodes) const int idof = madof[i]+dir-2; if (idof < madof[i+1]-1) if (msc[idof] < 0 && -msc[idof] <= (int)rf.size()) - retVal += rf(-msc[idof]); + retVal += rf[-msc[idof]-1]; } return retVal; } -bool SAM::getNodalReactions (int inod, const Vector& rf, Vector& nrf) const +bool SAM::getNodalReactions (int inod, const RealArray& rf, Vector& nrf) const { if (inod < 1 || inod > nnod) { @@ -1000,7 +1001,7 @@ bool SAM::getNodalReactions (int inod, const Vector& rf, Vector& nrf) const if (msc[ip] < 0 && -msc[ip] <= (int)rf.size()) { haveRF = true; - nrf[i] = rf(-msc[ip]); + nrf[i] = rf[-msc[ip]-1]; } return haveRF; diff --git a/src/LinAlg/SAM.h b/src/LinAlg/SAM.h index 3c1777fa..d88bd711 100644 --- a/src/LinAlg/SAM.h +++ b/src/LinAlg/SAM.h @@ -81,9 +81,9 @@ public: //! \return \e false if no free DOFs in the system, otherwise \e true //! //! \details This method must be called once before the first call to - //! \a assembleSystem for a given load case or time step. + //! assembleSystem() for a given load case or time step. bool initForAssembly(SystemMatrix& sysK, SystemVector& sysRHS, - Vector* reactionForces = nullptr, + RealArray* reactionForces = nullptr, bool dontLockSP = false) const; //! \brief Initializes a system matrix prior to the element assembly. @@ -91,7 +91,7 @@ public: //! \return \e false if no free DOFs in the system, otherwise \e true //! //! \details This method must be called once before the first call to - //! \a assembleSystem for a given load case or time step. + //! assembleSystem() for a given load case or time step. bool initForAssembly(SystemMatrix& sysM) const; //! \brief Initializes the system load vector prior to the element assembly. @@ -99,7 +99,7 @@ public: //! \param reactionForces Pointer to vector of nodal reaction forces //! \return \e false if no free DOFs in the system, otherwise \e true bool initForAssembly(SystemVector& sysRHS, - Vector* reactionForces = nullptr) const; + RealArray* reactionForces = nullptr) const; //! \brief Adds an element stiffness matrix into the system stiffness matrix. //! \param sysK The left-hand-side system stiffness matrix @@ -113,7 +113,7 @@ public: //! these are also added into the right-hand-side system load vector. bool assembleSystem(SystemMatrix& sysK, SystemVector& sysRHS, const Matrix& eK, int iel = 0, - Vector* reactionForces = nullptr) const; + RealArray* reactionForces = nullptr) const; //! \brief Adds an element matrix into the corresponding system matrix. //! \param sysM The left-hand-side system matrix @@ -133,7 +133,7 @@ public: //! these are added into the right-hand-side system load vector. bool assembleSystem(SystemVector& sysRHS, const Matrix& eK, int iel = 0, - Vector* reactionForces = nullptr) const; + RealArray* reactionForces = nullptr) const; //! \brief Adds an element load vector into the system load vector. //! \param sysRHS The right-hand-side system load vector @@ -143,7 +143,7 @@ public: //! \return \e true on successful assembly, otherwise \e false bool assembleSystem(SystemVector& sysRHS, const RealArray& eS, int iel = 0, - Vector* reactionForces = nullptr) const; + RealArray* reactionForces = nullptr) const; //! \brief Adds a node load vector into the system load vector. //! \param sysRHS The right-hand-side system load vector @@ -153,7 +153,7 @@ public: //! \return \e true on successful assembly, otherwise \e false bool assembleSystem(SystemVector& sysRHS, const Real* nS, int inod = 0, - Vector* reactionForces = nullptr) const; + RealArray* reactionForces = nullptr) const; //! \brief Adds a node load vector into the system load vector. //! \param sysRHS The right-hand-side system load vector @@ -258,13 +258,13 @@ public: //! \brief Computes the energy norm contributions from nodal reaction forces. //! \param[in] u The (incremental) nodal displacement vector - //! \param[in] rf Compressed reaction force vector for the entire model - virtual Real normReact(const Vector& u, const Vector& rf) const; + //! \param[in] rf Reaction force container for the entire model + virtual Real normReact(const RealArray& u, const RealArray& rf) const; //! \brief Returns the total reaction force in the given coordinate direction. //! \param[in] dir 1-based coordinate direction index - //! \param[in] rf Compressed reaction force vector for the entire model + //! \param[in] rf Reaction force container for the entire model //! \param[in] nodes The set of nodes to consider (null means all nodes) - Real getReaction(int dir, const Vector& rf, + Real getReaction(int dir, const RealArray& rf, const IntVec* nodes = nullptr) const; //! \brief Checks for total reaction force in the given coordinate direction. //! \param[in] dir 1-based coordinate direction index @@ -272,10 +272,10 @@ public: bool haveReaction(int dir, const IntVec* nodes = nullptr) const; //! \brief Returns a vector of reaction forces for a given node. //! \param[in] inod Identifier for the node to get the reaction forces for - //! \param[in] rf Compressed reaction force vector for the entire model + //! \param[in] rf Reaction force container for the entire model //! \param[out] nrf Nodal reaction forces //! \return \e true if the specified node has reaction forces - bool getNodalReactions(int inod, const Vector& rf, Vector& nrf) const; + bool getNodalReactions(int inod, const RealArray& rf, Vector& nrf) const; //! \brief Merges the assembly data from another %SIM with this. virtual bool merge(const SAM*, const std::map*) { return false; } @@ -295,10 +295,10 @@ protected: void assembleRHS(Real* RHS, Real value, int ieq) const; //! \brief Assembles reaction forces for the fixed and prescribed DOFs. - //! \param reac The vector of reaction forces - //! \param[in] eS The element load vector + //! \param rf Reaction force container for the entire model + //! \param[in] eS The element load vector //! \param[in] iel Identifier for the element that \a eS belongs to - void assembleReactions(Vector& reac, const RealArray& eS, int iel) const; + void assembleReactions(RealArray& rf, const RealArray& eS, int iel) const; //! \brief Expands a solution vector from equation-ordering to DOF-ordering. //! \param[in] solVec Pointer to solution vector, length = NEQ @@ -325,7 +325,7 @@ protected: int& nmmnpc; //!< Number of elements in MMNPC int& nmmceq; //!< Number of elements in MMCEQ - // The standard SAM arrays (see K. Bell's reports for detailed explanation). + // The standard %SAM arrays (see K. Bell's reports for detailed explanation). // We are using plane C-pointers for these items such that they more easily // can be passed directly as arguments to FORTRAN subroutines. int* mpmnpc; //!< Matrix of pointers to MNPCs in MMNPC diff --git a/src/SIM/NewmarkSIM.C b/src/SIM/NewmarkSIM.C index 04a0c9d1..8c10e09f 100644 --- a/src/SIM/NewmarkSIM.C +++ b/src/SIM/NewmarkSIM.C @@ -13,7 +13,7 @@ #include "NewmarkSIM.h" #include "SIMoutput.h" -#include "AlgEqSystem.h" +#include "SystemMatrix.h" #include "TimeStep.h" #include "IFEM.h" #include "Profiler.h" diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 9460747d..0c231302 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -22,6 +22,7 @@ #endif #include "IntegrandBase.h" #include "AlgEqSystem.h" +#include "SystemMatrix.h" #include "LinSolParams.h" #include "EigSolver.h" #include "GlbNorm.h" @@ -1717,20 +1718,21 @@ bool SIMbase::solutionNorms (const TimeDomain& time, double SIMbase::externalEnergy (const Vectors& psol, const TimeDomain&) const { - const Vector* reactionFrcs = this->getReactionForces(); - if (!reactionFrcs || !mySam || psol.empty()) return 0.0; + const RealArray* reactionForces = this->getReactionForces(); + if (!reactionForces || !mySam || psol.empty()) return 0.0; // Add norm contributions due to inhomogeneous Dirichlet boundary conditions. // That is, the path integral of the total solution vector times the // reaction forces at the prescribed DOFs. if (psol.size() == 1) - return mySam->normReact(psol.front(),*reactionFrcs); + return mySam->normReact(psol.front(),*reactionForces); if (prevForces.empty()) - extEnergy += 0.5*mySam->normReact(psol[0]-psol[1],*reactionFrcs); + extEnergy += 0.5*mySam->normReact(psol[0]-psol[1],*reactionForces); else - extEnergy += 0.5*mySam->normReact(psol[0]-psol[1],*reactionFrcs+prevForces); - prevForces = *reactionFrcs; + extEnergy += 0.5*mySam->normReact(psol[0]-psol[1], + prevForces.add(*reactionForces)); + prevForces = *reactionForces; return extEnergy; } @@ -1748,7 +1750,7 @@ double SIMbase::externalEnergy (const Vectors& psol, const TimeDomain&) const bool SIMbase::getCurrentReactions (RealArray& RF, const Vector& psol) const { - const Vector* reactionForces = this->getReactionForces(); + const RealArray* reactionForces = this->getReactionForces(); if (!reactionForces || !mySam) return false; RF.resize(1+nsd); @@ -1762,7 +1764,7 @@ bool SIMbase::getCurrentReactions (RealArray& RF, const Vector& psol) const bool SIMbase::getCurrentReactions (RealArray& RF, int pcode) const { - const Vector* reactionForces = this->getReactionForces(); + const RealArray* reactionForces = this->getReactionForces(); if (!reactionForces || !mySam) return false; IntVec glbNodes; @@ -1790,7 +1792,7 @@ bool SIMbase::haveReactions (int pcode) const } -const Vector* SIMbase::getReactionForces () const +const RealArray* SIMbase::getReactionForces () const { return myEqSys ? myEqSys->getReactions() : nullptr; } diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index 5737865c..b454bb2e 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -451,8 +451,8 @@ public: //! \brief Checks for total reaction forces associated with a boundary. //! \param[in] pcode Property code identifying the boundary (0 = all) bool haveReactions(int pcode = 0) const; - //! \brief Returns current reaction force vector. - virtual const Vector* getReactionForces() const; + //! \brief Returns current reaction force container. + virtual const RealArray* getReactionForces() const; //! \brief Computes the total external load of current time/load step. virtual bool getExtLoad(RealArray&, const TimeDomain&) const { return false; } diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C index b99b829c..b6f325ad 100644 --- a/src/SIM/SIMoutput.C +++ b/src/SIM/SIMoutput.C @@ -1493,7 +1493,7 @@ bool SIMoutput::dumpResults (const Vector& psol, double time, if (adm.dd.isPartitioned() && adm.getProcId() != 0) return true; // dump only for procId=0 when domain decomposition - const Vector* reactionForces = this->getReactionForces(); + const RealArray* reactionForces = this->getReactionForces(); RealFunc* psolScl = mySol ? mySol->getScalarSol() : nullptr; VecFunc* psolVec = mySol ? mySol->getVectorSol() : nullptr; VecFunc* ssolScl = mySol ? mySol->getScalarSecSol() : nullptr;