root/src/grid_deformation.cpp

DEFINITIONS

1. identity
2. logistic_growth
3. dalquist
4. mcf_sphere
5. set_timeProvider
6. evaluate

   1 /*! \file grid_deformation.cpp
   2     \brief Implementation of Deformation class
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power April 2016
   8           Originally written by Christian Power
   9                (power22c@gmail.com) Januar 2016
  10 
  11      Idea
  12      --------------------------------------------------
  13 
  14      Explicit flow functions are coded as inline functions.
  15 
  16     \author Christian Power
  17     \date 23. April 2016
  18     \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  19  */
  20 
  21 #include <cmath>
  22 #include <numeric>
  23 #include "grid.h"
  24 #include "grid_GridAndTime_impl.h"
  25 
  26 using namespace std;
  27 using Esfem::Grid::Deformation;
  28 //! \f$\R^3\f$
  29 using Domain = Esfem::Grid::Deformation::Domain;
  30 //! \f$\R^3\f$
  31 using Range = Esfem::Grid::Deformation::Range;
  32 
  33 static_assert(Esfem::Grid::Deformation::Domain::dimension == 3,
  34               "Bad domain dimension.");
  35 static_assert(Esfem::Grid::Deformation::Range::dimension == 3,
  36               "Bad range dimension.");
  37 
  38 //! \f$id\colon \R^3 \to \R^3\f$
  39 static inline void identity(const Domain& x, Range& y) noexcept{  
  40   y[0] = x[0]; 
  41   y[1] = x[1]; 
  42   y[2] = x[2]; 
  43 }
  44 
  45 
  46 //! \f$r(t) = \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1 - e^{-kt})}\f$
  47 /*! \param t Current time
  48   \param x Point from the initial surface
  49   \retval y \f$ y = r(t) x\f$
  50   \pre Initial surface is a sphere.
  51  */
  52 static inline void logistic_growth(const double t, const Domain& x, Range& y) noexcept{
  53   const double r_end = 2., r0 = 1., k = .5; // logistic function parameter
  54   const double r = r_end * r0 / (r_end*exp(-k*t) + r0*(1-exp(-k*t)));
  55   y = x;
  56   y *= r;
  57 }
  58 
  59 //! Dalquist test equation with \f$\lambda=1\f$
  60 static inline void dalquist(const double t, const Domain& x, Range& y){
  61   const double factor = exp(-t);
  62   y = x;
  63   y *= factor;
  64 }
  65 
  66 //! \f$ R(t) = \sqrt{ R_0^2 - 2nt}\f$
  67 static inline void mcf_sphere(const double t, const Domain& x, Range& y){
  68   const double norm_square
  69     // = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
  70     // = 1.;
  71     = std::inner_product(&x[0], &x[0] + Domain::dimension, &x[0], 0.);
  72   const double factor = sqrt( norm_square - 2 * Esfem::Grid::grid_dim() * t);
  73   y = x;
  74   y *= factor;
  75 }
  76 
  77 // ----------------------------------------------------------------------
  78 // Implementaion of Deformation
  79 
  80 struct Esfem::Grid::Deformation::Data{
  81   // Old hash map
  82   // Impl::Evolving_grid eg_ptr;
  83   //! New hash map
  84   Impl::hash::grid hg;
  85   //! I do not assume ownership
  86   const Dune::Fem::TimeProviderBase* tp_ptr {nullptr};
  87   //! All pointer to zero
  88   Data() = default;
  89   //! Use hash map
  90   Data(const std::string& fname) :hg {fname} {}
  91 };
  92 
  93 Esfem::Grid::Deformation::Deformation() :d_ptr {std::make_unique<Data>()} {}
  94 Esfem::Grid::Deformation::Deformation(const std::string& fname)
  95   :d_ptr {std::make_unique<Data>(fname)} {}
  96 Esfem::Grid::Deformation::~Deformation() = default;
  97 
  98 void Esfem::Grid::Deformation::
  99 set_timeProvider(const Dune::Fem::TimeProviderBase& tp){
 100   d_ptr -> tp_ptr = &tp;
 101 }
 102 void Esfem::Grid::Deformation::evaluate(const Domain& x, Range& y) const{
 103   // const double t = d_ptr -> tp_ptr->time();
 104 
 105   // mcf_sphere(t, x, y);
 106 
 107   // dalquist(t, x, y);
 108 
 109   // logistic_growth(t, x, y);
 110 
 111   // identity(x, y);
 112   
 113   y = d_ptr->hg[x];
 114 
 115   // const auto eg = *(d_ptr -> eg_ptr);
 116   // y = eg[x];
 117 }
 118 Deformation& Deformation::operator=(const Vec_FEfun& rhs){
 119   d_ptr->hg = rhs;
 120   return *this;
 121 }