This source file includes following definitions.
- r_div_rEnd
- eoc_velocity
- Explicit_initial_data
- Random_initial_data
- Random_initial_data
- clone
- interpolate
- evaluate
- clone
- interpolate
- evaluate
- clone
- interpolate
- evaluate
- interpolate
- evaluate
- interpolate
- evaluate
- interpolate
- evaluate
- evaluate
- hom_value
- pertubation
- print_configuration
- dof_filename
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 #include <iostream>
23 #include <sstream>
24 #include <fstream>
25 #include <cmath>
26 #include <numeric>
27 #include <config.h>
28 #include <dune/fem/io/parameter.hh>
29 #include <dune/fem/operator/lagrangeinterpolation.hh>
30 #include "secOrd_op_initData_impl.h"
31 #include "io_parameter.h"
32 #include "esfem_error.h"
33
34
35 using Esfem::Impl::Explicit_initial_data;
36 using Esfem::Impl::Random_initial_data;
37 using Esfem::Impl::Analytic_velocity;
38 using Esfem::Impl::sphere_1EF;
39 using Esfem::Impl::sphere_2EF;
40 using Esfem::Impl::sphere_3EF;
41 using Esfem::Impl::sphere_eigenFun;
42 using Esfem::Impl::sphere_mcf_sol;
43 using Esfem::Impl::sls_iData;
44 using Esfem::Impl::sls_v_iData;
45 using Esfem::Impl::sd_iData;
46 using Dune::Fem::Parameter;
47 using namespace std;
48
49 using Vec_domain = Analytic_velocity::Domain;
50
51 using Vec_range = Analytic_velocity::Range;
52
53
54
55
56
57
58
59
60 static inline double r_div_rEnd(const double t){
61 const double r0 = 1., r_end = 2., kt = .5 * t;
62 return r0 / (r_end * exp(-kt) + r0 * (1 - exp(-kt)) );
63 }
64
65
66
67
68
69
70
71 static inline void eoc_velocity(const double t, const Vec_domain& d, Vec_range& r){
72 const double k = .5,
73 r_d_rEnd = r_div_rEnd(t),
74 factor = k * (1 - r_d_rEnd);
75 r[0] = d[0] * factor;
76 r[1] = d[1] * factor;
77 r[2] = d[2] * factor;
78 }
79
80
81
82
83 Explicit_initial_data::
84 Explicit_initial_data(const Esfem::Grid::Grid_and_time& gt,
85 const Esfem::Growth type)
86 :tp {gt.time_provider()}
87 {
88 switch(type){
89 case Growth::promoting:
90 fun_impl =
91 [&tp = tp](const Domain& d, Range& r){
92 const double x = d[0];
93 const double y = d[1];
94 const double t = tp.time();
95 r = x * y * std::exp(-6. * t);
96 };
97 break;
98 case Growth::inhibiting:
99 fun_impl = [&tp = tp](const Domain& d, Range& r){
100 const double y = d[1];
101 const double z = d[2];
102 const double t = tp.time();
103 r = y * z * std::exp(-6. * t);
104 };
105 break;
106 default:
107 throw InitData_error {Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
108 break;
109 };
110 }
111
112
113
114
115 Random_initial_data::
116 Random_initial_data(const Esfem::Io::Parameter& p,
117 const Esfem::Growth type)
118 :Random_initial_data {hom_value(p, type), pertubation(p, type)}
119 {
120 std::cout << print_configuration(p, type) << std::endl;
121 }
122 Random_initial_data::
123 Random_initial_data(const double hom_value,
124 const double pertubation)
125 :random_fun {std::bind(Random_dist {hom_value, hom_value + pertubation},
126 Random_engine {})}
127 {}
128
129
130
131
132 sphere_1EF::sphere_1EF(const Grid::Grid_and_time& gt)
133 :tp {gt.time_provider()} {}
134 sphere_1EF* sphere_1EF::clone(){
135 return new sphere_1EF {*this};
136 }
137 void sphere_1EF::interpolate(Grid::Scal_FEfun& fef) const{
138 using Fef = Grid::Scal_FEfun::Dune_FEfun;
139 Dune::LagrangeInterpolation<Fef>::interpolateFunction(*this, fef);
140 }
141 void sphere_1EF::evaluate(const domain& x, range& y)const{
142 y = x[0]*x[1]*exp(-6*tp.time());
143 }
144
145
146
147
148 sphere_eigenFun::sphere_eigenFun(const Grid::Grid_and_time& gt)
149 :tp {gt.time_provider()} {}
150 sphere_eigenFun* sphere_eigenFun::clone(){
151 return new sphere_eigenFun {*this};
152 }
153 void sphere_eigenFun::interpolate(Grid::Vec_FEfun& vfef) const{
154 using Vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
155 Dune::LagrangeInterpolation<Vfef>::interpolateFunction(*this, vfef);
156 }
157 void sphere_eigenFun::evaluate(const Domain& x, Range& y) const{
158 y[0] = x[0]*x[1];
159 y[1] = x[1]*x[2];
160 y[2] = x[0]*x[2];
161 y *= exp(-6.*tp.time());
162 }
163
164
165
166
167 sphere_mcf_sol::sphere_mcf_sol(const Grid::Grid_and_time& gt)
168 :tp {gt.time_provider()} {}
169 sphere_mcf_sol* sphere_mcf_sol::clone(){
170 return new sphere_mcf_sol {*this};
171 }
172 void sphere_mcf_sol::interpolate(Grid::Vec_FEfun& rhs) const{
173 using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
174 Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
175 }
176 void sphere_mcf_sol::evaluate(const Domain& x, Range& y) const{
177 auto norm = 0.;
178 for(int i = 0; i < Domain::dimension; ++i) norm += x[i]*x[i];
179 norm = sqrt(norm);
180 const auto rt = sqrt(1. - 2 * 2 * tp.time());
181 y = x;
182 y *= rt / norm;
183 }
184
185
186
187
188 sls_iData::sls_iData(const Grid::Grid_and_time& gt)
189 :tp {gt.time_provider()},
190 rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
191 rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
192 k {Parameter::getValue<double>("logistic_growth.steepness", .5)}
193 {}
194 void sls_iData::interpolate(Grid::Vec_FEfun& rhs) const{
195 using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
196 Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
197 }
198 void sls_iData::evaluate(const Domain& x, Range& y) const{
199 const auto
200 norm = sqrt(inner_product(&x[0], &x[0]+Domain::dimension, &x[0], 0.)),
201 ekt = exp(-k*tp.time()),
202 lgf = rE*rA/(rE * ekt + rA * (1 - ekt));
203 y = x;
204 y *= lgf / norm;
205 }
206
207 sls_v_iData::sls_v_iData(const Grid::Grid_and_time& gt)
208 :tp {gt.time_provider()},
209 rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
210 rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
211 k {Parameter::getValue<double>("logistic_growth.steepness", .5)}
212 {}
213 void sls_v_iData::interpolate(Grid::Vec_FEfun& rhs) const{
214 using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
215 Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
216 }
217 void sls_v_iData::evaluate(const Domain& x, Range& y) const{
218 const auto
219 ekt = exp(-k*tp.time()),
220 rt = rE*rA/(rE * ekt + rA * (1. - ekt));
221 y = x;
222 y *= k * (1. - rt/ rE);
223 }
224
225 sd_iData::sd_iData(const Grid::Grid_and_time& gt) :tp {gt.time_provider()} {}
226 void sd_iData::interpolate(Grid::Vec_FEfun& rhs) const{
227 using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
228 Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
229 }
230 void sd_iData::evaluate(const Domain& x, Range& y) const{
231 const auto
232 norm = sqrt(inner_product(&x[0], &x[0]+Domain::dimension, &x[0], 0.)),
233 fac = exp(tp.time()) / norm;
234 y = x;
235 y *= fac;
236 }
237
238
239
240 Analytic_velocity::Analytic_velocity(const Esfem::Grid::Grid_and_time& gt)
241 :tp {gt.time_provider()}
242 {}
243
244 void Analytic_velocity::evaluate(const Domain& d, Range& r) const{
245 const double t = tp.time();
246 eoc_velocity(t, d, r);
247 }
248
249
250
251
252 double Esfem::Impl::
253 hom_value(const Esfem::Io::Parameter& p, const Esfem::Growth type){
254 double rv = 0.;
255 switch(type){
256 case Esfem::Growth::promoting:
257 rv = p.u_hom_value();
258 break;
259 case Esfem::Growth::inhibiting:
260 rv = p.w_hom_value();
261 break;
262 default:
263 throw InitData_error{Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
264 break;
265 };
266 return rv;
267 }
268
269 double Esfem::Impl::
270 pertubation(const Esfem::Io::Parameter& p, const Esfem::Growth type){
271 double rv = 0.;
272 switch(type){
273 case Esfem::Growth::promoting:
274 rv = p.u_pertubation();
275 break;
276 case Esfem::Growth::inhibiting:
277 rv = p.w_pertubation();
278 break;
279 default:
280 throw InitData_error{Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
281 break;
282 };
283 return rv;
284 }
285
286 std::string Esfem::Impl::print_configuration(const Esfem::Io::Parameter& p,
287 const Esfem::Growth type){
288 std::ostringstream oss;
289 oss << "Using Random_initial_data():\n"
290 "Random distribution: uniform_real_distribution<> with\n"
291 << "hom_value: " << hom_value(p, type) << '\n'
292 << "pertubation: " << pertubation(p, type);
293 return oss.str();
294 }
295
296 std::string Esfem::Impl::dof_filename(const Io::Parameter& p, const Growth type){
297 std::string rv {};
298 switch(type){
299 case Growth::promoting:
300 rv = p.u_init_dof();
301 break;
302 case Growth::inhibiting:
303 rv = p.w_init_dof();
304 break;
305 default:
306 throw InitData_error {Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
307 break;
308 };
309 return rv;
310 }
311
312
313
314
315 Esfem::SecOrd_op::Init_data::Data::Data(const Grid::Grid_and_time& gt,
316 const Growth type)
317 :eid_ptr {std::make_unique<Explicit_initial_data>(gt,type)}
318 {}
319
320 Esfem::SecOrd_op::Init_data::Data::Data(const Io::Parameter& p, const Growth type)
321 :dof_io_filename {Impl::dof_filename(p, type)},
322 rid_ptr {std::make_unique<Random_initial_data>(p, type)}
323 {}
324
325
326
327
328 Esfem::SecOrd_op::Exact_velocity::Data::Data(const Grid::Grid_and_time& gt)
329 :v_fun {gt}
330 {}