intrinsics_propagation.cpp 2.9 KB
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
/***************************************
	Auteur : Pierre Aubert
	Mail : aubertp7@gmail.com
	Licence : CeCILL-C
****************************************/

#include "phoenix_intrinsics.h"
#include <algorithm>
#include "intrinsics_propagation.h"


///Propagate the U and V species in the matVecVecU and matVecV
/**	@param[out] outMatVecVecU : updated matrix U version (with vectorial neighbours)
 * 	@param[out] outMatVecV : updated matrix V version (with vectorial neighbours)
 * 	@param matVecVecU : input of matrix U (with vectorial neighbours)
 * 	@param matVecV : input of matrix V (with vectorial neighbours)
 * 	@param nbRow : number of rows of the matrices
 * 	@param nbCol : number of columns of the matrices
 * 	@param padding : padding of the columns of all matrices
 * 	@param matBroadcastDeltaSquare : matrix of the delta square values (with broadcast neighbours)
 * 	@param nbStencilRow : number of rows of the matrix matBroadcastDeltaSquare
 * 	@param nbStencilCol : number of columns of the matrix matBroadcastDeltaSquare
 * 	@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 * outMatVecVecU, float * outMatVecV, const float * matVecVecU, const float * matVecV, long nbRow, long nbCol, long padding,
		       const float * matBroadcastDeltaSquare, long nbStencilRow, long nbStencilCol,
		       float diffudionRateU, float diffusionRateV, float feedRate, float killRate, float dt)
{
	long offsetStencilRow((nbStencilRow - 1l)/2l);
	long offsetStencilCol((nbStencilCol - 1l)/2l);
	
	size_t nbVecCol(nbCol/PLIB_VECTOR_SIZE_FLOAT);
	
	
	for(long i(0l); i < nbRow; ++i){
		long firstRowStencil(std::max(i - offsetStencilRow, 0l));
		long lastRowStencil(std::min(i + offsetStencilRow + 1l, nbRow));
		for(long j(0l); j < nbCol; ++j){
			long firstColStencil(std::max(j - offsetStencilCol, 0l));
			long lastColStencil(std::min(j + offsetStencilCol + 1l, nbCol));
			
			long stencilIndexRow(0l);
			float u(matVecVecU[i*nbCol + j]), v(matVecV[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(matBroadcastDeltaSquare[stencilIndexRow*nbStencilCol + stencilIndexCol]);
					fullU += (matVecVecU[k*nbCol + l] - u)*deltaSquare;
					fullV += (matVecV[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);
			
			outMatVecVecU[i*nbCol + j] = u + du*dt;
			outMatVecV[i*nbCol + j] = v + dv*dt;
		}
	}
}