Skip to content

Commit 0f691a7

Browse files
kmokstadkmokstad
committed
Added: Print node and local dof number of detected singular equation.
Added: Protected method SAM::printMEQN(). Changed: Reduce scope of some temporaries in various SAM methods.
1 parent c47e847 commit 0f691a7

File tree

4 files changed

+105
-89
lines changed

4 files changed

+105
-89
lines changed

‎src/LinAlg/SAM.C‎

Lines changed: 62 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -98,13 +98,21 @@ SAM::~SAM ()
9898
}
9999

100100

101-
std::pair<int,int> SAM::getNodeAndLocalDof (int idof) const
101+
std::pair<int,int> SAM::getNodeAndLocalDof (int idof, bool eqno) const
102102
{
103+
if (eqno)
104+
{
105+
if (int* eq = std::find(meqn,meqn+ndof,idof); eq != meqn+ndof)
106+
idof = std::distance(meqn,eq);
107+
else
108+
return { 0, 0 };
109+
}
110+
103111
for (int n = 1; n <= nnod; n++)
104112
if (madof[n] > idof)
105-
return std::make_pair(minex ? minex[n-1] : n, idof-madof[n-1]+1);
113+
return { minex ? minex[n-1] : n, idof-madof[n-1]+1 };
106114

107-
return std::make_pair(0,0);
115+
return { 0, 0 };
108116
}
109117

110118

@@ -138,7 +146,6 @@ void SAM::print (std::ostream& os) const
138146
os <<"\n\nSAM::mpar: "<< mpar[0];
139147
for (int i = 1; i < 30; i++)
140148
os << ((i%10) ? " " : "\n ") << mpar[i];
141-
os << std::endl;
142149

143150
if (mmnpc && mpmnpc && nel > 0)
144151
{
@@ -153,7 +160,6 @@ void SAM::print (std::ostream& os) const
153160
else if (mmnpc[i-1] < 0)
154161
os <<' '<< (minex ? -minex[-mmnpc[i-1]-1] : mmnpc[i-1]);
155162
}
156-
os << std::endl;
157163
}
158164

159165
if (ttcc && mmceq && mpmceq && nceq > 0)
@@ -164,13 +170,19 @@ void SAM::print (std::ostream& os) const
164170
os <<'\n'<< std::setw(4) << i+1 <<": ";
165171
this->printCEQ(os,i);
166172
}
167-
os << std::endl;
168173
}
174+
os << std::endl;
175+
176+
this->printMEQN(os);
177+
}
169178

