Split into compilation units

9 files changed, 281 insertions, 240 deletions

diff --git a/src/boundary_conditions.cc b/src/boundary_conditions.cc
new file mode 100644
index 0000000..3488463
--- /dev/null
+++ b/src/boundary_conditions.cc
@@ -0,0 +1,45 @@
+#include "boundary_conditions.h"
+
+std::pair<Vector<int>, Vector<int>> neighbors(Vector<int> v) {
+	if ( v[0] == 0 ) {
+		return {
+			{ -1, v[1] },
+			{  1, v[1] }
+		};
+	} else if ( v[1] == 0 ) {
+		return {
+			{ v[0], -1 },
+			{ v[0],  1 }
+		};
+	} else {
+		return {
+			{    0, v[1] },
+			{ v[0], 0    }
+		};
+	}
+}
+
+void computeWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal) {
+	const auto [neighborA, neighborB] = neighbors(normal);
+
+	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA);
+	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
+	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB);
+}
+
+void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, double vX) {
+	const auto [neighborA, neighborB] = neighbors(normal);
+
+	const double rho = pop.curr(cell).get(-1,0) + pop.curr(cell).get(0,0) + pop.curr(cell).get(1,0)
+		+ 2.*(
+			pop.curr(cell).get(-neighborA) +
+			pop.curr(cell).get(-normal   ) +
+			pop.curr(cell).get(-neighborB)
+		);
+
+	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA)
+	                              + 0.5*( pop.curr(cell).get( 1,0) - pop.curr(cell).get(-1,0) - vX*rho );
+	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
+	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB)
+	                              + 0.5*( pop.curr(cell).get(-1,0) - pop.curr(cell).get( 1,0) + vX*rho );
+}
diff --git a/src/boundary_conditions.h b/src/boundary_conditions.h
index c825c32..6b7de19 100644
--- a/src/boundary_conditions.h
+++ b/src/boundary_conditions.h
@@ -3,46 +3,6 @@
 #include "vector.h"
 #include "data_cell_buffer.h"
 
-std::pair<Vector<int>, Vector<int>> neighbors(Vector<int> v) {
-	if ( v[0] == 0 ) {
-		return {
-			{ -1, v[1] },
-			{  1, v[1] }
-		};
-	} else if ( v[1] == 0 ) {
-		return {
-			{ v[0], -1 },
-			{ v[0],  1 }
-		};
-	} else {
-		return {
-			{    0, v[1] },
-			{ v[0], 0    }
-		};
-	}
-}
+void computeWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal);
 
