From c80fffc4c6a8c1f5b27607105b4dd98340d483b8 Mon Sep 17 00:00:00 2001
From: Stefan Girke <stefan.girke@wwu.de>
Date: Thu, 18 Feb 2016 14:47:48 +0100
Subject: [PATCH] add ldg option to algorithm creator, rnamed fluxes

---
 dune/fem-dg/misc/algorithmcreatorselector.hh  | 86 +++++++++++++++----
 dune/fem-dg/misc/configurator.hh              |  5 +-
 dune/fem-dg/operator/dg/fluxoperator.hh       | 67 +++++++--------
 .../operator/fluxes/diffusion/averageflux.hh  | 23 ++---
 4 files changed, 110 insertions(+), 71 deletions(-)

diff --git a/dune/fem-dg/misc/algorithmcreatorselector.hh b/dune/fem-dg/misc/algorithmcreatorselector.hh
index a28cf5e8..cc6994c4 100644
--- a/dune/fem-dg/misc/algorithmcreatorselector.hh
+++ b/dune/fem-dg/misc/algorithmcreatorselector.hh
@@ -47,6 +47,7 @@
 
 //include operators
 #include <dune/fem-dg/operator/dg/primaloperator.hh>
+#include <dune/fem-dg/operator/dg/fluxoperator.hh>
 
 #include <dune/fem-dg/assemble/primalmatrix.hh>
 #include <dune/fem/space/discontinuousgalerkin.hh>
@@ -203,6 +204,21 @@ namespace Fem
     };
   }
 
+  /**
+   *  \brief Namespace containing an Enum class to describe the formulation
+   */
+  namespace Formulation
+  {
+    /**
+     * \ingroup FemDGParameter
+     */
+    enum class Enum
+    {
+      primal,
+      dual
+    };
+  }
+
 ///////////////////////////////////////////////////////////////////////////
 // GridPartSelector
 ///////////////////////////////////////////////////////////////////////////
@@ -268,11 +284,11 @@ namespace Fem
 // AdvectionDiffusionOperatorSelector
 ///////////////////////////////////////////////////////////////////////////
 
-  template< class OperatorTraits, AdvectionLimiter::Enum op >
+  template< class OperatorTraits, Formulation::Enum form, AdvectionLimiter::Enum op >
   class AdvectionDiffusionOperatorSelector;
 
   template< class OperatorTraits >
-  class AdvectionDiffusionOperatorSelector< OperatorTraits, AdvectionLimiter::Enum::unlimited >
+  class AdvectionDiffusionOperatorSelector< OperatorTraits, Formulation::Enum::primal, AdvectionLimiter::Enum::unlimited >
   {
     static const int advection = OperatorTraits::ModelType::hasAdvection;
     static const int diffusion = OperatorTraits::ModelType::hasDiffusion;
@@ -289,7 +305,7 @@ namespace Fem
   };
 
   template< class OperatorTraits >
-  class AdvectionDiffusionOperatorSelector< OperatorTraits, AdvectionLimiter::Enum::limited >
+  class AdvectionDiffusionOperatorSelector< OperatorTraits, Formulation::Enum::primal, AdvectionLimiter::Enum::limited >
   {
     static const int advection = OperatorTraits::ModelType::hasAdvection;
     static const int diffusion = OperatorTraits::ModelType::hasDiffusion;
@@ -304,44 +320,78 @@ namespace Fem
     typedef typename ImplExplOperatorSelectorType::ExplicitOperatorType      ExplicitOperatorType;
   };
 
