Changed: Use a lambda function in assembleSystem to reduce code duplication.

Changed: solveMatrixSystem is renamed to solveSystem (overloaded method).
Added: Support for dual solution field, for goal-oriented error estimates.

src/SIM/AdaptiveSIM.C

@@ -44,6 +44,9 @@ bool AdaptiveSIM::initAdaptor (size_t normGroup)
     return false;
 
   projs.resize(opt.project.size());
+  if (model.haveDualSol())
+    projd.resize(opt.project.size());
+
   if (opt.format >= 0)
   {
     prefix.reserve(opt.project.size());
@@ -101,7 +104,9 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF,
   };
 
   gNorm.clear();
+  dNorm.clear();
   eNorm.clear();
+  fNorm.clear();
 
   model.getProcessAdm().cout <<"\nAdaptive step "<< iStep << std::endl;
   if (iStep > 1)
@@ -127,17 +132,17 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF,
   if (!this->assembleAndSolveSystem())
     return failure();
 
-  eNorm.clear();
-  gNorm.clear();
-
   // Project the secondary solution onto the splines basis
   size_t idx = 0;
   model.setMode(SIM::RECOVERY);
   for (const SIMoptions::ProjectionMap::value_type& prj : opt.project)
     if (prj.first <= SIMoptions::NONE)
-      projs[idx++].clear(); // No projection for this norm group
+      idx++; // No projection for this norm group
     else if (!model.project(projs[idx++],solution.front(),prj.first))
       return failure();
+    else if (idx == adaptor && idx <= projd.size() && solution.size() > 1)
+      if (!model.project(projd[idx-1],solution[1],prj.first))
+        return failure();
 
   if (msgLevel > 1 && !projs.empty())
     model.getProcessAdm().cout << std::endl;
@@ -147,6 +152,9 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF,
   model.setQuadratureRule(opt.nGauss[1]);
   if (!model.solutionNorms(solution.front(),projs,eNorm,gNorm))
     return failure();
