root/src/maxH_main.cpp

DEFINITIONS

1. init_nodelist
2. init_elemlist
3. evolve_nodes
4. all_combinations
5. distance
6. max_h
7. print
8. print_what
9. main

   1 /*! \file maxH_main.cpp
   2     \brief Calculate the maximum h of a finite element grid
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power July 2016
   8           Originally written by Christian Power
   9                (power22c@gmail.com) July 2016
  10 
  11     \author Christian Power 
  12     \date 07. July 2016
  13     \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  14  */
  15 
  16 #include <iostream>
  17 #include <valarray>
  18 #include <vector>
  19 #include <fstream>
  20 #include <sstream>
  21 #include <dassert.h>
  22 
  23 using namespace std;
  24 
  25 //! File as a vector of lines
  26 struct file : vector<string>{
  27   //! Report an error
  28   struct bad : runtime_error{
  29     using runtime_error::runtime_error;
  30   };
  31 
  32   //! For error report
  33   const string filename;
  34 
  35   //! Read file
  36   explicit file(const string& fname)
  37     :filename {fname}
  38 
  39   {
  40     ifstream ifs {filename};
  41     if(!ifs) throw bad {Assert::compose(__FILE__, __LINE__, fname)};
  42     for(string s; getline(ifs, s); ) push_back(s);
  43     ifs.close();
  44   }
  45 };
  46 
  47 //! Nodes as a vector of valarrays of size 3
  48 template<class T>
  49 struct nodes : vector<valarray<T> >{
  50   //! Report an error
  51   struct bad : runtime_error{
  52     using runtime_error::runtime_error;
  53   };
  54   //! Read nodes from file
  55   /*! \warning Potentially narrowing double to int! */
  56   explicit nodes(const file& f){
  57     istringstream iss;
  58     iss.exceptions(ios_base::badbit | ios_base::failbit);
  59     size_t lineno {1};
  60     try{
  61       for(const auto& line : f){
  62         iss.str(line);
  63         double d1, d2, d3;
  64         iss >> d1 >> d2 >> d3;
  65         T t1(d1), t2(d2), t3(d3);  // this may narrow
  66         this->push_back({t1, t2, t3});
  67         iss.clear();
  68         ++lineno;
  69       }
  70     }
  71     catch(...){
  72       ostringstream oss;
  73       oss << f.filename << ", line " << lineno;
  74       throw_with_nested(bad {Assert::compose(__FILE__, __LINE__, oss.str())});
  75     }
  76   }
  77 };
  78 //! Balazs node list
  79 using nodelist = nodes<double>;
  80 //! Balazs element list
  81 using elemlist = nodes<int>;
  82 
  83 //! Nodelist with corresponding elemlist
  84 struct fempair{
  85   //! Report an error
  86   struct bad : runtime_error{
  87     using runtime_error::runtime_error;
  88   };
  89 
  90   //! The node list
  91   nodelist n;
  92   //! Zero based element list
  93   elemlist e;
  94   
  95   explicit fempair(const int& no) try
  96     :n {init_nodelist(no)}, e {init_elemlist(no)} {}
  97   catch(...){
  98     throw_with_nested(bad {Assert::compose(__FILE__, __LINE__ , "Constructor")});
  99   }
 100 private:
 101   //! Constructor helper
 102   nodelist init_nodelist(const int& no){
 103     ostringstream oss;
 104     oss << SPHEREPATH << "/raw_data/Sphere_nodes" << no << ".txt";
 105     file f {oss.str()};
 106     return nodelist {f};
 107   }
 108   //! Constructor helper
 109   elemlist init_elemlist(const int& no){
 110     ostringstream oss;
 111     oss << SPHEREPATH << "/raw_data/Sphere_elements" << no << ".txt";
 112     file f {oss.str()};
 113     elemlist e {f};
 114     // make e zero based
 115     for(auto& vi : e) vi -= 1;
 116     return e;
 117   }
 118 };
 119 
 120 template<class F>
 121 void evolve_nodes(fempair& fm, F f){
 122   for(auto& vd : fm.n) f(vd);
 123 }
 124 
 125 /*! \pre valarray has size of 3 */
 126 vector<pair<int, int> > all_combinations(const valarray<int>& v){
 127   Assert::dynamic(v.size() == 3, Assert::compose(__FILE__, __LINE__, "!= 3"));
 128   return {{v[0], v[1]}, {v[0], v[2]}, {v[1], v[2]}}; // has size 3
 129 }
 130 
 131 /*! \pre valarray has size of 3*/
 132 double distance(const valarray<double>& vd1, const valarray<double>& vd2){
 133   Assert::dynamic(vd1.size() == 3 && vd2.size() == 3,
 134                   Assert::compose(__FILE__, __LINE__, "!=3"));
 135   double rv {};
 136   for(int i = 0; i < 3; ++i) rv += pow(vd1[i] - vd2[i],2);
 137   return sqrt(rv);
 138 }
 139 
 140 double max_h(const fempair& fm){
 141   constexpr int no_edges {3};
 142   valarray<double> all_h(fm.e.size() * no_edges);
 143   auto ptr = begin(all_h);
 144   for(const auto& elem : fm.e){
 145     const auto comb = all_combinations(elem);
 146     for(const auto& p : comb){
 147       *ptr = distance(fm.n.at(p.first), fm.n.at(p.second));
 148       ++ptr;
 149     }
 150   }
 151   return all_h.max();
 152 }
 153 
 154 class evolution{
 155   double t {2.};
 156   double rE {2.};
 157   double rA {1.};
 158   double k {.5};
 159   double rt {(rE * rA)/( rE * exp(-k*t) + rA *( 1 - exp(-k*t)))};
 160 public:
 161   template<class C>
 162   void operator()(C& c){ c *= rt; }
 163 };
 164 
 165 //! For testing nodes<T>
 166 template<class T>
 167 void print(const valarray<T>& v){
 168   cout << "(";
 169   for(size_t i = 0; i < v.size() - 1; ++i) cout << v[i] << ", ";
 170   cout << v[v.size()-1] << ")" << endl;
 171 }
 172 
 173 void print_what(const exception& e){
 174   cerr << e.what() << '\n';
 175   try{
 176     rethrow_if_nested(e);
 177   }
 178   catch(const exception& e){
 179     print_what(e);
 180   }
 181   catch(...){
 182     cerr << "Unknown nested error\n";
 183   }
 184 }
 185 
 186 int main() try{
 187   for(int i = 1; i <= 5; ++i){
 188     fempair fp {i};
 189     evolve_nodes(fp, evolution {});
 190     const auto h = max_h(fp);
 191     cout << "h_max " << i << ": " << h << endl;
 192   }
 193  }
 194  catch(const exception& e){
 195    print_what(e);
 196    return 1;
 197  }
 198  catch(...){
 199    cerr << "Unknown error\n";
 200    return 2;
 201  }