-  template< class OperatorTraits, AdvectionLimiter::Enum op, OperatorSplit::Enum split, Matrix::Enum ass >
+  template< class OperatorTraits >
+  class AdvectionDiffusionOperatorSelector< OperatorTraits, Formulation::Enum::dual, AdvectionLimiter::Enum::unlimited >
+  {
+    static const int advection = OperatorTraits::ModelType::hasAdvection;
+    static const int diffusion = OperatorTraits::ModelType::hasDiffusion;
+    typedef LDGAdvectionDiffusionOperator< OperatorTraits >                  DgType;
+    typedef LDGAdvectionDiffusionOperator< OperatorTraits >                  DgAdvectionType;
+    typedef LDGAdvectionDiffusionOperator< OperatorTraits >                  DgDiffusionType;
+    typedef ImplExplOperatorSelector< DgType, DgAdvectionType, DgDiffusionType, advection, diffusion >
+                                                                             ImplExplOperatorSelectorType;
+  public:
+    typedef typename ImplExplOperatorSelectorType::FullOperatorType          FullOperatorType;
+    typedef typename ImplExplOperatorSelectorType::ImplicitOperatorType      ImplicitOperatorType;
+    typedef typename ImplExplOperatorSelectorType::ExplicitOperatorType      ExplicitOperatorType;
+
+  };
+
+  template< class OperatorTraits >
+  class AdvectionDiffusionOperatorSelector< OperatorTraits, Formulation::Enum::dual, AdvectionLimiter::Enum::limited >
+  {
+    static const int advection = OperatorTraits::ModelType::hasAdvection;
+    static const int diffusion = OperatorTraits::ModelType::hasDiffusion;
+    typedef LDGLimitedAdvectionDiffusionOperator< OperatorTraits >           DgType;
+    typedef LDGLimitedAdvectionDiffusionOperator< OperatorTraits >           DgAdvectionType;
+    typedef LDGLimitedAdvectionDiffusionOperator< OperatorTraits >           DgDiffusionType;
+    typedef ImplExplOperatorSelector< DgType, DgAdvectionType, DgDiffusionType, advection, diffusion >
+                                                                             ImplExplOperatorSelectorType;
+  public:
+    typedef typename ImplExplOperatorSelectorType::FullOperatorType          FullOperatorType;
+    typedef typename ImplExplOperatorSelectorType::ImplicitOperatorType      ImplicitOperatorType;
+    typedef typename ImplExplOperatorSelectorType::ExplicitOperatorType      ExplicitOperatorType;
+  };
+
+
+  template< class OperatorTraits, Formulation::Enum form, AdvectionLimiter::Enum op, OperatorSplit::Enum split, Matrix::Enum ass >
   class OperatorSelector;
 
   //matrixfree
-  template< class OperatorTraits, AdvectionLimiter::Enum op >
-  struct OperatorSelector< OperatorTraits, op, OperatorSplit::Enum::expl, Matrix::Enum::matrixfree >
+  template< class OperatorTraits, Formulation::Enum form, AdvectionLimiter::Enum op >
+  struct OperatorSelector< OperatorTraits, form, op, OperatorSplit::Enum::expl, Matrix::Enum::matrixfree >
   {
-    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, op >::ExplicitOperatorType type;
+    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, form, op >::ExplicitOperatorType type;
   };
 
-  template< class OperatorTraits, AdvectionLimiter::Enum op >
-  struct OperatorSelector< OperatorTraits, op, OperatorSplit::Enum::impl, Matrix::Enum::matrixfree >
+  template< class OperatorTraits, Formulation::Enum form, AdvectionLimiter::Enum op >
+  struct OperatorSelector< OperatorTraits, form, op, OperatorSplit::Enum::impl, Matrix::Enum::matrixfree >
   {
-    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, op >::ImplicitOperatorType type;
+    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, form, op >::ImplicitOperatorType type;
   };
 
-  template< class OperatorTraits, AdvectionLimiter::Enum op >
-  struct OperatorSelector< OperatorTraits, op, OperatorSplit::Enum::full, Matrix::Enum::matrixfree >
+  template< class OperatorTraits, Formulation::Enum form, AdvectionLimiter::Enum op >
+  struct OperatorSelector< OperatorTraits, form, op, OperatorSplit::Enum::full, Matrix::Enum::matrixfree >
   {
-    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, op >::FullOperatorType type;
+    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, form, op >::FullOperatorType type;
   };
 
-  template< class OperatorTraits, AdvectionLimiter::Enum op >
-  struct OperatorSelector< OperatorTraits, op, OperatorSplit::Enum::rhs, Matrix::Enum::matrixfree >
+  template< class OperatorTraits, Formulation::Enum form, AdvectionLimiter::Enum op >
+  struct OperatorSelector< OperatorTraits, form, op, OperatorSplit::Enum::rhs, Matrix::Enum::matrixfree >
   {
-    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, op >::FullOperatorType type;
+    typedef typename AdvectionDiffusionOperatorSelector< OperatorTraits, form, op >::FullOperatorType type;
   };
 
   //assembled
   //TODO improve DGPrimalMatrixAssembly for correct splitting
   template< class AssemblerTraitsImp, AdvectionLimiter::Enum op >