-void computeWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal) {
-	const auto [neighborA, neighborB] = neighbors(normal);
-
-	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA);
-	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
-	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB);
-}
-
-void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, double vX) {
-	const auto [neighborA, neighborB] = neighbors(normal);
-
-	const double rho = pop.curr(cell).get(-1,0) + pop.curr(cell).get(0,0) + pop.curr(cell).get(1,0)
-		+ 2.*(
-			pop.curr(cell).get(-neighborA) +
-			pop.curr(cell).get(-normal   ) +
-			pop.curr(cell).get(-neighborB)
-		);
-
-	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA)
-	                              + 0.5*( pop.curr(cell).get( 1,0) - pop.curr(cell).get(-1,0) - vX*rho );
-	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
-	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB)
-	                              + 0.5*( pop.curr(cell).get(-1,0) - pop.curr(cell).get( 1,0) + vX*rho );
-}
+void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, double vX);
diff --git a/src/box_obstacle.cc b/src/box_obstacle.cc
new file mode 100644
index 0000000..9cb6ad0
--- /dev/null
+++ b/src/box_obstacle.cc
@@ -0,0 +1,28 @@
+#include "box_obstacle.h"
+
+#include "boundary_conditions.h"
+
+BoxObstacle::BoxObstacle(std::size_t lX, std::size_t lY, std::size_t uX, std::size_t uY):
+	lower_x_(lX), lower_y_(lY), upper_x_(uX), upper_y_(uY) { }
+
+bool BoxObstacle::isInside(std::size_t x, std::size_t y) const {
+	return x > lower_x_
+		&& x < upper_x_
+		&& y > lower_y_
+		&& y < upper_y_;
+}
+
+void BoxObstacle::applyBoundary(DataCellBuffer& pop) const {
+	for ( std::size_t x = lower_x_+1; x < upper_x_; ++x ) {
+		computeWallCell(pop, {x, lower_y_}, { 0,-1});
+		computeWallCell(pop, {x, upper_y_}, { 0, 1});
+	}
+	for ( std::size_t y = lower_y_+1; y < upper_y_; ++y ) {
+		computeWallCell(pop, {lower_x_, y}, {-1, 0});
+		computeWallCell(pop, {upper_x_, y}, { 1, 0});
+	}
+	computeWallCell(pop, {lower_x_, lower_y_}, {-1,-1});
+	computeWallCell(pop, {upper_x_, lower_y_}, { 1,-1});
+	computeWallCell(pop, {upper_x_, upper_y_}, { 1, 1});
+	computeWallCell(pop, {lower_x_, upper_y_}, {-1, 1});
+}
diff --git a/src/box_obstacle.h b/src/box_obstacle.h
index 46660a7..149e8e1 100644
--- a/src/box_obstacle.h
+++ b/src/box_obstacle.h
@@ -1,7 +1,6 @@
 #pragma once
 
 #include "data_cell_buffer.h"
