added: support K-refinement cycles

CMakeLists.txt

@@ -41,6 +41,9 @@ else()
   ifem_add_test(Lshape_p2_b5.reg Darcy)
   ifem_add_test(Wavefront_k10_p2_b20.reg Darcy)
   ifem_add_vtf_test(Wavefront_k10_p2_b20.vreg Darcy)
+  if(PETSC_FOUND)
+  ifem_add_test(Wavefront_k10_p2_b20_kref.reg Darcy)
+  endif()
 endif()
 list(APPEND TEST_APPS Darcy)
 

Darcy.C

@@ -258,6 +258,17 @@ Vec3 Darcy::getBodyForce(const Vec3& X) const
 }
 
 
+void Darcy::setMode (SIM::SolutionMode mode)
+{
+  m_mode = mode;
+
+  if (mode >= SIM::RECOVERY)
+    primsol.resize(1);
+  else
+    primsol.clear();
+}
+
+
 DarcyNorm::DarcyNorm (Darcy& p, VecFunc* a) : NormBase(p), anasol(a)
 {
   nrcmp = myProblem.getNoFields(2);
@@ -366,12 +377,6 @@ bool DarcyNorm::finalizeElement (LocalIntegral& elmInt,
 }
 
 
-void DarcyNorm::addBoundaryTerms (Vectors& gNorm, double energy) const
-{
-  gNorm.front()[1] += energy;
-}
-
-
 size_t DarcyNorm::getNoFields (int group) const
 {
   if (group < 1)

Darcy.h

@@ -65,44 +65,48 @@ class Darcy : public IntegrandBase
   //! \brief Evaluates the potential source (if any) at specified point.
   double getPotential(const Vec3& X) const;
 
+  //! \brief Defines the solution mode before the element assembly is started.
+  //! \param[in] mode The solution mode to use
+  void setMode(SIM::SolutionMode mode) override;
+
   using IntegrandBase::getLocalIntegral;
   //! \brief Returns a local integral contribution object for given element.
   //! \param[in] nen Number of nodes on element
   //! \param[in] neumann Whether or not we are assembling Neumann BCs
-  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
-                                          bool neumann) const;
+  LocalIntegral* getLocalIntegral(size_t nen, size_t,
+                                  bool neumann) const override;
 
   using IntegrandBase::evalInt;
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
-                       const Vec3& X) const;
+  bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
+               const Vec3& X) const override;
 
   using IntegrandBase::evalBou;
   //! \brief Evaluates the integrand at a boundary point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[in] normal Boundary normal vector at current integration point
-  virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
-                       const Vec3& X, const Vec3& normal) const;
+  bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
+               const Vec3& X, const Vec3& normal) const override;
 
   //! \brief Sets up the permeability matrix
   //! \param[out] K \f$ nsd\times nsd\f$-matrix or its inverse
   //! \param[in] X Cartesian coordinates of current point
   //! \param[in] inverse If \e true, set up the inverse matrix instead
-  virtual bool formKmatrix(Matrix& K, const Vec3& X, bool inverse = false) const;
+  bool formKmatrix(Matrix& K, const Vec3& X, bool inverse = false) const;
 
   using IntegrandBase::evalSol;
   //! \brief Evaluated the secondary solution at a result point
   //! param[out] s Array of solution field values at current point
   //! param[in] fe Finite element data at current point
   //! param[in] X Cartesian coordinates of current point
   //! param(in) MNPC Nodal point correspondence for the basis function values
-  virtual bool evalSol(Vector& s, const FiniteElement& fe,
-                       const Vec3& X, const std::vector<int>& MNPC) const;
+  bool evalSol(Vector& s, const FiniteElement& fe,
+               const Vec3& X, const std::vector<int>& MNPC) const override;
 
   //! \brief Evaluates the finite element (FE) solution at an integration point
   //! \param[out] s The FE solution values at current point
@@ -114,23 +118,23 @@ class Darcy : public IntegrandBase
 
   //! \brief Returns the number of primary/secondary solution field components.
   //! \param[in] fld which field set to consider (1=primary, 2=secondary)
-  virtual size_t getNoFields(int fld = 2) const { return fld > 1 ? nsd : 1; }
+  size_t getNoFields(int fld = 2) const override { return fld > 1 ? nsd : 1; }
 
   //! \brief Returns the name of the primary solution field.
   //! \param[in] prefix Name prefix
