Exchange troubled cell

Some more features for LimitPass.

dune/fem-dg/operator/limiter/limitpass.hh

@@ -39,12 +39,6 @@
 #include <dune/fem-dg/operator/limiter/smoothness.hh>
 #include <dune/fem-dg/operator/limiter/indicatorbase.hh>
 
-#if HAVE_DUNE_OPTIM
-#define WANT_DUNE_OPTIM 1
-//#define WANT_DUNE_OPTIM 0
-#endif
-
-
 //*************************************************************
 namespace Dune
 {
@@ -257,7 +251,7 @@ namespace Fem
     typedef AdaptationMethod<GridType> AdaptationMethodType;
 
     //! id for choosing admissible linear functions
-    enum AdmissibleFunctions { DGFunctions = 0, ReconstructedFunctions = 1 , BothFunctions = 2 };
+    enum AdmissibleFunctions { DGFunctions = 0, ReconstructedFunctions = 1 , DGAndReconFunctions = 2, LPFunctions = 3 };
 
     //! returns true if pass is currently active in the pass tree
     using BaseType :: active ;
@@ -269,7 +263,6 @@ namespace Fem
     typedef TroubledCellIndicatorBase<ArgumentFunctionType> TroubledCellIndicatorType;
 
   protected:
-#if WANT_DUNE_OPTIM
     typedef typename GridPartType :: GridViewType  GridViewType ;
 
     struct BoundaryValue
@@ -287,8 +280,6 @@ namespace Fem
     };
 
     typedef Dune::FV::LPReconstruction< GridViewType, RangeType, BoundaryValue > LinearProgramming;
-#endif
-
 
   public:
     //- Public methods
@@ -333,11 +324,7 @@ namespace Fem
       dest_(0),
       spc_(spc),
       gridPart_(spc_.gridPart()),
-#if WANT_DUNE_OPTIM
-      linProg_( static_cast< GridViewType > (gridPart_), BoundaryValue( *this ),
-                Dune::Fem::Parameter::getValue<double>("finitevolume.linearprogramming.tol", 1e-8 )
-              ),
-#endif
+      linProg_(),
       indexSet_( gridPart_.indexSet() ),
       localIdSet_( gridPart_.grid().localIdSet()),
       cornerPointSetContainer_(),
@@ -369,6 +356,12 @@ namespace Fem
            ? new ModalSmoothnessIndicator< ArgumentFunctionType >() : nullptr ),
       counter_( 0 )
     {
+      if( admissibleFunctions_ == LPFunctions )
+      {
+        linProg_.reset( new LinearProgramming( static_cast< GridViewType > (gridPart_), BoundaryValue( *this ),
+                      Dune::Fem::Parameter::getValue<double>("finitevolume.linearprogramming.tol", 1e-8 )) );
+      }
+
       if( Parameter :: verbose () )
       {
         std::cout << "LimitDGPass (Grid is ";
@@ -443,9 +436,10 @@ namespace Fem
     AdmissibleFunctions getAdmissibleFunctions( const Dune::Fem::ParameterReader &parameter ) const
     {
       // default value
-      int val = 1;
+      //int val = 1; // reconstruction
+      int val = 3; // lp-reconstruction
       val = parameter.getValue("femdg.limiter.admissiblefunctions", val);
-      assert( val >= DGFunctions || val <= BothFunctions );
+      assert( val >= DGFunctions || val <= LPFunctions );
       return (AdmissibleFunctions) val;
     }
 
@@ -661,21 +655,22 @@ namespace Fem
         matrixCacheVec_.resize( numLevels );
       }
 
-#if WANT_DUNE_OPTIM
-      // helper class for evaluation of average value of discrete function
-      EvalAverage average( *this, U, discreteModel_);
-
-      values_.resize( size );
-      // evaluate average values on all cells
-      for( const auto& element : Dune::elements( gridPart_ ) )
+      if( linProg_ )
       {
-        average.evaluate( element, values_[ gridPart_.indexSet().index( element ) ] );
-      }
+        // helper class for evaluation of average value of discrete function
+        EvalAverage average( *this, U, discreteModel_);
 
-      gradients_.resize( size );
-      // get reconstructions
-      linProg_( gridPart_.indexSet(), values_, gradients_ );
-#endif
+        values_.resize( size );
+        // evaluate average values on all cells
+        for( const auto& element : Dune::elements( gridPart_ ) )
+        {
+          average.evaluate( element, values_[ gridPart_.indexSet().index( element ) ] );
+        }
+
+        gradients_.resize( size );
+        // get reconstructions
+        (*linProg_)( gridPart_.indexSet(), values_, gradients_ );
+      }
     }
 
     //! Some management (interface version)
@@ -794,8 +789,8 @@ namespace Fem
       // cache geometry type
       const GeometryType geomType = geo.type();
       // get bary center of element
-      const LocalDomainType& x = geoInfo_.localCenter( geomType );
-      const DomainType enBary = geo.global( x );
+      const LocalDomainType& wLocal = geoInfo_.localCenter( geomType );
+      const DomainType enBary = geomType.isNone() ? geo.center() : geo.global( wLocal );
 
       const IntersectionIteratorType endnit = gridPart_.iend( en );
       IntersectionIteratorType niter = gridPart_.ibegin(en);
@@ -829,9 +824,8 @@ namespace Fem
       // otherwise everything may be corrupted
       {
         // get barycenter of entity
-        const DomainType& enBaryLocal = (int(dimGrid) == int(dimDomain)) ?
-          geoInfo_.localCenter( geomType ) :
-          geo.local( enBary ) ;
+        const DomainType& enBaryLocal = (int(dimGrid) == int(dimDomain) && ! geomType.isNone() ) ?
+                wLocal : geo.local( enBary ) ;
 
         // check that average value is physical
         if(  discreteModel_.hasPhysical() &&
@@ -927,14 +921,12 @@ namespace Fem
       ++limitedElements_;
 
       const unsigned int enIndex = indexSet_.index( en );
-#if WANT_DUNE_OPTIM
-      bool useLinProg = true ;
-      if( useLinProg )
+      // use linear function from LP reconstruction
+      if( admissibleFunctions_ == LPFunctions )
       {
         deoMod_ = gradients_[ enIndex ];
       }
       else
-#endif
       {
         // obtain combination set
         ComboSetType& comboSet = storedComboSets_[ nbVals_.size() ];
@@ -1698,9 +1690,7 @@ namespace Fem
     const DiscreteFunctionSpaceType& spc_;
     GridPartType& gridPart_;
 
-#if WANT_DUNE_OPTIM
-    mutable LinearProgramming linProg_;
-#endif
+    std::unique_ptr< LinearProgramming > linProg_;
 
     const IndexSetType& indexSet_;
     const LocalIdSetType& localIdSet_;
@@ -1758,8 +1748,8 @@ namespace Fem
     mutable bool reconstruct_;
 
     // choice of admissible linear functions
-    const AdmissibleFunctions admissibleFunctions_;
-    mutable AdmissibleFunctions usedAdmissibleFunctions_ ;
+    const AdmissibleFunctions    admissibleFunctions_;
+    mutable AdmissibleFunctions  usedAdmissibleFunctions_ ;
 
     mutable std::vector< RangeType  > values_;
     mutable std::vector< GradientType > gradients_;