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 template<typename In, typename Integrator>
39 // requires Input_iterator<In>()
40 void evolve(In p_t_0, In p_t_end, Integrator integrator);
41
42 const std::array<double,3>& operator[] (const std::array<double,3>& vec) const;
43 private:
44 std::unordered_map<std::string, std::array<double, 3> > evoMap;
45 // evoMap: Γʰ(0) → Γʰ(n), evoMap[node⁰] = nodeⁿ;
46 std::string filename;
47 size_t precision;
48 };
49
50 void ie_heat_algorithm(); // BDF1 algorithm
51 void bdf2_heat_algorithm(); // debugging
52 void nonlinear_algorithm(); // the actuall BDF algorithm;
53 // void algorithm();
54 // void poisson_problem(); // for testing reasons
55
56 // All implementations details
57
58 // Implementation details for 'EvoMapType's helper functions
59
60 std::vector<std::string> getVertexList(const std::string& filename)
61 // opens the file 'filename' and return a vector of lines, where each line
62 // represents a vertex of the grid. The lines are whitespace cleaned.
63 {
64 // local helper functions
65 auto error = [](const std::string& msg){
66 throw
67 std::runtime_error("ERROR: BDF::getVertexList()\n"
68 + msg );
69 };
70 auto clean_whitespace = [](const std::string& line)
71 // A 'line' may begin with whitespace. Clean it!
72 {
73 std::string clean_line;
74 std::istringstream iss {line};
75 for(std::string s; iss >> s; )
76 clean_line += s + ' ';
77 clean_line.pop_back(); // delete last whitespace
78 return clean_line;
79 };
80 // end local helper functions
81
82 std::ifstream ifs {filename};
83 if(!ifs) error("File \"" + filename + "\" does not exist!");
84
85 std::vector<std::string> lines;
86 for(std::string line; getline(ifs, line); )
87 lines.push_back(clean_whitespace(line));
88 ifs.close();
89
90 // get list of vertices in terms of iterators
91 std::vector<std::string>::iterator first =
92 std::find(lines.begin(), lines.end(), "VERTEX");
93 if( first == lines.end() ) error("Missing keyword \"VERTEX\" in file \""
94 + filename + "\"!");
95 ++first; // "VERTEX" is not part of the list
96 std::vector<std::string>::iterator last =
97 std::find(first, lines.end(), "#");
98 if( last == lines.end() ) error("Missing ending symbol \"#\" in file \""
99 + filename + "\"!");
100
101 return std::vector<std::string> {first, last};
102
103 }
104
105 size_t get_precision(const std::vector<std::string>& vertex_list){
106 if(vertex_list.size() == 0) return 0;
107
108 std::vector<size_t> precision_vec;
109 for(const std::string& line : vertex_list){
110 std::istringstream interpreter {line};
111 for(std::string s; interpreter >> s; ){
112 precision_vec.push_back(s.size());
113 }
114 }
115
116 std::vector<size_t>::iterator p_precision =
117 std::max_element(precision_vec.begin(), precision_vec.end());
118 return *p_precision;
119 }
120
121 template<typename T, size_t N >
122 std::string array_to_string(const std::array<T, N> ad,
123 const unsigned int precision = 8){
124 std::ostringstream oss;
125 oss << std::scientific;
126 // setprecision and scientific is essential! For setprecision it is mandatory
127 // that the number be greater then precision in the dgf file.
128 for(std::size_t s = 0; s < ad.size() -1 ; ++s)
129 oss << std::setprecision(precision) << ad[s] << ' ';
130 oss << ad.back(); // the last one without ' '
131 return oss.str();
132 }
133
134 template<typename T, size_t N >
135 std::array<T, N> string_to_array(const std::string s){
136 auto error = [] (const std::string& msg) {
137 throw std::runtime_error("ERROR: BDF::string_to_array()\n" + msg);
138 };
139 std::istringstream interpreter {s};
140 std::array<T, N> return_array;
141 for(auto& it : return_array)
142 if(!(interpreter >> it)) error("operator>>() failed!");
143 if(!(interpreter >> std::ws).eof()) // stuff left in stream?
144 error("whitespace left!");
145 return return_array;
146 }
147
148 // Implementation details for 'EvoMapType'
149
150 BDF::EvoMapType::EvoMapType( const std::string& input_filename ) :
151 filename {input_filename}
152 {
153 std::vector<std::string> lines = getVertexList(filename);
154 precision = get_precision(lines);
155 //std::cout << "precision = " << precision << std::endl; // debugging
156 for(auto& line : lines)
157 evoMap[ array_to_string( string_to_array<double, 3>(line), precision) ] =
158 string_to_array<double, 3>(line);
159 }
160
161 const std::array<double, 3>&
162 BDF::EvoMapType::operator[] (const std::array<double, 3>& node) const {
163 return evoMap.at(array_to_string(node, precision));
164 }
165
166 template<typename In, typename Integrator>
167 // requires Input_iterator<In>() &&
168 // Vector3d Integrator.integrate(In first, In last, Vector3d x_0)
169 // Vector3d from the 'Eigen' library
170 void BDF::EvoMapType::evolve(In p_t_0, In p_t_end, Integrator integrator){
171 for(auto& it : evoMap){
172 Eigen::Map<Eigen::Vector3d> x_0 {it.second.data()}; // get vector
173 x_0 = integrator.integrate(p_t_0, p_t_end, x_0); // perform action
174 std::array<double, 3> it_new; // translate it back
175 std::copy(x_0.data(), x_0.data() + x_0.size(), it_new.begin());
176 it.second = it_new;
177
178 // nice testing idea by Balázs
179 // std::cout << "surface error: "
180 // << (x_0(0)*x_0(0) / (1. + (sin(2. * M_PI * time_vec.back() ) / 4.)
181 // ) ) + x_0(1) * x_0(1) + x_0(2) * x_0(2) - 1
182 // << std::endl;
183 }
184 }
185
186 }
187
188
189
190
191 // #include "dune_bdf.hxx"
192 // CAPG remark: (CAPG convention) usually we include in *.hxx files only
193 // template related implementation, but because the dune make files are
194 // complicated, we put all the *.cpp stuff there
195
196 #endif // #ifndef DUNE_BDF_HPP