Exchange troubled cell

Some more features for LimitPass.

dune/fem-dg/models/modelwrapper.hh

@@ -16,6 +16,8 @@
 
 #include <dune/fem/schemes/diffusionmodel.hh>
 
+#include <dune/fem/misc/threads/threadsafevalue.hh>
+
 // fem-dg includes
 #include <dune/fem-dg/models/defaultmodel.hh>
 #include <dune/fem-dg/models/defaultprobleminterfaces.hh>
@@ -198,7 +200,7 @@ namespace Fem
         limitedRange_[ i ] = i;
 
       // if method has been filled then modified will be set differently
-      advection_.limitedRange( limitedRange_ );
+      advection().limitedRange( limitedRange_ );
     }
 
 #ifdef EULER_WRAPPER_TEST
@@ -213,13 +215,13 @@ namespace Fem
 #if HAVE_DUNE_FEMPY
       // update model times (only if time method is available on these models)
       //! TODO problem without virtualization advection_.setTime(t);
-      //! TODO problem without virtualization diffusion_.setTime(t);
+      //! TODO problem without virtualization diffusion().setTime(t);
       ::detail::CallSetTime< AdvectionModelType,
                              ::detail::CheckTimeMethod< AdvectionModelType >::value >
-        ::setTime( const_cast< AdvectionModelType& > (advection_), t );
+        ::setTime( const_cast< AdvectionModelType& > (advection()), t );
       ::detail::CallSetTime< DiffusionModelType,
                              ::detail::CheckTimeMethod< DiffusionModelType >::value >
-        ::setTime( const_cast< DiffusionModelType& > (diffusion_), t );
+        ::setTime( const_cast< DiffusionModelType& > (diffusion()), t );
 #endif
     }
 
@@ -229,9 +231,9 @@ namespace Fem
     void setEntity( const Entity& entity ) const
     {
       if( hasAdvection )
-        advection_.init( entity );
+        advection().init( entity );
       if( hasDiffusion )
-        diffusion_.init( entity );
+        diffusion().init( entity );
     }
 
     inline bool hasStiffSource() const { return AdditionalType::hasStiffSource; }
@@ -259,7 +261,7 @@ namespace Fem
                                RangeType & s) const
     {
       assert( hasDiffusion );
-      diffusion_.source( local.quadraturePoint(), u, du, s );
+      diffusion().source( local.quadraturePoint(), u, du, s );
       return 0;
     }
 
@@ -271,7 +273,7 @@ namespace Fem
                                   RangeType& s) const
     {
       assert( hasAdvection );
-      advection_.source( local.quadraturePoint(), u, du, s );
+      advection().source( local.quadraturePoint(), u, du, s );
       return 0;
     }
 
@@ -282,7 +284,7 @@ namespace Fem
                            JacobianRangeType& result ) const
     {
       assert( hasAdvection );
-      advection_.diffusiveFlux( local.quadraturePoint(), u, du, result );
+      advection().diffusiveFlux( local.quadraturePoint(), u, du, result );
     }
 
     template <class LocalEvaluation>
@@ -292,7 +294,7 @@ namespace Fem
     {
       if( hasDiffusion )
       {
-        maxValue = diffusion_.diffusionTimeStep( local.entity(), local.quadraturePoint(), 0.0, u );
+        maxValue = diffusion().diffusionTimeStep( local.entity(), local.quadraturePoint(), 0.0, u );
       }
     }
 
@@ -301,7 +303,7 @@ namespace Fem
                                      const T& circumEstimate,
                                      const RangeType& u ) const
     {
-      return diffusion_.diffusionTimeStep( local.entity(), local.quadraturePoint(), circumEstimate, u );
+      return diffusion().diffusionTimeStep( local.entity(), local.quadraturePoint(), circumEstimate, u );
     }
 
     // is not used
@@ -314,7 +316,7 @@ namespace Fem
       assert( hasAdvection );
       assert( 0 ); // if this is used we have to check if this is correct
       // TODO: u != ubar and du != dubar
-      advection_.linDiffusiveFlux( u, du, local.quadraturePoint(), u, du, A);
+      advection().linDiffusiveFlux( u, du, local.quadraturePoint(), u, du, A);
     }
 
     template <class LocalEvaluation>
@@ -329,7 +331,7 @@ namespace Fem
       RangeType u; // fake return variable
       const int id = getBoundaryId( local );
       // the following fails since is is called with
-      return advection_.hasBoundaryValue(id, time(), local.entity(), local.position(), u, u);
+      return advection().hasBoundaryValue(id, time(), local.entity(), local.position(), u, u);
     }
 
     // return uRight for insertion into the numerical flux
