diff --git a/src/SIM/SIM1D.C b/src/SIM/SIM1D.C index 87e66f6a..d5108399 100644 --- a/src/SIM/SIM1D.C +++ b/src/SIM/SIM1D.C @@ -17,7 +17,6 @@ #include "ASMs1DSpec.h" #include "Functions.h" #include "Utilities.h" - #include "tinyxml.h" @@ -411,6 +410,8 @@ void SIM1D::setQuadratureRule (size_t ng) for (size_t i = 0; i < myModel.size(); i++) if (!myModel.empty()) static_cast(myModel[i])->setGauss(ng); + + this->initIntegrationBuffers(); } diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C index cdf6814f..421d9ce9 100644 --- a/src/SIM/SIM2D.C +++ b/src/SIM/SIM2D.C @@ -15,9 +15,7 @@ #include "ASMs2DC1.h" #include "Functions.h" #include "Utilities.h" - #include "tinyxml.h" - #ifdef USE_OPENMP #include #endif @@ -203,7 +201,7 @@ bool SIM2D::parseGeometryTag(const TiXmlElement* elem) std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch << std::endl; static_cast(myModel[patch-1])->closeEdges(pedir); - // cannot do multi-threaded assembly with periodicities + // Cannot do multi-threaded assembly with periodicities #ifdef USE_OPENMP omp_set_num_threads(1); #endif @@ -390,11 +388,11 @@ bool SIM2D::parse (char* keyWord, std::istream& is) std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch << std::endl; static_cast(myModel[patch-1])->closeEdges(pedir); - // cannot do multi-threaded assembly with periodicities + } + // Cannot do multi-threaded assembly with periodicities #ifdef USE_OPENMP omp_set_num_threads(1); #endif - } } else if (!strncasecmp(keyWord,"CONSTRAINTS",11)) @@ -552,6 +550,8 @@ void SIM2D::setQuadratureRule (size_t ng) for (size_t i = 0; i < myModel.size(); i++) if (!myModel.empty()) static_cast(myModel[i])->setGauss(ng); + + this->initIntegrationBuffers(); } diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C index b1bd87c3..c74973d6 100644 --- a/src/SIM/SIM3D.C +++ b/src/SIM/SIM3D.C @@ -17,9 +17,11 @@ #include "ASMs3DSpec.h" #include "Functions.h" #include "Utilities.h" -#include - #include "tinyxml.h" +#include +#ifdef USE_OPENMP +#include +#endif SIM3D::SIM3D (bool checkRHS, unsigned char n1, unsigned char n2) @@ -211,11 +213,12 @@ bool SIM3D::parse (char* keyWord, std::istream& is) std::cout <<"\tPeriodic "<< char('H'+pfdir) <<"-direction P"<< patch << std::endl; static_cast(myModel[patch-1])->closeFaces(pfdir); - // cannot do multi-threaded assembly with periodicities + } + + // Cannot do multi-threaded assembly with periodicities #ifdef USE_OPENMP omp_set_num_threads(1); #endif - } } else if (!strncasecmp(keyWord,"CONSTRAINTS",11)) @@ -565,6 +568,8 @@ void SIM3D::setQuadratureRule (size_t ng) for (size_t i = 0; i < myModel.size(); i++) if (!myModel.empty()) static_cast(myModel[i])->setGauss(ng); + + this->initIntegrationBuffers(); } diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 7678d6f8..2e740b9c 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -58,6 +58,7 @@ SIMbase::SIMbase () : g2l(&myGlb2Loc) vizIncr = 1; format = 1; mixedFEM = false; + nIntGP = nBouGP = 0; MPCLess::compareSlaveDofOnly = true; // to avoid multiple slave definitions @@ -658,13 +659,9 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) if (renum > 0) std::cout <<"\nRenumbered "<< renum <<" nodes."<< std::endl; - size_t nIntP = 0, nBouP = 0; ASMs2DC1::renumberNodes(*g2l); for (mit = myModel.begin(); mit != myModel.end(); mit++) - { (*mit)->renumberNodes(*g2l); - (*mit)->getNoIntPoints(nIntP); // count number of global integration points - } // Process the specified Dirichlet boundary conditions std::cout <<"\nResolving Dirichlet boundary conditions"<< std::endl; @@ -687,12 +684,6 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) ok = false; break; - case Property::NEUMANN: - // Count number of global boundary integration points - if (p->patch > 0 && p->patch <= myModel.size()) - myModel[p->patch-1]->getNoBouPoints(nBouP,p->ldim,p->lindx); - break; - default: break; } @@ -709,9 +700,6 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) // Resolve possibly chaining of the MPC equations if (!allMPCs.empty()) ASMbase::resolveMPCchains(allMPCs); - // Let the integrand know how many integration points in total do we have - if (myProblem) myProblem->initIntegration(nIntP,nBouP); - // Preprocess the result points for (ResPointVec::iterator p = myPoints.begin(); p != myPoints.end();) { @@ -753,6 +741,23 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) } +void SIMbase::initIntegrationBuffers () +{ + nIntGP = nBouGP = 0; + for (size_t i = 0; i < myModel.size(); i++) + myModel[i]->getNoIntPoints(nIntGP); // count the interior integration points + + for (PropertyVec::const_iterator p = myProps.begin(); p != myProps.end(); p++) + if (p->pcode == Property::NEUMANN) + // Count the boundary integration points + if (p->patch > 0 && p->patch <= myModel.size()) + myModel[p->patch-1]->getNoBouPoints(nBouGP,p->ldim,p->lindx); + + // Let the integrand know how many integration points in total do we have + if (myProblem) myProblem->initIntegration(nIntGP,nBouGP); +} + + bool SIMbase::setPropertyType (int code, Property::Type ptype, int pindex) { if (code < 0) @@ -762,12 +767,11 @@ bool SIMbase::setPropertyType (int code, Property::Type ptype, int pindex) return false; } - for (size_t j = 0; j < myProps.size(); j++) - if (myProps[j].pindx == (size_t)code && - myProps[j].pcode == Property::UNDEFINED) + for (PropertyVec::iterator p = myProps.begin(); p != myProps.end(); p++) + if (p->pindx == (size_t)code && p->pcode == Property::UNDEFINED) { - myProps[j].pcode = ptype; - if (pindex >= 0) myProps[j].pindx = pindex; + p->pcode = ptype; + if (pindex >= 0) p->pindx = pindex; } return true; @@ -948,7 +952,7 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol, else ok = false; - if (j == 0 && ok) + if (j == 0) // All patches are referring to the same material, and we assume it has // been initialized during input processing (thus no initMaterial call here) for (i = 0; i < myModel.size() && ok; i++) @@ -1001,6 +1005,8 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol, else ok = false; + if (!ok) std::cerr <<" *** SIMbase::assembleSystem: Failure.\n"<< std::endl; + return ok && this->finalizeAssembly(newLHSmatrix); } @@ -1162,6 +1168,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time, } myProblem->initIntegration(time); + norm->initIntegration(nIntGP,nBouGP); const Vector& primsol = psol.front(); size_t i, j, k; @@ -1214,7 +1221,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time, else ok = false; - if (j == 0 && ok) + if (j == 0) // All patches are referring to the same material, and we assume it has // been initialized during input processing (thus no initMaterial call here) for (i = 0; i < myModel.size() && ok; i++) @@ -1257,6 +1264,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time, else ok = false; + if (!ok) std::cerr <<" *** SIMbase::solutionNorms: Failure.\n"<< std::endl; + // Clean up the dynamically allocated norm objects. This will also perform // the actual global norm assembly, in case the element norms are stored, // and always when multi-threading is used. diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index 90ad6369..939e8be3 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -470,6 +470,9 @@ public: NormBase* getNormIntegrand() const; protected: + //! \brief Allocates the problem-dependent integration point buffers, if any. + void initIntegrationBuffers(); + //! \brief Defines the type of a property set. //! \param[in] code The property code to be associated with the property type //! \param[in] ptype The property type to be associated with the given code @@ -614,7 +617,11 @@ protected: AlgEqSystem* myEqSys; //!< The actual linear equation system SAMpatch* mySam; //!< Auxiliary data for FE assembly management LinSolParams* mySolParams; //!< Input parameters for PETSc + private: + size_t nIntGP; //!< Number of interior integration points in the whole model + size_t nBouGP; //!< Number of boundary integration points in the whole model + //! \brief Parses a subelement of the tag from input file //! \param[in] elem The XML element to parse bool parseGeometryTag(const TiXmlElement* elem);