From 2fa77a2e83263a9df9ed3210125596f7ac7d986c Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 25 Nov 2025 16:42:30 +0100 Subject: [PATCH 1/5] Added: Command-line option -validateGroups --- src/SIM/SIMoptions.C | 11 ++++++----- src/SIM/SIMoptions.h | 3 ++- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/SIM/SIMoptions.C b/src/SIM/SIMoptions.C index 9fa949751..f32ab7775 100644 --- a/src/SIM/SIMoptions.C +++ b/src/SIM/SIMoptions.C @@ -41,6 +41,7 @@ SIMoptions::SIMoptions () #else num_threads_SLU = 1; #endif + validateGroups = false; eig = 0; nev = 10; @@ -129,17 +130,15 @@ bool SIMoptions::parseDiscretizationTag (const tinyxml2::XMLElement* elem) } else if (!strcasecmp(elem->Value(),"geometry")) { - std::string type; const tinyxml2::XMLElement* child = elem->FirstChildElement(); for (; child; child = child->NextSiblingElement()) if (!strcasecmp(child->Value(),"patchfile")) - if (utl::getAttribute(child,"type",type) && type == "lrspline") - discretization = ASM::LRSpline; + if (std::string type; utl::getAttribute(child,"type",type)) + if (type == "lrspline") discretization = ASM::LRSpline; } else if (!strcasecmp(elem->Value(),"nGauss")) { - int defaultG = 0; - if (utl::getAttribute(elem,"default",defaultG)) + if (int defaultG = 0; utl::getAttribute(elem,"default",defaultG)) nGauss[0] = nGauss[1] = defaultG > 0 ? 10+defaultG : defaultG; else if (elem->FirstChild()) { std::string value(elem->FirstChild()->Value()); @@ -316,6 +315,8 @@ bool SIMoptions::parseOldOptions (int argc, char** argv, int& i) solver = LinAlg::PETSC; else if (!strcmp(argv[i],"-istl")) solver = LinAlg::ISTL; + else if (!strcmp(argv[i],"-validateGroups")) + validateGroups = true; else if (!strncmp(argv[i],"-lag",4)) discretization = ASM::Lagrange; else if (!strncmp(argv[i],"-tri",4)) diff --git a/src/SIM/SIMoptions.h b/src/SIM/SIMoptions.h index ea9fa9e1c..cdec999d6 100644 --- a/src/SIM/SIMoptions.h +++ b/src/SIM/SIMoptions.h @@ -73,6 +73,7 @@ class SIMoptions ASM::Discretization discretization; //!< Spatial discretization option LinAlg::MatrixType solver; //!< The linear equation solver to use + bool validateGroups; //!< Flag validation of the generated threading groups int num_threads_SLU; //!< Number of threads for SuperLU_MT // Eigenvalue solver options @@ -106,7 +107,7 @@ class SIMoptions enum ProjectionMethod { NONE, GLOBAL, VDSA, QUASI, LEASTSQ, DGL2, CGL2, CGL2_INT, SCR }; //! \brief Projection method name mapping. - typedef std::map ProjectionMap; + using ProjectionMap = std::map; ProjectionMap project; //!< The projection methods to use }; From 38b441cea42392fb8e441342c795afe1109f6a0f Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 25 Nov 2025 17:00:40 +0100 Subject: [PATCH 2/5] Added: Method DiagMatrix::assembleStruct() for validation of groups --- src/LinAlg/DiagMatrix.C | 69 ++++++++++++++++++++++++++++++++++------- src/LinAlg/DiagMatrix.h | 8 +++++ 2 files changed, 65 insertions(+), 12 deletions(-) diff --git a/src/LinAlg/DiagMatrix.C b/src/LinAlg/DiagMatrix.C index b603fcd93..fc6f93557 100644 --- a/src/LinAlg/DiagMatrix.C +++ b/src/LinAlg/DiagMatrix.C @@ -62,7 +62,7 @@ void DiagMatrix::dump (std::ostream& os, LinAlg::StorageFormat format, bool DiagMatrix::assemble (const Matrix& eM, const SAM& sam, int e) { - if (myMat.size() != (size_t)sam.neq) + if (static_cast(myMat.size()) != sam.neq) return false; StdVector B; @@ -74,22 +74,67 @@ bool DiagMatrix::assemble (const Matrix& eM, const SAM& sam, int e) bool DiagMatrix::assemble (const Matrix& eM, const SAM& sam, SystemVector&, const IntVec& meq) { - if (meq.size() != 1 || eM.size() != 1) + const size_t nedof = meq.size(); + if (eM.rows() < nedof || eM.cols() < nedof) { - std::cerr <<" *** DiagMatrix::assemble: Only for one-dof elements, nedof = " - << meq.size() <<" size(eM) = " << eM.size() << std::endl; + std::cerr <<" *** DiagMatrix::assemble: Invalid element matrix, nedof = " + << nedof <<" size(eM) = "<< eM.rows() <<","<< eM.cols() + << std::endl; return false; } - int ieq = meq.front(); - if (ieq < 1 || ieq > sam.neq) - { - std::cerr <<" *** DiagMatrix::assemble: ieq="<< ieq <<" is out or range [1," - << sam.neq <<"]"<< std::endl; - return false; - } + // Add elements corresponding to free dofs in eM, + // and (appropriately weighted) elements corresponding to constrained dofs, + // into the diagonal system matrix + for (size_t j = 1; j <= nedof; j++) + if (int jeq = meq[j-1]; jeq > 0) + myMat(jeq) += eM(j,j); + else if (int jceq = -meq[j-1]; jceq > 0) + for (int jp = sam.mpmceq[jceq-1]; jp < sam.mpmceq[jceq]-1; jp++) + if (int jeq = sam.mmceq[jp] > 0 ? sam.meqn[sam.mmceq[jp]-1]:0; jeq > 0) + for (size_t i = 1; i <= nedof; i++) + if (meq[i-1] == jeq) + myMat(jeq) += sam.ttcc[jp]*(eM(i,j) + eM(j,i)); + else if (int iceq = -meq[i-1]; iceq > 0) + for (int ip = sam.mpmceq[iceq-1]; ip < sam.mpmceq[iceq]-1; ip++) + if (sam.mmceq[ip] > 0 && sam.meqn[sam.mmceq[ip]-1] == jeq) + myMat(jeq) += sam.ttcc[ip]*sam.ttcc[jp]*eM(i,j); + + return this->flagNonZeroEqs(meq); +} + + +/*! + This method is doing the same as the above assemble() method, + but where all element matrices are assumed identity matrices. + The purpose is to verify that each equation only gets a single contribution + within each threading group. +*/ - myMat(ieq) += eM(1,1); +bool DiagMatrix::assembleStruct (int val, const SAM& sam, const IntVec& meq) +{ + // Lambda function marking matrix contributions due to linear couplings. + auto&& addEq = [&diag=myMat,val](int ieq) + { + int current = diag(ieq)/1000.0; + if (current != val) + diag(ieq) += 1000.0*val; + }; + + const size_t nedof = meq.size(); + for (size_t j = 1; j <= nedof; j++) + if (int jeq = meq[j-1]; jeq > 0) + myMat(jeq) += 1.0; + else if (int jceq = -meq[j-1]; jceq > 0) + for (int jp = sam.mpmceq[jceq-1]; jp < sam.mpmceq[jceq]-1; jp++) + if (int jeq = sam.mmceq[jp] > 0 ? sam.meqn[sam.mmceq[jp]-1]:0; jeq > 0) + for (size_t i = 1; i <= nedof; i++) + if (meq[i-1] == jeq) + addEq(jeq); + else if (int iceq = -meq[i-1]; iceq > 0) + for (int ip = sam.mpmceq[iceq-1]; ip < sam.mpmceq[iceq]-1; ip++) + if (sam.mmceq[ip] > 0 && sam.meqn[sam.mmceq[ip]-1] == jeq) + addEq(jeq); return this->flagNonZeroEqs(meq); } diff --git a/src/LinAlg/DiagMatrix.h b/src/LinAlg/DiagMatrix.h index 4b65a92ac..7b230cdc0 100644 --- a/src/LinAlg/DiagMatrix.h +++ b/src/LinAlg/DiagMatrix.h @@ -101,6 +101,14 @@ class DiagMatrix : public SystemMatrix //! \details To be used when there is no underlying SAM virtual bool assemble(const Matrix& eM, const IntVec& meq); + //! \brief Assembles a validation matrix assuming unit element contributions. + //! \param[in] val Value identifying current threading group + //! \param[in] sam Auxiliary data for FE assembly management + //! \param[in] meq Matrix of element equation numbers + //! + //! \details For validation of threading groups. + bool assembleStruct(int val, const SAM& sam, const IntVec& meq); + //! \brief Multiplication with a scalar. virtual void mult(Real alpha) { myMat *= alpha; } From aed8d7aedede0e33a8220692d54e6f884bf5936f Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 25 Nov 2025 17:05:49 +0100 Subject: [PATCH 3/5] Added: Optional validation of threading groups based on SAM data --- src/ASM/ASMbase.h | 4 ++++ src/ASM/ASMs2D.C | 44 ++++++++++++++++++++++++++++++++++++++++++++ src/ASM/ASMs2D.h | 3 +++ src/SIM/SIMbase.C | 32 ++++++++++++++++++++------------ 4 files changed, 71 insertions(+), 12 deletions(-) diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 604fdd234..46e454181 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -45,6 +45,7 @@ class RealFunc; class VecFunc; class Vec3; class Tensor; +class SAM; namespace ASM { class InterfaceChecker; } namespace utl { template class Function; } @@ -515,10 +516,13 @@ class ASMbase virtual void generateThreadGroupsFromElms(const IntVec&) {} //! \brief Generates element groups for multi-threading based on a partition. virtual void generateProjThreadGroupsFromElms(const IntVec&) {} + //! \brief Validates the threading groups based on the assembly data in %SAM. + virtual bool validateThreadGroups(const SAM*) const { return true; } //! \brief Hook for changing number of threads. virtual void changeNumThreads() {} + // Methods for integration of finite element quantities. // These are the main computational methods of the ASM class hierarchy. // ==================================================================== diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index f9fbff213..500dbd7b4 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -36,6 +36,8 @@ #include "Point.h" #include "Tensor.h" #include "MPC.h" +#include "SAM.h" +#include "DiagMatrix.h" #include "IFEM.h" #include #include @@ -3356,6 +3358,48 @@ void ASMs2D::generateProjThreadGroupsFromElms (const IntVec& elms) } +bool ASMs2D::validateThreadGroups (const SAM* sam) const +{ + IFEM::cout <<"\nValidating element groups for multi-threaded assembly." + << std::endl; + + if (threadGroups[0].size() == 1) + return true; // Only one group (no multi-threading) + + size_t nErr = 0; + int iTGroup = 0; + for (size_t g = 0; g < threadGroups.size(); g++) + for (const IntVec& group : threadGroups[g]) + { + ++iTGroup; + size_t lErr = nErr; + DiagMatrix sysMat(sam->getNoEquations()); + for (int iel : group) + if (IntVec meen; !sam->getElmEqns(meen,MLGE[iel])) + ++nErr; + else if (!sysMat.assembleStruct(iTGroup,*sam,meen)) + ++nErr; + + for (size_t ieq = 1; ieq <= sysMat.dim(0); ieq++) + if (int count = sysMat(ieq); count > 1 && count/1000 != iTGroup) + { + std::cerr <<" *** Threading group "<< iTGroup <<" has "<< count + <<" contributors to equation "<< ieq; + std::pair dof = sam->getNodeAndLocalDof(ieq,true); + std::cerr <<" (node "<< dof.first <<" local dof "<< dof.second <<")" + << std::endl; + nErr += count; + } + + if (lErr == nErr) + IFEM::cout <<" * Thread group "<< iTGroup <<" (size "<< group.size() + <<") is OK"<< std::endl; + } + + return nErr == 0; +} + + bool ASMs2D::addRigidCpl (int lindx, int ldim, int basis, int& gMaster, const Vec3& Xmaster, bool extraPt) { diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h index cff90c254..47350f5a8 100644 --- a/src/ASM/ASMs2D.h +++ b/src/ASM/ASMs2D.h @@ -698,6 +698,9 @@ class ASMs2D : public ASMstruct, public ASM2D //! \brief Generate element groups from a partition. virtual void generateProjThreadGroupsFromElms(const IntVec& elms); + //! \brief Validates the threading groups based on the assembly data in %SAM. + virtual bool validateThreadGroups(const SAM* sam) const; + //! \brief Returns 0-based index of first node on integration basis. virtual int getFirstItgElmNode() const { return 0; } //! \brief Returns 0-based index of last node on integration basis. diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index cc93d6cc5..5899a984c 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -362,18 +362,6 @@ bool SIMbase::preprocessC (const IntVec& ignored, bool fixDup, double time0) if (timeDependent && !this->initDirichlet(time0)) return false; - // Generate element groups for multi-threading - bool silence = msgLevel < 1 || (msgLevel < 3 && nGlPatches > 1); - for (ASMbase* pch : myModel) - if (!pch->empty()) - pch->generateThreadGroups(*myProblem,silence,lagMTOK); - - for (const Property& p : myProps) - if (p.pcode == Property::NEUMANN || - p.pcode == Property::NEUMANN_GENERIC || - p.pcode == Property::ROBIN) - this->generateThreadGroups(p,silence); - // Preprocess the result points this->preprocessResultPoints(); @@ -402,6 +390,26 @@ bool SIMbase::preprocessC (const IntVec& ignored, bool fixDup, double time0) return false; } + // Generate element groups for multi-threading + bool silence = msgLevel < 1 || (msgLevel < 3 && nGlPatches > 1); + if (opt.validateGroups) silence = false; + for (ASMbase* pch : myModel) + if (!pch->empty()) + pch->generateThreadGroups(*myProblem,silence,lagMTOK); + + for (const Property& p : myProps) + if (p.pcode == Property::NEUMANN || + p.pcode == Property::NEUMANN_GENERIC || + p.pcode == Property::ROBIN) + if (ASMbase* pch = this->getPatch(p.patch); pch && !pch->empty()) + if (abs(p.ldim)+1 == pch->getNoParamDim()) + pch->generateThreadGroups(p.lindx%10,silence,lagMTOK); + + if (opt.validateGroups) + for (ASMbase* pch : myModel) + if (!pch->empty() && !pch->validateThreadGroups(mySam)) + return false; + if (mdFlag > 0) nDofS = mySam->getNoDOFs(); else if (!adm.dd.setup(adm,*this)) From 387858baa4f39e2c72eb12366b582eed7ce85dc4 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 26 Nov 2025 06:38:28 +0100 Subject: [PATCH 4/5] 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. --- src/LinAlg/SAM.C | 145 ++++++++++++++++++----------------------- src/LinAlg/SAM.h | 5 +- src/LinAlg/SPRMatrix.C | 38 +++++++++-- src/LinAlg/SPRMatrix.h | 7 +- 4 files changed, 105 insertions(+), 90 deletions(-) diff --git a/src/LinAlg/SAM.C b/src/LinAlg/SAM.C index d0430f906..53f2cac87 100644 --- a/src/LinAlg/SAM.C +++ b/src/LinAlg/SAM.C @@ -98,13 +98,21 @@ SAM::~SAM () } -std::pair SAM::getNodeAndLocalDof (int idof) const +std::pair SAM::getNodeAndLocalDof (int idof, bool eqno) const { + if (eqno) + { + if (int* eq = std::find(meqn,meqn+ndof,idof); eq != meqn+ndof) + idof = std::distance(meqn,eq); + else + return { 0, 0 }; + } + for (int n = 1; n <= nnod; n++) if (madof[n] > idof) - return std::make_pair(minex ? minex[n-1] : n, idof-madof[n-1]+1); + return { minex ? minex[n-1] : n, idof-madof[n-1]+1 }; - return std::make_pair(0,0); + return { 0, 0 }; } @@ -138,7 +146,6 @@ void SAM::print (std::ostream& os) const os <<"\n\nSAM::mpar: "<< mpar[0]; for (int i = 1; i < 30; i++) os << ((i%10) ? " " : "\n ") << mpar[i]; - os << std::endl; if (mmnpc && mpmnpc && nel > 0) { @@ -153,7 +160,6 @@ void SAM::print (std::ostream& os) const else if (mmnpc[i-1] < 0) os <<' '<< (minex ? -minex[-mmnpc[i-1]-1] : mmnpc[i-1]); } - os << std::endl; } if (ttcc && mmceq && mpmceq && nceq > 0) @@ -164,13 +170,19 @@ void SAM::print (std::ostream& os) const os <<'\n'<< std::setw(4) << i+1 <<": "; this->printCEQ(os,i); } - os << std::endl; } + os << std::endl; + + this->printMEQN(os); +} + +void SAM::printMEQN (std::ostream& os) const +{ if (meqn && madof && nnod > 0 && neq > 0) { char code[7], dofType[7]; - os <<"\n\nNode --> Equations"; + os <<"\nNode --> Equations"; for (int n = 0; n < nnod; n++) { int dof, i = 0; @@ -195,22 +207,25 @@ void SAM::print (std::ostream& os) const void SAM::printStatusCodes (std::ostream& os) const { - os <<"\nNode\tDOFs\t MSC"; - for (int i = 0; i < nnod; i++) + if (msc && madof && nnod > 0) { - os <<"\n"<< std::setw(4) << 1+i - <<"\t"<< madof[i] <<" - "<< madof[i+1]-1 <<"\t"; - for (int j = madof[i]; j < madof[i+1]; j++) - os <<" "<< msc[j-1]; + os <<"\nNode\tDOFs\t MSC"; + for (int i = 0; i < nnod; i++) + { + os <<"\n"<< std::setw(4) << 1+i + <<"\t"<< madof[i] <<" - "<< madof[i+1]-1 <<"\t"; + for (int j = madof[i]; j < madof[i+1]; j++) + os <<" "<< msc[j-1]; + } + os << std::endl; } - os << std::endl; } bool SAM::initSystemEquations () { #ifdef SP_DEBUG - std::cout <<"SAM::initSystemEquations()"<< std::endl; + std::cout <<"\nSAM::initSystemEquations()"<< std::endl; #endif if (!msc && ndof > 0) return false; if ((!msc || !mpmceq) && nceq > 0) return false; @@ -307,6 +322,9 @@ bool SAM::initSystemEquations () else if (msc[idof] == 2) meqn[idof] = j++; #endif +#if SP_DEBUG > 1 + this->printMEQN(std::cout); +#endif if (ierr == 0) return true; @@ -365,20 +383,15 @@ void SAM::getDofCouplings (std::vector& dofc) const this->getElmEqns(meen,iel); for (size_t j = 0; j < meen.size(); j++) - { - int jeq = meen[j]; - if (jeq > 0) + if (int jeq = meen[j]; jeq > 0) { dofc[jeq-1].insert(jeq); for (size_t i = 0; i < j; i++) - { - int ieq = meen[i]; - if (ieq > 0) + if (int ieq = meen[i]; ieq > 0) { dofc[ieq-1].insert(jeq); dofc[jeq-1].insert(ieq); } - } } else if (jeq < 0) { @@ -410,7 +423,6 @@ void SAM::getDofCouplings (std::vector& dofc) const } } } - } } } @@ -514,7 +526,7 @@ bool SAM::assembleSystem (SystemVector& sysRHS, #ifdef USE_F77SAM int* work = new int[eS.size()]; addev2_(&eS.front(), ttcc, mpar, madof, meqn, mpmnpc, mmnpc, mpmceq, mmceq, - iel, eS.size(), 6, sysRHS.getPtr(), work, ierr); + iel, eS.size(), 6, sysRHS.getPtr(), work, ierr); delete[] work; #else IntVec meen; @@ -543,13 +555,10 @@ bool SAM::assembleSystem (SystemVector& sysRHS, const Real* nS, int inod, if (reactionForces) { - int k = 0; + int k = 0, nReact = reactionForces->size(); for (int j = madof[inod-1]; j < madof[inod]; j++, k++) - { - int ipR = -msc[j-1]; - if (ipR > 0 && ipR <= static_cast(reactionForces->size())) + if (int ipR = -msc[j-1]; ipR > 0 && ipR <= nReact) (*reactionForces)[ipR-1] += nS[k]; - } } return true; @@ -614,18 +623,12 @@ void SAM::assembleReactions (RealArray& rf, const RealArray& eS, int iel) const int ip = mpmnpc[iel-1]; int nenod = mpmnpc[iel] - ip; for (int i = 0; i < nenod; i++, ip++) - { - int node = mmnpc[ip-1]; - if (node < 0) + if (int node = mmnpc[ip-1]; node < 0) k += madof[-node] - madof[-node-1]; else if (node > 0) for (int j = madof[node-1]; j < madof[node] && k < eS.size(); j++, k++) - { - int ipR = -msc[j-1]; - if (ipR > 0 && ipR <= static_cast(rf.size())) + if (int ipR = -msc[j-1]; ipR > 0 && ipR <= static_cast(rf.size())) rf[ipR-1] += eS[k]; - } - } } @@ -678,13 +681,10 @@ bool SAM::getElmEqns (IntVec& meen, int iel, size_t nedof) const #else meen.reserve(nedof > 0 ? nedof : neldof); for (int i = 0; i < nenod; i++, ip++) - { - int node = mmnpc[ip]; - if (node > 0) + if (int node = mmnpc[ip]; node > 0) meen.insert(meen.end(),meqn+madof[node-1]-1,meqn+madof[node]-1); else if (node < 0) meen.insert(meen.end(),madof[-node]-madof[-node-1],0); - } neldof = meen.size(); #endif if (nedof == 0 || neldof == static_cast(nedof)) return true; @@ -789,27 +789,20 @@ bool SAM::expandVector (const Real* solVec, Vector& dofVec, Real scaleSD) const dofVec.resize(ndof,true); #ifdef USE_F77SAM expand_(solVec, ttcc, mpmceq, mmceq, meqn, Real(1), scaleSD, ndof, neq, - dofVec.ptr()); + dofVec.ptr()); #else for (int idof = 0; idof < ndof; idof++) - { - int ieq = meqn[idof]; - int iceq = -ieq; - if (ieq > 0) + if (int ieq = meqn[idof]; ieq > 0) dofVec[idof] += solVec[ieq-1]; - else if (iceq > 0) + else if (int iceq = -meqn[idof]; iceq > 0) { int ip = mpmceq[iceq-1]; dofVec[idof] += scaleSD*ttcc[ip-1]; for (; ip < mpmceq[iceq]-1; ip++) - if (mmceq[ip] > 0) - { - ieq = meqn[mmceq[ip]-1]; - if (ieq > 0 && ieq <= neq) - dofVec[idof] += ttcc[ip]*solVec[ieq-1]; - } + if (mmceq[ip] > 0) + if (int ieq = meqn[mmceq[ip]-1]; ieq > 0 && ieq <= neq) + dofVec[idof] += ttcc[ip]*solVec[ieq-1]; } - } #endif return true; @@ -822,11 +815,8 @@ bool SAM::getSolVec (RealArray& solVec, const RealArray& dofVec) const solVec.resize(neq,Real(0)); for (int idof = 0; idof < ndof; idof++) - { - int ieq = meqn[idof]; - if (ieq > 0) + if (int ieq = meqn[idof]; ieq > 0) solVec[ieq-1] = dofVec[idof]; - } return true; } @@ -837,16 +827,13 @@ bool SAM::applyDirichlet (Vector& dofVec) const if (!meqn) return false; for (int idof = 0; idof < ndof; idof++) - { - int iceq = -meqn[idof]; - if (iceq > 0) + if (int iceq = -meqn[idof]; iceq > 0) { int ip = mpmceq[iceq-1]; dofVec[idof] = ttcc[ip-1]; } else if (iceq == 0) dofVec[idof] = Real(0); - } return true; } @@ -900,8 +887,7 @@ Real SAM::normInf (const Vector& x, size_t& comp, char dofType) const else if (nodeType.empty() && dof_type.empty()) { // All DOFs are of the same type and the number of nodal DOFs is constant - int nndof = madof[1] - madof[0]; - if (static_cast(comp) <= nndof) + if (int nndof = madof[1] - madof[0]; static_cast(comp) <= nndof) return x.normInf(--comp,nndof); else return Real(0); @@ -923,15 +909,12 @@ Real SAM::normInf (const Vector& x, size_t& comp, char dofType) const } } else if (nodeType[i] == dofType) - { - int idof = madof[i] + icmp - 2; - if (idof >= 0 && idof < madof[i+1]-1) + if (int idof = madof[i] + icmp - 2; idof >= 0 && idof < madof[i+1]-1) if (fabs(x[idof]) > retVal) { comp = i+1; retVal = fabs(x[idof]); } - } return retVal; } @@ -962,11 +945,9 @@ bool SAM::haveReaction (int dir, const IntVec* nodes) const if (dir > 0 && nspdof > 0) for (int i = 0; i < nnod; i++) if (!nodes || std::find(nodes->begin(),nodes->end(),i+1) != nodes->end()) - { - int idof = madof[i]+dir-2; - if (idof < madof[i+1]-1 && msc[idof] < 0 && -msc[idof] <= nspdof) - return true; - } + if (int idof = madof[i]+dir-2; idof < madof[i+1]-1) + if (msc[idof] < 0 && -msc[idof] <= nspdof) + return true; return false; } @@ -979,16 +960,14 @@ Real SAM::getReaction (int dir, const RealArray& rf, const IntVec* nodes) const if (dir > 0 && nspdof > 0) for (int i = 0; i < nnod; i++) if (!nodes || std::find(nodes->begin(),nodes->end(),i+1) != nodes->end()) - { - int idof = madof[i]+dir-2; - int ipr = idof < madof[i+1]-1 ? -msc[idof]-1 : -1; - if (ipr >= 0 && ipr < static_cast(rf.size())) - { - retVal += rf[ipr]; - if (nodes) // clear this force component to avoid it being added twice - const_cast(rf)[ipr] = 0.0; - } - } + if (int idof = madof[i]+dir-2; idof < madof[i+1]-1) + if (int ipr = -msc[idof]-1; ipr >= 0) + if (ipr < static_cast(rf.size())) + { + retVal += rf[ipr]; + if (nodes) // clear force component to avoid it being added twice + const_cast(rf)[ipr] = 0.0; + } return retVal; } diff --git a/src/LinAlg/SAM.h b/src/LinAlg/SAM.h index e5651bf5f..32b64b5a8 100644 --- a/src/LinAlg/SAM.h +++ b/src/LinAlg/SAM.h @@ -191,9 +191,10 @@ class SAM //! \brief Returns the internal node number and local index for a global DOF. //! \param[in] idof Global DOF-number in the range [1,NDOF] + //! \param[in] eqno If \e true, \a idof is the 1-based global equation number //! \return first = internal node number in the range [1,NNOD] //! \return second = local DOF number in the range [1,NNDOF] - std::pair getNodeAndLocalDof(int idof) const; + std::pair getNodeAndLocalDof(int idof, bool eqno = false) const; //! \brief Finds the equation number corresponding to a local nodal DOF. //! \param[in] inod Identifier for the node to get the equation number for @@ -318,6 +319,8 @@ class SAM void printStatusCodes(std::ostream& os) const; //! \brief Prints out a constraint equation to the given stream. void printCEQ(std::ostream& os, int iceq) const; + //! \brief Prints out the equation number mapping to the given stream. + void printMEQN(std::ostream& os) const; private: int mpar[50]; //!< Matrix of parameters diff --git a/src/LinAlg/SPRMatrix.C b/src/LinAlg/SPRMatrix.C index 037049e47..d0193413c 100644 --- a/src/LinAlg/SPRMatrix.C +++ b/src/LinAlg/SPRMatrix.C @@ -239,6 +239,7 @@ SPRMatrix::SAM64::SAM64 (const SAM& sam) mpmceq = copyTo64(sam.mpmceq,sam.nceq+1); mmceq = copyTo64(sam.mmceq,sam.nmmceq); meqn = copyTo64(sam.meqn,sam.ndof); + minex = sam.minex; } @@ -252,6 +253,23 @@ SPRMatrix::SAM64::~SAM64 () delete[] mmceq; delete[] meqn; } + + +std::pair SPRMatrix::SAM64::getNodeAndLocalDof (int ieq, bool) const +{ + const Int_ nnod = mpar[0]; + const Int_ ndof = mpar[2]; + if (Int_* idof = std::find(meqn,meqn+ndof,ieq); idof != meqn+ndof) + ieq = std::distance(meqn,idof); + else + return { 0, 0 }; + + for (int n = 1; n <= nnod; n++) + if (madof[n] > static_cast(ieq)) + return { minex ? minex[n-1] : n, ieq-madof[n-1]+1 }; + + return { 0, 0 }; +} #endif @@ -677,11 +695,23 @@ bool SPRMatrix::solve (SystemVector& B, Real*) Int_ iop = mpar[0] < 5 ? 3 : 4; sprsol_(iop, mpar, mtrees, msifa, values, B.getPtr(), B.dim(), 1, tol, iWork.data(), rWork.data(), 6, ierr); - if (!ierr) return true; - - std::cerr <<"SAM::SPRSOL: Failure "<< ierr << std::endl; + if (ierr < 0) + std::cerr <<"SAM::SPRSOL: Failure ("<< ierr <<")"; + else if (ierr == 3) + std::cerr <<"SAM::SPRSOL: Negative pivot"; + else + return true; + + std::pair dof = mySam->getNodeAndLocalDof(mpar[26],true); + if (dof.first > 0) + std::cerr <<" detected in node "<< dof.first <<", local DOF "<< dof.second; + else if (mpar[26] > 0) + std::cerr <<", for equation "<< mpar[26] <<" (unknown) "; + std::cerr << std::endl; +#else + ierr = -99; #endif - return false; + return ierr >= 0; } diff --git a/src/LinAlg/SPRMatrix.h b/src/LinAlg/SPRMatrix.h index d2aa594fd..a224b6361 100644 --- a/src/LinAlg/SPRMatrix.h +++ b/src/LinAlg/SPRMatrix.h @@ -79,7 +79,7 @@ class SPRMatrix : public SystemMatrix //! \details When multi-point constraints are present, contributions from //! these are also added into the system right-hand-side vector, \a B. virtual bool assemble(const Matrix& eM, const SAM& sam, - SystemVector& B, int e); + SystemVector& B, int e); //! \brief Adds an element matrix into the associated system matrix. //! \param[in] eM The element matrix //! \param[in] sam Auxiliary data for FE assembly management @@ -127,7 +127,7 @@ class SPRMatrix : public SystemMatrix //! \param[in] shift Eigenvalue shift //! \param[in] iop Option telling whether to factorize matrix \a A or \b B. bool solveEig(SPRMatrix& B, RealArray& val, Matrix& vec, int nev, - Real shift = 0.0, int iop = 1); + Real shift = 0.0, int iop = 1); //! \brief Returns the L-infinity norm of the matrix. virtual Real Linfnorm() const; @@ -178,11 +178,14 @@ class SPRMatrix : public SystemMatrix Int_* mpmceq; //!< Matrix of pointers to MCEQs in MMCEQ Int_* mmceq; //!< Matrix of matrices of constraint equation definitions Int_* meqn; //!< Matrix of equation numbers + int* minex; //!< Matrix of internal to external node number mapping //! \brief The constructor copies the 32-bit arrays into 64-bit versions. SAM64(const SAM& sam); //! \brief The destructor frees the dynamically allocated arrays. ~SAM64(); + //! \brief Returns the node and local DOF number for the equation \a ieq. + std::pair getNodeAndLocalDof(int ieq, bool) const; }; bool sharedSam = false; //!< True if SAM object is shared with another matrix From 572a0e4a590b98bec364685776e7ab6db91b0887 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 3 Dec 2025 15:09:22 +0100 Subject: [PATCH 5/5] Fixed: Properly account for the MPC-couplings when coloring the elements --- src/ASM/ASMu2DLag.C | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/src/ASM/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index 81370a5ec..bb45f551a 100644 --- a/src/ASM/ASMu2DLag.C +++ b/src/ASM/ASMu2DLag.C @@ -463,6 +463,7 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence, IntVec status(nel,0); using IntSet = std::set; + IntSet slaveNodes; std::vector nodeConn(nnod); // node-to-element connectivity std::vector nodeNode(nnod); // node-to-node connectivity via MPCs @@ -471,16 +472,17 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence, for (size_t i = 0; i < mpc->getNoMaster(); i++) if (int mastr = this->getNodeIndex(mpc->getMaster(i).node); mastr > 0) { - nodeNode[slave-1].insert(mastr-1); + slaveNodes.insert(slave-1); nodeNode[mastr-1].insert(slave-1); } size_t fixedElements = 0; for (size_t iel = 0; iel < nel; iel++) if (separateGroup1noded && MNPC[iel].size() == 1 && - nodeNode[MNPC[iel].front()].empty()) + slaveNodes.find(MNPC[iel].front()) == slaveNodes.end()) { - status[iel] = 1; // Separate group for all single-noded elements + // Separate color for all single-noded elements + status[iel] = 1; // whose element node is not a slave if (++fixedElements == 1) threadGroups[0].resize(1, { static_cast(iel) }); else @@ -488,11 +490,30 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence, } else for (int node : MNPC[iel]) - { nodeConn[node].insert(iel); - for (int master : nodeNode[node]) - nodeConn[master].insert(iel); - } + + // Second pass - account for the MPC couplings + for (size_t master = 0; master < nodeNode.size(); master++) + if (!nodeNode[master].empty()) + { + IntSet commonElms(nodeConn[master]); + for (int slave : nodeNode[master]) + commonElms.insert(nodeConn[slave].begin(),nodeConn[slave].end()); + nodeConn[master].insert(commonElms.begin(),commonElms.end()); + for (int slave : nodeNode[master]) + nodeConn[slave] = commonElms; + } + +#if SP_DEBUG > 1 + std::cout <<"\nNode to element connectivity (incl. MPC couplings):"; + for (size_t inod = 0; inod < nnod; inod++) + { + std::cout <<"\n"<< inod <<" ->"; + for (int e : nodeConn[inod]) + std::cout <<" "<< e; + } + std::cout << std::endl; +#endif for (size_t nColors = fixedElements > 0; fixedElements < nel; ++nColors) { @@ -501,7 +522,7 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence, [](int& s) { if (s < 0) s = 0; }); // Look for available elements - IntVec thisColor; + IntVec& thisColor = threadGroups[0].emplace_back(); for (size_t i = 0; i < nel; ++i) if (status[i] == 0) { @@ -514,8 +535,6 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence, if (status[j] == 0) // if not assigned a color yet status[j] = -1; // set as unavailable (with current color) } - - threadGroups[0].emplace_back(thisColor); } if (!silence)