From 0f33e80af7a918390b1e5e504a64f5060e83ba12 Mon Sep 17 00:00:00 2001 From: kmo Date: Sun, 5 Aug 2012 18:10:17 +0000 Subject: [PATCH] Added: Print of maximum values of some secondary solution quantities. Fixed: Reset rigid body positions and the reference norm on cut-back. git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1839 e10b68d5-8a6e-419e-a041-bce267b0401d --- Apps/FiniteDefElasticity/NonlinearDriver.C | 22 +++++++++-- .../NonlinearElasticityUL.C | 5 +++ src/Integrands/Elasticity.C | 39 +++++++++++++++++++ src/Integrands/Elasticity.h | 12 +++++- 4 files changed, 73 insertions(+), 5 deletions(-) diff --git a/Apps/FiniteDefElasticity/NonlinearDriver.C b/Apps/FiniteDefElasticity/NonlinearDriver.C index c90fcc57..15f2a69a 100644 --- a/Apps/FiniteDefElasticity/NonlinearDriver.C +++ b/Apps/FiniteDefElasticity/NonlinearDriver.C @@ -13,6 +13,7 @@ #include "NonlinearDriver.h" #include "SIMbase.h" +#include "Elasticity.h" #include "DataExporter.h" #include "tinyxml.h" @@ -138,6 +139,7 @@ int NonlinearDriver::solveProblem (bool skip2nd, bool energyNorm, if (dtDump <= 0.0) dtDump = params.stopTime + 1.0; double nextDump = params.time.t + dtDump; double nextSave = params.time.t + model->opt.dtSave; + const Elasticity* elp = dynamic_cast(model->getProblem()); int iStep = 0; // Save initial state to VTF if (model->opt.format >= 0 && params.multiSteps()) @@ -156,10 +158,11 @@ int NonlinearDriver::solveProblem (bool skip2nd, bool energyNorm, if (stat == NonLinSIM::DIVERGED) { // Try cut-back with a smaller time step when diverging - if (params.cutback()) - std::copy(solution[1].begin(),solution[1].end(),solution[0].begin()); - else - break; + if (!params.cutback()) break; + + std::copy(solution[1].begin(),solution[1].end(),solution[0].begin()); + model->updateConfiguration(solution.front()); + refNorm = 1.0; // Reset the reference norm } // Solve the nonlinear FE problem at this load step @@ -188,15 +191,26 @@ int NonlinearDriver::solveProblem (bool skip2nd, bool energyNorm, { // Save solution variables to VTF for visualization if (model->opt.format >= 0) + { if (!this->saveStep(++iStep,params.time.t,skip2nd)) return 6; + // Print out the maximum values + if (elp) + { + elp->printMaxVals(std::cout,outPrec,7); + elp->printMaxVals(std::cout,outPrec,8); + } + } + // Save solution variables to HDF5 if (writer) if (!writer->dumpTimeLevel()) return 7; nextSave = params.time.t + model->opt.dtSave; + if (nextSave > params.stopTime) + nextSave = params.stopTime; // Always save the final step } } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C index 1b71ba03..9c4e8fef 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C @@ -68,6 +68,11 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode) case SIM::BUCKLING: eKm = 1; eKg = 2; + break; + + case SIM::RECOVERY: + maxVal.clear(); + break; default: ; diff --git a/src/Integrands/Elasticity.C b/src/Integrands/Elasticity.C index ec99e005..ba66d44a 100644 --- a/src/Integrands/Elasticity.C +++ b/src/Integrands/Elasticity.C @@ -21,6 +21,7 @@ #include "Vec3Oper.h" #include "AnaSol.h" #include "VTF.h" +#include #ifndef epsR //! \brief Zero tolerance for the radial coordinate. @@ -118,6 +119,10 @@ void Elasticity::setMode (SIM::SolutionMode mode) eS = 1; break; + case SIM::RECOVERY: + maxVal.clear(); + break; + default: ; } @@ -523,6 +528,15 @@ bool Elasticity::evalSol (Vector& s, for (int i = 1; i <= material->getNoIntVariables(); i++) s.push_back(material->getInternalVariable(i)); + // Find the maximum values for each quantity + if (maxVal.empty()) + for (size_t j = 0; j < s.size(); j++) + maxVal.push_back(std::make_pair(X,s[j])); + else + for (size_t j = 0; j < s.size(); j++) + if (fabs(s[j]) > fabs(maxVal[j].second)) + maxVal[j] = std::make_pair(X,s[j]); + return true; } @@ -631,6 +645,31 @@ const char* Elasticity::getField2Name (size_t i, const char* prefix) const } +void Elasticity::printMaxVals (std::ostream& os, + std::streamsize precision, size_t comp) const +{ + size_t i1 = 1, i2 = maxVal.size(); + if (comp > i2) + return; + else if (comp > 0) + i1 = i2 = comp; + + for (size_t i = i1; i <= i2; i++) + { + const char* name = this->getField2Name(i-1); + os <<" Max "<< name <<":"; + for (size_t j = strlen(name); j < 16; j++) std::cout <<' '; + std::streamsize flWidth = 8 + precision; + std::streamsize oldPrec = os.precision(precision); + std::ios::fmtflags oldF = os.flags(std::ios::scientific | std::ios::right); + os << std::setw(flWidth) << maxVal[i-1].second; + os.precision(oldPrec); + os.flags(oldF); + std::cout <<" X = "<< maxVal[i-1].first << std::endl; + } +} + + NormBase* Elasticity::getNormIntegrand (AnaSol* asol) const { if (asol) diff --git a/src/Integrands/Elasticity.h b/src/Integrands/Elasticity.h index a45875a6..6860297c 100644 --- a/src/Integrands/Elasticity.h +++ b/src/Integrands/Elasticity.h @@ -142,6 +142,13 @@ public: //! \param[in] prefix Name prefix for all components virtual const char* getField2Name(size_t i, const char* prefix = 0) const; + //! \brief Prints out the maximum secondary solution values. + //! \param os Output stream to write the values to + //! \param[in] precision Number of digits after the decimal point + //! \param[in] comp Which component to print (0 means all) + void printMaxVals(std::ostream& os, std::streamsize precision, + size_t comp = 0) const; + protected: //! \brief Calculates some kinematic quantities at current point. //! \param[in] eV Element solution vector @@ -222,7 +229,10 @@ protected: VecFunc* fluxFld; //!< Pointer to explicit boundary traction field VecFunc* bodyFld; //!< Pointer to body force field - mutable std::vector tracVal; //!< Traction field point values + typedef std::pair PointValue; //!< Convenience type + + mutable std::vector maxVal; //!< Maximum result values + mutable std::vector tracVal; //!< Traction field point values unsigned short int nsd; //!< Number of space dimensions (1, 2 or 3) unsigned short int nDF; //!< Dimension on deformation gradient (2 or 3)