179+
180+
void SAM::printMEQN (std::ostream& os) const
181+
{
170182
if (meqn && madof && nnod > 0 && neq > 0)
171183
{
172184
char code[7], dofType[7];
173-
os <<"\n\nNode --> Equations";
185+
os <<"\nNode --> Equations";
174186
for (int n = 0; n < nnod; n++)
175187
{
176188
int dof, i = 0;
@@ -195,22 +207,25 @@ void SAM::print (std::ostream& os) const
195207

196208
void SAM::printStatusCodes (std::ostream& os) const
197209
{
198-
os <<"\nNode\tDOFs\t MSC";
199-
for (int i = 0; i < nnod; i++)
210+
if (msc && madof && nnod > 0)
200211
{
201-
os <<"\n"<< std::setw(4) << 1+i
202-
<<"\t"<< madof[i] <<" - "<< madof[i+1]-1 <<"\t";
203-
for (int j = madof[i]; j < madof[i+1]; j++)
204-
os <<" "<< msc[j-1];
212+
os <<"\nNode\tDOFs\t MSC";
213+
for (int i = 0; i < nnod; i++)
214+
{
215+
os <<"\n"<< std::setw(4) << 1+i
216+
<<"\t"<< madof[i] <<" - "<< madof[i+1]-1 <<"\t";
217+
for (int j = madof[i]; j < madof[i+1]; j++)
218+
os <<" "<< msc[j-1];
219+
}
220+
os << std::endl;
205221
}
206-
os << std::endl;
207222
}
208223

209224

210225
bool SAM::initSystemEquations ()
211226
{
212227
#ifdef SP_DEBUG
213-
std::cout <<"SAM::initSystemEquations()"<< std::endl;
228+
std::cout <<"\nSAM::initSystemEquations()"<< std::endl;
214229
#endif
215230
if (!msc && ndof > 0) return false;
216231
if ((!msc || !mpmceq) && nceq > 0) return false;
@@ -307,6 +322,9 @@ bool SAM::initSystemEquations ()
307322
else if (msc[idof] == 2)
308323
meqn[idof] = j++;
309324
#endif
325+
#if SP_DEBUG > 1
326+
this->printMEQN(std::cout);
327+
#endif
310328

311329
if (ierr == 0) return true;
312330

@@ -365,20 +383,15 @@ void SAM::getDofCouplings (std::vector<IntSet>& dofc) const
365383
this->getElmEqns(meen,iel);
366384

367385
for (size_t j = 0; j < meen.size(); j++)
368-
{
369-
int jeq = meen[j];
370-
if (jeq > 0)
386+
if (int jeq = meen[j]; jeq > 0)
371387
{
372388
dofc[jeq-1].insert(jeq);
373389
for (size_t i = 0; i < j; i++)
374-
{
375-
int ieq = meen[i];
376-
if (ieq > 0)
390+
if (int ieq = meen[i]; ieq > 0)
377391
{
378392
dofc[ieq-1].insert(jeq);
379393
dofc[jeq-1].insert(ieq);
380394
}
381-
}
382395
}
383396
else if (jeq < 0)
384397
{
@@ -410,7 +423,6 @@ void SAM::getDofCouplings (std::vector<IntSet>& dofc) const
410423
}
411424
}
412425
}
413-
}
414426
}
415427
}
416428

@@ -514,7 +526,7 @@ bool SAM::assembleSystem (SystemVector& sysRHS,
514526
#ifdef USE_F77SAM
515527
int* work = new int[eS.size()];
516528
addev2_(&eS.front(), ttcc, mpar, madof, meqn, mpmnpc, mmnpc, mpmceq, mmceq,
517-
iel, eS.size(), 6, sysRHS.getPtr(), work, ierr);
529+
iel, eS.size(), 6, sysRHS.getPtr(), work, ierr);
518530
delete[] work;
519531
#else
520532
IntVec meen;
@@ -543,13 +555,10 @@ bool SAM::assembleSystem (SystemVector& sysRHS, const Real* nS, int inod,
543555

544556
if (reactionForces)
545557
{
546-
int k = 0;
558+
int k = 0, nReact = reactionForces->size();
547559
for (int j = madof[inod-1]; j < madof[inod]; j++, k++)
548-
{
549-
int ipR = -msc[j-1];
550-
if (ipR > 0 && ipR <= static_cast<int>(reactionForces->size()))
560+
if (int ipR = -msc[j-1]; ipR > 0 && ipR <= nReact)
551561
(*reactionForces)[ipR-1] += nS[k];
552-
}
553562
}
554563

555564
return true;
@@ -614,18 +623,12 @@ void SAM::assembleReactions (RealArray& rf, const RealArray& eS, int iel) const
614623
int ip = mpmnpc[iel-1];
615624
int nenod = mpmnpc[iel] - ip;
616625
for (int i = 0; i < nenod; i++, ip++)
617-
{
618-
int node = mmnpc[ip-1];
619-
if (node < 0)
626+
if (int node = mmnpc[ip-1]; node < 0)
620627
k += madof[-node] - madof[-node-1];
621628
else if (node > 0)
622629
for (int j = madof[node-1]; j < madof[node] && k < eS.size(); j++, k++)
623-
{
624-
int ipR = -msc[j-1];
625-
if (ipR > 0 && ipR <= static_cast<int>(rf.size()))
630+
if (int ipR = -msc[j-1]; ipR > 0 && ipR <= static_cast<int>(rf.size()))
626631
rf[ipR-1] += eS[k];
627-
}
628-
}
629632
}
630633

