root/old_code/tumor_growth/tumor_growth.h

INCLUDED FROM

DEFINITIONS

1. evaluate
2. evaluate
3. jac_ran_to_jac_ran_lhs
4. jac_ran_to_jac_ran_rhs
5. ran_to_ran
6. ran_u_normal_rhs
7. Brusselator_model
8. ran_to_ran
9. ran_to_ran_ones
10. ran_to_ran_lhs
11. jac_ran_to_jac_ran_lhs
12. tri_ran_to_ran_lhs
13. generate_rhs
14. generate_rhs_old_surface
15. generate_rhs_new_surface

   1 /*! \file tumor_growth_operator.h
   2 
   3     \brief Tumor growth model (2012-Elliott+Styles) + our regularizing term
   4 
   5      Revision history:
   6 
   7           Revised by Christian Power dd.mm.yyyy
   8           Originally written by Christian Power
   9                (power22c@gmail.com) 22.11.2015
  10 
  11      This header provides model classes and operator classes to solve 
  12      a tumor growth model proposed by Elliott and Styles via the ESFEM.  
  13      Parameters: 
  14      - a,b,Dᶜ,γ ∈ R for the scalar equation with  γ ∼ vol(Γ).
  15      - ε,δ,α ∈ R for the surface equation, where ε and α are small 
  16        regularization parameter.
  17      PDE: 
  18      Find u:Γ→R (growth-promoting), w:Γ→R (growth-inhibiting) and
  19      X: Γ⁰×[0,T]→Rᵐ⁺¹ such that
  20                 ∂•u + u div(v) - Δu   = f₁(u,w),
  21                 ∂•w + w div(v) - DᶜΔw = f₂(u,w),
  22      where for f₁,f₂ we use the Brusselator model
  23                 f₁(u,w) = γ(a-u+u²w) ∧ f₂(u,w) = γ(b-u²w),
  24      and for the surface 
  25                 v - αΔv = (ε (-H) + δu) n = εΔX + δun,
  26                     ∂ᵗX = v(X).
  27      FEM discretization:
  28      Find u:I→Rᴺ (growth-promoting nodal values), w:I→Rᴺ (growth-inhibiting 
  29      nodal values) and X:I→R³ᴺ (surface nodal values) such that
  30                 (M(X) + α A(X))∂ᵗX = ε A(X) X + δ M(u,n)
  31                 ∂ᵗ(M(X)u) + A(X)u = γ(b M(X)1 + M(X;u,w)u)
  32                 ∂ᵗ(M(X)w) + Dᶜ A(X)w = γ(b M(X)1 - M(X;u,u)w)
  33 
  34      Pseudocode (Elliott+Styles discretization):
  35      - Given Xⁿ,uⁿ,wⁿ, solve for Xⁿ⁺¹
  36   
  37        (M₃ⁿ + (α+ετ)A₃ⁿ)Xⁿ⁺¹ = (M₃ⁿ+αA₃ⁿ)Xⁿ + τδ M₃ⁿ(uⁿ,nⁿ),
  38   
  39        where nⁿ is elementwise normal.  For the moment we choose an alternative
  40        normal namely nⁿ = A₃ⁿ id / ‖A₃ⁿ id‖.
  41      - Given Xⁿ⁺¹,uⁿ,wⁿ, solve for uⁿ⁺¹
  42   
  43        (1+τγ)(Mu)ⁿ⁺¹ + τ (Au)ⁿ⁺¹ - τγ Mⁿ⁺¹(uⁿ,wⁿ)uⁿ⁺¹ = (Mu)ⁿ + τγa Mⁿ⁺¹1
  44   
  45        where M(a,b) a 4 tensor is, namely Mⁱʲᵏˡ = ∫χⁱ̧χʲχᵏχˡ and 1 means the 
  46        constant 1 finite element function.
  47      - Given Xⁿ⁺¹,uⁿ⁺¹,wⁿ, solve for wⁿ⁺¹
  48   
  49        (Mw)ⁿ⁺¹ + τDᶜ (Aw)ⁿ⁺¹ + τγ Mⁿ⁺¹(uⁿ⁺¹,uⁿ⁺¹)wⁿ⁺¹ = (Mw)ⁿ + τγb Mⁿ⁺¹1
  50 
  51 
  52      Created by Christian Power on 22.11.2015
  53      Copyright (c) 2015 Christian Power.  All rights reserved.
  54  */
  55 
  56 #ifndef TUMOR_GROWTH_OPERATOR_H
  57 #define TUMOR_GROWTH_OPERATOR_H
  58 
  59 #include <cmath>
  60 #include <random>
  61 #include <functional>
  62 
  63 #include <dune/fem/function/common/function.hh>
  64 #include <dune/fem/solver/timeprovider.hh>
  65 #include <dune/fem/quadrature/quadrature.hh>
  66 
  67 #include <dune/common/fmatrix.hh>
  68 
  69 #include <dune/fem/quadrature/cachingquadrature.hh>
  70 #include <dune/fem/operator/common/operator.hh>
  71 
  72 #include <heat.hh>
  73 
  74 namespace Tumor_growth{
  75   // ----------------------------------------------------------------------
  76   // Initial data for the Brusselator model
  77 
  78   template<typename Function_space>
  79   struct Initial_data_u :
  80     Dune::Fem::Function<Function_space, Initial_data_u<Function_space> >{
  81 
  82     explicit Initial_data_u(const double mean_value, const double variance)
  83       : random_fun (std::bind(Random_dist {mean_value, variance},
  84                           Random_engine {}))
  85     {}
  86 
  87     using Domain = typename Function_space::DomainType;
  88     using Range = typename Function_space::RangeType;
  89 
  90     void evaluate(const Domain& p, Range& q) const{
  91       q = random_fun();
  92     }
  93   private:
  94     using Random_dist = std::normal_distribution<>;
  95     using Random_engine = std::default_random_engine;
  96     std::function<double()> random_fun;
  97   };
  98 
  99   template<typename Function_space>
 100   struct Initial_data_w :
 101     Dune::Fem::Function<Function_space, Initial_data_w<Function_space> >{
 102 
 103     explicit Initial_data_w(const double mean_value, const double variance)
 104       : base {mean_value}, pertubation {variance}
 105     {}
 106     
 107     using Domain = typename Function_space::DomainType;
 108     using Range = typename Function_space::RangeType;
 109 
 110     void evaluate(const Domain& p, Range& q) const{
 111       // const double  x = p[0], y = p[1], z = p[2];
 112       q = Range {base};
 113     }
 114   private:
 115     double base;
 116     double pertubation;
 117   };
 118   
 119   // ----------------------------------------------------------------------
 120   // Surface and Brusselator Models
 121 
 122   template <typename Vector_function_space, typename Scalar_function_space>
 123   struct Surface_evolution_model{
 124     using Vector = Vector_function_space;
 125     using Scalar = Scalar_function_space;
 126 
 127     template<typename T>
 128     using Domain = typename T::DomainType;
 129     template<typename T>
 130     using Range = typename T::RangeType;
 131     template<typename T>
 132     using JacobianRange = typename T::JacobianRangeType;
 133     template<typename T>
 134     using DomainField = typename T::DomainFieldType;
 135     template<typename T>
 136     using RangeField = typename T::RangeFieldType;
 137     static const int dim_vec_domain = Domain<Vector>::dimension;
 138     static const int dim_vec_range = Range<Vector>::dimension;
 139     static const int dim_scalar_domain = Domain<Scalar>::dimension;
 140     static const int dim_scalar_range = Range<Scalar>::dimension;
 141 
 142     explicit
 143     Surface_evolution_model (const Dune::Fem::TimeProviderBase& time_provider,
 144                              const double alpha_regularization_parameter,
 145                              const double epsilon_regularization_parameter,
 146                              const double delta_parameter) :
 147       tp {time_provider}, reg_alpha {alpha_regularization_parameter},
 148       reg_epsilon {epsilon_regularization_parameter}, delta {delta_parameter}
 149     {}
 150 
 151     // Pseudocode:
 152     // (M₃ⁿ + (α+ετ)A₃ⁿ)Xⁿ⁺¹ = (M₃ⁿ+αA₃ⁿ)Xⁿ + τδ M₃ⁿ(uⁿ,nⁿ)
 153     // alpha_regularization_parameter == α,
 154     // epsilon_regularization_parameter == ε,
 155     // delta_parameter == δ,
 156     // time_provider.deltaT() == τ.
 157 
 158     template <typename Entity, typename Point>
 159     void jac_ran_to_jac_ran_lhs(const Entity& entity, const Point& xx, 
 160                                 const JacobianRange<Vector>& dX_p,
 161                                 JacobianRange<Vector>& a_dX_p) const;
 162     // a_dX_p = (α+ετ)A₃ⁿdX_p
 163     template <typename Entity, typename Point>
 164     void jac_ran_to_jac_ran_rhs(const Entity& entity, const Point& xx, 
 165                                 const JacobianRange<Vector>& dX_p,
 166                                 JacobianRange<Vector>& a_dX_p) const;
 167     // a_dX_p = αA₃ⁿdX_p
 168     template <typename Entity, typename Point>
 169     void ran_to_ran(const Entity& entity, const Point& xx, 
 170                     const Range<Vector>& X_p, Range<Vector>& X_p_chi) const;
 171     // X_p_chi = M₃ⁿX_p
 172     template <typename Entity, typename Point>
 173     void ran_u_normal_rhs(const Entity& entity, const Point& xx, 
 174                           const Range<Scalar>& u_p,
 175                           Range<Vector>& n_p_chi) const;
 176     // n_p_chi = τδM₃ⁿ(u_p,nⁿ)
 177   private:
 178     const Dune::Fem::TimeProviderBase& tp;
 179     const double reg_alpha;
 180     const double reg_epsilon;
 181     const double delta;
 182   };
 183 
 184   enum class Growth {promoting, inhibiting};
 185 
 186   template <typename Scalar_function_space>
 187   class Brusselator_model{
 188   public:
 189     using Scalar = Scalar_function_space;
 190     using Range = typename Scalar::RangeType;
 191     using JacobianRange = typename Scalar::JacobianRangeType;
 192 
 193     explicit
 194     Brusselator_model(const double uniform_time_step,
 195                       const double gamma, 
 196                       const double a_or_b, 
 197                       const Growth scalar_type,
 198                       const double diffusion = 1.);
 199     // Pseudocode:
 200     // growth-promoting
 201     // (1+τγ)(Mu)ⁿ⁺¹ + τ (Au)ⁿ⁺¹ - τγ Mⁿ⁺¹(uⁿ,wⁿ)uⁿ⁺¹ = (Mu)ⁿ + τγa Mⁿ⁺¹1
 202     // growth-inhibiting
 203     // (Mw)ⁿ⁺¹ + τDᶜ (Aw)ⁿ⁺¹ + τγ Mⁿ⁺¹(uⁿ⁺¹,uⁿ⁺¹)wⁿ⁺¹ = (Mw)ⁿ + τγb Mⁿ⁺¹1
 204 
 205     template <typename Entity, typename Point>
 206     void ran_to_ran(const Entity& entity, const Point& xx, 
 207                     const Range& u_p, Range& u_p_chi) const;
 208     // u_p_chi = (Mu)ⁿ || u_p_chi = (Mw)ⁿ 
 209 
 210     template <typename Entity, typename Point>
 211     void ran_to_ran_ones(const Entity& entity, const Point& xx, 
 212                          Range& ones_chi) const;
 213     // ones_chi = τγa Mⁿ⁺¹1 || ones_chi = τγb Mⁿ⁺¹1 
 214 
 215     template <typename Entity, typename Point>
 216     void ran_to_ran_lhs(const Entity& entity, const Point& xx, 
 217                          const Range& u_p, Range& u_p_chi) const;
 218     // u_p_chi = (1+τγ)(Mu)ⁿ⁺¹ || u_p_chi = (Mw)ⁿ⁺¹ 
 219 
 220     template <typename Entity, typename Point>
 221     void jac_ran_to_jac_ran_lhs(const Entity& entity, const Point& xx, 
 222                                 const JacobianRange& du_p,
 223                                 JacobianRange& du_p_dchi) const;
 224     // du_p_dchi = τ(Au)ⁿ⁺¹ || du_p_dchi = τDᶜ(Aw)ⁿ⁺¹ 
 225 
 226     template <typename Entity, typename Point>
 227     void tri_ran_to_ran_lhs(const Entity& entity, const Point& xx,
 228                             const Range& first, const Range& second,
 229                             const Range& u_p, Range& uwu_p_chi) const;
 230     // uwu_p_chi = - τγ Mⁿ⁺¹(uⁿ,wⁿ)uⁿ⁺¹ || uwu_p_chi = τγ Mⁿ⁺¹(uⁿ⁺¹,uⁿ⁺¹)wⁿ⁺¹
 231 
 232   private:
 233     using Domain = typename Scalar::DomainType;
 234 
 235     // ----------
 236     // data members
 237     double mass_matrix_lhs {0.};
 238     double stiffness_matrix_lhs {0.};
 239     double quadrilinear_mass_matrix_lhs {0.};
 240     double mass_matrix_rhs {0.};
 241   };
 242 
 243   // ----------------------------------------------------------------------
 244   // Parabolic operators
 245 
 246   template <typename Vec_discrete_function, typename Scalar_discrete_function, 
 247             typename Model>
 248   /* requires Model:
 249      template <typename Entity, typename Point>
 250      void jac_ran_to_jac_ran_lhs(const Entity& entity, const Point& xx, 
 251      const JacobianRange<Vector>& dX_p,
 252      JacobianRange<Vector>& a_dX_p) const;
 253      template <typename Entity, typename Point>
 254      void jac_ran_to_jac_ran_rhs(const Entity& entity, const Point& xx, 
 255      const JacobianRange<Vector>& dX_p,
 256      JacobianRange<Vector>& a_dX_p) const;
 257      template <typename Entity, typename Point>
 258      void ran_to_ran(const Entity& entity, const Point& xx, 
 259      const Range<Vector>& X_p, Range<Vector>& X_p_chi) const;
 260      template <typename Entity, typename Point>
 261      void ran_u_normal_rhs(const Entity& entity, const Point& xx, 
 262      const Range<Scalar>& u_p,
 263      Range<Vector>& n_p_chi) const;
 264   */
 265   struct Surface_parabolic_operator
 266     : Dune::Fem::Operator<Vec_discrete_function>{
 267     using Vector = Vec_discrete_function;
 268     using Scalar = Scalar_discrete_function;
 269   protected:
 270     using Discrete_function_space = typename Vector::DiscreteFunctionSpaceType;
 271     using Iterator = typename Discrete_function_space::IteratorType;
 272     using Entity = typename Iterator::Entity;
 273     using Geometry = typename Entity::Geometry;
 274     using GridPart = typename Discrete_function_space::GridPartType;
 275     using QuadratureType = Dune::Fem::CachingQuadrature<GridPart, 0>;
 276 
 277     template<typename T>
 278     using Local_function = typename T::LocalFunctionType;
 279     template<typename T>
 280     using Range = typename Local_function<T>::RangeType;
 281     template<typename T>
 282     using Jacobian_range = typename Local_function<T>::JacobianRangeType;
 283 
 284   public:
 285     explicit
 286     Surface_parabolic_operator(const Model& model_input,
 287                                const Scalar_discrete_function& u_previous,
 288                                const double bdf_alpha_lead_coeff = 1.)
 289       : model (model_input), u_approx (u_previous), bdf_factor {bdf_alpha_lead_coeff}
 290     {}
 291   
 292     virtual void operator() (const Vec_discrete_function& u, 
 293                              Vec_discrete_function& w) const;
 294     // Pseudocode:
 295     // w = (M₃ⁿ + (α+ετ)A₃ⁿ)u
 296 
 297     void generate_rhs(const Vec_discrete_function& u, Vec_discrete_function& w) const;
 298     // Pseudocode:
 299     // w = (M₃ⁿ+αA₃ⁿ)u + τδ M₃ⁿ(u_previous,nⁿ)
 300 
 301   private:
 302     const Model model;
 303     const Scalar_discrete_function& u_approx;
 304     const double bdf_factor;
 305   };
 306 
 307   template <typename Scalar_discrete_function, typename Model>
 308   /* requires Model:
 309      template <typename Entity, typename Point>
 310      void ran_to_ran(const Entity& entity, const Point& xx, 
 311                      const Range& u_p, Range& u_p_chi) const;
 312      template <typename Entity, typename Point>
 313      void ran_to_ran_ones(const Entity& entity, const Point& xx, 
 314                           Range& ones_chi) const;
 315      template <typename Entity, typename Point>
 316      void ran_to_ran_lhs(const Entity& entity, const Point& xx, 
 317                          const Range& u_p, Range& u_p_chi) const;
 318      template <typename Entity, typename Point>
 319      void jac_ran_to_jac_ran_lhs(const Entity& entity, const Point& xx, 
 320                                  const JacobianRange& du_p,
 321                                  JacobianRange& du_p_dchi) const;
 322      template <typename Entity, typename Point>
 323      void tri_ran_to_ran_lhs(const Entity& entity, const Point& xx,
 324                              const Range& first, const Range& second,
 325                              const Range& u_p, Range& uwu_p_chi) con st;
 326   */
 327   struct Scalar_parabolic_operator
 328     : Dune::Fem::Operator<Scalar_discrete_function>{
 329   public:
 330     using FE_function = Scalar_discrete_function;
 331 
 332     explicit
 333     Scalar_parabolic_operator(const Model& model_input,
 334                               const FE_function& u_previous,
 335                               const FE_function& u_or_w_previous,
 336                               const double bdf_alpha_lead_coeff = 1.)
 337       : model (model_input), u_approx (u_previous), u_or_w_approx (u_or_w_previous),
 338         bdf_factor {bdf_alpha_lead_coeff}
 339     {}
 340 
 341     virtual void operator() (const FE_function& u, FE_function& w) const;
 342     // Pseudocode:
 343     // w = (1+τγ)(Mu)ⁿ⁺¹ + τ (Au)ⁿ⁺¹ - τγ Mⁿ⁺¹(uⁿ,wⁿ)uⁿ⁺¹  ||
 344     // w = (Mw)ⁿ⁺¹ + τDᶜ (Aw)ⁿ⁺¹ + τγ Mⁿ⁺¹(uⁿ⁺¹,uⁿ⁺¹)wⁿ⁺¹
 345 
 346     void generate_rhs_old_surface(const FE_function& u, FE_function& w) const;
 347     // Pseudocode: w = (Mu)ⁿ || w = (Mw)ⁿ
 348 
 349     void generate_rhs_new_surface(FE_function& w) const;
 350     // Pseudocode: w = τγa Mⁿ⁺¹1 || w = τγb Mⁿ⁺¹1
 351   private:
 352     using Discrete_function_space = typename FE_function::DiscreteFunctionSpaceType;
 353     using Iterator = typename Discrete_function_space::IteratorType;
 354     using Entity = typename Iterator::Entity;
 355     using Geometry = typename Entity::Geometry;
 356     using GridPart = typename Discrete_function_space::GridPartType;
 357     using QuadratureType = Dune::Fem::CachingQuadrature<GridPart, 0>;
 358 
 359     using Local_function = typename FE_function::LocalFunctionType;
 360     using Range = typename Local_function::RangeType;
 361     using Jacobian_range = typename Local_function::JacobianRangeType;
 362     
 363     // ----------
 364     // Data members
 365     const Model model;
 366     const FE_function& u_approx;
 367     const FE_function& u_or_w_approx;
 368     const double bdf_factor;
 369   };
 370 
 371   // ======================================================================
 372   // All implementation details
 373 
 374   // ----------------------------------------------------------------------
 375   // Implementation details for 'Surface_evolution_model'
 376   template <typename Vector_function_space, typename Scalar_function_space>
 377   template <typename Entity, typename Point>
 378   inline void Surface_evolution_model<Vector_function_space, Scalar_function_space>::
 379   jac_ran_to_jac_ran_lhs(const Entity& entity, const Point& xx, 
 380                          const JacobianRange<Vector>& dX_p,
 381                          JacobianRange<Vector>& a_dX_p) const{
 382     // DomainType xGlobal = entity.geometry().global(Dune::coordinate(xx));
 383     // xGlobal.size() == dim_domain == 3;
 384     // DomainFieldType x = xGlobal[0], y = xGlobal[1], z = xGlobal[2];
 385     // double t = tp.time(), dT = tp.deltaT();
 386     // 'dX_p' is a R²ˣ³ Matrix, so it has two indices.               
 387     a_dX_p = dX_p;
 388     a_dX_p *= (reg_alpha + reg_epsilon * tp.deltaT());
 389   }
 390   template <typename Vector_function_space, typename Scalar_function_space>
 391   template <typename Entity, typename Point>
 392   inline void Surface_evolution_model<Vector_function_space, Scalar_function_space>::
 393   jac_ran_to_jac_ran_rhs(const Entity& entity, const Point& xx, 
 394                          const JacobianRange<Vector>& dX_p,
 395                          JacobianRange<Vector>& a_dX_p) const{
 396     // DomainType xGlobal = entity.geometry().global(Dune::coordinate(xx));
 397     // xGlobal.size() == dim_domain == 3;
 398     // DomainFieldType x = xGlobal[0], y = xGlobal[1], z = xGlobal[2];
 399     // double t = tp.time(), dT = tp.deltaT();
 400     // 'dX_p' is a R²ˣ³ Matrix, so it has two indices.               
 401     a_dX_p = dX_p;
 402     a_dX_p *= reg_alpha;
 403   }
 404   template <typename Vector_function_space, typename Scalar_function_space>
 405   template <typename Entity, typename Point>
 406   inline void Surface_evolution_model<Vector_function_space, Scalar_function_space>::
 407   ran_to_ran(const Entity& entity, const Point& xx, 
 408              const Range<Vector>& X_p, Range<Vector>& X_p_chi) const{
 409     // DomainType xGlobal = entity.geometry().global(Dune::coordinate(xx));
 410     // xGlobal.size() == dim_domain == 3;
 411     // DomainFieldType x = xGlobal[0], y = xGlobal[1], z = xGlobal[2];
 412     // double t = tp.time(), dT = tp.deltaT();
 413     X_p_chi = X_p;
 414   }
 415   template <typename Vector_function_space, typename Scalar_function_space>
 416   template <typename Entity, typename Point>
 417   void Surface_evolution_model<Vector_function_space, Scalar_function_space>::
 418   ran_u_normal_rhs(const Entity& entity, const Point& xx, const Range<Scalar>& u_p,
 419                    Range<Vector>& n_p_chi) const{
 420     // Domain<Vector> xGlobal = entity.geometry().global(Dune::coordinate(xx));
 421     // xGlobal.size() == dim_domain == 3;
 422     // DomainFieldType x = xGlobal[0], y = xGlobal[1], z = xGlobal[2];
 423     // double t = tp.time(), dT = tp.deltaT();
 424 
 425     // Consistency check.  Comment this later out.
 426     // if(entity.geometry().corners() != 3) throw std::runtime_error
 427     //   {"Error in Surface_evolution_model::ran_u_normal_rhs().\n"
 428     //    "Number of corners unequal to 3."};
 429     Domain<Vector> p0 = entity.geometry().corner(0);
 430     Domain<Vector> p1 = entity.geometry().corner(1);
 431     Domain<Vector> p2 = entity.geometry().corner(2);
 432     
 433     // Calculating normal vector
 434 
 435     // Tangent vector components
 436     Domain<Vector> v = p2 - p0;
 437     Domain<Vector> w = p1 - p0;
 438 
 439     auto eucl_norm = [](const Range<Vector>& vv){
 440       return std::sqrt(vv[0] * vv[0] + vv[1] * vv[1] + vv[2] * vv[2]);
 441     };
 442 
 443     // Formula for cross product
 444     // n = v × w / ‖v × w‖
 445     n_p_chi[0] = v[1] * w[2] - v[2] * w[1];
 446     n_p_chi[1] = - v[0] * w[2] + v[2] * w[0];
 447     n_p_chi[2] = v[0] * w[1] - v[1] * w[0];
 448     n_p_chi /= eucl_norm(n_p_chi);
 449 
 450     // The final result
 451     n_p_chi *= u_p * tp.deltaT() * delta; 
 452   }
 453 
 454   // ----------------------------------------------------------------------
 455   // Implementation details for 'Brusselator_model'
 456 
 457   template <typename Scalar_function_space>
 458   Brusselator_model<Scalar_function_space>::
 459   Brusselator_model(const double uniform_time_step,
 460                     const double gamma, const double a_or_b, const Growth scalar_type,
 461                     const double diffusion)
 462   {
 463     switch(scalar_type){
 464     case Growth::promoting:
 465       mass_matrix_lhs = 1. + uniform_time_step * gamma;
 466       stiffness_matrix_lhs = uniform_time_step;
 467       quadrilinear_mass_matrix_lhs = (-1.) * uniform_time_step * gamma;
 468       break;
 469     case Growth::inhibiting:
 470       mass_matrix_lhs = 1.;
 471       stiffness_matrix_lhs = uniform_time_step * diffusion;
 472       quadrilinear_mass_matrix_lhs = uniform_time_step * gamma;
 473       break;
 474     default:
 475       throw std::runtime_error {"Error in constructor of Brusselator_model.\n"
 476           "Your scalar_type is not implemented."};
 477       break;
 478     };
 479     mass_matrix_rhs = uniform_time_step * gamma * a_or_b;
 480   }
 481 
 482   template <typename Scalar_function_space>
 483   template <typename Entity, typename Point>
 484   inline void Brusselator_model<Scalar_function_space>::
 485   ran_to_ran(const Entity& entity, const Point& xx, 
 486              const Range& u_p, Range& u_p_chi) const{
 487     u_p_chi = u_p;
 488   }
 489   template <typename Scalar_function_space>
 490   template <typename Entity, typename Point>
 491   inline void Brusselator_model<Scalar_function_space>::
 492   ran_to_ran_ones(const Entity& entity, const Point& xx, 
 493                   Range& ones_chi) const{
 494     // ones_chi = τγa Mⁿ⁺¹1 || ones_chi = τγb Mⁿ⁺¹1     
 495     ones_chi = Range {mass_matrix_rhs};
 496   }
 497   template <typename Scalar_function_space>
 498   template <typename Entity, typename Point>
 499   inline void Brusselator_model<Scalar_function_space>::
 500   ran_to_ran_lhs(const Entity& entity, const Point& xx, 
 501                  const Range& u_p, Range& u_p_chi) const{
 502     // u_p_chi = (1+τγ)(Mu)ⁿ⁺¹ || u_p_chi = (Mw)ⁿ⁺¹
 503     u_p_chi = u_p;
 504     u_p_chi *= mass_matrix_lhs;
 505   }
 506   template <typename Scalar_function_space>
 507   template <typename Entity, typename Point>
 508   inline void Brusselator_model<Scalar_function_space>::
 509   jac_ran_to_jac_ran_lhs(const Entity& entity, const Point& xx, 
 510                          const JacobianRange& du_p,
 511                          JacobianRange& du_p_dchi) const{
 512     // du_p_dchi = τ(Au)ⁿ⁺¹ || du_p_dchi = τDᶜ(Aw)ⁿ⁺¹
 513     du_p_dchi = du_p;
 514     du_p_dchi *= stiffness_matrix_lhs;
 515   }
 516   template <typename Scalar_function_space>
 517   template <typename Entity, typename Point>
 518   inline void Brusselator_model<Scalar_function_space>::
 519   tri_ran_to_ran_lhs(const Entity& entity, const Point& xx,
 520                      const Range& first, const Range& second,
 521                      const Range& u_p, Range& uwu_p_chi) const{
 522     // uwu_p_chi = - τγ Mⁿ⁺¹(uⁿ,wⁿ)uⁿ⁺¹ || uwu_p_chi = τγ Mⁿ⁺¹(uⁿ⁺¹,uⁿ⁺¹)wⁿ⁺¹
 523     uwu_p_chi = u_p;
 524     uwu_p_chi *= first * second * quadrilinear_mass_matrix_lhs;
 525   }
 526   
 527   // ----------------------------------------------------------------------
 528   // Implementation details for 'Surface_parabolic_operator'
 529 
 530   template <typename Vec_discrete_function, typename Scalar_discrete_function, 
 531             typename Model>
 532   void Surface_parabolic_operator<Vec_discrete_function, 
 533                                   Scalar_discrete_function, Model>::
 534   operator() (const Vec_discrete_function& u, Vec_discrete_function& w) const{
 535     w.clear();
 536 
 537     const Discrete_function_space& dfSpace = w.space();
 538     for(const Entity& entity : dfSpace){
 539       // auto == tmp class Iterator<Vector>
 540       // auto& == tmp class Entity<Vector> convertible to Entity<Vector>
 541       // The constructor of 'Entity<Vector>' is protected, hence auto&
 542       // Remark: Entity<Vector> == Entity<Scalar>!
 543   
 544       const Geometry& geometry = entity.geometry();
 545       
 546       const Local_function<Vector> uLocal = u.localFunction(entity);
 547       Local_function<Vector> wLocal = w.localFunction(entity);
 548     
 549       QuadratureType quadrature(entity, uLocal.order() + wLocal.order() + 1);
 550     
 551       for(size_t pt = 0; pt < quadrature.nop(); ++pt){
 552     
 553         Range<Vector> u_x;      // u_x == X(x);
 554         uLocal.evaluate(quadrature[pt], u_x);   // initialize
 555     
 556         Jacobian_range<Vector> du_x;    // du_x == ∇X(x) (jacobian matrix)
 557         uLocal.jacobian(quadrature[pt], du_x);  // initialize
 558     
 559         Range<Vector> u_x_psi;
 560         // u_x_psi ≈ ∫ X Ψ
 561         model.ran_to_ran(entity, quadrature[pt], u_x, u_x_psi); 
 562         // initialize
 563         // u_x_chi *= bdf_factor;
 564   
 565         Jacobian_range<Vector> du_x_dpsi;
 566         // du_x_dpsi ≈ (α+ετ)∫ ∇X · ∇Ψ  (∇ stands for gradient matrix)
 567         model.jac_ran_to_jac_ran_lhs(entity, quadrature[pt], du_x, du_x_dpsi);
 568         // initialize 
 569     
 570         const typename QuadratureType::CoordinateType& x = quadrature.point(pt);
 571         const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 572         // weigth = ω_pt (det(dψᵀ·dψ))^½ 
 573                 
 574         u_x_psi *= weight;
 575         du_x_dpsi *= weight;
 576         wLocal.axpy(quadrature[pt], u_x_psi, du_x_dpsi);
 577         // calculating matrix * vector without assembling the matrices
 578       }
 579     }
 580     w.communicate();
 581   }
 582   template <typename Vec_discrete_function, typename Scalar_discrete_function, 
 583             typename Model>
 584   void Surface_parabolic_operator<Vec_discrete_function, 
 585                                   Scalar_discrete_function, Model>::
 586   generate_rhs(const Vec_discrete_function& u, Vec_discrete_function& w) const{
 587     w.clear();
 588 
 589     const Discrete_function_space& dfSpace = w.space();
 590     for(const Entity& entity : dfSpace){ 
 591       // auto& == tmp class Entity<Vector> convertible to Entity<vector>&
 592       // Entity<Vector> == Entity<Scalar>!
 593 
 594       const Geometry& geometry = entity.geometry();
 595 
 596       const Local_function<Vector> uLocal = u.localFunction(entity);
 597       const Local_function<Scalar> u_approx_local 
 598         = u_approx.localFunction(entity);
 599       Local_function<Vector> wLocal = w.localFunction(entity);
 600 
 601       QuadratureType quadrature(entity, uLocal.order() + wLocal.order() + 1);
 602 
 603       for(size_t pt = 0; pt < quadrature.nop(); ++pt){
 604         Range<Vector> u_x;      // u_x == X(x);
 605         uLocal.evaluate(quadrature[pt], u_x);   // initialize
 606 
 607         Range<Scalar> u_approx_x;       // u_approx_x = u(x);
 608         u_approx_local.evaluate(quadrature[pt], u_approx_x);
 609 
 610         Jacobian_range<Vector> du_x;    // du_x == ∇X(x) (jacobian matrix)
 611         uLocal.jacobian(quadrature[pt], du_x);  // initialize
 612 
 613         Range<Vector> u_x_psi;
 614         // u_x_psi ≈ ∫ X Ψ
 615         model.ran_to_ran(entity, quadrature[pt], u_x, u_x_psi); // initialize
 616         // u_x_chi *= bdf_factor;
 617 
 618         Jacobian_range<Vector> du_x_dpsi;
 619         // du_x_dpsi ≈ α ∫ ∇X · ∇Ψ  (∇ stands for gradient matrix)
 620         model.jac_ran_to_jac_ran_rhs(entity, quadrature[pt], du_x, du_x_dpsi);
 621         // initialize 
 622 
 623         Range<Vector> un_x_psi;
 624         // un_x_psi ≈ τδ ∫ u_approx n Ψ
 625         model.ran_u_normal_rhs(entity, quadrature[pt], u_approx_x, un_x_psi);
 626 
 627         const typename QuadratureType::CoordinateType& x = quadrature.point(pt);
 628         const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 629         // weigth = ω_pt (det(dψᵀ·dψ))^½ 
 630                 
 631         u_x_psi *= weight;
 632         un_x_psi *= weight;
 633         du_x_dpsi *= weight;
 634 
 635         wLocal.axpy(quadrature[pt], u_x_psi, du_x_dpsi);
 636         wLocal.axpy(quadrature[pt], un_x_psi, 0);
 637         // calculating matrix * vector without assembling the matrices
 638       }
 639     }
 640     w.communicate();
 641   }
 642 
 643   // ----------------------------------------------------------------------
 644   // Implementation details for 'Scalar_parabolic_operator'
 645 
 646   template <typename Scalar_discrete_function, typename Model>
 647   void Scalar_parabolic_operator<Scalar_discrete_function, Model>::
 648   operator() (const FE_function& u, FE_function& w) const{
 649     w.clear();
 650 
 651     const Discrete_function_space& dfSpace = w.space();
 652     for(const Entity& entity : dfSpace){  
 653       const Geometry& geometry = entity.geometry();
 654 
 655       const Local_function uLocal = u.localFunction(entity);
 656       const auto u_approx_local = u_approx.localFunction(entity);
 657       const auto u_or_w_approx_local = u_or_w_approx.localFunction(entity);
 658       Local_function wLocal = w.localFunction(entity);
 659 
 660       QuadratureType quadrature(entity, uLocal.order() + wLocal.order() + 1);
 661       for(size_t pt = 0; pt < quadrature.nop(); ++pt){
 662         Range u_x;      // u_x == u(x);
 663         uLocal.evaluate(quadrature[pt], u_x);   // initialize
 664     
 665         Jacobian_range du_x;    // du_x == ∇u(x) (jacobian matrix)
 666         uLocal.jacobian(quadrature[pt], du_x);  // initialize
 667     
 668         Range u_x_chi;
 669         // u_x_chi ≈ (1+τγ) ∫ u χ || u_x_chi ≈ ∫ u χ
 670         model.ran_to_ran_lhs(entity, quadrature[pt], u_x, u_x_chi);     
 671         // initialize
 672         // u_x_chi *= bdf_factor;
 673   
 674         Jacobian_range du_x_dchi;
 675         // du_x_dchi ≈ τ ∫ ∇u · ∇χ  (∇ stands for gradient matrix)
 676         // || du_x_dchi ≈ τDᶜ ∫ ∇u · ∇χ
 677         model.jac_ran_to_jac_ran_lhs(entity, quadrature[pt], du_x, du_x_dchi);
 678         // initialize 
 679     
 680         const typename QuadratureType::CoordinateType& x = quadrature.point(pt);
 681         const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 682         // weigth = ω_pt (det(dψᵀ·dψ))^½ 
 683                 
 684         u_x_chi *= weight;
 685         du_x_dchi *= weight;
 686         wLocal.axpy(quadrature[pt], u_x_chi, du_x_dchi);
 687         // calculating matrix * vector without assembling the matrices
 688       }
 689 
 690       QuadratureType 
 691         bigger_quadrature(entity, uLocal.order() + u_approx_local.order() 
 692                           + u_or_w_approx_local.order() + wLocal.order() + 1);
 693       for(size_t pt = 0; pt < bigger_quadrature.nop(); ++pt){
 694         Range u_x;      // u_x == u(x);
 695         uLocal.evaluate(bigger_quadrature[pt], u_x);    // initialize
 696         Range u_app;    // u_app == u_approx(x);
 697         u_approx_local.evaluate(bigger_quadrature[pt], u_app);  // initialize
 698         Range u_or_w_app;       // u_or_w_app == u_or_w_approx(x);
 699         u_or_w_approx_local.evaluate(bigger_quadrature[pt], u_or_w_app);
 700         // initialize
 701 
 702         Range uwu_x_chi;        
 703         // uwu_x_chi ≈ -τγ ∫ uⁿ wⁿ uⁿ⁺¹ χ || uwu_x_chi ≈ τγ ∫ uⁿ⁺¹uⁿ⁺¹wⁿ⁺¹χ
 704         model.tri_ran_to_ran_lhs(entity, bigger_quadrature[pt],
 705                                  u_app, u_or_w_app, u_x, uwu_x_chi);
 706 
 707         const typename QuadratureType::CoordinateType& x = bigger_quadrature.point(pt);
 708         const double weight
 709           = bigger_quadrature.weight(pt) * geometry.integrationElement(x);
 710         // weigth = ω_pt (det(dψᵀ·dψ))^½ 
 711 
 712         wLocal.axpy(bigger_quadrature[pt], uwu_x_chi, 0);
 713       }
 714     }
 715     w.communicate();
 716   }  
 717   template <typename Scalar_discrete_function, typename Model>
 718   void Scalar_parabolic_operator<Scalar_discrete_function, Model>::
 719   generate_rhs_old_surface (const FE_function& u, FE_function& w) const{
 720     w.clear();
 721   
 722     const Discrete_function_space& dfSpace = w.space();
 723     for(const Entity& entity : dfSpace){
 724   
 725       const Geometry& geometry = entity.geometry();
 726       
 727       const Local_function uLocal = u.localFunction(entity);
 728       Local_function wLocal = w.localFunction(entity);
 729     
 730       QuadratureType quadrature(entity, uLocal.order() + wLocal.order() + 1);
 731     
 732       for(size_t pt = 0; pt < quadrature.nop(); ++pt){
 733     
 734         Range u_x;      // u_x == X(x);
 735         uLocal.evaluate(quadrature[pt], u_x);   // initialize
 736     
 737         Range u_x_chi;
 738         // u_x_chi ≈ ∫ X χ
 739         model.ran_to_ran(entity, quadrature[pt], u_x, u_x_chi); 
 740         // initialize
 741   
 742         const typename QuadratureType::CoordinateType& x = quadrature.point(pt);
 743         const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 744         // weigth = ω_pt (det(dψᵀ·dψ))^½ 
 745                 
 746         u_x_chi *= weight;
 747         wLocal.axpy(quadrature[pt], u_x_chi, 0);
 748         // calculating matrix * vector without assembling the matrices
 749       }
 750     }
 751     w.communicate();
 752   }
 753   template <typename Scalar_discrete_function, typename Model>
 754   void Scalar_parabolic_operator<Scalar_discrete_function, Model>::
 755   generate_rhs_new_surface (FE_function& w) const{
 756     w.clear();
 757   
 758     const Discrete_function_space& dfSpace = w.space();
 759     for(const Entity& entity : dfSpace){
 760       const Geometry& geometry = entity.geometry();      
 761 
 762       Local_function wLocal = w.localFunction(entity);
 763    
 764       QuadratureType quadrature(entity, wLocal.order() + 1);
 765     
 766       for(size_t pt = 0; pt < quadrature.nop(); ++pt){
 767 
 768         Range ones_chi;
 769         // ones_chi ≈ τγa ∫ 1 χ || ones_chi ≈ τγb ∫ 1 χ
 770         model.ran_to_ran_ones(entity, quadrature[pt], ones_chi);        
 771         // initialize
 772 
 773         const typename QuadratureType::CoordinateType& x = quadrature.point(pt);
 774         const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 775         // weigth = ω_pt (det(dψᵀ·dψ))^½ 
 776                 
 777         ones_chi *= weight;
 778         wLocal.axpy(quadrature[pt], ones_chi, 0);
 779         // calculating matrix * vector without assembling the matrices
 780       }
 781     }
 782     w.communicate();
 783   }
 784 }       // namespace Tumor_growth
 785 
 786 #endif // TUMOR_GROWTH_OPERATOR_H
 787 
 788 /*! Log:
 789     2015.11.27: Save a promising tumor_growth_code()
 790 
 791                 I have run the experiment with u_base_value = 1, w_base_value = 0, 
 792                 no perturbation values
 793 
 794     2015.11.26: Crush a bug in the Scalar_parabolic_operator::operator() function
 795     
 796                 I received a assertion error from an dune array class.  The reason was 
 797                 that a bigger_quadrature object had to replace a quadrature object
 798                 in the second for-loop.  
 799 
 800     2015.11.25: Add Brusselator_operator class without testing it
 801     2015.11.25: Add Brusselator_model class without testing it
 802     2015.11.25: Corrected mistake in typedef/using declaration of Surface_operator
 803     2015.11.24: Give operator() and generate_rhs()
 804     2015.11.21: Save tested version of the discrete normal vector
 805  */