root/old_code/tumor_growth/tumor_bdf.h

DEFINITIONS

1. getVertexList
2. get_relevant_digits
3. get_precision
4. array_to_string
5. string_to_array
6. evolve
7. save_original_vertices
8. evolve

   1 #ifndef DUNE_BDF_HPP
   2 #define DUNE_BDF_HPP
   3 
   4 /*********************************************************************
   5  *  dune_bdf.hpp, dune_bdf.hxx                                       *
   6  *                                                                   *
   7  *  Defines auxiliar classes like a string vector with the correct   *
   8  *  file names.                                                      *
   9  *                                                                   *
  10  *  Revision history:                                                *
  11  *  none (this is experimental and error phrone)                     *
  12  *                                                                   *
  13  *                                                                   *
  14  *  Created by Christian Power on 09.03.15.                          *
  15  *  Copyright (c) 2015 Christian Power. All rights reserved.         *
  16  *                                                                   *
  17  *********************************************************************/
  18 
  19 #include <iostream>             // debugging
  20 #include <vector>
  21 #include <string>
  22 #include <fstream>
  23 #include <sstream>
  24 #include <cmath>
  25 #include <exception>
  26 #include <unordered_map>
  27 //#include <map>
  28 #include <array>
  29 #include <algorithm>
  30 #include <Eigen/Dense>
  31 
  32 namespace BDF
  33 {
  34   class EvoMapType{
  35   public:
  36     EvoMapType(const std::string& input_filename);
  37 
  38     void save_original_vertices();
  39 
  40     template<typename In, typename Integrator>
  41         // requires Input_iterator<In>()
  42     void evolve(In p_t_0, In p_t_end, Integrator integrator);
  43 
  44     template<typename DiscreteFunction>
  45     void evolve(const DiscreteFunction& new_vertices);
  46     
  47     const std::array<double,3>& operator[] (const std::array<double,3>& vec) const;
  48   private:
  49     std::unordered_map<std::string, std::array<double, 3> > evoMap;
  50     // evoMap: Γʰ(0) → Γʰ(n), evoMap[node⁰] = nodeⁿ;
  51     std::vector<std::string> original_vertices;
  52     std::string filename;
  53     size_t precision;
  54   };
  55 
  56   void ie_heat_algorithm();             // BDF1 algorithm
  57   void bdf2_heat_algorithm();           // debugging
  58   void nonlinear_algorithm();           // the actuall BDF algorithm;
  59   // void algorithm();
  60   // void poisson_problem();    // for testing reasons
  61 
  62   // All implementations details
  63 
  64   // Implementation details for 'EvoMapType's helper functions
  65 
  66   std::vector<std::string> getVertexList(const std::string& filename)
  67   // opens the file 'filename' and return a vector of lines, where each line
  68   // represents a vertex of the grid.  The lines are whitespace cleaned.
  69   {
  70     // local helper functions
  71     auto error = [](const std::string& msg){ 
  72       throw 
  73       std::runtime_error("ERROR: BDF::getVertexList()\n"
  74                          + msg ); 
  75     };
  76     auto clean_whitespace = [](const std::string& line)
  77     // A 'line' may begin with whitespace. Clean it!
  78       {
  79         std::string clean_line;
  80         std::istringstream iss {line};
  81         for(std::string s; iss >> s; )
  82           clean_line += s + ' ';
  83         clean_line.pop_back();  // delete last whitespace
  84         return clean_line;
  85       };
  86     // end local helper functions
  87 
  88     std::ifstream ifs {filename};
  89     if(!ifs) error("File \"" + filename + "\" does not exist!");
  90 
  91     std::vector<std::string> lines;
  92     for(std::string line; getline(ifs, line); )
  93       lines.push_back(clean_whitespace(line));
  94     ifs.close();
  95 
  96     // get list of vertices in terms of iterators
  97     std::vector<std::string>::iterator first =
  98       std::find(lines.begin(), lines.end(), "VERTEX");
  99     if( first == lines.end() ) error("Missing keyword \"VERTEX\" in file \""
 100                                      + filename + "\"!");
 101     ++first;    // "VERTEX" is not part of the list
 102     std::vector<std::string>::iterator last =
 103       std::find(first, lines.end(), "#");
 104     if( last == lines.end() ) error("Missing ending symbol \"#\" in file \""
 105                                     + filename + "\"!");
 106         
 107     return std::vector<std::string> {first, last};
 108         
 109   }
 110 
 111   size_t get_relevant_digits(const std::string& number){
 112     // assume scientific format 0.9298e-10
 113     // want in this example the number 4 back
 114     const auto dot_iter = number.find(".");
 115     const auto exp_iter = number.find("e");
 116     const auto dist = exp_iter - dot_iter ;
 117     return dist > 0 ? dist : 0;
 118   }
 119 
 120   size_t get_precision(const std::vector<std::string>& vertex_list){
 121     if(vertex_list.size() == 0) return 0;
 122 
 123     std::vector<size_t> precision_vec;
 124     for(const std::string& line : vertex_list){
 125       std::istringstream interpreter {line};
 126       for(std::string s; interpreter >> s; ){
 127         const auto precision = get_relevant_digits(s);
 128         precision_vec.push_back(precision);
 129       }
 130     }
 131 
 132     // std::cout << "precision_vec:" << std::endl;
 133     // for(const auto it : precision_vec)
 134     //   std::cout << it << ' ';
 135     // std::cout << std::endl;
 136 
 137     std::vector<size_t>::iterator  p_precision = 
 138       std::max_element(precision_vec.begin(), precision_vec.end());
 139     return *p_precision;
 140   }
 141 
 142   template<typename T, size_t N >
 143   std::string array_to_string(const std::array<T, N> ad,
 144                               const unsigned int precision = 8){
 145     std::ostringstream oss;
 146     oss << std::scientific;
 147     // setprecision and scientific is essential! For setprecision it is mandatory
 148     // that the number  be greater then precision in the dgf file.
 149     for(std::size_t s = 0; s < ad.size() -1 ; ++s)
 150       oss << std::setprecision(precision) << ad[s] << ' ';
 151     oss << ad.back();   // the last one without ' '
 152     return oss.str();
 153   }
 154 
 155   template<typename T, size_t N >
 156   std::array<T, N> string_to_array(const std::string s){
 157     auto error = [] (const std::string& msg) { 
 158       throw std::runtime_error("ERROR: BDF::string_to_array()\n" + msg); 
 159     };
 160     std::istringstream interpreter {s};
 161     std::array<T, N> return_array;
 162     for(auto& it : return_array)
 163       if(!(interpreter >> it)) error("operator>>() failed!");
 164     if(!(interpreter >> std::ws).eof()) // stuff left in stream?
 165       error("whitespace left!");
 166     return return_array;
 167   }
 168 
 169   // Implementation details for 'EvoMapType'
 170 
 171   BDF::EvoMapType::EvoMapType( const std::string& input_filename ) :
 172     filename {input_filename}
 173         {
 174           std::vector<std::string> lines = getVertexList(filename);
 175           precision = get_precision(lines);
 176           //std::cout << "precision = " << precision << std::endl; // debugging
 177           for(auto& line : lines)
 178             evoMap[ array_to_string( string_to_array<double, 3>(line), precision) ] = 
 179               string_to_array<double, 3>(line);
 180         }
 181 
 182   const std::array<double, 3>&
 183   BDF::EvoMapType::operator[] (const std::array<double, 3>& node) const { 
 184     const std::array<double, 3>* cpad3 = nullptr;
 185     const auto node_as_string = array_to_string(node, precision);
 186     try{
 187       cpad3 = &evoMap.at(node_as_string);
 188     }
 189     catch(std::exception& e){
 190       // std::cerr << "Error in EvoMapType::operator[].\n"
 191       //                << e.what() << '\n';
 192       std::cout << "Your node\n"
 193         + node_as_string 
 194         + "\nWas not good.  Full node list." << std::endl;
 195       for(const auto& it : evoMap)
 196         std::cout << it.first << std::endl;
 197     }
 198     // if(cpad3) throw std::runtime_error {"We have a nullptr.  Aborting."};
 199     return *cpad3;
 200   }
 201 
 202   template<typename In, typename Integrator>
 203   // requires Input_iterator<In>() && 
 204   // Vector3d Integrator.integrate(In first, In last, Vector3d x_0)
 205   // Vector3d from the 'Eigen' library
 206   void BDF::EvoMapType::evolve(In p_t_0, In p_t_end, Integrator integrator){
 207     for(auto& it : evoMap){
 208       Eigen::Map<Eigen::Vector3d> x_0 {it.second.data()};       // get vector
 209       x_0 = integrator.integrate(p_t_0, p_t_end, x_0);          // perform action
 210       std::array<double, 3> it_new;                             // translate it back
 211       std::copy(x_0.data(), x_0.data() + x_0.size(), it_new.begin());
 212       it.second = it_new;
 213 
 214       // nice testing idea by Balázs
 215       // std::cout << "surface error: " 
 216       //           << (x_0(0)*x_0(0) / (1. + (sin(2. * M_PI * time_vec.back() ) / 4.)
 217       //                                ) ) + x_0(1) * x_0(1) + x_0(2) * x_0(2) - 1
 218       //           << std::endl;
 219     }
 220   }
 221 
 222   void EvoMapType::save_original_vertices(){
 223     // Important: I assume that the sorting of the nodal values of a 
 224     // finite element function is the same like the the vertex list from the 
 225     // original file
 226 
 227     original_vertices = getVertexList(filename);
 228     precision = get_precision(original_vertices);
 229     for(auto& vertex : original_vertices)
 230       vertex = array_to_string( string_to_array<double, 3>(vertex), precision);
 231 
 232   }
 233 
 234   template<typename DiscreteFunction>
 235   void BDF::EvoMapType::evolve(const DiscreteFunction& new_vertices){
 236     // I assume original_vertices are already constructed.  
 237     // C.f. save_original_vertices()
 238 
 239     constexpr auto dim = 3;
 240     const auto no_vertices = original_vertices.size();
 241     const auto dof_times_dim = new_vertices.size();
 242 
 243     if(no_vertices * dim != dof_times_dim) throw std::runtime_error
 244       {"Error in EvoMapType::evolve().\nPerhaps forgotten to use"
 245           " EvoMapType::save_original_vertices()?"};
 246 
 247     auto p_dof_component = new_vertices.dbegin();
 248     for(size_t it_vertex = 0; it_vertex < no_vertices; ++it_vertex){
 249       // initialize an array containing the new vertex
 250       std::array<double, dim> array_new_vertex;
 251       for(size_t  it_component = 0; it_component < dim; ++it_component){
 252         array_new_vertex.at(it_component) = *p_dof_component;
 253         ++p_dof_component;
 254       }
 255       // evolve
 256       const auto str_original_vertex = original_vertices.at(it_vertex);
 257       evoMap[str_original_vertex] = array_new_vertex;
 258     }
 259   }
 260 }       // namespace BDF
 261 
 262 
 263 
 264 
 265 // #include "dune_bdf.hxx"
 266 // CAPG remark: (CAPG convention) usually we include in *.hxx files only
 267 // template related implementation, but because the dune make files are
 268 // complicated, we put all the *.cpp stuff there
 269 
 270 #endif // #ifndef DUNE_BDF_HPP