added support for K refinement cycles

CMakeLists.txt

@@ -75,6 +75,10 @@ IF(NOT MPI_FOUND OR IFEM_SERIAL_TESTS_IN_PARALLEL)
     ifem_add_vtf_test(round_stair_3patch_3D_LR_adap.vreg Poisson)
   endif()
   ifem_add_vtf_test(Square.vreg Poisson)
+  if(PETSC_FOUND)
+    ifem_add_test(Cube-kref.reg Poisson)
+    ifem_add_test(Square-kref.reg Poisson)
+  endif()
 endif()
 list(APPEND TEST_APPS Poisson)
 

SIMPoisson.h

@@ -16,15 +16,22 @@
 
 #include "SIMMultiPatchModelGen.h"
 #include "Poisson.h"
+#include "AlgEqSystem.h"
 #include "AnaSol.h"
+#include "ASMbase.h"
 #include "SIM1D.h"
 #include "SIM2D.h"
 #include "SIM3D.h"
 #include "Utilities.h"
 #include "DataExporter.h"
 #include "IFEM.h"
+#include "LinSolParams.h"
 #include "tinyxml.h"
+#include "TimeStep.h"
 #include <fstream>
+#ifdef HAS_PETSC
+#include "PETScMatrix.h"
+#endif
 
 
 /*!
@@ -34,6 +41,8 @@
 template<class Dim> class SIMPoisson : public SIMMultiPatchModelGen<Dim>
 {
 public:
+  using SetupProps = bool; //!< A bool is the only setup property
+
   //! \brief Default constructor.
   explicit SIMPoisson(bool checkRHS = false)
     : SIMMultiPatchModelGen<Dim>(1,checkRHS),
@@ -137,8 +146,8 @@ template<class Dim> class SIMPoisson : public SIMMultiPatchModelGen<Dim>
     return this->writeGlvBC(nBlock);
   }
 
-  //! \brief Assembles and solves the linear system.
-  bool solveStep(TimeStep&)
+  //! \brief Initialize simulator.
+  bool init()
   {
     if (!this->setMode(Dim::opt.eig == 0 ? SIM::STATIC : SIM::VIBRATION))
       return false;
@@ -149,6 +158,12 @@ template<class Dim> class SIMPoisson : public SIMMultiPatchModelGen<Dim>
       this->initSystem(Dim::opt.solver,1,0);
     this->setQuadratureRule(Dim::opt.nGauss[0],true);
 
+    return true;
+  }
+
+  //! \brief Assembles and solves the linear system.
+  bool solveStep(TimeStep&)
+  {
     if (!this->assembleSystem())
       return false;
     else if (vizRHS && Dim::opt.eig == 0)
@@ -187,7 +202,7 @@ template<class Dim> class SIMPoisson : public SIMMultiPatchModelGen<Dim>
   }
 
   //! \brief Saves solution-dependent quantities to file for postprocessing.
-  bool saveStep(TimeStep&, int& nBlock)
+  bool saveStep(TimeStep& tp, int& nBlock)
   {
     if (!asciiFile.empty())
     {
@@ -253,7 +268,10 @@ template<class Dim> class SIMPoisson : public SIMMultiPatchModelGen<Dim>
     if (!this->writeGlvN(myNorm,1,nBlock,prefix.data()))
       return false;
 
-    return this->writeGlvStep(1,0.0,1);
+    if (tp.iter == 0)
+      return this->writeGlvStep(1,0.0,1);
+    else
+      return this->writeGlvStep(tp.iter,tp.iter,1);
   }
 
   //! \brief Prints a norm group to the log stream.

Test/Cube-kref.reg

@@ -0,0 +1,44 @@
+Cube-kref.xinp -kref
+
+	Doing 3 k-refinement cycles.
+	K-solver tolerance 1e-07
+Number of elements    64
+Number of nodes       125
+Number of dofs        125
+Number of unknowns    27
+K-refinement cycle 1
+Number of elements    64
+Number of nodes       216
+Number of dofs        216
+Number of unknowns    64
+L2-norm            : 0.242676
+Max displacement   : 0.994383
+Energy norm           a(u^h,u^h)^0.5 : 1.92322
+External energy          (h,u^h)^0.5 : 1.92322
+Exact norm                a(u,u)^0.5 : 1.92382
+Exact error      a(e,e)^0.5, e=u-u^h : 0.0483303
+Exact relative error (%) : 2.5122
+K-refinement cycle 2
+Number of elements    64
+Number of nodes       343
+Number of dofs        343
+Number of unknowns    125
+L2-norm            : 0.225482
+Max displacement   : 1.38805
+Energy norm           a(u^h,u^h)^0.5 : 1.92374
+External energy          (h,u^h)^0.5 : 1.92374
+Exact norm                a(u,u)^0.5 : 1.92382
+Exact error      a(e,e)^0.5, e=u-u^h : 0.0180016
+Exact relative error (%) : 0.935719
+K-refinement cycle 3
+Number of elements    64
+Number of nodes       512
+Number of dofs        512
+Number of unknowns    216
+L2-norm            : 0.227619
+Max displacement   : 1.16167
+Energy norm           a(u^h,u^h)^0.5 : 1.92382
+External energy          (h,u^h)^0.5 : 1.92382
+Exact norm                a(u,u)^0.5 : 1.92382
+Exact error      a(e,e)^0.5, e=u-u^h : 0.000660832
+Exact relative error (%) : 0.0343499

Test/Cube-kref.xinp

@@ -0,0 +1,40 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<simulation>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <raiseorder patch="1" u="0" v="0" w="0"/>
+    <refine type="uniform" patch="1" u="3" v="3" w="3"/>
+    <topologysets>
+      <set name="Boundary" type="face">
+        <item patch="1">1 2 3 4 5 6</item>
+      </set>
+    </topologysets>
+  </geometry>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Boundary" comp="1"/>
+  </boundaryconditions>
+
+  <!-- Problem-specific block !-->
+  <poisson>
+    <source type="expression">3*PI*PI*sin(x*PI)*sin(y*PI)*sin(z*PI)</source>
+    <anasol type="expression">
+      <primary>sin(PI*x)*sin(PI*y)*sin(z*PI)</primary>
+      <secondary>-PI*cos(PI*x)*sin(PI*y)*sin(z*PI) |
+                 -PI*sin(PI*x)*cos(PI*y)*sin(PI*z) |
+                 -PI*sin(PI*x)*sin(PI*y)*cos(PI*z)
+      </secondary>
+    </anasol>
+  </poisson>
+
+  <krefinement cycles="3" rtol="1e-7"/>
+
+  <linearsolver>
+    <type>cg</type>
+    <rtol>1e-9</rtol>
+  </linearsolver>
+
+</simulation>

Test/Square-kref.reg

@@ -0,0 +1,44 @@
+Square-kref.xinp -kref
+
+	Doing 3 k-refinement cycles.
+	K-solver tolerance 1e-07
+Number of elements    64
+Number of nodes       81
+Number of dofs        81
+Number of unknowns    72
+K-refinement cycle 1
+Number of elements    64
+Number of nodes       100
+Number of dofs        100
+Number of unknowns    90
+L2-norm            : 0.986916
+Max displacement   : 2.00012
+Energy norm           a(u^h,u^h)^0.5 : 5.32009
+External energy          (h,u^h)^0.5 : 5.32009
+Exact norm                a(u,u)^0.5 : 5.32155
+Exact error      a(e,e)^0.5, e=u-u^h : 0.125387
+Exact relative error (%) : 2.3562
+K-refinement cycle 2
+Number of elements    64
+Number of nodes       121
+Number of dofs        121
+Number of unknowns    110
+L2-norm            : 1.02204
+Max displacement   : 2.2147
+Energy norm           a(u^h,u^h)^0.5 : 5.32152
+External energy          (h,u^h)^0.5 : 5.32152
+Exact norm                a(u,u)^0.5 : 5.32155
+Exact error      a(e,e)^0.5, e=u-u^h : 0.0180157
+Exact relative error (%) : 0.338542
+K-refinement cycle 3
+Number of elements    64
+Number of nodes       144
+Number of dofs        144
+Number of unknowns    132
+L2-norm            : 1.04048
+Max displacement   : 2.10254
+Energy norm           a(u^h,u^h)^0.5 : 5.32155
+External energy          (h,u^h)^0.5 : 5.32155
+Exact norm                a(u,u)^0.5 : 5.32155
+Exact error      a(e,e)^0.5, e=u-u^h : 0.00224356
+Exact relative error (%) : 0.0421598

Test/Square-kref.xinp

@@ -0,0 +1,46 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<simulation>
+
+  <!-- General - geometry definitions !-->
+  <geometry dim="2" scale="2.0">
+    <raiseorder patch="1" u="0" v="0"/>
+    <refine patch="1" u="7" v="7"/>
+    <topologysets>
+      <set name="Homogenous" type="edge">
+        <item patch="1">4</item>
+      </set>
+      <set name="Neumann" type="edge">
+        <item patch="1">3</item>
+      </set>
+    </topologysets>
+  </geometry>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Homogenous" comp="1"/>
+    <neumann set="Neumann" type="anasol"/>
+  </boundaryconditions>
+
+  <!-- Problem-specific block !-->
+  <poisson>
+    <source type="expression">PI*PI*cos(PI*x)*(2-y)</source>
+    <anasol type="expression">
+      <primary>cos(PI*x)*(2-y)</primary>
+      <secondary>PI*sin(PI*x)*(2-y)|cos(PI*x)</secondary>
+    </anasol>
+    <galerkin>1.0</galerkin>
+  </poisson>
+
+  <krefinement cycles="3" rtol="1e-7"/>
+
+  <linearsolver>
+    <type>cg</type>
+    <rtol>1e-9</rtol>
+  </linearsolver>
+
+  <discretization>
+    <nGauss default="0"/>
+  </discretization>
+
+</simulation>

Test/Square.vreg

@@ -16,3 +16,7 @@ Step 1
     Scalar: 'FE q_y'
     Scalar: 'Exact q_x'
     Scalar: 'Exact q_y'
+  Element block 2
+    32 nodes
+    32 elements
+    Vector: 'Heat flux'

main_Poisson.C

@@ -15,6 +15,7 @@
 #include "SIMPoisson.h"
 #include "SIMargsBase.h"
 #include "SIMSolverAdap.h"
+#include "SIMSolverKRef.h"
 #include "Utilities.h"
 #include "Profiler.h"
 
@@ -67,10 +68,31 @@ int runSimulator(char* infile, bool checkRHS,
   if (model.opt.dumpHDF5(infile))
     solver.handleDataOutput(model.opt.hdf5);
 
+  model.init();
   return solver.solveProblem(infile,"Solving the Poisson problem");
 }
 
 
+#ifdef HAS_PETSC
+template<class Dim>
+int runSimulatorK(char* infile, bool checkRHS,
+                  const std::vector<int>& ignoredPatches,
+                  bool fixDup, bool vizRHS, bool dumpASCII)
+{
+  SIMKRefWrap<SIMPoisson<Dim>> model(checkRHS), model2(checkRHS);
+  SIMSolverKRef<SIMPoisson<Dim>> solver(model, model2);
+
+  if (!solver.read(infile))
+    return 1;
+
+  if (model.opt.dumpHDF5(infile))
+    solver.handleDataOutput(model.opt.hdf5);
+
+  return solver.solveProblem(infile,"Solving the Poisson problem");
+}
+#endif
+
+
 /*!
   \brief Main program for the NURBS-based isogeometric Poisson equation solver.
 
@@ -119,6 +141,7 @@ int main (int argc, char** argv)
   bool vizRHS = false;
   bool fixDup = false;
   bool dumpASCII = false;
+  bool kref = false;
   char* infile = nullptr;
   SIMargsBase args("poisson");
 
@@ -137,6 +160,8 @@ int main (int argc, char** argv)
       SIMbase::ignoreDirichlet = true;
     else if (!strcmp(argv[i],"-checkRHS"))
       checkRHS = true;
+    else if (!strcmp(argv[i],"-kref"))
+      kref = true;
     else if (!strcmp(argv[i],"-vizRHS"))
       vizRHS = true;
     else if (!strcmp(argv[i],"-fixDup"))
@@ -197,8 +222,25 @@ int main (int argc, char** argv)
     case 3:
       return runSimulator<SIM3D,SIMSolverAdap>(infile, checkRHS, ignoredPatches,
                                                fixDup, vizRHS, dumpASCII);
-  }
-  else
+    }
+  else if (kref) {
+#ifndef HAS_PETSC
+    std::cerr << "K-refinement requested, but built without PETSc support. Bailing.." << std::endl;
+    return 1;
+#else
+    switch (args.dim) {
+    case 1:
+      return runSimulatorK<SIM1D>(infile, checkRHS, ignoredPatches,
+                                  fixDup, vizRHS, dumpASCII);
+    case 2:
+      return runSimulatorK<SIM2D>(infile, checkRHS, ignoredPatches,
+                                  fixDup, vizRHS, dumpASCII);
+    case 3:
+      return runSimulatorK<SIM3D>(infile, checkRHS, ignoredPatches,
+                                  fixDup, vizRHS, dumpASCII);
+    }
+#endif
+  } else
     switch (args.dim) {
     case 1:
       return runSimulator<SIM1D>(infile, checkRHS, ignoredPatches,