-#include "boundary_conditions.h"
 
 struct BoxObstacle {
 	const std::size_t lower_x_;
@@ -9,28 +8,9 @@ struct BoxObstacle {
 	const std::size_t upper_x_;
 	const std::size_t upper_y_;
 
-	BoxObstacle(std::size_t lX, std::size_t lY, std::size_t uX, std::size_t uY):
-		lower_x_(lX), lower_y_(lY), upper_x_(uX), upper_y_(uY) { }
+	BoxObstacle(std::size_t lX, std::size_t lY, std::size_t uX, std::size_t uY);
 
-	bool isInside(std::size_t x, std::size_t y) const {
-		return x > lower_x_
-		    && x < upper_x_
-		    && y > lower_y_
-		    && y < upper_y_;
-	}
+	bool isInside(std::size_t x, std::size_t y) const;
 
-	void applyBoundary(DataCellBuffer& pop) const {
-		for ( std::size_t x = lower_x_+1; x < upper_x_; ++x ) {
-			computeWallCell(pop, {x, lower_y_}, { 0,-1});
-			computeWallCell(pop, {x, upper_y_}, { 0, 1});
-		}
-		for ( std::size_t y = lower_y_+1; y < upper_y_; ++y ) {
-			computeWallCell(pop, {lower_x_, y}, {-1, 0});
-			computeWallCell(pop, {upper_x_, y}, { 1, 0});
-		}
-		computeWallCell(pop, {lower_x_, lower_y_}, {-1,-1});
-		computeWallCell(pop, {upper_x_, lower_y_}, { 1,-1});
-		computeWallCell(pop, {upper_x_, upper_y_}, { 1, 1});
-		computeWallCell(pop, {lower_x_, upper_y_}, {-1, 1});
-	}
+	void applyBoundary(DataCellBuffer& pop) const;
 };
diff --git a/src/data_cell_buffer.h b/src/data_cell_buffer.h
index 06c2644..033fb24 100644
--- a/src/data_cell_buffer.h
+++ b/src/data_cell_buffer.h
@@ -5,46 +5,46 @@
 #include "data_cell.h"
 
 class DataCellBuffer {
-	private:
-		const std::size_t dim_x_;
-		const std::size_t dim_y_;
-
-		std::unique_ptr<DataCell[]> curr_;
-		std::unique_ptr<DataCell[]> prev_;
-
-	public:
-		DataCellBuffer(std::size_t dimX, std::size_t dimY):
-			dim_x_(dimX),
-			dim_y_(dimY),
-			curr_(new DataCell[dimX*dimY]),
-			prev_(new DataCell[dimX*dimY]) { }
-
-		std::size_t dimX() const {
-			return dim_x_;
-		}
-
-		std::size_t dimY() const {
-			return dim_y_;
-		}
-
-		void swap() {
-			curr_.swap(prev_);
-		}
-
-		inline DataCell& curr(std::size_t x, std::size_t y) {
-			return curr_[y*dim_x_ + x];
-		}
-
-		inline DataCell& curr(Vector<std::size_t> pos) {
-			return curr(pos[0], pos[1]);
-		}
-
-		inline DataCell& prev(std::size_t x, std::size_t y) {
-			return prev_[y*dim_x_ + x];
-		}
-
-		inline DataCell& prev(Vector<std::size_t> pos) {
-			return prev(pos[0], pos[1]);
-		}
+private:
+	const std::size_t dim_x_;
+	const std::size_t dim_y_;
+
+	std::unique_ptr<DataCell[]> curr_;
+	std::unique_ptr<DataCell[]> prev_;
+
+public:
+	DataCellBuffer(std::size_t dimX, std::size_t dimY):
+		dim_x_(dimX),
+		dim_y_(dimY),
+		curr_(new DataCell[dimX*dimY]),
+		prev_(new DataCell[dimX*dimY]) { }
+
+	std::size_t dimX() const {
+		return dim_x_;
+	}
+
+	std::size_t dimY() const {
+		return dim_y_;
+	}
+
+	void swap() {
+		curr_.swap(prev_);
+	}
+
+	inline DataCell& curr(std::size_t x, std::size_t y) {
+		return curr_[y*dim_x_ + x];
+	}
+
+	inline DataCell& curr(Vector<std::size_t> pos) {
+		return curr(pos[0], pos[1]);
+	}
+
+	inline DataCell& prev(std::size_t x, std::size_t y) {
+		return prev_[y*dim_x_ + x];
+	}
+
+	inline DataCell& prev(Vector<std::size_t> pos) {
+		return prev(pos[0], pos[1]);
+	}
 };
 
diff --git a/src/fluid_buffer.cc b/src/fluid_buffer.cc
new file mode 100644
index 0000000..e73c818
--- /dev/null
+++ b/src/fluid_buffer.cc
@@ -0,0 +1,75 @@
+#include "fluid_buffer.h"
+
+#include <fstream>
+
+FluidBuffer::FluidBuffer(std::size_t dimX, std::size_t dimY):
+	dim_x_(dimX),
+	dim_y_(dimY),
+	density_(new Density[dimX*dimY]),
+	velocity_(new Velocity[dimX*dimY]) { }
+
+std::size_t FluidBuffer::dimX() const {
+	return dim_x_;
+}
+
+std::size_t FluidBuffer::dimY() const {
+	return dim_y_;
+}
+
+Density& FluidBuffer::density(std::size_t x, std::size_t y) {
+	return density_[y*dim_x_ + x];
+}
+
+Density& FluidBuffer::density(Vector<std::size_t> pos) {
+	return density(pos[0], pos[1]);
+}
+
+Velocity& FluidBuffer::velocity(std::size_t x, std::size_t y) {
+	return velocity_[y*dim_x_ + x];
+}
+
+Velocity& FluidBuffer::velocity(Vector<std::size_t> pos) {
+	return velocity(pos[0], pos[1]);
+}
+
+void FluidBuffer::writeAsVTK(int time) {
+	std::ofstream fout;
+	fout.open(("result/data_t" + std::to_string(time) + ".vtk").c_str());
+
+	fout << "# vtk DataFile Version 3.0\n";
+	fout << "lbm_output\n";
+	fout << "ASCII\n";
+	fout << "DATASET RECTILINEAR_GRID\n";
+	fout << "DIMENSIONS " << dimX() << " " << dimY() << " 1" << "\n";
+
+	fout << "X_COORDINATES " << dimX() << " float\n";
+	for( std::size_t x = 0; x < dimX(); ++x ) {
+		fout << x << " ";
+	}
+
+	fout << "\nY_COORDINATES " << dimY() << " float\n";
+	for( std::size_t y = 0; y < dimY(); ++y ) {
+		fout << y << " ";
+	}
+
+	fout << "\nZ_COORDINATES " << 1 << " float\n";
+	fout << 0 << "\n";
+	fout << "POINT_DATA " << dimX() * dimY() << "\n";
+
+	fout << "VECTORS velocity float\n";
+	for ( std::size_t y = 0; y < dimY(); ++y ) {
+		for ( std::size_t x = 0; x < dimX(); ++x ) {
+			fout << velocity(x,y)[0] << " " << velocity(x,y)[1] << " 0\n";
+		}
+	}
+
+	fout << "SCALARS density float 1\n";
+	fout << "LOOKUP_TABLE default\n";
+	for ( std::size_t y = 0; y < dimY(); ++y ) {
+		for ( std::size_t x = 0; x < dimX(); ++x ) {
+			fout << density(x,y) << "\n";
+		}
+	}
+
+	fout.close();
+}
diff --git a/src/fluid_buffer.h b/src/fluid_buffer.h
index a6c41ba..1d802c0 100644
--- a/src/fluid_buffer.h
+++ b/src/fluid_buffer.h
@@ -1,6 +1,5 @@
 #pragma once
 
-#include <fstream>
 #include <memory>
 
 #include "vector.h"
@@ -9,83 +8,24 @@ using Velocity = Vector<double>;
 using Density  = double;
 
 class FluidBuffer {
-	private:
-		const std::size_t dim_x_;
-		const std::size_t dim_y_;
+private:
+	const std::size_t dim_x_;
+	const std::size_t dim_y_;
 
-		std::unique_ptr<Density[]> density_;
-		std::unique_ptr<Velocity[]> velocity_;
+	std::unique_ptr<Density[]> density_;
+	std::unique_ptr<Velocity[]> velocity_;
 
-	public:
-		FluidBuffer(std::size_t dimX, std::size_t dimY):
-			dim_x_(dimX),
-			dim_y_(dimY),
-			density_(new Density[dimX*dimY]),
-			velocity_(new Velocity[dimX*dimY]) { }
+public:
+	FluidBuffer(std::size_t dimX, std::size_t dimY);
 
-		std::size_t dimX() const {
-			return dim_x_;
-		}
+	std::size_t dimX() const;
+	std::size_t dimY() const;
 
-		std::size_t dimY() const {
-			return dim_y_;
-		}
+	Density& density(std::size_t x, std::size_t y);
+	Density& density(Vector<std::size_t> pos);
 
-		Density& density(std::size_t x, std::size_t y) {
-			return density_[y*dim_x_ + x];
-		}
+	Velocity& velocity(std::size_t x, std::size_t y);
+	Velocity& velocity(Vector<std::size_t> pos);
 
-		Density& curr(Vector<std::size_t> pos) {
-			return density(pos[0], pos[1]);
-		}
-
-		Velocity& velocity(std::size_t x, std::size_t y) {
-			return velocity_[y*dim_x_ + x];
-		}
-
-		Velocity& velocity(Vector<std::size_t> pos) {
-			return velocity(pos[0], pos[1]);
-		}
-
-		void writeAsVTK(int time) {
-			std::ofstream fout;
-			fout.open(("result/data_t" + std::to_string(time) + ".vtk").c_str());
-
-			fout << "# vtk DataFile Version 3.0\n";
-			fout << "lbm_output\n";
-			fout << "ASCII\n";
-			fout << "DATASET RECTILINEAR_GRID\n";
-			fout << "DIMENSIONS " << dimX() << " " << dimY() << " 1" << "\n";
-
-			fout << "X_COORDINATES " << dimX() << " float\n";
-			for( std::size_t x = 0; x < dimX(); ++x ) {
-				fout << x << " ";
-			}
-
-			fout << "\nY_COORDINATES " << dimY() << " float\n";
-			for( std::size_t y = 0; y < dimY(); ++y ) {
-				fout << y << " ";
-			}
-
-			fout << "\nZ_COORDINATES " << 1 << " float\n";
-			fout << 0 << "\n";
-			fout << "POINT_DATA " << dimX() * dimY() << "\n";
-
-			fout << "VECTORS velocity float\n";
-			for ( std::size_t y = 0; y < dimY(); ++y ) {
-				for ( std::size_t x = 0; x < dimX(); ++x ) {
-					fout << velocity(x,y)[0] << " " << velocity(x,y)[1] << " 0\n";
-				}
-			}
-
-			fout << "SCALARS density float 1\n";
-			fout << "LOOKUP_TABLE default\n";
-			for ( std::size_t y = 0; y < dimY(); ++y ) {
-				for ( std::size_t x = 0; x < dimX(); ++x ) {
-					fout << density(x,y) << "\n";
-				}
-			}
-
-			fout.close();
-		}
+	void writeAsVTK(int time);
 };
diff --git a/src/lbm.cc b/src/lbm.cc
new file mode 100644
index 0000000..ca84cee
--- /dev/null
+++ b/src/lbm.cc
@@ -0,0 +1,68 @@
+#include "lbm.h"
+
+void Cell::equilibrize(Density d, Velocity v) {
+	for ( int i = -1; i <= 1; ++i ) {
+		for ( int j = -1; j <= 1; ++j ) {
+			get(i,j) = weight.get(i,j) * d * (1 + 3*v.comp(i,j) + 4.5*sq(v.comp(i,j)) - 1.5*sq(v.norm()));
+		}
+	}
+}
+
+double Cell::sum() const {
+	return get(-1, 1) + get( 0, 1) + get( 1, 1) + get(-1, 0) + get( 0, 0) + get( 1, 0) + get(-1,-1) + get( 0,-1) + get( 1,-1);
+}
+
+Velocity Cell::velocity(Density d) const {
+	return 1./d * Velocity{
+		get( 1, 0) - get(-1, 0) + get( 1, 1) - get(-1,-1) + get( 1,-1) - get(-1,1),
+		get( 0, 1) - get( 0,-1) + get( 1, 1) - get(-1,-1) - get( 1,-1) + get(-1,1)
+	};
+}
+
+void streamFluidCell(DataCellBuffer& pop, std::size_t x, std::size_t y) {
+	if ( x != 0 && x != pop.dimX()-1 && y != 0 && y != pop.dimY()-1 ) {
+		// stream internal cells
+		for ( int i = -1; i <= 1; ++i ) {
+			for ( int j = -1; j <= 1; ++j ) {
+				pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
+			}
+		}
+	} else {
+		// stream boundary cells,
+		// missing populations to be determined by boundary conditions
+		for ( int i = -1; i <= 1; ++i ) {
+			for ( int j = -1; j <= 1; ++j ) {
+				if ( (x > 0 || i >= 0) && x+i <= pop.dimX()-1 && (y > 0 || j >= 0) && y+j <= pop.dimY()-1 ) {
+					pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
+				}
+			}
+		}
+	}
+}
+
+void collideFluidCell(
+	double omega,
+	DataCellBuffer& pop, FluidBuffer& fluid,
+	std::size_t x, std::size_t y) {
+	// compute equilibrium
+	Cell eq;
+	eq.equilibrize(fluid.density(x,y), fluid.velocity(x,y));
+
+	// collide (BGK, relax towards equilibrium)
+	if ( x != 0 && x != pop.dimX()-1 && y != 0 && y != pop.dimY()-1 ) {
+		for ( int i = -1; i <= 1; ++i ) {
+			for ( int j = -1; j <= 1; ++j ) {
+				pop.curr(x,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j));
+			}
+		}
+	} else {
+		// partial collide for boundary cells
+		for ( int i = -1; i <= 1; ++i ) {
+			for ( int j = -1; j <= 1; ++j ) {
+				if ( (x > 0 || i >= 0) && x+i <= pop.dimX()-1 && (y > 0 || j >= 0) && y+j <= pop.dimY()-1 ) {
+					pop.curr(x,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j));
+				}
+			}
+		}
+	}
+}
diff --git a/src/lbm.h b/src/lbm.h
index 685fa08..96348c1 100644
--- a/src/lbm.h
+++ b/src/lbm.h
@@ -11,71 +11,16 @@ constexpr DataCell weight{
 };
 
 struct Cell : DataCell {
-	void equilibrize(Density d, Velocity v) {
-		for ( int i = -1; i <= 1; ++i ) {
-			for ( int j = -1; j <= 1; ++j ) {
-				get(i,j) = weight.get(i,j) * d * (1 + 3*v.comp(i,j) + 4.5*sq(v.comp(i,j)) - 1.5*sq(v.norm()));
-			}
-		}
-	}
+	void equilibrize(Density d, Velocity v);
 
-	double sum() const {
-		return get(-1, 1) + get( 0, 1) + get( 1, 1) + get(-1, 0) + get( 0, 0) + get( 1, 0) + get(-1,-1) + get( 0,-1) + get( 1,-1);
-	}
+	double sum() const;
 
-	Velocity velocity(Density d) const {
-		return 1./d * Velocity{
-			get( 1, 0) - get(-1, 0) + get( 1, 1) - get(-1,-1) + get( 1,-1) - get(-1,1),
-			get( 0, 1) - get( 0,-1) + get( 1, 1) - get(-1,-1) - get( 1,-1) + get(-1,1)
-		};
-	}
+	Velocity velocity(Density d) const;
 };
 