-  virtual std::string getField1Name(size_t, const char* prefix = 0) const;
+  std::string getField1Name(size_t, const char* prefix = 0) const override;
 
   //! \brief Returns the name of a secondary solution field component
   //! \param[in] i Field component index
   //! \param[in] prefix Name prefix for all components
-  virtual std::string getField2Name(size_t i, const char* prefix = 0) const;
+  std::string getField2Name(size_t i, const char* prefix = 0) const override;
 
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
   //! \note The Integrand object is allocated dynamically and has to be deleted
   //! manually when leaving the scope of the pointer variable receiving the
   //! returned pointer value.
   //! \param[in] asol Pointer to analytical solution (optional)
-  virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
+  NormBase* getNormIntegrand(AnaSol* asol = 0) const override;
 
 private:
   VecFunc* permvalues;    //!< Permeability function.
@@ -167,45 +171,40 @@ class DarcyNorm : public NormBase
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
-                       const Vec3& X) const;
+  bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
+               const Vec3& X) const override;
 
   using NormBase::evalBou;
   //! \brief Evaluates the integrand at a boundary point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite Element quantities
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[in] normal Boundary normal vector at current integration point
-  virtual bool evalBou(LocalIntegral& elmInt,  const FiniteElement& fe,
-                       const Vec3& X, const Vec3& normal) const;
+  bool evalBou(LocalIntegral& elmInt,  const FiniteElement& fe,
+               const Vec3& X, const Vec3& normal) const override;
 
   using NormBase::finalizeElement;
   //! \brief Finalizes the element norms after the numerical integration.
   //! \details This method is used to compute effectivity indices.
   //! \param elmInt The local integral object to receive the contributions
-  virtual bool finalizeElement(LocalIntegral& elmInt, const TimeDomain&,size_t);
+  bool finalizeElement(LocalIntegral& elmInt, const TimeDomain&,size_t) override;
 
   //! \brief Returns whether this norm has explicit boundary contributions.
-  virtual bool hasBoundaryTerms() const { return true; }
-
-  //! \brief Adds external energy terms to relevant norms.
-  //! \param gNorm Global norm quantities
-  //! \param[in] energy Global external energy
-  virtual void addBoundaryTerms(Vectors& gNorm, double energy) const;
+  bool hasBoundaryTerms() const override { return true; }
 
   //! \brief Returns the number of norm groups or size of a specified group.
   //! \param[in] group The norm group to return the size of
   //! (if zero, return the number of groups)
-  virtual size_t getNoFields(int group = 0) const;
+  size_t getNoFields(int group = 0) const override;
 
   //! \brief Returns the name of a norm quantity.
   //! \param[in] i The norm group (one-based index)
   //! \param[in] j The norm number (one-based index)
   //! \param[in] prefix Common prefix for all norm names
-  virtual std::string getName(size_t i, size_t j, const char* prefix) const;
+  std::string getName(size_t i, size_t j, const char* prefix) const override;
 
   //! \brief Returns whether a norm quantity stores element contributions.
-  virtual bool hasElementContributions(size_t i, size_t j) const;
+  bool hasElementContributions(size_t i, size_t j) const override;
 
 private:
   VecFunc* anasol; //!< Analytical heat flux

SIMDarcy.h

@@ -24,6 +24,7 @@
 #include "tinyxml.h"
 #include "Property.h"
 #include "DataExporter.h"
