4.1.1.3 : La source naive_propagation.cpp
Écrivons le fichier naive_propagation.cpp :
Nous aurons besoin des fonctions
std::min et
std::max qui sont disponible dans
algorithm :
Ensuite, nous incluons notre header :
1
|
#include "naive_propagation.h"
|
Et nous commençons l'implementation de notre fonction :
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
|
///Propagate the U and V species in the matU and matV
/** @param[out] outMatU : updated matrix U version
* @param[out] outMatV : updated matrix V version
* @param matU : input of matrix U
* @param matV : input of matrix V
* @param nbRow : number of rows of the matrices
* @param nbCol : number of columns of the matrices
* @param matDeltaSquare : matrix of the delta square values
* @param nbStencilRow : number of rows of the matrix matDeltaSquare
* @param nbStencilCol : number of columns of the matrix matDeltaSquare
* @param diffudionRateU : diffusion rate of the U specie
* @param diffudionRateV : diffusion rate of the V specie
* @param feedRate : rate of the process which feeds U and drains U, V and P
* @param killRate : rate of the process which converts V into P
* @param dt : time interval between two steps
*/
void grayscott_propagation(float * outMatU, float * outMatV, const float * matU, const float * matV, long nbRow, long nbCol,
const float * matDeltaSquare, long nbStencilRow, long nbStencilCol,
float diffudionRateU, float diffusionRateV, float feedRate, float killRate, float dt)
{
|
Nous déterminons les
offset de notre
stencil (le nombre de couches à partir de la cellule centrale) :
1
2
|
long offsetStencilRow((nbStencilRow - 1l)/2l);
long offsetStencilCol((nbStencilCol - 1l)/2l);
|
Nous bouclons sur les lignes et les colonnes de nos matrices pour mettre à jour toutes nos cellules :
1
2
|
for(long i(0l); i < nbRow; ++i){
for(long j(0l); j < nbCol; ++j){
|
Il faut maintenant déterminer les bornes de nos calculs (voir section
4.1.1.1) :
1
2
3
4
5
|
long firstRowStencil(std::max(i - offsetStencilRow, 0l));
long firstColStencil(std::max(j - offsetStencilCol, 0l));
long lastRowStencil(std::min(i + offsetStencilRow + 1l, nbRow));
long lastColStencil(std::min(j + offsetStencilCol + 1l, nbCol));
|
Définissons quelques variables temporaires :
1
2
3
|
long stencilIndexRow(0l);
float u(matU[i*nbCol + j]), v(matV[i*nbCol + j]);
float fullU(0.0f), fullV(0.0f);
|
Nous devons maintenant boucler sur les lignes et les colonnes de notre
stencil :
1
2
3
|
for(long k(firstRowStencil); k < lastRowStencil; ++k){
long stencilIndexCol(0l);
for(long l(firstColStencil); l < lastColStencil; ++l){
|
Nous pouvons enfin calculer notre gradient :
1
2
3
|
float deltaSquare(matDeltaSquare[stencilIndexRow*nbStencilCol + stencilIndexCol]);
fullU += (matU[k*nbCol + l] - u)*deltaSquare;
fullV += (matV[k*nbCol + l] - v)*deltaSquare;
|
Il ne faut pas oublier d'incrémenter les indices qui nous permettent de parcourir la matrice
matDeltaSquare] :
1
2
3
4
|
++stencilIndexCol;
}
++stencilIndexRow;
}
|
On finalise le calcul :
1
2
3
|
float uvSquare(u*v*v);
float du(diffudionRateU*fullU - uvSquare + feedRate*(1.0f - u));
float dv(diffusionRateV*fullV + uvSquare - (feedRate + killRate)*v);
|
Et on sauvegarde le résultat :
1
2
|
outMatU[i*nbCol + j] = u + du*dt;
outMatV[i*nbCol + j] = v + dv*dt;
|
Fin des deux boucles sur les lignes et les colonnes :
Fin de la fonction :
Le fichier
naive_propagation.cpp complet :
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
|
#include <algorithm>
#include "naive_propagation.h"
///Propagate the U and V species in the matU and matV
/** @param[out] outMatU : updated matrix U version
* @param[out] outMatV : updated matrix V version
* @param matU : input of matrix U
* @param matV : input of matrix V
* @param nbRow : number of rows of the matrices
* @param nbCol : number of columns of the matrices
* @param matDeltaSquare : matrix of the delta square values
* @param nbStencilRow : number of rows of the matrix matDeltaSquare
* @param nbStencilCol : number of columns of the matrix matDeltaSquare
* @param diffudionRateU : diffusion rate of the U specie
* @param diffudionRateV : diffusion rate of the V specie
* @param feedRate : rate of the process which feeds U and drains U, V and P
* @param killRate : rate of the process which converts V into P
* @param dt : time interval between two steps
*/
void grayscott_propagation(float * outMatU, float * outMatV, const float * matU, const float * matV, long nbRow, long nbCol,
const float * matDeltaSquare, long nbStencilRow, long nbStencilCol,
float diffudionRateU, float diffusionRateV, float feedRate, float killRate, float dt)
{
long offsetStencilRow((nbStencilRow - 1l)/2l);
long offsetStencilCol((nbStencilCol - 1l)/2l);
for(long i(0l); i < nbRow; ++i){
for(long j(0l); j < nbCol; ++j){
long firstRowStencil(std::max(i - offsetStencilRow, 0l));
long firstColStencil(std::max(j - offsetStencilCol, 0l));
long lastRowStencil(std::min(i + offsetStencilRow + 1l, nbRow));
long lastColStencil(std::min(j + offsetStencilCol + 1l, nbCol));
long stencilIndexRow(0l);
float u(matU[i*nbCol + j]), v(matV[i*nbCol + j]);
float fullU(0.0f), fullV(0.0f);
for(long k(firstRowStencil); k < lastRowStencil; ++k){
long stencilIndexCol(0l);
for(long l(firstColStencil); l < lastColStencil; ++l){
float deltaSquare(matDeltaSquare[stencilIndexRow*nbStencilCol + stencilIndexCol]);
fullU += (matU[k*nbCol + l] - u)*deltaSquare;
fullV += (matV[k*nbCol + l] - v)*deltaSquare;
++stencilIndexCol;
}
++stencilIndexRow;
}
float uvSquare(u*v*v);
float du(diffudionRateU*fullU - uvSquare + feedRate*(1.0f - u));
float dv(diffusionRateV*fullV + uvSquare - (feedRate + killRate)*v);
outMatU[i*nbCol + j] = u + du*dt;
outMatV[i*nbCol + j] = v + dv*dt;
}
}
}
|
Vous pouvez le télécharger
ici.