631634

@@ -678,13 +681,10 @@ bool SAM::getElmEqns (IntVec& meen, int iel, size_t nedof) const
678681
#else
679682
meen.reserve(nedof > 0 ? nedof : neldof);
680683
for (int i = 0; i < nenod; i++, ip++)
681-
{
682-
int node = mmnpc[ip];
683-
if (node > 0)
684+
if (int node = mmnpc[ip]; node > 0)
684685
meen.insert(meen.end(),meqn+madof[node-1]-1,meqn+madof[node]-1);
685686
else if (node < 0)
686687
meen.insert(meen.end(),madof[-node]-madof[-node-1],0);
687-
}
688688
neldof = meen.size();
689689
#endif
690690
if (nedof == 0 || neldof == static_cast<int>(nedof)) return true;
@@ -789,25 +789,19 @@ bool SAM::expandVector (const Real* solVec, Vector& dofVec, Real scaleSD) const
789789
dofVec.resize(ndof,true);
790790
#ifdef USE_F77SAM
791791
expand_(solVec, ttcc, mpmceq, mmceq, meqn, Real(1), scaleSD, ndof, neq,
792-
dofVec.ptr());
792+
dofVec.ptr());
793793
#else
794794
for (int idof = 0; idof < ndof; idof++)
795-
{
796-
int ieq = meqn[idof];
797-
int iceq = -ieq;
798-
if (ieq > 0)
795+
if (int ieq = meqn[idof]; ieq > 0)
799796
dofVec[idof] += solVec[ieq-1];
800-
else if (iceq > 0)
797+
else if (int iceq = -meqn[idof]; iceq > 0)
801798
{
802799
int ip = mpmceq[iceq-1];
803800
dofVec[idof] += scaleSD*ttcc[ip-1];
804801
for (; ip < mpmceq[iceq]-1; ip++)
805-
if (mmceq[ip] > 0)
806-
{
807-
ieq = meqn[mmceq[ip]-1];
808-
if (ieq > 0 && ieq <= neq)
809-
dofVec[idof] += ttcc[ip]*solVec[ieq-1];
810-
}
802+
if (mmceq[ip] > 0)
803+
if (int ieq = meqn[mmceq[ip]-1]; ieq > 0 && ieq <= neq)
804+
dofVec[idof] += ttcc[ip]*solVec[ieq-1];
811805
}
812806
}
813807
#endif
@@ -822,11 +816,8 @@ bool SAM::getSolVec (RealArray& solVec, const RealArray& dofVec) const
822816

823817
solVec.resize(neq,Real(0));
824818
for (int idof = 0; idof < ndof; idof++)
825-
{
826-
int ieq = meqn[idof];
827-
if (ieq > 0)
819+
if (int ieq = meqn[idof]; ieq > 0)
828820
solVec[ieq-1] = dofVec[idof];
829-
}
830821

831822
return true;
832823
}
@@ -837,16 +828,13 @@ bool SAM::applyDirichlet (Vector& dofVec) const
837828
if (!meqn) return false;
838829

839830
for (int idof = 0; idof < ndof; idof++)
840-
{
841-
int iceq = -meqn[idof];
842-
if (iceq > 0)
831+
if (int iceq = -meqn[idof]; iceq > 0)
843832
{
844833
int ip = mpmceq[iceq-1];
845834
dofVec[idof] = ttcc[ip-1];
846835
}
847836
else if (iceq == 0)
848837
dofVec[idof] = Real(0);
849-
}
850838

