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

Added: Assembly of scalar quantities together with the equation system.

Browse Source 
Added: Debug print of extracted element solution vectors.
Changed: Default value of argument witRF in SIMbase::initSystem to false.
Changed: Removed three spaces in the summary print for scalar solutions.
Fixed: Index error in the second getCurrentReactions method.
This commit is contained in:
Knut Morten Okstad
2017-01-26 14:30:42 +01:00
committed by Arne Morten Kvarving
parent 2250ca5d45
commit 81edcadbe8
11 changed files with 86 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
src/ASM/AlgEqSystem.C
View File

@@ -17,8 +17,9 @@
bool AlgEqSystem::init (SystemMatrix::Type mtype, const LinSolParams* spar,
size_t nmat, size_t nvec, bool withReactions,
LinAlg::LinearSystemType ltype, int num_threads_SLU)
size_t nmat, size_t nvec, size_t nscl,
bool withReactions, LinAlg::LinearSystemType ltype,
int num_threads_SLU)
{
// Using the sign of the num_threads_SLU argument to flag this (convenience)
bool dontLockSparsityPattern = num_threads_SLU < 0;
@@ -32,6 +33,7 @@ bool AlgEqSystem::init (SystemMatrix::Type mtype, const LinSolParams* spar,
A.resize(nmat);
b.resize(nvec,nullptr);
c.resize(nscl,0.0);
R.clear();
for (i = 0; i < A.size(); i++)
@@ -85,6 +87,7 @@ void AlgEqSystem::clear ()
A.clear();
b.clear();
c.clear();
R.clear();
}
@@ -111,6 +114,9 @@ void AlgEqSystem::initialize (bool initLHS)
for (i = 0; i < b.size(); i++)
b[i]->init();
for (i = 0; i < c.size(); i++)
c[i] = 0.0;
R.fill(0.0);
}
@@ -177,6 +183,10 @@ bool AlgEqSystem::assemble (const LocalIntegral* elmObj, int elmId)
status = sam.assembleSystem(*A[i]._A, elMat->A[i], elmId);
}
// Assembly of scalar qunatities
for (i = 0; i < c.size() && i < elMat->c.size(); i++)
c[i] += elMat->c[i];
if (!status)
std::cerr <<" *** AlgEqSystem::assemble: Failure for element "<< elmId
<<": size(A)="<< A.size() <<","<< elMat->A.size()
@@ -207,6 +217,14 @@ bool AlgEqSystem::finalize (bool newLHS)
#if SP_DEBUG > 2
else
std::cout <<"\nSystem right-hand-side vector:"<< *b[i];
if (!c.empty())
{
std::cout <<"\nScalar quantities:";
for (size_t i = 0; i < c.size(); i++)
std::cout <<" "<< c[i];
std::cout << std::endl;
}
#endif
return true;

6
src/ASM/AlgEqSystem.h
View File

@@ -42,11 +42,12 @@ public:
//! \param[in] spar Input parameters for the linear equation solver
//! \param[in] nmat Number of system matrices to allocate
//! \param[in] nvec Number of system vectors to allocate
//! \param[in] nscl Number of scalar quantities to allocate
//! \param[in] withReactions If \e false, no reaction forces will be computed
//! \param[in] ltype Linear system type
//! \param[in] num_threads_SLU Number of threads for SuperLU_MT
bool init(SystemMatrix::Type mtype, const LinSolParams* spar,
size_t nmat, size_t nvec, bool withReactions,
size_t nmat, size_t nvec, size_t nscl, bool withReactions,
LinAlg::LinearSystemType ltype, int num_threads_SLU = 1);
//! \brief Erases the system matrices and frees dynamically allocated storage.
@@ -81,6 +82,8 @@ public:
SystemMatrix* getMatrix(size_t i = 0) { return i < A.size() ? A[i]._A : 0; }
//! \brief Returns the \a i'th right-hand-side vector of the equation system.
SystemVector* getVector(size_t i = 0) { return i < b.size() ? b[i] : 0; }
//! \brief Returns the \a i'th scalar quantity.
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; }
@@ -97,6 +100,7 @@ private:
std::vector<SysMatrixPair> A; //!< The actual coefficient matrices
std::vector<SystemVector*> b; //!< The actual right-hand-side vectors
std::vector<double> c; //!< Global scalar quantities
Vector R; //!< Nodal reaction forces
const SAM& sam; //!< Data for FE assembly management

