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 :

1

#include <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 :

1
2

		}
	}


Fin de la fonction :

1

}


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

/***************************************
	Auteur : Pierre Aubert
	Mail : pierre.aubert@lapp.in2p3.fr
	Licence : CeCILL-C
****************************************/

#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.