Add new options for symmetry/penalty in local operator

New options are OOP (no penalty term) and IIP (no symmetry term)
Also adapted error estimator (for OOP scheme)
Moved error estimator into Operator namespace
Changed initialization orders in both operators to suppress warnings

dune/dorie/interface/adaptivity.hh

@@ -42,7 +42,7 @@ private:
 	/// Adaptivity GFS
 	using AGFS = typename AGFSHelper::Type;
 	/// Error estimator local operator
-	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,Boundary>;
+	using ESTLOP = Dune::Dorie::Operator::FluxErrorEstimator<Traits,Param,Boundary>;
 	/// Empty constraints container
 	using NoTrafo = Dune::PDELab::EmptyTransformation;
 	/// Solution vector type

dune/dorie/solver/operator_DG.hh

@@ -19,7 +19,9 @@ namespace Operator {
 struct RichardsDGMethod
 {
   enum Type { NIPG, //!< Non-symmetric interior penalty
-              SIPG //!< Symmetric interior penalty (default)
+              SIPG, //!< Symmetric interior penalty (default)
+              OOB, //!< Oden, Babuska, Baumann (no penalty term)
+              IIP //!< Incomplete interior penalty (no symmetry)
             };
 };
 
@@ -85,13 +87,13 @@ protected:
   const SourceTerm& sourceTerm; //!< class providing source term information
   const RichardsDGMethod::Type method; //!< Interior penalty type
   const RichardsDGWeights::Type weights; //!< DG method weights switch
-  const RF penalty_factor; //!< Interior penalty factor
   const Mapper mapper; //!< maps grid elements (cells) to natural numbers
 
   const int intorderadd; //!< integration order addend
   const int quadrature_factor; //!< factor to FEM integration order
-  RF time; //!< operator internal time for BC, parameter, and source term queries
-  RF theta; //!< Penalty factor
+  RF penalty_factor; //!< Interior penalty factor
+  RF theta; //!< Symmetry factor
+  RF time;  //!< operator internal time for BC, parameter, and source term queries
   //! Cache for fixing the BC type over one entire time step
   mutable std::map<typename std::size_t,BoundaryCondition::Type> bc_type_cache;
 
@@ -115,28 +117,33 @@ public:
    *  \param view_ Grid view
    *  \param boundary_ Object prividing boundary conditions
    *  \param sourceTerm_ Object providing source term information
-   *  \param method_ DG Penalty type \see RichardsDGMethod::Type
+   *  \param method_ DG method type \see RichardsDGMethod::Type
    *  \param weights_ DG weigths switch \see RichardsDGWeights::Type
    *  \param intorderadd_ Addend to integration order
    *  \param quadrature_factor_ Factor to FEM integration order
    */
   RichardsDGSpatialOperator (const Dune::ParameterTree& config, const Parameter& param_, const GridView view_,
       Boundary& boundary_, const SourceTerm& sourceTerm_,
-      //RichardsDGMethod::Type method_ = RichardsDGMethod::NIPG,
       RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
       RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
       int intorderadd_ = 2, int quadrature_factor_ = 2)
     : Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
     Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
     Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
-    param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), weights(weights_),
-    penalty_factor(config.get<RF>("dg.penaltyFactor")), mapper(view_),
+    param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), weights(weights_), mapper(view_),
     intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
+    penalty_factor(config.get<RF>("dg.penaltyFactor")),
+    time(0.0),
     cache(20)
   {
-    theta = 1.;
-    if (method == RichardsDGMethod::SIPG) theta = -1.;
-    //theta = 0.;
+    theta = 0.; // IIP
+    if (method == RichardsDGMethod::NIPG || method == RichardsDGMethod::OOB)
+      theta = 1.;
+    else if (method == RichardsDGMethod::SIPG)
+      theta = -1.;
+    
+    if (method == RichardsDGMethod::OOB)
+      penalty_factor = 0.;
   }
 
   /// Volume integral depending on test and ansatz functions

dune/dorie/solver/operator_error_indicator.hh

@@ -12,6 +12,7 @@
 
 namespace Dune {
   namespace Dorie {
+  namespace Operator { 
 
     /// Local operator for residual-based error estimation.
     /*
@@ -48,18 +49,23 @@ namespace Dune {
 
       const Parameter& param;
       const Boundary& boundary;
-      const RF penalty_factor;
       const int intorderadd;
       const int quadrature_factor;
+      RF penalty_factor;
       RF time;
 
       FluxErrorEstimator (const Dune::ParameterTree& config,
         const Parameter& param_, const Boundary& boundary_,
+        RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
         int intorderadd_ = 2, int quadrature_factor_ = 2)
       : param(param_), boundary(boundary_),
+        intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
         penalty_factor(config.get<RF>("dg.penaltyFactor")),
-        intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
-      { }
+        time(0.0)
+      {
+        if (method_ == RichardsDGMethod::OOB)
+          penalty_factor = 0.0;
+      }
 
       /// skeleton integral depending on test and ansatz functions
       /** each face is only visited ONCE!
@@ -330,7 +336,8 @@ namespace Dune {
 
     };
 
-  }
-}
+  } // namespace Operator
+  } // namespace Dorie
+} // namespace Dune
 
 #endif // DUNE_DORIE_ERROR_INDICATOR_HH
\ No newline at end of file