-  struct OperatorSelector< AssemblerTraitsImp, op, OperatorSplit::Enum::full, Matrix::Enum::assembled >
+  struct OperatorSelector< AssemblerTraitsImp, Formulation::Enum::primal, op, OperatorSplit::Enum::full, Matrix::Enum::assembled >
   {
     typedef DGPrimalMatrixAssembly< AssemblerTraitsImp > type;
   };
 
   template< class AssemblerTraitsImp, AdvectionLimiter::Enum op >
-  struct OperatorSelector< AssemblerTraitsImp, op, OperatorSplit::Enum::rhs, Matrix::Enum::assembled >
+  struct OperatorSelector< AssemblerTraitsImp, Formulation::Enum::primal, op, OperatorSplit::Enum::rhs, Matrix::Enum::assembled >
   {
     typedef DGPrimalMatrixAssembly< AssemblerTraitsImp > type;
   };
diff --git a/dune/fem-dg/misc/configurator.hh b/dune/fem-dg/misc/configurator.hh
index f47bf002..0f7b4d54 100644
--- a/dune/fem-dg/misc/configurator.hh
+++ b/dune/fem-dg/misc/configurator.hh
@@ -36,7 +36,8 @@ namespace Fem
             AdvectionLimiter::Enum advLimitId,
             Matrix::Enum matrixId,
             AdvectionFlux::Enum advFluxId,
