Bring in some more structure

8 files changed, 368 insertions, 285 deletions

diff --git a/boltz.cc b/boltz.cc
index e36e337..6f8788e 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -1,117 +1,10 @@
 #include <iostream>
-#include <numeric>
-#include <iostream>
-#include <fstream>
-#include <memory>
 #include <vector>
 #include <algorithm>
 
-#include "vector.h"
-
-struct DataCell {
-	double data[3][3];
-
-	inline double& get(int x, int y) {
-		return data[1+x][1-y];
-	}
-
-	inline double& get(Vector<int> v) {
-		return get(v[0], v[1]);
-	}
-
-	inline double get(int x, int y) const {
-		return data[1+x][1-y];
-	}
-
-	inline double get(Vector<int> v) const {
-		return get(v[0], v[1]);
-	}
-};
-
-using Velocity = Vector<double>;
-using Density = double;
-
-constexpr DataCell weight{
-	1./36., 1./9., 1./36.,
-	1./9.,  4./9., 1./9.,
-	1./36., 1./9., 1./36
-};
-
-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()));
-			}
-		}
-	}
-
-	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);
-	}
-
-	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)
-		};
-	}
-};
-
-class CellBuffer {
-	private:
-		const std::size_t dim_x_;
-		const std::size_t dim_y_;
-
-		std::unique_ptr<Cell[]> curr_;
-		std::unique_ptr<Cell[]> prev_;
-
-	public:
-		CellBuffer(std::size_t dimX, std::size_t dimY):
-			dim_x_(dimX),
-			dim_y_(dimY),
-			curr_(new Cell[dimX*dimY]),
-			prev_(new Cell[dimX*dimY]) { }
-
-		void swap() {
-			curr_.swap(prev_);
-		}
-
-		inline Cell& curr(std::size_t x, std::size_t y) {
-			return curr_[y*dim_x_ + x];
-		}
-
-		inline Cell& curr(Vector<std::size_t> pos) {
-			return curr(pos[0], pos[1]);
-		}
-
-		inline Cell& prev(std::size_t x, std::size_t y) {
-			return prev_[y*dim_x_ + x];
-		}
-
-		inline Cell& prev(Vector<std::size_t> pos) {
-			return prev(pos[0], pos[1]);
-		}
-};
-
-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    }
-		};
-	}
-}
+#include "lbm.h"
+#include "boundary_conditions.h"
+#include "box_obstacle.h"
 
 constexpr std::size_t dimX = 120;
 constexpr std::size_t dimY = dimX;
@@ -122,110 +15,8 @@ constexpr double reynolds = 1000;
 constexpr double tau   = 3 * uLid * (dimX-1) / reynolds + 0.5;
 constexpr double omega = 1. / tau;
 