851839
return true;
852840
}
@@ -900,8 +888,7 @@ Real SAM::normInf (const Vector& x, size_t& comp, char dofType) const
900888
else if (nodeType.empty() && dof_type.empty())
901889
{
902890
// All DOFs are of the same type and the number of nodal DOFs is constant
903-
int nndof = madof[1] - madof[0];
904-
if (static_cast<int>(comp) <= nndof)
891+
if (int nndof = madof[1] - madof[0]; static_cast<int>(comp) <= nndof)
905892
return x.normInf(--comp,nndof);
906893
else
907894
return Real(0);
@@ -923,15 +910,12 @@ Real SAM::normInf (const Vector& x, size_t& comp, char dofType) const
923910
}
924911
}
925912
else if (nodeType[i] == dofType)
926-
{
927-
int idof = madof[i] + icmp - 2;
928-
if (idof >= 0 && idof < madof[i+1]-1)
913+
if (int idof = madof[i] + icmp - 2; idof >= 0 && idof < madof[i+1]-1)
929914
if (fabs(x[idof]) > retVal)
930915
{
931916
comp = i+1;
932917
retVal = fabs(x[idof]);
933918
}
934-
}
935919

936920
return retVal;
937921
}
@@ -962,11 +946,9 @@ bool SAM::haveReaction (int dir, const IntVec* nodes) const
962946
if (dir > 0 && nspdof > 0)
963947
for (int i = 0; i < nnod; i++)
964948
if (!nodes || std::find(nodes->begin(),nodes->end(),i+1) != nodes->end())
965-
{
966-
int idof = madof[i]+dir-2;
967-
if (idof < madof[i+1]-1 && msc[idof] < 0 && -msc[idof] <= nspdof)
968-
return true;
969-
}
949+
if (int idof = madof[i]+dir-2; idof < madof[i+1]-1)
950+
if (msc[idof] < 0 && -msc[idof] <= nspdof)
951+
return true;
970952

971953
return false;
972954
}
@@ -979,16 +961,14 @@ Real SAM::getReaction (int dir, const RealArray& rf, const IntVec* nodes) const
979961
if (dir > 0 && nspdof > 0)
980962
for (int i = 0; i < nnod; i++)
981963
if (!nodes || std::find(nodes->begin(),nodes->end(),i+1) != nodes->end())
982-
{
983-
int idof = madof[i]+dir-2;
984-
int ipr = idof < madof[i+1]-1 ? -msc[idof]-1 : -1;
985-
if (ipr >= 0 && ipr < static_cast<int>(rf.size()))
986-
{
987-
retVal += rf[ipr];
988-
if (nodes) // clear this force component to avoid it being added twice
989-
const_cast<RealArray&>(rf)[ipr] = 0.0;
990-
}
991-
}
964+
if (int idof = madof[i]+dir-2; idof < madof[i+1]-1)
965+
if (int ipr = -msc[idof]-1; ipr >= 0)
966+
if (ipr < static_cast<int>(rf.size()))
967+
{
968+
retVal += rf[ipr];
969+
if (nodes) // clear force component to avoid it being added twice
970+
const_cast<RealArray&>(rf)[ipr] = 0.0;
971+
}
992972

993973
return retVal;
994974
}

‎src/LinAlg/SAM.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -191,9 +191,10 @@ class SAM
191191

192192
//! \brief Returns the internal node number and local index for a global DOF.
193193
//! \param[in] idof Global DOF-number in the range [1,NDOF]
194+
//! \param[in] eqno If \e true, \a idof is the 1-based global equation number
194195
//! \return first = internal node number in the range [1,NNOD]
195196
//! \return second = local DOF number in the range [1,NNDOF]
196-
std::pair<int,int> getNodeAndLocalDof(int idof) const;
197+
std::pair<int,int> getNodeAndLocalDof(int idof, bool eqno = false) const;
197198

198199
//! \brief Finds the equation number corresponding to a local nodal DOF.
199200
//! \param[in] inod Identifier for the node to get the equation number for
@@ -318,6 +319,8 @@ class SAM
318319
void printStatusCodes(std::ostream& os) const;
319320
//! \brief Prints out a constraint equation to the given stream.
320321
void printCEQ(std::ostream& os, int iceq) const;
322+
//! \brief Prints out the equation number mapping to the given stream.
323+
void printMEQN(std::ostream& os) const;
321324

322325
private:
323326
int mpar[50]; //!< Matrix of parameters

0 commit comments

Comments
 (0)