22
src/ASM/ElmMats.C
View File

@@ -14,6 +14,14 @@
#include "ElmMats.h"
void ElmMats::resize (size_t nA, size_t nB, size_t nC)
{
A.resize(nA);
b.resize(nB);
c.resize(nC);
}
void ElmMats::redim (size_t ndim)
{
for (std::vector<Matrix>::iterator ait = A.begin(); ait != A.end(); ++ait)
@@ -27,10 +35,13 @@ void ElmMats::redim (size_t ndim)
const Matrix& ElmMats::getNewtonMatrix () const
{
if (A.empty())
{
std::cerr <<" *** ElMats::getNewtonMatrix: No element matrices"<< std::endl;
static Matrix empty;
return empty;
}
#if SP_DEBUG > 2
else
std::cout <<"\nElement coefficient matrix"<< A.front();
std::cout <<"\nElement coefficient matrix"<< A.front();
#endif
return A.front();
}
@@ -39,10 +50,13 @@ const Matrix& ElmMats::getNewtonMatrix () const
const Vector& ElmMats::getRHSVector () const
{
if (b.empty())
{
std::cerr <<" *** ElMats::getRHSVector: No element vectors"<< std::endl;
static Vector empty;
return empty;
}
#if SP_DEBUG > 2
else
std::cout <<"\nElement right-hand-side vector"<< b.front();
std::cout <<"\nElement right-hand-side vector"<< b.front();
#endif
return b.front();
}

4
src/ASM/ElmMats.h
View File

@@ -41,7 +41,8 @@ public:
//! \brief Defines the number of element matrices and vectors.
//! \param[in] nA Number of element matrices
//! \param[in] nB Number of element vectors
void resize(size_t nA, size_t nB) { A.resize(nA); b.resize(nB); }
//! \param[in] nC Number of scalar quantities
void resize(size_t nA, size_t nB, size_t nC = 0);
//! \brief Sets the dimension of the element matrices and vectors.
//! \param[in] ndim Number of rows and columns in the matrices/vectors
@@ -59,6 +60,7 @@ public:
std::vector<Matrix> A; //!< The element coefficient matrices
std::vector<Vector> b; //!< The element right-hand-side vectors
std::vector<double> c; //!< The scalar quantities
bool rhsOnly; //!< If \e true, only the right-hand-sides are assembled
bool withLHS; //!< If \e true, left-hand-side element matrices are present

9
src/ASM/IntegrandBase.C
View File

@@ -51,6 +51,11 @@ bool IntegrandBase::initElement (const std::vector<int>& MNPC,
if (!primsol[i].empty())
ierr = utl::gather(MNPC,npv,primsol[i],elmInt.vec[i]);
#if SP_DEBUG > 2
for (size_t j = 0; j < elmInt.vec.size(); j++)
std::cout <<"Element solution vector "<< j+1 << elmInt.vec[j];
#endif
if (ierr == 0) return true;
std::cerr <<" *** IntegrandBase::initElement: Detected "
@@ -97,6 +102,10 @@ bool IntegrandBase::initElementBou (const std::vector<int>& MNPC,
if (!primsol.empty() && !primsol.front().empty())
ierr = utl::gather(MNPC,npv,primsol.front(),elmInt.vec.front());
#if SP_DEBUG > 2
std::cout <<"Element solution vector"<< elmInt.vec.front();
#endif
if (ierr == 0) return true;
std::cerr <<" *** IntegrandBase::initElementBou: Detected "

2
src/SIM/AdaptiveSIM.C
View File

@@ -204,7 +204,7 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep)
// Assemble the linear FE equation system
model.setMode(SIM::STATIC,true);
model.initSystem(opt.solver, 1, model.getNoRHS());
model.initSystem(opt.solver,1,model.getNoRHS(),0,true);
model.setQuadratureRule(opt.nGauss[0],true);
if (!model.assembleSystem())

2
src/SIM/EigenModeSIM.C
View File

@@ -91,7 +91,7 @@ bool EigenModeSIM::initSol (size_t nSol)
// Solve the eigenvalue problem giving the natural eigenfrequencies
model.setMode(SIM::VIBRATION);
model.initSystem(opt.solver,2,0,false);
model.initSystem(opt.solver,2,0);
model.setQuadratureRule(opt.nGauss[0],true);
if (!model.assembleSystem())
return false;

4
src/SIM/MultiStepSIM.C
View File

@@ -64,9 +64,9 @@ bool MultiStepSIM::initSol (size_t nSol)
}
bool MultiStepSIM::initEqSystem (bool withRF)
bool MultiStepSIM::initEqSystem (bool withRF, size_t nScl)
{
return model.initSystem(opt.solver,1,nRHSvec,withRF);
return model.initSystem(opt.solver,1,nRHSvec,nScl,withRF);
}

3
src/SIM/MultiStepSIM.h
View File

@@ -56,7 +56,8 @@ public:
//! \brief Allocates the FE system matrices.
//! \param[in] withRF Whether nodal reaction forces should be computed or not
bool initEqSystem(bool withRF = true);
//! \param[in] nScl Number of global scalar quantities to integrate
bool initEqSystem(bool withRF = true, size_t nScl = 0);
//! \brief Advances the time/load step one step forward.
//! \param param Time stepping parameters

21
src/SIM/SIMbase.C
View File

@@ -397,7 +397,8 @@ RealFunc* SIMbase::getSclFunc (int code) const
}
bool SIMbase::initSystem (int mType, size_t nMats, size_t nVec, bool withRF)
bool SIMbase::initSystem (int mType, size_t nMats, size_t nVec, size_t nScl,
bool withRF)
{
if (!mySam) return true; // Silently ignore when no algebraic system
@@ -439,7 +440,7 @@ bool SIMbase::initSystem (int mType, size_t nMats, size_t nVec, bool withRF)
}
return myEqSys->init(static_cast<SystemMatrix::Type>(mType),
mySolParams, nMats, nVec, withRF,
mySolParams, nMats, nVec, nScl, withRF,
myProblem->getLinearSystemType(), opt.num_threads_SLU);
}
@@ -890,6 +891,12 @@ bool SIMbase::extractLoadVec (Vector& loadVec) const
}
double SIMbase::extractScalar (size_t i) const
{
return myEqSys ? myEqSys->getScalar(i) : 0.0;
}
bool SIMbase::applyDirichlet (Vector& glbVec) const
{
return mySam->applyDirichlet(glbVec);
@@ -998,8 +1005,10 @@ void SIMbase::printSolutionSummary (const Vector& solution, int printSol,
if (compName)
adm.cout <<"\n >>> Solution summary <<<\n\nL2-norm : ";
else
else if (nf > 1)
adm.cout <<" Primary solution summary: L2-norm : ";
else
adm.cout <<" Primary solution summary: L2-norm : ";
adm.cout << utl::trunc(dNorm);
if (nf == 1 && utl::trunc(dMax[0]) != 0.0)
@@ -1007,7 +1016,7 @@ void SIMbase::printSolutionSummary (const Vector& solution, int printSol,
if (compName)
adm.cout <<"\nMax "<< compName <<" : ";
else
adm.cout <<"\n Max value : ";
adm.cout <<"\n Max value : ";
adm.cout << dMax[0] <<" node "<< iMax[0];
}
else if (nf > 1)
@@ -1386,7 +1395,7 @@ bool SIMbase::getCurrentReactions (RealArray& RF, const Vector& psol) const
RF.resize(1+nsd);
RF.front() = 2.0*mySam->normReact(psol,*reactionForces);
for (size_t dir = 1; dir < RF.size(); dir++)
for (unsigned char dir = 1; dir <= nsd; dir++)
RF[dir] = mySam->getReaction(dir,*reactionForces);
return true;
@@ -1406,7 +1415,7 @@ bool SIMbase::getCurrentReactions (RealArray& RF, int pcode) const
RF.resize(nsd);
for (unsigned char dir = 1; dir <= nsd; dir++)
RF[dir] = mySam->getReaction(dir,*reactionForces,&glbNodes);
RF[dir-1] = mySam->getReaction(dir,*reactionForces,&glbNodes);
return true;
}

17
src/SIM/SIMbase.h
View File

@@ -91,9 +91,10 @@ public:
//! \param[in] mType The matrix format to use
//! \param[in] nMats Number of system matrices
//! \param[in] nVec Number of system right-hand-side vectors
//! \param[in] nScl Number of global scalar quantities
//! \param[in] withRF Whether nodal reaction forces should be computed or not
bool initSystem(int mType, size_t nMats = 1, size_t nVec = 1,
bool withRF = true);
bool initSystem(int mType, size_t nMats = 1, size_t nVec = 1, size_t nScl = 0,
bool withRF = false);
//! \brief Associates a system vector to a system matrix.
//! \sa AlgEqSystem::setAssociatedVector
@@ -240,6 +241,8 @@ public:
//! \brief Extracts the assembled load vector for inspection/visualization.
//! \param[out] loadVec Global load vector in DOF-order
bool extractLoadVec(Vector& loadVec) const;
//! \brief Extracts an assembled global scalar quantity.
double extractScalar(size_t i = 0) const;
//! \brief Applies the Dirichlet conditions to given vector.
//! \param[out] glbVec Global vector in DOF-order
@@ -280,7 +283,7 @@ public:
//! \param[in] eps Only record the energies larger than this tolerance
//! \param[out] worst Node and local DOF number and values of the worst DOFs
void getWorstDofs(const Vector& x, const Vector& r, size_t nWorst, double eps,
std::map<std::pair<int,int>,RealArray>& worst) const;
std::map<std::pair<int,int>,RealArray>& worst) const;
//! \brief Evaluates some iteration norms for convergence assessment.
//! \param[in] x Global primary solution vector
@@ -289,7 +292,7 @@ public:
//! \param[out] rNorm Residual norm of solution increment
//! \param[out] dNorm Displacement norm of solution increment
void iterationNorms(const Vector& x, const Vector& r,
double& eNorm, double& rNorm, double& dNorm) const;
double& eNorm, double& rNorm, double& dNorm) const;
//! \brief Evaluates some norms of the primary solution vector
//! \param[in] x Global primary solution vector
@@ -335,7 +338,7 @@ public:
//!
//! \details Use this version for linear/stationary problems only.
bool solutionNorms(const Vectors& psol, const Vectors& ssol,
Matrix& eNorm, Vectors& gNorm)
Matrix& eNorm, Vectors& gNorm)
{ return this->solutionNorms(TimeDomain(),psol,ssol,gNorm,&eNorm); }
//! \brief Integrates some solution norm quantities.
//! \param[in] psol Primary solution vectors
@@ -381,8 +384,8 @@ public:
//! \param[in] iA Index of system matrix \b A in \a myEqSys->A
//! \param[in] iB Index of system matrix \b B in \a myEqSys->A
bool systemModes(std::vector<Mode>& solution,
int nev, int ncv, int iop, double shift,
size_t iA = 0, size_t iB = 1);
int nev, int ncv, int iop, double shift,
size_t iA = 0, size_t iB = 1);
//! \brief Performs a generalized eigenvalue analysis of the assembled system.
//! \param[out] solution Computed eigenvalues and associated eigenvectors
//! \param[in] iA Index of system matrix \b A in \a myEqSys->A