+#include "TimeStep.h"
 
 
 /*!
@@ -33,8 +34,9 @@
 template<class Dim> class SIMDarcy : public Dim
 {
 public:
+  using SetupProps = bool; //!< Only a bool setup property.
   //! \brief Default constructor.
-  SIMDarcy() : Dim(1), drc(Dim::dimension), solVec(&sol)
+  SIMDarcy(bool checkRHS = false) : Dim(1,checkRHS), drc(Dim::dimension), solVec(&sol)
   {
     Dim::myProblem = &drc;
     aCode[0] = aCode[1] = 0;
@@ -160,7 +162,7 @@ template<class Dim> class SIMDarcy : public Dim
   void setSol(const Vector* sol) { solVec = sol; }
 
   //! \brief Return solution vector.
-  const Vector& getSolution(int=0) { return *solVec; }
+  Vector& getSolution(int=0) { return sol; }
 
   //! \brief Register fields for data export
   void registerFields(DataExporter& exporter)
@@ -178,13 +180,13 @@ template<class Dim> class SIMDarcy : public Dim
     if (Dim::opt.format < 0)
       return true;
 
-    nBlock = 0;
     return this->writeGlvG(geoBlk,fileName);
   }
 
   //! \brief Saves the converged results to VTF file of a given time step.
+  //! \param[in] tp Time stepping information
   //! \param[in] nBlock Running VTF block counter
-  bool saveStep(const TimeStep&, int& nBlock)
+  bool saveStep(const TimeStep& tp, int& nBlock)
   {
     if (Dim::opt.format < 0)
       return true;
@@ -202,18 +204,27 @@ template<class Dim> class SIMDarcy : public Dim
         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 Computes the solution for the current time step.
-  bool solveStep(TimeStep&)
+  //! \brief Initialize simulator.
+  bool init()
   {
-    if (!this->setMode(SIM::DYNAMIC))
+    if (!this->setMode(SIM::STATIC))
       return false;
 
     this->initSystem(Dim::opt.solver,1,1,0,true);
     this->setQuadratureRule(Dim::opt.nGauss[0],true);
 
+    return true;
+  }
+
+  //! \brief Computes the solution for the current time step.
+  bool solveStep(TimeStep&)
+  {
     if (!this->assembleSystem())
       return false;
 
@@ -273,6 +284,7 @@ template<class Dim> class SIMDarcy : public Dim
     // Evaluate solution norms
     Matrix eNorm;
     Vectors gNorm;
+    this->setMode(SIM::RECOVERY);
     this->setQuadratureRule(Dim::opt.nGauss[1]);
     if (this->solutionNorms(this->getSolution(),eNorm,gNorm))
       this->printNorms(gNorm);

Test/Wavefront_k10_p2_b20_kref.reg

@@ -0,0 +1,46 @@
+Wavefront_k10_p2_b20_kref.xinp -kref
+
+K-refinement cycle 1
+Number of elements    16
+Number of nodes       25
+Number of dofs        25
+Number of constraints 16
+Number of unknowns    9
+Number of elements    16
+Number of nodes       36
+Number of dofs        36
+Number of constraints 20
+Number of unknowns    16
+L2-norm            : 3.41662
+Max pressure       : 6.63686
+Energy norm           a(u^h,u^h)^0.5 : 32.9142
+External energy          (h,u^h)^0.5 : 38.6113
+Exact norm                a(u,u)^0.5 : 30.1213
+Exact error      a(e,e)^0.5, e=u-u^h : 48.0936
+Exact relative error (%) : 159.666
+K-refinement cycle 2
+Number of elements    16
+Number of nodes       49
+Number of dofs        49
+Number of constraints 24
+Number of unknowns    25
+L2-norm            : 3.01658
+Max pressure       : 13.4129
+Energy norm           a(u^h,u^h)^0.5 : 48.3129
+External energy          (h,u^h)^0.5 : 43.9144
+Exact norm                a(u,u)^0.5 : 27.5448
+Exact error      a(e,e)^0.5, e=u-u^h : 41.3086
+Exact relative error (%) : 149.969
+K-refinement cycle 3
+Number of elements    16
+Number of nodes       64
+Number of dofs        64
+Number of constraints 28
+Number of unknowns    36
+L2-norm            : 2.00667
+Max pressure       : 4.96432
+Energy norm           a(u^h,u^h)^0.5 : 22.374
+External energy          (h,u^h)^0.5 : 21.0557
+Exact norm                a(u,u)^0.5 : 24.2404
+Exact error      a(e,e)^0.5, e=u-u^h : 17.7332
+Exact relative error (%) : 73.1558

Test/Wavefront_k10_p2_b20_kref.xinp

@@ -0,0 +1,34 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+<simulation>
+  <!-- General - geometry definitions !-->
+  <geometry dim="2" sets="true">
+    <raiseorder patch="1" u="0" v="0"/>
+    <refine type="uniform" patch="1" u="3" v="3"/>
+  </geometry>
+
+  <!-- General - numerical integration scheme !-->
+  <discretization>
+    <nGauss default="0"/>
+  </discretization>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Boundary" comp="1" type="anasol"/>
+  </boundaryconditions>
+
+  <!-- Problem-specific block !-->
+  <darcy>
+    <permvalues>98.1|9.81</permvalues>
+    <bodyforce>0|0</bodyforce>
+    <source type="Wavefront"/>
+    <anasol type="Wavefront"/>
+  </darcy>
+
+  <krefinement cycles="3" rtol="1e-7"/>
+
+  <linearsolver>
+    <type>gmres</type>
+    <rtol>1e-9</rtol>
+  </linearsolver>
+
+</simulation>