root/src/secOrd_op_rhs_impl.cpp

DEFINITIONS

1. dassert
2. update_cache
3. evaluate
4. addScaled_to
5. addScaled_to
6. addScaled_to

   1 /*! \file secOrd_op_rhs_impl.cpp
   2     \brief Implementation for `secOrd_op_rhs_impl.h`
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power April 2016
   8           Originally written by Christian Power
   9                (power22c@gmail.com) April 2016
  10 
  11      Idea
  12      --------------------------------------------------
  13 
  14      Insert the sage-syntax solution of the right-hand side here.
  15 
  16      \author Christian Power 
  17      \date 26. April 2016
  18      \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  19  */
  20 
  21 #include <cmath>
  22 #include <numeric>
  23 #include "config.h"
  24 #include <dune/fem/io/parameter.hh>
  25 #include <dassert.h>
  26 #include "esfem_error.h"
  27 #include "secOrd_op_rhs_impl.h"
  28 #include "esfem_error.h"
  29 
  30 // // Scalar valued dune finite element function
  31 // using Scal_fef = Esfem::Grid::Scal_FEfun::Dune_FEfun;
  32 // // Vector valued dune finite element function
  33 // using Vec_fef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
  34 // // Geometry for an element
  35 // using Geometry
  36 // = Scal_fef::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
  37 // // Technicality for quadrature 
  38 // using Grid_part = Scal_fef::GridPartType;
  39 // // Quadrature on a straight element
  40 // using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
  41 using Esfem::Impl::Rhs_fun;
  42 using Esfem::Impl::Vec_rhs_fun;
  43 using Esfem::Impl::sls_rhs;
  44 using Esfem::Impl::sd_rhs;
  45 using Esfem::Impl::sdp_u_rhs;
  46 using Dune::Fem::Parameter;
  47 using namespace std;
  48 
  49 // ----------------------------------------------------------------------
  50 // Implementation of Rhs_fun
  51 
  52 Rhs_fun::Rhs_fun(const Dune::Fem::TimeProviderBase& tpb, const Growth type)
  53   :tp {tpb}
  54 {
  55   const double rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)};
  56   const double rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)};
  57   const double k {Parameter::getValue<double>("logistic_growth.steepness", .5)};
  58   const double Dc {Parameter::getValue<double>("tumor_growth.heat.Dc", 10.)};
  59   const double ep {Parameter::getValue<double>("tumor_growth.surface.epsilon", .01)};
  60   const double al {Parameter::getValue<double>("tumor_growth.surface.alpha", 1e-3)};
  61   const double delta {Parameter::getValue<double>("tumor_growth.surface.delta", .4)};
  62 
  63   dassert(rE > rA, Assert::compose(__FILE__, __LINE__, "r_end <= r_start"));
  64   dassert(rA > 0, Assert::compose(__FILE__, __LINE__, "r_start < 0"));
  65   dassert(k > 0, Assert::compose(__FILE__, __LINE__, "Steepness non-positive"));
  66   // Other parameter are tested in Io::Parameter
  67   
  68   switch(type){
  69   case Growth::promoting:
  70     fun_impl = [&tp = tp, rA, k, rE]
  71       (const Domain& d, Range& r){
  72       const double x = d[0];
  73       const double y = d[1];
  74       const double z = d[2];
  75       const double t = tp.time();
  76       r = 2*k*x*y*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1)*exp(-6*t) - (k*(pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1) + k*(pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1) + k*(pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1))*x*y*exp(-6*t) + 4*pow(x,3)*y*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),3) + 4*x*pow(y,3)*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),3) - 6*x*y*exp(-6*t) + 2*(pow(x,3)*y*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) + (pow(x,3)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x/(pow(x,2) + pow(y,2) + pow(z,2)))*y*exp(-6*t) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1) + (4*pow(x,2)*pow(y,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - pow(x,2)*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)) - (pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*exp(-6*t))*x*y/(pow(x,2) + pow(y,2) + pow(z,2)) + (4*pow(x,2)*pow(y,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - pow(y,2)*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)) - (pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*exp(-6*t))*x*y/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(x*pow(y,3)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) + x*(pow(y,3)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - y/(pow(x,2) + pow(y,2) + pow(z,2)))*exp(-6*t) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1) + 2*(2*pow(x,2)*y*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - y*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*x*z/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(2*x*pow(y,2)*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*y*z/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(2*x*y*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1);
  77     };
  78     break;
  79   case Growth::inhibiting:
  80     fun_impl = [&tp = tp, rA, k, rE, al, Dc]
  81       (const Domain& d, Range& r){
  82       const double x = d[0];
  83       const double y = d[1];
  84       const double z = d[2];
  85       const double t = tp.time();
  86       r = 2*k*y*z*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1)*exp(-6*t) - (k*(pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1) + k*(pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1) + k*(pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1))*y*z*exp(-6*t) - 6*y*z*exp(-6*t) + (4*pow(x,2)*pow(y,3)*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),3) + 4*pow(x,2)*y*pow(z,3)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),3) + 2*(2*pow(x,2)*y*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - y*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1) + 2*(2*x*pow(y,2)*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*x*y/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(pow(y,3)*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) + (pow(y,3)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - y/(pow(x,2) + pow(y,2) + pow(z,2)))*z*exp(-6*t) - y*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1) + 2*(2*x*y*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*x*z/(pow(x,2) + pow(y,2) + pow(z,2)) + (4*pow(y,2)*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - pow(y,2)*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)) - (pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*exp(-6*t))*y*z/(pow(x,2) + pow(y,2) + pow(z,2)) + (4*pow(y,2)*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - pow(z,2)*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)) - (pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*exp(-6*t))*y*z/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(y*pow(z,3)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) + y*(pow(z,3)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - z/(pow(x,2) + pow(y,2) + pow(z,2)))*exp(-6*t) - y*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1))*Dc;
  87     };
  88     break;
  89   default:
  90     throw Rhs_error {Assert::compose(__FILE__, __LINE__, "Rhs_fun()")};
  91     break;
  92   };
  93 }
  94 
  95 void Rhs_fun::dassert(const bool assertion, const std::string& msg){
  96   Assert::dynamic<Assert::level(Assert::default_level), Esfem::Parameter_error>
  97     (assertion, msg);
  98 }
  99 
 100 // ----------------------------------------------------------------------
 101 // Implementation of Vec_rhs_fun
 102 
 103 Vec_rhs_fun::Vec_rhs_fun(const Dune::Fem::TimeProviderBase& tpb)
 104   : tp {tpb},
 105   alpha {Parameter::getValue<double>("tumor_growth.surface.alpha", 1e-3)},
 106   epsilon {Parameter::getValue<double>("tumor_growth.surface.epsilon", .01)},
 107   r_start {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
 108   r_end {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
 109   k {Parameter::getValue<double>("logistic_growth.steepness", .5)},
 110   delta {Parameter::getValue<double>("tumor_growth.surface.delta", .4)}
 111 { 
 112   // no checking has to be done, since other classes already do this
 113   // cache[0] = tp.time() - 1; // cache[0] != tp.time()
 114   // update_cache();
 115 }
 116 
 117 void Vec_rhs_fun::update_cache() const{
 118   if(cache[0] == tp.time()) return;
 119   const double t = tp.time();
 120   const double e_kt = exp(-k * t);
 121   const double r_t = r_end * r_start / (r_end * e_kt + r_start * (1 - e_kt) );
 122   cache[0] = t;
 123   cache[1] = k * ( 1 - r_t / r_end);
 124   cache[2] = 2 * ( alpha * cache[1] + epsilon);
 125   cache[3] = delta * exp(-6 * t);
 126 }
 127 
 128 void Vec_rhs_fun::evaluate(const Domain& d, Range& r) const{
 129   r = 0;
 130   /*
 131   const double x = d[0];
 132   const double y = d[1];
 133   const double z = d[2];
 134   const double norm_square = x*x + y*y + z*z;
 135   const double factor = -(1. + (alpha + epsilon) * 2. / norm_square);
 136   r[0] = factor * d[0];
 137   r[1] = factor * d[1];
 138   r[2] = factor * d[2];
 139   */
 140   /*
 141   update_cache();
 142   const double factor = cache[2] * abs_d + cache[3] / abs_d - cache[4] * x * y;
 143   r[0] = factor * x / abs_d;
 144   r[1] = factor * y / abs_d;
 145   r[2] = factor * z / abs_d;
 146   */
 147 }
 148 Vec_rhs_fun::Range Vec_rhs_fun::operator()(const Domain& d) const{
 149   Range r {0};
 150   evaluate(d,r);
 151   return r;
 152 }
 153 
 154 sls_rhs::sls_rhs(const Grid::Grid_and_time& gt) 
 155   :tp {gt.time_provider()},
 156    lvec {"lvec", gt.vec_fe_space()},
 157    rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
 158    rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
 159    a {Parameter::getValue<double>("tumor_growth.surface.alpha", 1e-3)},
 160    e {Parameter::getValue<double>("tumor_growth.surface.epsilon", .01)},
 161    k {Parameter::getValue<double>("logistic_growth.steepness", .5)},
 162    delta {Parameter::getValue<double>("tumor_growth.surface.delta", .4)}
 163 {}
 164 auto sls_rhs::operator()(const dom& d) const -> ran{
 165   auto u_fun = [&tp = tp](auto x, auto y){
 166     return x * y * exp(-6.*tp.time());
 167   };
 168   const auto 
 169     norm = sqrt(inner_product(&d[0], &d[0]+ dom::dimension, &d[0], 0.)),
 170     a_til = k * (1 - norm/rE),
 171     mc = dim / norm,
 172     factor = a_til + ((a * a_til + e ) * mc - delta *u_fun(d[0], d[1]) )/ norm;
 173   ran r = d;
 174   r *= factor;
 175   return r;
 176 }
 177 void sls_rhs::addScaled_to(Grid::Vec_FEfun& rhs){
 178   assemble_RHS(*this, lvec);
 179   Grid::Vec_FEfun::Dune_FEfun& fef = rhs;
 180   fef.axpy(tp.deltaT(), lvec);
 181 }
 182 
 183 sd_rhs::sd_rhs(const Grid::Grid_and_time& gt)
 184   :tp {gt.time_provider()},
 185    lvec {"lvec", gt.vec_fe_space()},
 186    a {Parameter::getValue<double>("tumor_growth.surface.alpha", 1e-3)},
 187    e {Parameter::getValue<double>("tumor_growth.surface.epsilon", .01)}
 188 {}
 189 auto sd_rhs::operator()(const dom& d) const -> ran{
 190   const auto fac = (1. + (a+e)*dim*exp(-2. * tp.time()));
 191   ran r {d};
 192   r *= fac;
 193   return r;
 194 }
 195 void sd_rhs::addScaled_to(Grid::Vec_FEfun& rhs){
 196   assemble_RHS(*this, lvec);
 197   Grid::Vec_FEfun::Dune_FEfun& fef = rhs;
 198   fef.axpy(tp.deltaT(), lvec);
 199 }
 200 
 201 sdp_u_rhs::sdp_u_rhs(const Grid::Grid_and_time& gt)
 202   :tp {gt.time_provider()},
 203    lscal {"lscal", gt.fe_space()},
 204    rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
 205    rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
 206    k {Parameter::getValue<double>("logistic_growth.steepness", .5)}
 207 {}
 208 auto sdp_u_rhs::operator()(const dom& d) const -> ran{
 209   const auto x = d[0], y = d[1], z = d[2], t = tp.time();
 210   return 2*k*x*y*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1)*exp(-6*t) - (k*(pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1) + k*(pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1) + k*(pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*(rA/(rA*(exp(-k*t) - 1) - rE*exp(-k*t)) + 1))*x*y*exp(-6*t) + 4*pow(x,3)*y*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),3) + 4*x*pow(y,3)*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),3) - 6*x*y*exp(-6*t) + 2*(pow(x,3)*y*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) + (pow(x,3)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x/(pow(x,2) + pow(y,2) + pow(z,2)))*y*exp(-6*t) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1) + (4*pow(x,2)*pow(y,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - pow(x,2)*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)) - (pow(x,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*exp(-6*t))*x*y/(pow(x,2) + pow(y,2) + pow(z,2)) + (4*pow(x,2)*pow(y,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - pow(y,2)*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)) - (pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1)*exp(-6*t))*x*y/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(x*pow(y,3)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) + x*(pow(y,3)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - y/(pow(x,2) + pow(y,2) + pow(z,2)))*exp(-6*t) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(y,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1) + 2*(2*pow(x,2)*y*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - y*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*x*z/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(2*x*pow(y,2)*z*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x*z*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*y*z/(pow(x,2) + pow(y,2) + pow(z,2)) + 2*(2*x*y*pow(z,2)*exp(-6*t)/pow(pow(x,2) + pow(y,2) + pow(z,2),2) - x*y*exp(-6*t)/(pow(x,2) + pow(y,2) + pow(z,2)))*(pow(z,2)/(pow(x,2) + pow(y,2) + pow(z,2)) - 1);
 211 }
 212 void sdp_u_rhs::addScaled_to(Grid::Scal_FEfun& rhs){
 213   assemble_RHS(*this, lscal);
 214   Grid::Scal_FEfun::Dune_FEfun& fef = rhs;
 215   fef.axpy(tp.deltaT(), lscal);
 216 }