added: add support for adaptive simulators using the ISolver interface

in particular, this allows adaptivity for stationary, nonlinear systems
such as NavierStokes or Boussinesq.

Apps/Common/SIMSolverAdap.h

@@ -16,6 +16,39 @@
 
 #include "SIMSolver.h"
 #include "AdaptiveSIM.h"
+#include "TimeStep.h"
+
+
+/*!
+  \brief Adaptive simulator driver using the ISolver interface.
+*/
+
+template<class T1>
+class AdaptiveISolver : public AdaptiveSIM
+{
+public:
+  //! \brief The constructor initializes default adaptation parameters.
+  //! \param sim The FE model
+  //! \param sa If \e true, this is a stand-alone driver
+  AdaptiveISolver(T1& sim, bool sa = true) : AdaptiveSIM(sim, sa), model(sim) {}
+
+protected:
+  //! \brief Assemble and solve the equation system.
+  virtual bool assembleAndSolveSystem()
+  {
+    TimeStep dummy;
+    model.init(dummy);
+    bool result = model.solveStep(dummy);
+    if (result)
+      solution = model.getSolutions();
+
+    rCond = 1.0;
+
+    return result;
+  }
+
+  T1& model; //!< Reference to the actual sim
+};
 
 
 /*!
@@ -24,17 +57,18 @@
   ISolver interface. It provides an adaptive loop with data output.
 */
 
-template<class T1> class SIMSolverAdap : public SIMSolverStat<T1>
+template<class T1, class AdapSim>
+class SIMSolverAdapImpl : public SIMSolverStat<T1>
 {
 public:
   //! \brief The constructor forwards to the parent class constructor.
-  explicit SIMSolverAdap(T1& s1) : SIMSolverStat<T1>(s1), aSim(s1,false)
+  explicit SIMSolverAdapImpl(T1& s1) : SIMSolverStat<T1>(s1), aSim(s1,false)
   {
     this->S1.setSol(&aSim.getSolution());
   }
 
   //! \brief Empty destructor.
-  virtual ~SIMSolverAdap() {}
+  virtual ~SIMSolverAdapImpl() {}
 
   //! \brief Reads solver data from the specified input file.
   virtual bool read(const char* file) { return this->SIMadmin::read(file); }
@@ -68,7 +102,10 @@ protected:
     return aSim.parse(elem);
   }
 
-  AdaptiveSIM aSim; //!< Adaptive simulation driver
+  AdapSim aSim; //!< Adaptive simulation driver
 };
 
+template<class T1>
+using SIMSolverAdap = SIMSolverAdapImpl<T1, AdaptiveSIM>; //!< Convenience alias template
+
 #endif

src/SIM/AdaptiveSIM.C

@@ -222,6 +222,28 @@ bool AdaptiveSIM::initAdaptor (size_t normGroup)
 }
 
 
+bool AdaptiveSIM::assembleAndSolveSystem()
+{
+  // Assemble the linear FE equation system
+  if (!model.setMode(SIM::STATIC,true))
+    return false;
+  model.setQuadratureRule(opt.nGauss[0],true);
+  if (!model.assembleSystem())
+    return false;
+
+  // Solve the linear system of equations
+  int printSol = 1;
+  solution.resize(model.getNoRHS());
+  for (size_t i = 0; i < solution.size(); i++)
+    if (!model.solveSystem(solution[i],printSol,&rCond,"displacement",i==0,i))
+      return false;
+    else if (solution.size() > 2)
+      printSol = 0; // Print summary only for the first two solutions
+
+  return true;
+}
+
+
 bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF,
                              std::streamsize precision)
 {
@@ -242,21 +264,11 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF,
     // Output the initial grid to eps-file
     model.refine(LR::RefineData(),"mesh_001.eps");
 
-  // Assemble the linear FE equation system
-  model.setMode(SIM::STATIC,true);
+  // Assemble and solve the FE equation system
   model.initSystem(opt.solver,1,model.getNoRHS(),0,withRF);
-  model.setQuadratureRule(opt.nGauss[0],true);
-  if (!model.assembleSystem())
-    return false;
 
-  // Solve the linear system of equations
-  int printSol = 1;
-  solution.resize(model.getNoRHS());
-  for (size_t i = 0; i < solution.size(); i++)
-    if (!model.solveSystem(solution[i],printSol,&rCond,"displacement",i==0,i))
-      return false;
-    else if (solution.size() > 2)
-      printSol = 0; // Print summary only for the first two solutions
+  if (!this->assembleAndSolveSystem())
+    return false;
 
   eNorm.clear();
   gNorm.clear();

src/SIM/AdaptiveSIM.h

@@ -76,6 +76,13 @@ public:
   //! \param[in] elem The XML element to parse
   virtual bool parse(const TiXmlElement* elem);
 
+protected:
+  Vectors solution; //!< All solutions (galerkin projections)
+  double rCond;     //!< Actual reciprocal condition number
+
+  //! \brief Assemble and solve the equation system.
+  virtual bool assembleAndSolveSystem();
+
 private:
   SIMoutput& model; //!< The isogeometric FE model
   bool       alone; //!< If \e false, this class is wrapped by SIMSolver
@@ -85,7 +92,6 @@ private:
   double beta;         //!< Refinement percentage in each step
   double errTol;       //!< Global error stop tolerance
   double condLimit;    //!< Upper limit on condition number
-  double rCond;        //!< Actual reciprocal condition number
   int    maxStep;      //!< Maximum number of adaptive refinements
   int    maxDOFs;      //!< Maximum number of degrees of freedom
   int    symmetry;     //!< Always refine a multiplum of this
@@ -103,7 +109,6 @@ private:
 
   size_t  adaptor;  //!< Norm group to base the mesh adaptation on
   size_t  adNorm;   //!< Which norm to base the mesh adaptation on
-  Vectors solution; //!< All solutions (galerkin projections)
   Vectors gNorm;    //!< Global norms
   Matrix  eNorm;    //!< Element norms