aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--CMakeLists.txt4
-rw-r--r--src/boundary_conditions.cc45
-rw-r--r--src/boundary_conditions.h44
-rw-r--r--src/box_obstacle.cc28
-rw-r--r--src/box_obstacle.h26
-rw-r--r--src/data_cell_buffer.h82
-rw-r--r--src/fluid_buffer.cc75
-rw-r--r--src/fluid_buffer.h88
-rw-r--r--src/lbm.cc68
-rw-r--r--src/lbm.h65
10 files changed, 285 insertions, 240 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a5b185d..8818966 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,4 +13,8 @@ include_directories(
add_executable(
boltz
boltz.cc
+ src/lbm.cc
+ src/fluid_buffer.cc
+ src/boundary_conditions.cc
+ src/box_obstacle.cc
)
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);