@@ -342,7 +344,7 @@ namespace Fem
 #ifndef NDEBUG
       const bool isDirichlet =
 #endif
-      advection_.boundaryValue(id, time(), local.entity(), local.quadraturePoint(),
+      advection().boundaryValue(id, time(), local.entity(), local.quadraturePoint(),
                                     local.intersection().unitOuterNormal( local.localPosition() ),
                                     uLeft, uRight);
       assert( isDirichlet );
@@ -363,7 +365,7 @@ namespace Fem
 #ifndef NDEBUG
       const bool isFluxBnd =
 #endif
-      advection_.boundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, gLeft);
+      advection().boundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, gLeft);
       gLeft *= len;
       assert( isFluxBnd );
       return 0; // TODO: do something better here
@@ -376,7 +378,7 @@ namespace Fem
                     JacobianRangeType& diff ) const
     {
       assert( hasDiffusion );
-      diffusion_.diffusiveFlux( local.quadraturePoint(), u, du, diff);
+      diffusion().diffusiveFlux( local.quadraturePoint(), u, du, diff);
     }
 
 
@@ -395,7 +397,7 @@ namespace Fem
 #ifndef NDEBUG
       const bool isFluxBnd =
 #endif
-      diffusion_.diffusionBoundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, jacLeft, gLeft);
+      diffusion().diffusionBoundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, jacLeft, gLeft);
       gLeft *= len;
       assert( isFluxBnd );
       return 0; // QUESTION: do something better here? Yes, return time step restriction if possible
@@ -411,7 +413,7 @@ namespace Fem
       // TODO: add a max speed for the diffusion time step control
       // this needs to be added in diffusionTimeStep
       assert( hasAdvection );
-      advspeed = advection_.maxSpeed( time(), local.entity(), local.quadraturePoint(), unitNormal, u );
+      advspeed = advection().maxSpeed( time(), local.entity(), local.quadraturePoint(), unitNormal, u );
       totalspeed = advspeed;
     }
 
@@ -425,7 +427,7 @@ namespace Fem
     inline DomainType velocity (const LocalEvaluation& local,
                                 RangeType& u) const
     {
-      return advection_.velocity( time(), local.entity(), local.quadraturePoint(), u );
+      return advection().velocity( time(), local.entity(), local.quadraturePoint(), u );
     }
 
     /////////////////////////////////////////////////////////////////
@@ -437,7 +439,7 @@ namespace Fem
                           const RangeType& u,
                           DomainType& velocity) const
     {
-      velocity = advection_.velocity( time(), en, x, u );
+      velocity = advection().velocity( time(), en, x, u );
     }
 
     // we have physical check for this model
@@ -452,7 +454,7 @@ namespace Fem
                          const DomainType& x,
                          const RangeType& u) const
     {
-      return advection_.physical( entity, x, u ) > 0;
+      return advection().physical( entity, x, u ) > 0;
     }
 
     // adjust average value if necessary
@@ -463,7 +465,7 @@ namespace Fem
                              RangeType& u ) const
     {
       // nothing to be done here for this test case
-      advection_.adjustAverageValue( entity, xLocal, u );
+      advection().adjustAverageValue( entity, xLocal, u );
     }
 
     // calculate jump between left and right value
@@ -474,7 +476,7 @@ namespace Fem
                      const RangeType& uRight,
                      RangeType& jump) const
     {
-      jump = advection_.jump( it, x, uLeft, uRight );
+      jump = advection().jump( it, x, uLeft, uRight );
     }
 
     // calculate jump between left and right value
@@ -485,7 +487,7 @@ namespace Fem
                                      const RangeType& uRight,
                                      RangeType& indicator) const
     {
-      indicator = advection_.jump( it, x, uLeft, uRight );
+      indicator = advection().jump( it, x, uLeft, uRight );
     }
 
     // misc for eoc calculation
@@ -495,8 +497,14 @@ namespace Fem
     }
 
   protected:
-    const AdvectionModelType& advection_;
-    const DiffusionModelType& diffusion_;
+    const AdvectionModelType& advection () const { return *advection_; }
+    const DiffusionModelType& diffusion () const { return *diffusion_; }
+
+    Dune::Fem::ThreadSafeValue< AdvectionModelType > advection_;
+    Dune::Fem::ThreadSafeValue< DiffusionModelType > diffusion_;
+
+    //const AdvectionModelType& advection_;
+    //const DiffusionModelType& diffusion_;
 
     const ProblemType& problem_;
     LimitedRangeType limitedRange_;