Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#include <vector>
#include <iostream>
#include <iomanip>
#include <cassert>
#include <numeric>
using real = double ;
const real KILL_RATE { 0.062f };
const real FEED_RATE { 0.03f };
const real DT { 1.0f };
const real DIFFUSION_RATE_U { 0.1f };
const real DIFFUSION_RATE_V { 0.05f };
template<typename Collection>
void drop( Collection & data, std::size_t nb_rows, std::size_t nb_cols, real value ) {
const std::size_t center_row { nb_rows/2ul };
const std::size_t center_col { nb_cols/2ul };
for (std::size_t i = (center_row-2ul) ; i < (center_row+2ul); i++) {
for (std::size_t j = (center_col-3ul); j < (center_col+3ul); j++) {
data[i*nb_cols+j] = value ;
}
}
}
void process(
std::vector<real> const & iu, std::vector<real> const & iv,
std::vector<real> & ou, std::vector<real> & ov,
std::size_t nb_rows, std::size_t nb_cols ) {
for (std::size_t r = 0 ; r < (nb_rows-2); r++) {
for (std::size_t c = 0; c < (nb_cols-2); c++) {
real u = iu[(r+1)*nb_cols+c+1];
real v = iv[(r+1)*nb_cols+c+1];
real uvv = u*v*v;
real full_u = 0.0f;
real full_v = 0.0f;
for(std::size_t k = 0ul; k < 3ul; ++k){
for(std::size_t l = 0ul; l < 3ul; ++l){
full_u += iu[(r+k)*nb_cols+c+l] - u;
full_v += iv[(r+k)*nb_cols+c+l] - v;
}
}
real du = DIFFUSION_RATE_U*full_u - uvv + FEED_RATE*(1.0f - u);
real dv = DIFFUSION_RATE_V*full_v + uvv - (FEED_RATE + KILL_RATE)*v;
ou[(r+1)*nb_cols+c+1] = u + du*DT;
ov[(r+1)*nb_cols+c+1] = v + dv*DT;
}
}
}
int main( int argc, char * argv[] ) {
// runtime parameters
assert(argc==5) ;
std::size_t nb_rows {std::stoul(argv[1])};
std::size_t nb_cols {std::stoul(argv[2])};
std::size_t nb_images {std::stoul(argv[3])};
std::size_t nb_iterations {std::stoul(argv[4])};
std::cout << std::fixed << std::setprecision(6);
assert(nb_images<=1000);
try {
// Temporary work images
std::vector<real> u1(nb_rows*nb_cols,1.f);
std::vector<real> v1(nb_rows*nb_cols,0.f);
std::vector<real> u2(nb_rows*nb_cols,1.f);
std::vector<real> v2(nb_rows*nb_cols,0.f);
drop(u1, nb_rows, nb_cols, 0.f);
drop(v1, nb_rows, nb_cols, 1.f);
// Final images sequence
std::vector<std::vector<real>> images;
images.push_back(u1);
// Iterations
std::cout<<std::endl ;
for ( std::size_t image = 0 ; image < nb_images ; ++image ) {
for ( std::size_t iter = 0 ; iter < nb_iterations ; ++iter ) {
process( u1, v1, u2, v2, nb_rows, nb_cols );
std::swap(u1,u2) ; std::swap (v1,v2) ;
}
std::cout << '.' << std::flush;
images.push_back(u1);
}
std::cout << std::endl;
// Final product
std::cout << std::endl;
real product = std::inner_product(u1.begin(), u1.end(), v1.begin(),0.f);
std::cout
<< "product : " << product
<< std::endl;
}
catch (std::exception & e) {
std::cout << e.what() << std::endl;
}
catch (const char * e) {
std::cout << e << std::endl;
}
// End
return 0;
}