+  if (!projd.empty() && solution.size() > 1)
+    if (!model.solutionNorms(solution[1],projd,fNorm,dNorm))
+      return failure();
 
   model.setMode(SIM::RECOVERY);
   if (!model.dumpResults(solution.front(),0.0,
@@ -172,15 +180,21 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec)
   if (outPrec > 0)
   {
     std::streamsize oldPrec = IFEM::cout.precision(outPrec);
-    this->printNorms(gNorm,eNorm);
+    this->printNorms(gNorm,dNorm,eNorm);
     IFEM::cout.precision(oldPrec);
   }
   else
-    this->printNorms(gNorm,eNorm);
+    this->printNorms(gNorm,dNorm,eNorm);
+
+  // Calculate refinement indicators including the dual error estimates
+  Vector refIn = eNorm.getRow(eRow);
+  if (eRow <= fNorm.rows())
+    for (size_t i = 1; i <= refIn.size(); i++)
+      refIn(i) = sqrt(refIn(i)*fNorm(eRow,i));
 
   // Set up refinement parameters
   LR::RefineData prm;
-  if (this->calcRefinement(prm,iStep,gNorm,eNorm.getRow(eRow)) <= 0)
+  if (this->calcRefinement(prm,iStep,gNorm,refIn) <= 0)
     return false;
 
   // Now refine the mesh and write out resulting grid

src/SIM/AdaptiveSIM.h

@@ -79,12 +79,15 @@ protected:
 
 private:
   Vectors gNorm; //!< Global norms
+  Vectors dNorm; //!< Dual global norms
   Matrix  eNorm; //!< Element norms
+  Matrix  fNorm; //!< Dual element norms
 
   int geoBlk; //!< Running VTF geometry block counter
   int nBlock; //!< Running VTF result block counter
 
   std::vector<Vector>      projs;  //!< Projected secondary solutions
+  std::vector<Vector>      projd;  //!< Projected dual solutions
   std::vector<std::string> prefix; //!< Norm prefices for VTF-output
 
 protected:

src/SIM/AdaptiveSetup.C

@@ -438,7 +438,7 @@ int AdaptiveSetup::calcRefinement (LR::RefineData& prm, int iStep,
 }
 
 
-void AdaptiveSetup::printNorms (const Vectors& gNorm,
+void AdaptiveSetup::printNorms (const Vectors& gNorm, const Vectors& dNorm,
                                 const Matrix& eNorm, size_t w) const
 {
   model.printNorms(gNorm,w);
@@ -457,6 +457,14 @@ void AdaptiveSetup::printNorms (const Vectors& gNorm,
     if (model.haveAnaSol() && adaptor != 0)
       IFEM::cout <<"\nEffectivity index  : "
                  << model.getEffectivityIndex(gNorm,adaptor,adNorm);
+    if (adaptor < dNorm.size())
+    {
+      const Vector& bNorm = dNorm[adaptor];
+      IFEM::cout <<"\nError estimate, dual solution"
+                 << utl::adjustRight(w-29,"(z)") << bNorm(adNorm)
+                 <<"\nRelative error (%) : "
+                 << 100.0*bNorm(adNorm)/hypot(dNorm.front()(1),bNorm(2));
+    }
     IFEM::cout <<"\n"<< std::endl;
 
     delete norm;

src/SIM/AdaptiveSetup.h

@@ -69,9 +69,10 @@ public:
 
   //! \brief Prints out global norms to the log stream.
   //! \param[in] gNorm Global norms
+  //! \param[in] dNorm Global dual norms
   //! \param[in] eNorm Element norms
   //! \param[in] w Field width for global norm labels
-  void printNorms(const Vectors& gNorm,
+  void printNorms(const Vectors& gNorm, const Vectors& dNorm,
                   const Matrix& eNorm, size_t w = 36) const;
 
   //! \brief Dumps current mesh to external file(s) for inspection.

src/SIM/SIMbase.C

@@ -661,7 +661,7 @@ size_t SIMbase::getNoEquations () const
 
 size_t SIMbase::getNoRHS () const
 {
-  return myEqSys ? myEqSys->getNoRHS() : 1;
+  return myEqSys ? myEqSys->getNoRHS() : (this->haveDualSol() ? 2 : 1);
 }
 
 
@@ -809,6 +809,37 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
 {
   PROFILE1("Element assembly");
 
+  // Lambda function for assembling the interior terms for a given patch
+  auto&& assembleInterior = [this,time,prevSol](IntegrandBase* integrand,
+                                                GlobalIntegral& integral,
+                                                ASMbase* pch, int pidx)
+  {
+    if (msgLevel > 1)
+      IFEM::cout <<"\nAssembling interior matrix terms for P"<< pidx
+                 << std::endl;
+
+    if (!this->initBodyLoad(pidx))
+      return false;
+
+    if (!this->extractPatchSolution(integrand,prevSol,pidx-1))
+      return false;
+
+    if (dualField && myProblem->getExtractionField())
+    {
+      Matrix extrField(*myProblem->getExtractionField());
+      if (!pch->L2projection(extrField,dualField))
+        return false;
+    }
+
+    if (mySol)
+      mySol->initPatch(pch->idx);
+
+#if SP_DEBUG > 2
+    integrand->printSolution(std::cout,pch->idx+1);
+#endif
+    return pch->integrate(*integrand,integral,time);
+  };
+
   bool ok = true;
   bool isAssembling = (myProblem->getMode() > SIM::INIT &&
                        myProblem->getMode() < SIM::RECOVERY);
@@ -840,6 +871,7 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
       for (p = myProps.begin(); p != myProps.end() && ok; ++p)
         if (p->pcode == Property::MATERIAL &&
             (it->first == 0 || it->first == p->pindx))
+        {
           if (!(pch = this->getPatch(p->patch)))
           {
             std::cerr <<" *** SIMbase::assembleSystem: Patch index "<< p->patch
@@ -849,39 +881,17 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
           else if (this->initMaterial(p->pindx))
           {
             lp = p->patch;
-            if (msgLevel > 1)
+            ok = assembleInterior(it->second,sysQ,pch,lp);
-              IFEM::cout <<"\nAssembling interior matrix terms for P"<< lp
-                         << std::endl;
-            ok &= this->initBodyLoad(lp);
-            ok &= this->extractPatchSolution(it->second,prevSol,lp-1);
-            if (mySol)
-              mySol->initPatch(pch->idx);
-#if SP_DEBUG > 2
-            it->second->printSolution(std::cout,pch->idx+1);
-#endif
-            ok &= pch->integrate(*it->second,sysQ,time);
           }
           else
             ok = false;
+        }
 
       if (lp == 0 && it->first == 0)
         // All patches refer to the same material, and we assume it has been
         // initialized during input processing (thus no initMaterial call here)
-        for (size_t k = 0; k < myModel.size() && ok; k++)
+        for (lp = 0; lp < myModel.size() && ok; lp++)
-        {
+          ok = assembleInterior(it->second,sysQ,myModel[lp],1+lp);
-          lp = k+1;
-          if (msgLevel > 1)
-            IFEM::cout <<"\nAssembling interior matrix terms for P"<< lp
-                       << std::endl;
-          ok &= this->initBodyLoad(lp);
-          ok &= this->extractPatchSolution(it->second,prevSol,k);
-          if (mySol)
-            mySol->initPatch(myModel[k]->idx);
-#if SP_DEBUG > 2
-          it->second->printSolution(std::cout,myModel[k]->idx+1);
-#endif
-          ok &= myModel[k]->integrate(*it->second,sysQ,time);
-        }
     }
 
     // Assemble contributions from the Neumann boundary conditions
@@ -1087,17 +1097,17 @@ bool SIMbase::solveSystem (Vector& solution, int printSol, double* rCond,
 }
 
 
-bool SIMbase::solveMatrixSystem (Vectors& solution, int printSol,
+bool SIMbase::solveSystem (Vectors& solution, int printSol, const char* cmpName)
-                                 const char* compName)
 {
-  solution.resize(this->getNoRHS());
+  size_t nSol = myEqSys ? myEqSys->getNoRHS() : 0;
-  for (size_t i = 0; i < solution.size(); i++)
+  if (solution.size() < nSol)
-    if (!this->solveSystem(solution[i],printSol,nullptr,compName,i==0,i))
+    solution.resize(nSol);
-      return false;
-    else if (solution.size() > 2)
-      printSol = 0; // Print summary only for the first two solutions
 
-  return true;
+  bool status = nSol > 0;
+  for (size_t i = 0; i < nSol && status; i++)
+    status = this->solveSystem(solution[i],printSol,nullptr,cmpName,i==0,i);
+
+  return status;
 }
 
 

src/SIM/SIMbase.h

@@ -285,11 +285,11 @@ public:
   }
 
   //! \brief Solves a linear system of equations with multiple right-hand-sides.
-  //! \param[out] solution Global primary solution vector
+  //! \param[out] solution Global primary solution vectors
   //! \param[in] printSol Print solution if its size is less than \a printSol
-  //! \param[in] compName Solution name to be used in norm output
+  //! \param[in] cmpName Solution name to be used in norm output
-  bool solveMatrixSystem(Vectors& solution, int printSol = 0,
+  bool solveSystem(Vectors& solution, int printSol = 0,
-                         const char* compName = "displacement");
+                   const char* cmpName = "displacement");
 
   //! \brief Finds the DOFs showing the worst convergence behavior.
   //! \param[in] x Global primary solution vector