old_code/2015linfty/dune

/* [<][>][^][v][top][bottom][index][help] */
This source file includes following definitions.
evaluate
jacobian
evaluate
jacobian
evaluate
jacobian
evaluate
jacobian
   1 /*********************************************************************
   2  *  dune_bdf.hpp, dune_bdf.hxx                                       *
   3  *                                                                   *
   4  *  Defines auxiliar classes like a string vector with the correct   *
   5  *  file names.                                                      *
   6  *                                                                   *
   7  *  Revision history:                                                *
   8  *  none (this is experimental and error phrone)                     *
   9  *                                                                   *
  10  *                                                                   *
  11  *  Created by Christian Power on 19.06.14.                          *
  12  *  Copyright (c) 2014 Christian Power. All rights reserved.         *
  13  *                                                                   *
  14  *********************************************************************/
  15 
  16 #ifndef DUNE_BDF_HXX
  17 #define DUNE_BDF_HXX
  18 
  19 // all implementation details
  20 
  21 // helper function
  22 
  23 // Usually in only template definitions are put into a *.hxx file
  24 // (CAPG convention), but the dune makefiles are very complicated,
  25 // so I put also the function definition in there.
  26 
  27 // starting the interface for EvoMapType::evolve()
  28 struct SolverType : 
  29         NUMERIK::NewtonLinearSolverType< Eigen::Vector3d, Eigen::Matrix3d, SolverType>{
  30         Eigen::Vector3d operator() (const Eigen::Matrix3d& A, const Eigen::Vector3d& b) const {
  31                 return A.fullPivLu().solve( (-1) * b);
  32         }
  33 };
  34 
  35 struct CAPGNormType : NUMERIK::NormFunctorType<double, Eigen::Vector3d, CAPGNormType>{
  36         double operator() (Eigen::Vector3d v) const{
  37                 return v.norm();
  38         }
  39 };
  40 
  41 struct SurfaceVectorfieldType : 
  42         NUMERIK::C1SpaceTimeFunctionType
  43         <Eigen::Vector3d, Eigen::Matrix3d, SurfaceVectorfieldType>{
  44         typedef C1SpaceTimeFunctionType
  45         <Eigen::Vector3d, Eigen::Matrix3d, SurfaceVectorfieldType> BaseType;
  46         using BaseType::get_time;
  47         using BaseType::get_position;
  48         using BaseType::get_dT;
  49 
  50         SurfaceVectorfieldType(const double starting_time, const Eigen::Vector3d& starting_pos) : 
  51                 BaseType(starting_time, starting_pos) {}
  52         Eigen::Vector3d evaluate(Eigen::Vector3d vec) const{
  53                 double x = vec(0), y = vec(1), z = vec(2), t = get_time();
  54                 Eigen::Vector3d f_t_vec( 
  55                         16.*M_PI*pow(x,3)*cos(2*M_PI*t)/( (pow(y,2) + pow(z,2) + 16.*pow(x,2)/pow(sin(2.*M_PI*t) + 4.,2) )*pow(sin(2.*M_PI*t) + 4.,3) ),
  56                         4.*M_PI*pow(x,2)*y*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16.*pow(x,2)/pow(sin(2.*M_PI*t) + 4.,2) )*pow(sin(2.*M_PI*t) + 4.,2)),
  57                         4.*M_PI*pow(x,2)*z*cos(2.*M_PI*t)/((pow(y,2) + pow(z,2) + 16.*pow(x,2)/pow(sin(2*M_PI*t) + 4.,2) )*pow(sin(2.*M_PI*t) + 4,2) ) );
  58                 return  vec - get_dT() * f_t_vec - get_position();
  59                 // x_next - dT * f(t_next, x_next) - x_n
  60         }
  61         Eigen::Matrix3d jacobian(Eigen::Vector3d vec) const{
  62                 double x = vec(0), y = vec(1), z = vec(2), t = get_time();
  63                 Eigen::Matrix3d jac_f_t_vec;
  64                 jac_f_t_vec << 
  65                         // first row
  66                         48*M_PI*pow(x,2)*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,3)) - 512*M_PI*pow(x,4)*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,5)), 
  67                         -32*M_PI*pow(x,3)*y*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,3)), 
  68                         -32*M_PI*pow(x,3)*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,3)),
  69                         // second row
  70                         8*M_PI*x*y*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2)) - 128*M_PI*pow(x,3)*y*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,4)), 
  71                         -8*M_PI*pow(x,2)*pow(y,2)*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)) + 4*M_PI*pow(x,2)*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2)), 
  72                         -8*M_PI*pow(x,2)*y*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)), 
  73                         // third row
  74                         8*M_PI*x*z*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2)) - 128*M_PI*pow(x,3)*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,4)), 
  75                         -8*M_PI*pow(x,2)*y*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)), 
  76                         -8*M_PI*pow(x,2)*pow(z,2)*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)) + 4*M_PI*pow(x,2)*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2));
  77                 return Eigen::Matrix3d::Identity() - get_dT() * jac_f_t_vec;
  78                 // id_3x3 - dT * jac(f)(t_next, x_next)
  79         }
  80 };
  81 
  82 struct ALESurfaceVectorfieldType : 
  83         NUMERIK::C1SpaceTimeFunctionType
  84         <Eigen::Vector3d, Eigen::Matrix3d, ALESurfaceVectorfieldType>{
  85         typedef C1SpaceTimeFunctionType
  86         <Eigen::Vector3d, Eigen::Matrix3d, ALESurfaceVectorfieldType> BaseType;
  87         using BaseType::get_time;
  88         using BaseType::get_position;
  89         using BaseType::get_dT;
  90 
  91         ALESurfaceVectorfieldType(const double starting_time, 
  92                                                           const Eigen::Vector3d& starting_pos) : 
  93                 BaseType(starting_time, starting_pos) {}
  94 
  95         Eigen::Vector3d evaluate(Eigen::Vector3d vec) const{
  96                 double x = vec(0), y = vec(1), z = vec(2), t = get_time();
  97                 Eigen::Vector3d f_t_vec( 
  98                         .5 * M_PI * cos(2. * M_PI * t) / ( 2. * sqrt( 1. + .25 * sin(2. * M_PI * t) ) ) *x, 0, 0 );
  99                 return  vec - get_dT() * f_t_vec - get_position();
 100                 // x_next - dT * f(t_next, x_next) - x_n
 101         }
 102         Eigen::Matrix3d jacobian(Eigen::Vector3d vec) const{
 103                 double x = vec(0), y = vec(1), z = vec(2), t = get_time();
 104                 Eigen::Matrix3d jac_f_t_vec;
 105                 jac_f_t_vec << 
 106                         // first row
 107                         .5 * M_PI * cos(2. * M_PI * t) / ( 2. * sqrt( 1. + .25 * sin(2. * M_PI * t) ) ) , 0. , 0.,
 108                         // second row
 109                         0. , 0. , 0.,
 110                         // third row
 111                         0. , 0. , 0.;
 112                 return Eigen::Matrix3d::Identity() - get_dT() * jac_f_t_vec;
 113                 // id_3x3 - dT * jac(f)(t_next, x_next)
 114         }
 115 };
 116 
 117 struct StationaryType : 
 118         NUMERIK::C1SpaceTimeFunctionType <Eigen::Vector3d, Eigen::Matrix3d, StationaryType>{
 119         
 120         typedef C1SpaceTimeFunctionType
 121         <Eigen::Vector3d, Eigen::Matrix3d, StationaryType> BaseType;
 122         using BaseType::get_time;
 123         using BaseType::get_position;
 124         using BaseType::get_dT;
 125         
 126         StationaryType(const double starting_time, 
 127                                    const Eigen::Vector3d& starting_pos) : 
 128                 BaseType(starting_time, starting_pos) {}
 129 
 130         Eigen::Vector3d evaluate(Eigen::Vector3d vec) const{
 131                 return  vec - get_position();
 132         }
 133         Eigen::Matrix3d jacobian(Eigen::Vector3d vec) const{
 134                 return Eigen::Matrix3d::Identity();
 135         }
 136 };
 137 
 138 struct RotatingType : 
 139         NUMERIK::C1SpaceTimeFunctionType <Eigen::Vector3d, Eigen::Matrix3d, RotatingType>{ 
 140         
 141         typedef C1SpaceTimeFunctionType
 142         <Eigen::Vector3d, Eigen::Matrix3d, RotatingType> BaseType;
 143         using BaseType::get_position;
 144         
 145         RotatingType(const double starting_time, 
 146                                    const Eigen::Vector3d& starting_pos) : 
 147                 BaseType(starting_time, starting_pos) {}
 148 
 149         Eigen::Vector3d evaluate(Eigen::Vector3d vec) const{
 150                 return  vec - get_position();
 151         }
 152         Eigen::Matrix3d jacobian(Eigen::Vector3d vec) const{
 153                 return Eigen::Matrix3d::Identity();
 154         }
 155 };
 156 #endif // #ifndef DUNE_BDF_HXX
/* [<][>][^][v][top][bottom][index][help] */
root/old_code/2015linfty/dune_bdf.hxx

DEFINITIONS
