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

@@ -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<const Elasticity*>(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
     }
   }
 

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:
       ;

src/Integrands/Elasticity.C

@@ -21,6 +21,7 @@
 #include "Vec3Oper.h"
 #include "AnaSol.h"
 #include "VTF.h"
+#include <iomanip>
 
 #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)

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<Vec3Pair> tracVal; //!< Traction field point values
+  typedef std::pair<Vec3,double> PointValue; //!< Convenience type
+
+  mutable std::vector<PointValue> maxVal;  //!< Maximum result values
+  mutable std::vector<Vec3Pair>   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)