-void streamFluidCell(DataCellBuffer& pop, std::size_t x, std::size_t y) {
-	if ( x != 0 && x != pop.dimX()-1 && y != 0 && y != pop.dimY()-1 ) {
-		// stream internal cells
-		for ( int i = -1; i <= 1; ++i ) {
-			for ( int j = -1; j <= 1; ++j ) {
-				pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
-			}
-		}
-	} else {
-		// stream boundary cells,
-		// missing populations to be determined by boundary conditions
-		for ( int i = -1; i <= 1; ++i ) {
-			for ( int j = -1; j <= 1; ++j ) {
-				if ( (x > 0 || i >= 0) && x+i <= pop.dimX()-1 && (y > 0 || j >= 0) && y+j <= pop.dimY()-1 ) {
-					pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
-				}
-			}
-		}
-	}
-}
+void streamFluidCell(DataCellBuffer& pop, std::size_t x, std::size_t y);
 
 void collideFluidCell(
 	double omega,
 	DataCellBuffer& pop, FluidBuffer& fluid,
-	std::size_t x, std::size_t y) {
-	// compute equilibrium
-	Cell eq;
-	eq.equilibrize(fluid.density(x,y), fluid.velocity(x,y));
-
-	// collide (BGK, relax towards equilibrium)
-	if ( x != 0 && x != pop.dimX()-1 && y != 0 && y != pop.dimY()-1 ) {
-		for ( int i = -1; i <= 1; ++i ) {
-			for ( int j = -1; j <= 1; ++j ) {
-				pop.curr(x,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j));
-			}
-		}
-	} else {
-		// partial collide for boundary cells
-		for ( int i = -1; i <= 1; ++i ) {
-			for ( int j = -1; j <= 1; ++j ) {
-				if ( (x > 0 || i >= 0) && x+i <= pop.dimX()-1 && (y > 0 || j >= 0) && y+j <= pop.dimY()-1 ) {
-					pop.curr(x,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j));
-				}
-			}
-		}
-	}
-}
-
+	std::size_t x, std::size_t y);