Commit ac063d4e authored by Erik Strand's avatar Erik Strand

Switch to periodic boundary conditions

parent 832649ce
...@@ -78,7 +78,7 @@ void GlutGrapher::display() const { ...@@ -78,7 +78,7 @@ void GlutGrapher::display() const {
//.................................................................................................. //..................................................................................................
void GlutGrapher::idle() { void GlutGrapher::idle() {
bool const wait_to_start = true; bool const wait_to_start = false;
if (wait_to_start && frame_ == 0) { if (wait_to_start && frame_ == 0) {
std::cout << "Press enter to start a 5 second countdown" << std::endl; std::cout << "Press enter to start a 5 second countdown" << std::endl;
do {} while (std::cin.get() != '\n'); do {} while (std::cin.get() != '\n');
......
...@@ -12,12 +12,18 @@ void mouse(int, int, int, int) { exit(0); } ...@@ -12,12 +12,18 @@ void mouse(int, int, int, int) { exit(0); }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
SparseMatrix<Scalar> d_x_first_order(int32_t const dim) { SparseMatrix<Scalar> d_x_first_order(int32_t const dim) {
return MatrixBuilder(dim).stripe(0, -1).stripe(1, 1); // stripe(0, -1) means add a stripe of -1s on the diagonal
// stripe(1, 1) means add a stripe of 1s one above the diagonal
// element(dim - 1, 0, 1) means add a single 1 at row dim - 1 and column 0
return MatrixBuilder(dim).stripe(0, -1).stripe(1, 1).element(dim - 1, 0, 1);
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
SparseMatrix<Scalar> d_x2_first_order(int32_t const dim) { SparseMatrix<Scalar> d_x2_first_order(int32_t const dim) {
return MatrixBuilder(dim).stripe(-1, 1).stripe(0, -2).stripe(1, 1); return MatrixBuilder(dim)
.stripe(-1, 1).element(0, dim - 1, 1)
.stripe(0, -2)
.stripe(1, 1).element(dim - 1, 0, 1);
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
...@@ -49,7 +55,7 @@ int main(int, char**) { ...@@ -49,7 +55,7 @@ int main(int, char**) {
Scalar const alpha = advection_coefficient * h / d; Scalar const alpha = advection_coefficient * h / d;
Scalar const delta = diffusion_coefficient * h / (d * d); Scalar const delta = diffusion_coefficient * h / (d * d);
// transport matrices (assumes dirichlet boundary conditions identically zero) // transport matrices
SparseMatrix<Scalar> const advection_matrix = alpha * d_x_first_order(dim); SparseMatrix<Scalar> const advection_matrix = alpha * d_x_first_order(dim);
SparseMatrix<Scalar> const diffusion_matrix = delta * d_x2_first_order(dim); SparseMatrix<Scalar> const diffusion_matrix = delta * d_x2_first_order(dim);
SparseMatrix<Scalar> const reaction_matrix = h * diagonal_ones(dim) * reaction_coefficients.asDiagonal(); SparseMatrix<Scalar> const reaction_matrix = h * diagonal_ones(dim) * reaction_coefficients.asDiagonal();
...@@ -69,7 +75,7 @@ int main(int, char**) { ...@@ -69,7 +75,7 @@ int main(int, char**) {
std::cout << transport_matrix << '\n'; std::cout << transport_matrix << '\n';
*/ */
// define the initial state // define the initial state
VectorX<Scalar> state = VectorX<Scalar>::Zero(dim); VectorX<Scalar> state = VectorX<Scalar>::Zero(dim);
for (int32_t i = 0; i < dim; ++i) { for (int32_t i = 0; i < dim; ++i) {
Scalar const x = (i + static_cast<Scalar>(0.5)) * d; Scalar const x = (i + static_cast<Scalar>(0.5)) * d;
......
...@@ -15,6 +15,12 @@ MatrixBuilder& MatrixBuilder::stripe(int32_t k, Scalar value) { ...@@ -15,6 +15,12 @@ MatrixBuilder& MatrixBuilder::stripe(int32_t k, Scalar value) {
return *this; return *this;
} }
//..................................................................................................
MatrixBuilder& MatrixBuilder::element(int32_t row, int32_t col, Scalar value) {
triplets_.emplace_back(row, col, value);
return *this;
}
//.................................................................................................. //..................................................................................................
MatrixBuilder::operator SparseMatrix<Scalar>() const { MatrixBuilder::operator SparseMatrix<Scalar>() const {
SparseMatrix<Scalar> result(n_, n_); SparseMatrix<Scalar> result(n_, n_);
......
...@@ -11,8 +11,11 @@ namespace simucene { ...@@ -11,8 +11,11 @@ namespace simucene {
class MatrixBuilder { class MatrixBuilder {
public: public:
MatrixBuilder(int32_t n): n_(n) {} MatrixBuilder(int32_t n): n_(n) {}
// k = 0 means on the diagonal, 1 means one above the diagonal, -1 one below, etc. // k = 0 means on the diagonal, 1 means one above the diagonal, -1 one below, etc.
MatrixBuilder& stripe(int32_t k, Scalar value); MatrixBuilder& stripe(int32_t k, Scalar value);
MatrixBuilder& element(int32_t row, int32_t col, Scalar value);
operator SparseMatrix<Scalar>() const; operator SparseMatrix<Scalar>() const;
private: private:
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment