Exchange troubled cell

Some more features for LimitPass.

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

@@ -38,12 +38,6 @@
 
 #include <dune/fem-dg/operator/limiter/smoothness.hh>
 
-#if HAVE_DUNE_OPTIM
-#define WANT_DUNE_OPTIM 1
-//#define WANT_DUNE_OPTIM 0
-#endif
-
-
 //*************************************************************
 namespace Dune
 {
@@ -256,7 +250,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 ;
@@ -265,10 +259,9 @@ namespace Fem
     typedef CheckCartesian< GridPartType >  CheckCartesianType;
 
     //! function describing an external troubled cell indicator
-    typedef const std::function< double(const DestinationType&, const LocalFunctionType& ) > TroubledCellIndicatorType;
+    typedef const std::function< double(const ArgumentFunctionType&, const LocalFunctionType& ) > TroubledCellIndicatorType;
 
   protected:
-#if WANT_DUNE_OPTIM
     typedef typename GridPartType :: GridViewType  GridViewType ;
 
     struct BoundaryValue
@@ -286,8 +279,6 @@ namespace Fem
     };
 
     typedef Dune::FV::LPReconstruction< GridViewType, RangeType, BoundaryValue > LinearProgramming;
-#endif
-
 
   public:
     //- Public methods
@@ -332,11 +323,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_(),
@@ -365,9 +352,15 @@ namespace Fem
       admissibleFunctions_( getAdmissibleFunctions( parameter ) ),
       usedAdmissibleFunctions_( admissibleFunctions_ ),
       extTroubledCellIndicator_( indicator_ == 2
-           ? &ModalSmoothnessIndicator< DestinationType >::troubledCellIndicator : nullptr ),
+           ? &ModalSmoothnessIndicator< ArgumentFunctionType >::troubledCellIndicator : 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 ";
@@ -437,9 +430,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;
     }
 
@@ -655,21 +649,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)
@@ -788,8 +783,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);
@@ -823,9 +818,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() &&
@@ -921,14 +915,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() ];
@@ -1419,7 +1411,7 @@ namespace Fem
     }
 
     // calculate shock detector
-    bool calculateIndicator(const DestinationType& U,
+    bool calculateIndicator(const ArgumentFunctionType& U,
                             const LocalFunctionType& uEn,
                             const Geometry& geo,
                             const bool initLimiter,
@@ -1598,7 +1590,6 @@ namespace Fem
       // multiply h pol ord with circume
       const double circFactor = (circume > 0.0) ? (hPowPolOrder/(circume * tolFactor_ )) : 0.0;
 
-      double modalInd = -1;
       if( extTroubledCellIndicator_ )
       {
         shockIndicator = ( tol_1_ * extTroubledCellIndicator_( U, uEn ) );
@@ -1693,9 +1684,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_;
@@ -1753,8 +1742,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_;