-CellBuffer pop(dimX, dimY);
-
-Density  density [dimX][dimY];
-Velocity velocity[dimX][dimY];
-
-void streamFluidCell(std::size_t x, std::size_t y) {
-	if ( x != 0 && x != dimX - 1 && y != 0 && y != 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+i >= 0 && x+i <= dimX - 1 && y+j >= 0 && y+j <= dimY - 1 ) {
-					pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
-				}
-			}
-		}
-	}
-}
-
-void collideFluidCell(std::size_t x, std::size_t y) {
-	// compute equilibrium
-	Cell eq;
-	eq.equilibrize(density[x][y], velocity[x][y]);
-
-	// collide (BGK, relax towards equilibrium)
-	if ( x != 0 && x != dimX - 1 && y != 0 && y != 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+i >= 0 && x+i <= dimX - 1 && y+j >= 0 && y+j <= 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));
-				}
-			}
-		}
-	}
-}
-
-void computeWallCell(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(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);
-}
-
-struct BoxObstacle {
-	const std::size_t lower_x_;
-	const std::size_t lower_y_;
-	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) { }
-
-	bool isInside(std::size_t x, std::size_t y) const {
-		return x > lower_x_
-		    && x < upper_x_
-		    && y > lower_y_
-		    && y < upper_y_;
-	}
-
-	void applyBoundary() const {
-		for ( std::size_t x = lower_x_+1; x < upper_x_; ++x ) {
-			computeWallCell({x, lower_y_}, { 0,-1});
-			computeWallCell({x, upper_y_}, { 0, 1});
-		}
-		for ( std::size_t y = lower_y_+1; y < upper_y_; ++y ) {
-			computeWallCell({lower_x_, y}, {-1, 0});
-			computeWallCell({upper_x_, y}, { 1, 0});
-		}
-		computeWallCell({lower_x_, lower_y_}, {-1,-1});
-		computeWallCell({upper_x_, lower_y_}, { 1,-1});
-		computeWallCell({upper_x_, upper_y_}, { 1, 1});
-		computeWallCell({lower_x_, upper_y_}, {-1, 1});
-	}
-};
+DataCellBuffer pop(dimX, dimY);
+FluidBuffer fluid(dimX, dimY);
 
 std::vector<BoxObstacle> obstacles{
 	{20, 20,  40,  40},
@@ -239,108 +30,69 @@ std::vector<BoxObstacle> obstacles{
 	{80, 80, 100, 100},
 };
 
-void computeLbmStep(std::size_t t) {
+void init() {
+	for ( std::size_t x = 0; x < dimX; ++x ) {
+		for ( std::size_t y = 0; y < dimY; ++y ) {
+			fluid.density(x,y)  = 1.0;
+			fluid.velocity(x,y) = { 0.0, 0.0 };
+
+			static_cast<Cell&>(pop.curr(x,y)).equilibrize(
+				fluid.density(x,y), fluid.velocity(x,y));
+			static_cast<Cell&>(pop.prev(x,y)).equilibrize(
+				fluid.density(x,y), fluid.velocity(x,y));
+		}
+	}
+}
+
+void computeLbmStep() {
 	pop.swap();
 
 	for ( std::size_t x = 0; x < dimX; ++x ) {
 		for ( std::size_t y = 0; y < dimY; ++y ) {
-			if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) -> bool {
+			if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) {
 				return !o.isInside(x, y);
 			}) ) {
-				streamFluidCell(x, y);
+				streamFluidCell(pop, x, y);
 			}
 		}
 	}
 
 	// straight wall cell bounce back
 	for ( std::size_t x = 0; x < dimX; ++x ) {
-		computeZouHeVelocityWallCell({x, dimY-1}, { 0,-1}, uLid);
+		computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uLid);
 	}
 	for ( std::size_t x = 1; x < dimX-1; ++x ) {
-		computeWallCell({x, 0}, { 0, 1});
+		computeWallCell(pop, {x, 0}, { 0, 1});
 	}
 	for ( std::size_t y = 1; y < dimY-1; ++y ) {
-		computeWallCell({0,      y}, { 1, 0});
-		computeWallCell({dimX-1, y}, {-1, 0});
+		computeWallCell(pop, {0,      y}, { 1, 0});
+		computeWallCell(pop, {dimX-1, y}, {-1, 0});
 	}
 
 	// edge wall cell bounce back
-	computeWallCell({0,      0 }, { 1, 1});
-	computeWallCell({dimX-1, 0 }, {-1, 1});
+	computeWallCell(pop, {0,      0 }, { 1, 1});
+	computeWallCell(pop, {dimX-1, 0 }, {-1, 1});
 
 	// obstacles
-	for ( auto& o : obstacles ) {
-		o.applyBoundary();
+	for ( const auto& box : obstacles ) {
+		box.applyBoundary(pop);
 	}
 
 	for ( std::size_t x = 0; x < dimX; ++x ) {
 		for ( std::size_t y = 0; y < dimY; ++y ) {
-			density[x][y]  = pop.curr(x,y).sum();
-			velocity[x][y] = pop.curr(x,y).velocity(density[x][y]);
+			Cell& cell = static_cast<Cell&>(pop.curr(x,y));
+			fluid.density(x,y)  = cell.sum();
+			fluid.velocity(x,y) = cell.velocity(fluid.density(x,y));
 
-			if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) -> bool {
+			if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) {
 				return !o.isInside(x, y);
 			}) ) {
-				collideFluidCell(x, y);
+				collideFluidCell(omega, pop, fluid, x, y);
 			}
 		}
 	}
 }
 
