root/old_code/tumor_growth/tumor_deformation.h

DEFINITIONS

1. set_time_provider
2. evaluate
3. evaluate
4. evaluate

   1 #ifndef DEFORMATION_HH
   2 #define DEFORMATION_HH
   3 
   4 #include <dune/common/exceptions.hh>
   5 #include <dune/grid/geometrygrid/coordfunction.hh>
   6 #include <dune/fem/space/common/functionspace.hh>
   7 
   8 #include<cassert>
   9 #include<dune_bdf.hpp>
  10 
  11 #include <Eigen/Dense>
  12 
  13 // DeformationCoordFunction
  14 // ------------------------
  15 
  16 class DeformationCoordFunction
  17   : public Dune::AnalyticalCoordFunction<double, 3, 3, DeformationCoordFunction>{
  18   typedef Dune::AnalyticalCoordFunction<double, 3, 3, DeformationCoordFunction> 
  19   BaseType;
  20 public:
  21   typedef Dune::Fem::FunctionSpace<double, double, 3, 3> FunctionSpaceType;
  22 
  23   typedef BaseType::DomainVector DomainVector;
  24   typedef BaseType::RangeVector RangeVector;
  25 
  26   DeformationCoordFunction() {};
  27   DeformationCoordFunction (BDF::EvoMapType& evoMap)
  28     : evoMapPt {&evoMap} {}
  29   void set_time_provider(Dune::Fem::TimeProviderBase& time_provider){
  30     p_tp = &time_provider;
  31   }
  32   void evaluate (const DomainVector& x, RangeVector& y) const{
  33     // double t = p_tp->time();
  34     
  35     // '2015linfty' example
  36     // y[0] = sqrt(1. + sin(2. * M_PI * t )/4.) * x[0];
  37     // y[1] = x[1];
  38     // y[2] = x[2];             
  39     
  40     // ODE solver code or nonlinear evolution code
  41     // std:: cout <<  "Hello\n"
  42     //         << x[0] << ' ' << x[1] << ' ' << x[2] << std::endl;
  43 
  44     std::array<double, 3> domain_node { {x[0], x[1], x[2]} };
  45 
  46     const std::array<double, 3> range_node = (*evoMapPt)[domain_node];
  47     for(std::size_t s = 0; s < range_node.size(); ++s)
  48       y[s] = range_node[s];
  49 
  50     // 'baseball bat' example
  51     //const double newtime
  52     // = (time < 0.5 ? 0.0 : (time < 0.75 ? 4*(time - 0.5) : 1.0));
  53     //const double newtime = std::min( time, 1.0 );
  54     //const double r1 = std::abs( x[ 0 ] );
  55     //const double target 
  56     // = (1.0 - (r1*r1))*((r1*r1) + 0.05) + (r1*r1)*sqrt(1.0 - (r1*r1));
  57                         
  58     //const double r2 = std::sqrt( x[1]*x[1] + x[2]*x[2] );
  59     //const double factor 
  60     // = std::exp( -2*newtime )*r2 + (1.0 - std::exp( -2*newtime ))*target;
  61 
  62     //y[ 0 ] = 2 * x[ 0 ] + newtime*(x[ 0 ] > 0 ? 2.0 : -1.0 )*x[ 0 ];
  63     //y[ 1 ] = factor * x[ 1 ] / (r2 + 0.000001);
  64     //y[ 2 ] = factor * x[ 2 ] / (r2 + 0.000001);
  65     
  66     //y[ 0 ] = x[ 0 ];
  67     //y[ 1 ] = x[ 1 ]; //(1+ .05 * std::sin(4*M_PI*time));
  68     //y[ 2 ] = x[ 2 ]; //* (1+ .5 * std::cos(2*M_PI*time));       
  69   }
  70 private:
  71   BDF::EvoMapType* evoMapPt = nullptr;  
  72         // use this if your surface evolves via an ODE
  73   Dune::Fem::TimeProviderBase* p_tp = nullptr;
  74 };
  75 
  76 // ------------------------------------------------------------
  77 // CAPG: For what is this needed?????
  78 // ------------------------------------------------------------
  79 
  80 //! deformation depending on a discrete function 
  81 template <class DiscreteFunctionType>
  82 class DeformationDiscreteFunction
  83   : public Dune::DiscreteCoordFunction<double, 3, 
  84                                        DeformationDiscreteFunction<DiscreteFunctionType> 
  85                                        >{
  86   typedef Dune::DiscreteCoordFunction<double, 3, 
  87                                       DeformationDiscreteFunction<DiscreteFunctionType>
  88                                       > BaseType;
  89   typedef typename DiscreteFunctionType::GridType GridType;
  90   typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
  91   typedef typename DiscreteFunctionType::RangeType  RangeType;
  92 public:
  93   DeformationDiscreteFunction (const DiscreteFunctionType& vertices )
  94     : vertices_( vertices )
  95   {}
  96 
  97   template< class HostEntity , class RangeVector >
  98   void evaluate ( const HostEntity& hostEntity, unsigned int corner,
  99                   RangeVector& y ) const{
 100     DUNE_THROW(Dune::NotImplemented,"evaluate not implemented for codim > 0");
 101   }
 102 
 103   template <class RangeVector> 
 104   void evaluate ( const typename GridType:: template Codim<0>::Entity& entity, 
 105                   unsigned int corner, RangeVector& y ) const{
 106     typedef typename GridType::ctype  ctype;
 107     enum { dim = GridType::dimension };
 108     
 109     const Dune::GenericReferenceElement<ctype, dim>& refElement
 110       = Dune::GenericReferenceElements<ctype, dim>::general( entity.type() );
 111     
 112     LocalFunctionType localVertices = vertices_.localFunction( entity );        
 113     
 114     localVertices.evaluate( refElement.position( corner, dim ), y );
 115   }
 116   // void setTime ( const double time ){}
 117 protected:
 118   const DiscreteFunctionType& vertices_;
 119 };
 120 
 121 #endif // #ifndef DEFORMATION_HH