From ac063d4ed8af475f7cb65981307f4d7b48f4a301 Mon Sep 17 00:00:00 2001
From: Erik Strand <erik.strand@cba.mit.edu>
Date: Sun, 15 Sep 2019 16:27:26 -0400
Subject: [PATCH] Switch to periodic boundary conditions
---
simucene/glut_grapher.cpp | 2 +-
simucene/main.cpp | 14 ++++++++++----
simucene/matrix_builder.cpp | 6 ++++++
simucene/matrix_builder.h | 3 +++
4 files changed, 20 insertions(+), 5 deletions(-)
diff --git a/simucene/glut_grapher.cpp b/simucene/glut_grapher.cpp
index 5cf4f91..256a594 100644
--- a/simucene/glut_grapher.cpp
+++ b/simucene/glut_grapher.cpp
@@ -78,7 +78,7 @@ void GlutGrapher::display() const {
//..................................................................................................
void GlutGrapher::idle() {
- bool const wait_to_start = true;
+ bool const wait_to_start = false;
if (wait_to_start && frame_ == 0) {
std::cout << "Press enter to start a 5 second countdown" << std::endl;
do {} while (std::cin.get() != '\n');
diff --git a/simucene/main.cpp b/simucene/main.cpp
index 3b3b3ae..1489d5b 100644
--- a/simucene/main.cpp
+++ b/simucene/main.cpp
@@ -12,12 +12,18 @@ void mouse(int, int, int, int) { exit(0); }
//--------------------------------------------------------------------------------------------------
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) {
- 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**) {
Scalar const alpha = advection_coefficient * h / 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 diffusion_matrix = delta * d_x2_first_order(dim);
SparseMatrix<Scalar> const reaction_matrix = h * diagonal_ones(dim) * reaction_coefficients.asDiagonal();
@@ -69,7 +75,7 @@ int main(int, char**) {
std::cout << transport_matrix << '\n';
*/
- // define the initial state
+ // define the initial state
VectorX<Scalar> state = VectorX<Scalar>::Zero(dim);
for (int32_t i = 0; i < dim; ++i) {
Scalar const x = (i + static_cast<Scalar>(0.5)) * d;
diff --git a/simucene/matrix_builder.cpp b/simucene/matrix_builder.cpp
index 9586396..eba588d 100644
--- a/simucene/matrix_builder.cpp
+++ b/simucene/matrix_builder.cpp
@@ -15,6 +15,12 @@ MatrixBuilder& MatrixBuilder::stripe(int32_t k, Scalar value) {
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 {
SparseMatrix<Scalar> result(n_, n_);
diff --git a/simucene/matrix_builder.h b/simucene/matrix_builder.h
index 7c727c2..07b2550 100644
--- a/simucene/matrix_builder.h
+++ b/simucene/matrix_builder.h
@@ -11,8 +11,11 @@ namespace simucene {
class MatrixBuilder {
public:
MatrixBuilder(int32_t n): n_(n) {}
+
// k = 0 means on the diagonal, 1 means one above the diagonal, -1 one below, etc.
MatrixBuilder& stripe(int32_t k, Scalar value);
+ MatrixBuilder& element(int32_t row, int32_t col, Scalar value);
+
operator SparseMatrix<Scalar>() const;
private:
--
GitLab