-void init() {
-	for ( std::size_t x = 0; x < dimX; ++x ) {
-		for ( std::size_t y = 0; y < dimY; ++y ) {
-			velocity[x][y] = { 0.0, 0.0 };
-			density[x][y] = 1.0;
-
-			pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]);
-			pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]);
-		}
-	}
-}
-
-void writeCurrentStateAsVTK(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();
-}
-
 int main() {
 	init();
 
@@ -349,12 +101,12 @@ int main() {
 	std::cout << "tau:  " << tau      << std::endl;
 
 	for ( std::size_t t = 0; t <= 6000; ++t ) {
-		computeLbmStep(t);
+		computeLbmStep();
 
 		if ( t % 100 == 0 ) {
 			std::cout << ".";
 			std::cout.flush();
-			writeCurrentStateAsVTK(t);
+			fluid.writeAsVTK(t);
 		}
 	}
 
diff --git a/src/boundary_conditions.h b/src/boundary_conditions.h
new file mode 100644
index 0000000..c825c32
--- /dev/null
+++ b/src/boundary_conditions.h
@@ -0,0 +1,48 @@
+#pragma once
+
+#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) {
+	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/box_obstacle.h b/src/box_obstacle.h
new file mode 100644
index 0000000..46660a7
--- /dev/null
+++ b/src/box_obstacle.h
@@ -0,0 +1,36 @@
+#pragma once
+
+#include "data_cell_buffer.h"
+#include "boundary_conditions.h"
+
+struct BoxObstacle {
+	const std::size_t lower_x_;
+	const std::size_t lower_y_;
+	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) { }
+
+	bool isInside(std::size_t x, std::size_t y) const {
+		return x > lower_x_
+		    && x < upper_x_
+		    && y > lower_y_
+		    && y < upper_y_;
+	}
+
+	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});
+	}
+};
diff --git a/src/data_cell.h b/src/data_cell.h
new file mode 100644
index 0000000..fbb5dfb
--- /dev/null
+++ b/src/data_cell.h
@@ -0,0 +1,23 @@
+#pragma once
+
+#include "vector.h"
+
+struct DataCell {
+	double data[3][3];
+
+	inline double& get(int x, int y) {
+		return data[1+x][1-y];
+	}
+
+	inline double& get(Vector<int> v) {
+		return get(v[0], v[1]);
+	}
+
+	inline double get(int x, int y) const {
+		return data[1+x][1-y];
+	}
+
+	inline double get(Vector<int> v) const {
+		return get(v[0], v[1]);
+	}
+};
diff --git a/src/data_cell_buffer.h b/src/data_cell_buffer.h
new file mode 100644
index 0000000..06c2644
--- /dev/null
+++ b/src/data_cell_buffer.h
@@ -0,0 +1,50 @@
+#pragma once
+
+#include <memory>
+
+#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]);
+		}
+};
+
diff --git a/src/fluid_buffer.h b/src/fluid_buffer.h
new file mode 100644
index 0000000..a6c41ba
--- /dev/null
+++ b/src/fluid_buffer.h
@@ -0,0 +1,91 @@
+#pragma once
+
+#include <fstream>
+#include <memory>
+
+#include "vector.h"
+
+using Velocity = Vector<double>;
+using Density  = double;
+
+class FluidBuffer {
+	private:
+		const std::size_t dim_x_;
+		const std::size_t dim_y_;
+
+		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]) { }
+
+		std::size_t dimX() const {
+			return dim_x_;
+		}
+
+		std::size_t dimY() const {
+			return dim_y_;
+		}
+
+		Density& density(std::size_t x, std::size_t y) {
+			return density_[y*dim_x_ + x];
+		}
+
+		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();
+		}
+};
diff --git a/src/lbm.h b/src/lbm.h
new file mode 100644
index 0000000..685fa08
--- /dev/null
+++ b/src/lbm.h
@@ -0,0 +1,81 @@
+#pragma once
+
+#include "data_cell.h"
+#include "data_cell_buffer.h"
+#include "fluid_buffer.h"
+
+constexpr DataCell weight{
+	1./36., 1./9., 1./36.,
+	1./9.,  4./9., 1./9.,
+	1./36., 1./9., 1./36
+};
+
+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()));
+			}
+		}
+	}
+
+	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);
+	}
+
+	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)
+		};
+	}
+};
+
+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/vector.h b/src/vector.h
index 3a3ed5c..5e105c4 100644
--- a/src/vector.h
+++ b/src/vector.h
@@ -1,3 +1,5 @@
+#pragma once
+
 #include <cmath>
 
 inline double sq(double x) noexcept {