OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed: Use std::vector<double> as reaction force container

Browse Source
This commit is contained in:
Knut Morten Okstad 2021-11-03 11:11:15 +01:00
parent f4d4c7247a
commit ea0ad9a4fe
8 changed files with 67 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/ASM/AlgEqSystem.C
View File

@ -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<double>* 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++)

13
src/ASM/AlgEqSystem.h
View File

@ -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<double>* 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<SystemVector*> b; //!< The actual right-hand-side vectors
std::vector<double> c; //!< Global scalar quantities
std::vector<double>* d; //!< Multithreading buffer for the scalar values
Vector R; //!< Nodal reaction forces
std::vector<double> R; //!< Nodal reaction forces
const SAM& sam; //!< Data for FE assembly management
const ProcessAdm* adm; //!< Parallel process administrator

43
src/LinAlg/SAM.C
View File

@ -417,7 +417,7 @@ bool SAM::getDofCouplings (std::vector<IntSet>& 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;

36
src/LinAlg/SAM.h
View File

@ -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<int,int>*) { 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

2
src/SIM/NewmarkSIM.C
View File

@ -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"

20
src/SIM/SIMbase.C
View File

@ -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;
}

4
src/SIM/SIMbase.h
View File

@ -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; }

2
src/SIM/SIMoutput.C
View File

@ -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;