-            PrimalDiffusionFlux::Enum diffFluxId >
+            PrimalDiffusionFlux::Enum diffFluxId,
+            Formulation::Enum formId = Formulation::Enum::primal >
   class AlgorithmConfigurator
   {
     template< class AnalyticalTraitsImp, class NewModelImp >
@@ -100,7 +101,7 @@ namespace Fem
 
     //Operator/Assembler
     template< class OpTraits, OperatorSplit::Enum opSplit = OperatorSplit::Enum::full, Matrix::Enum matrix = matrixId >
-    using Operators = typename OperatorSelector< OpTraits, advLimitId, opSplit, matrix >::type;
+    using Operators = typename OperatorSelector< OpTraits, formId, advLimitId, opSplit, matrix >::type;
 
     //Matrix Containers
     template< class DomainDFSpace, class RangeDFSpace  >
diff --git a/dune/fem-dg/operator/dg/fluxoperator.hh b/dune/fem-dg/operator/dg/fluxoperator.hh
index 67bd22bf..21b80ea1 100644
--- a/dune/fem-dg/operator/dg/fluxoperator.hh
+++ b/dune/fem-dg/operator/dg/fluxoperator.hh
@@ -1,9 +1,5 @@
-#ifndef DUNE_FEM_DG_FLUXOPERATOR_HH
-#define DUNE_FEM_DG_FLUXOPERATOR_HH
-
-#ifdef HEADERCHECK
-#define FLUX 1
-#endif
+#ifndef DUNE_FEM_LDG_FLUXOPERATOR_HH
+#define DUNE_FEM_LDG_FLUXOPERATOR_HH
 
 #include <string>
 
@@ -27,11 +23,11 @@ namespace Dune
 namespace Fem
 {
 
-  // DGAdvectionDiffusionOperator
+  // LDGAdvectionDiffusionOperator
   //-----------------------------
 
   template< class Traits, bool advection = true , bool diffusion = true >
-  class DGAdvectionDiffusionOperator :
+  class LDGAdvectionDiffusionOperator :
     public Fem::SpaceOperatorInterface
       < typename Traits :: DestinationType >
   {
@@ -107,7 +103,7 @@ namespace Fem
     typedef typename Traits :: ExtraParameterTupleType ExtraParameterTupleType;
 
   public:
-    DGAdvectionDiffusionOperator( GridPartType& gridPart,
+    LDGAdvectionDiffusionOperator( GridPartType& gridPart,
                                   ProblemType& problem,
                                   ExtraParameterTupleType tuple =  ExtraParameterTupleType(),
                                   const std::string keyPrefix = "" ) :
@@ -182,10 +178,7 @@ namespace Fem
              <<", $\\eta = ";
       diffFlux_.diffusionFluxPenalty( stream );
       stream <<"$, {\\bf Adv. Flux:} ";
-      if (FLUX==1)
-        stream <<"LLF";
-      else if (FLUX==2)
-        stream <<"HLL";
+      //TODO has to be implemented
       stream <<",\\\\\n";
       return stream.str();
     }
@@ -230,7 +223,7 @@ namespace Fem
   //--------------------
 
   template< class OpTraits >
-  class DGAdvectionOperator : public
+  class LDGAdvectionOperator : public
     DGAdvectionDiffusionOperatorBase< LDGAdvectionTraits< OpTraits, true > >
   {
     typedef LDGAdvectionTraits< OpTraits, true> Traits;
@@ -240,7 +233,7 @@ namespace Fem
     typedef typename BaseType :: ProblemType   ProblemType;
     typedef typename BaseType :: ExtraParameterTupleType ExtraParameterTupleType;
 
-    DGAdvectionOperator( GridPartType& gridPart, ProblemType& problem,
+    LDGAdvectionOperator( GridPartType& gridPart, ProblemType& problem,
                          ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
                          const std::string keyPrefix = "" )
       : BaseType( gridPart, problem, tuple, keyPrefix )
@@ -251,33 +244,33 @@ namespace Fem
       std::stringstream stream;
       stream <<"{\\bf Adv. Op.}, flux formulation, order: " << Traits::polynomialOrder+1
              <<", {\\bf Adv. Flux:} ";
-      if (FLUX==1)
+      /*if (FLUX==1)
         stream <<"LLF";
       else if (FLUX==2)
-        stream <<"HLL";
+        stream <<"HLL";*/
       stream <<",\\\\\n";
       return stream.str();
     }
   };
 
 
-  // DGDiffusionOperator
+  // LDGDiffusionOperator
   //--------------------
 
   template< class Traits >
-  class DGDiffusionOperator : public
-    DGAdvectionDiffusionOperator< Traits, false >
+  class LDGDiffusionOperator : public
+    LDGAdvectionDiffusionOperator< Traits, false >
   {
-    typedef DGAdvectionDiffusionOperator< Traits, false >  BaseType;
+    typedef LDGAdvectionDiffusionOperator< Traits, false >  BaseType;
 
   public:
     typedef typename BaseType :: GridPartType  GridPartType;
     typedef typename BaseType :: ProblemType   ProblemType;
     typedef typename BaseType :: ExtraParameterTupleType ExtraParameterTupleType;
 
-    DGDiffusionOperator( GridPartType& gridPart, ProblemType& problem,
-                         ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
-                         const std::string keyPrefix = "" )
+    LDGDiffusionOperator( GridPartType& gridPart, ProblemType& problem,
+                          ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
+                          const std::string keyPrefix = "" )
       : BaseType( gridPart, problem, tuple, keyPrefix )
     {}
 
@@ -298,7 +291,7 @@ namespace Fem
   };
 
 
-  // DGLimitedAdvectionDiffusionOperator
+  // LDGLimitedAdvectionDiffusionOperator
   //------------------------------------
 
   /** \class DGLimitedAdvectionDiffusionOperator
@@ -312,7 +305,7 @@ namespace Fem
    */
   template< class Traits,
             bool advection = true >
-  class DGLimitedAdvectionDiffusionOperator :
+  class LDGLimitedAdvectionDiffusionOperator :
     public Fem::SpaceOperatorInterface
       < typename Traits :: DestinationType >
   {
@@ -403,9 +396,9 @@ namespace Fem
     };
 
   public:
-    DGLimitedAdvectionDiffusionOperator( GridPartType& gridPart, ProblemType& problem,
-                                         ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
-                                         const std::string keyPrefix = "" )
+    LDGLimitedAdvectionDiffusionOperator( GridPartType& gridPart, ProblemType& problem,
+                                          ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
+                                          const std::string keyPrefix = "" )
       : model_( problem )
       , numflux_( model_ )
       , gridPart_( gridPart )
@@ -427,7 +420,7 @@ namespace Fem
       problem1_.setIndicator( &indicator_ );
     }
 
-    ~DGLimitedAdvectionDiffusionOperator() { delete uTmp_; }
+    ~LDGLimitedAdvectionDiffusionOperator() { delete uTmp_; }
 
     void setTime(const double time) {
 	    pass3_.setTime( time );
@@ -484,10 +477,10 @@ namespace Fem
              <<", penalty: ";
       diffFlux_.diffusionFluxPenalty( stream );
       stream <<", {\\bf Adv. Flux:} ";
-      if (FLUX==1)
+      /*if (FLUX==1)
         stream <<"LLF";
       else if (FLUX==2)
-        stream <<"HLL";
+        stream <<"HLL";*/
       stream <<",\\\\\n";
       return stream.str();
     }
@@ -533,12 +526,12 @@ namespace Fem
       < Traits, u, advection, diffusion> DiscreteModelType;
   };
 
-  // DGAdaptationIndicatorOperator
+  // LDGAdaptationIndicatorOperator
   //------------------------------
 
   template< class OpTraits,
             bool advection, bool diffusion = false >
-  struct DGAdaptationIndicatorOperator : public
+  struct LDGAdaptationIndicatorOperator : public
     DGAdvectionDiffusionOperatorBase<
        AdaptationIndicatorTraits< OpTraits, advection, diffusion > >
   {
@@ -548,9 +541,9 @@ namespace Fem
     typedef typename BaseType :: ProblemType   ProblemType ;
     typedef typename BaseType :: ExtraParameterTupleType ExtraParameterTupleType;
 
-    DGAdaptationIndicatorOperator( GridPartType& gridPart, ProblemType& problem,
-                                   ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
-                                   const std::string keyPrefix = "" )
+    LDGAdaptationIndicatorOperator( GridPartType& gridPart, ProblemType& problem,
+                                    ExtraParameterTupleType  tuple =  ExtraParameterTupleType(),
+                                    const std::string keyPrefix = "" )
       : BaseType( gridPart, problem, tuple, keyPrefix )
     {
       if ( Fem::Parameter::verbose() )
diff --git a/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh b/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh
index 7c2fd226..04b9cb74 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh
@@ -1,5 +1,5 @@
-#ifndef DUNE_LDGFLUX_HH
-#define DUNE_LDGFLUX_HH
+#ifndef DUNE_FEM_DG_LDGAVERAGEFLUX_HH
+#define DUNE_FEM_DG_LDGAVERAGEFLUX_HH
 
 // Dune-Fem includes
 #include "fluxbase.hh"
@@ -17,7 +17,7 @@ namespace Fem
   template <class DiscreteFunctionSpaceImp,
             class ModelImp,
             class FluxParameterImp >
-  class LDGDiffusionFlux :
+  class LDGAverageDiffusionFlux :
     public DGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp >
   {
     typedef DGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp > BaseType;
@@ -45,14 +45,9 @@ namespace Fem
 
     typedef typename BaseType :: ParameterType  ParameterType;
 
-#if DUNE_VERSION_NEWER_REV(DUNE_FEM,1,1,0)
     // type of gradient space
     typedef typename DiscreteFunctionSpaceType ::
         template ToNewDimRange< dimGradRange > :: Type   DiscreteGradientSpaceType;
-#else
-    // type of gradient space
-    typedef CombinedSpace< DiscreteFunctionSpaceType, dimGradRange>  DiscreteGradientSpaceType;
-#endif
 
     typedef typename DiscreteGradientSpaceType :: RangeType GradientRangeType;
     typedef typename DiscreteGradientSpaceType :: JacobianRangeType GradientJacobianType;
@@ -62,7 +57,7 @@ namespace Fem
 
   private:
     // no copying
-    LDGDiffusionFlux(const LDGDiffusionFlux& other);
+    LDGAverageDiffusionFlux(const LDGAverageDiffusionFlux& other);
   protected:
     using BaseType :: determineDirection;
     using BaseType :: model_;
@@ -75,9 +70,9 @@ namespace Fem
     /**
      * \brief constructor
      */
-    LDGDiffusionFlux(GridPartType& gridPart,
-                     const ModelImp& mod,
-                     const ParameterType& param ) :
+    LDGAverageDiffusionFlux(GridPartType& gridPart,
+                            const ModelImp& mod,
+                            const ParameterType& param ) :
       BaseType( mod, true, param ),
       penalty_( parameter().penalty() ),
       // Set CFL number for penalty term (compare diffusion in first pass)
@@ -85,7 +80,7 @@ namespace Fem
     {
       if( Fem::Parameter::verbose () )
       {
-        std::cout << "LDGDiffusionFlux: penalty = " << penalty_ << std::endl;
+        std::cout << "LDGAverageDiffusionFlux: penalty = " << penalty_ << std::endl;
       }
     }
 
@@ -374,7 +369,7 @@ namespace Fem
     const double penalty_;
     const bool penaltyTerm_;
 
-  }; // end LDGDiffusionFlux
+  }; // end LDGAverageDiffusionFlux
 
 } // end namespace
 } // end namespace
-- 
GitLab