From ad6564eeb1581f65d4cb5bacfb66769ec9501037 Mon Sep 17 00:00:00 2001
From: Magnus Ulimoen <flymagnus@gmail.com>
Date: Mon, 15 Jun 2020 22:01:42 +0200
Subject: [PATCH] sparse upwind operators

---
 sbp/src/operators.rs           | 20 ++++++++++++++++++++
 sbp/src/operators/upwind4.rs   | 11 +++++++++++
 sbp/src/operators/upwind4h2.rs | 11 +++++++++++
 sbp/src/operators/upwind9.rs   | 11 +++++++++++
 sbp/src/operators/upwind9h2.rs | 11 +++++++++++
 5 files changed, 64 insertions(+)

diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs
index 2906c94..8479510 100644
--- a/sbp/src/operators.rs
+++ b/sbp/src/operators.rs
@@ -94,12 +94,18 @@ impl<SBP: SbpOperator1d + Copy> SbpOperator2d for SBP {
 pub trait UpwindOperator1d: SbpOperator1d + Send + Sync {
     fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
     fn as_sbp(&self) -> &dyn SbpOperator1d;
+
+    #[cfg(feature = "sparse")]
+    fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float>;
 }
 
 pub trait UpwindOperator2d: SbpOperator2d + Send + Sync {
     fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
     fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
     fn as_sbp(&self) -> &dyn SbpOperator2d;
+
+    fn op_xi(&self) -> &dyn UpwindOperator1d;
+    fn op_eta(&self) -> &dyn UpwindOperator1d;
 }
 
 impl<UOeta: UpwindOperator1d, UOxi: UpwindOperator1d> UpwindOperator2d for (&UOeta, &UOxi) {
@@ -116,6 +122,13 @@ impl<UOeta: UpwindOperator1d, UOxi: UpwindOperator1d> UpwindOperator2d for (&UOe
     fn as_sbp(&self) -> &dyn SbpOperator2d {
         self
     }
+
+    fn op_xi(&self) -> &dyn UpwindOperator1d {
+        self.1
+    }
+    fn op_eta(&self) -> &dyn UpwindOperator1d {
+        self.0
+    }
 }
 
 impl<UO: UpwindOperator1d + Copy> UpwindOperator2d for UO {
@@ -128,6 +141,13 @@ impl<UO: UpwindOperator1d + Copy> UpwindOperator2d for UO {
     fn as_sbp(&self) -> &dyn SbpOperator2d {
         self
     }
+
+    fn op_xi(&self) -> &dyn UpwindOperator1d {
+        self
+    }
+    fn op_eta(&self) -> &dyn UpwindOperator1d {
+        self
+    }
 }
 
 pub trait InterpolationOperator: Send + Sync {
diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs
index d943dfc..69fb97f 100644
--- a/sbp/src/operators/upwind4.rs
+++ b/sbp/src/operators/upwind4.rs
@@ -354,6 +354,17 @@ impl UpwindOperator1d for Upwind4 {
     fn as_sbp(&self) -> &dyn SbpOperator1d {
         self
     }
+
+    #[cfg(feature = "sparse")]
+    fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
+        super::sparse_from_block(
+            Self::DISS_BLOCK,
+            Self::DISS_DIAG,
+            super::Symmetry::Symmetric,
+            super::OperatorType::Normal,
+            n,
+        )
+    }
 }
 
 impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind4) {
diff --git a/sbp/src/operators/upwind4h2.rs b/sbp/src/operators/upwind4h2.rs
index 5c6ee61..6346313 100644
--- a/sbp/src/operators/upwind4h2.rs
+++ b/sbp/src/operators/upwind4h2.rs
@@ -175,4 +175,15 @@ impl UpwindOperator1d for Upwind4h2 {
     fn as_sbp(&self) -> &dyn SbpOperator1d {
         self
     }
+
+    #[cfg(feature = "sparse")]
+    fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
+        super::sparse_from_block(
+            Self::DISS_BLOCK,
+            Self::DISS_DIAG,
+            super::Symmetry::Symmetric,
+            super::OperatorType::H2,
+            n,
+        )
+    }
 }
diff --git a/sbp/src/operators/upwind9.rs b/sbp/src/operators/upwind9.rs
index 62722b4..e4e473a 100644
--- a/sbp/src/operators/upwind9.rs
+++ b/sbp/src/operators/upwind9.rs
@@ -119,6 +119,17 @@ impl UpwindOperator1d for Upwind9 {
     fn as_sbp(&self) -> &dyn SbpOperator1d {
         self
     }
+
+    #[cfg(feature = "sparse")]
+    fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
+        super::sparse_from_block(
+            Self::DISS_BLOCK,
+            Self::DISS_DIAG,
+            super::Symmetry::Symmetric,
+            super::OperatorType::Normal,
+            n,
+        )
+    }
 }
 
 impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind9) {
diff --git a/sbp/src/operators/upwind9h2.rs b/sbp/src/operators/upwind9h2.rs
index 15c3290..143af1f 100644
--- a/sbp/src/operators/upwind9h2.rs
+++ b/sbp/src/operators/upwind9h2.rs
@@ -146,6 +146,17 @@ impl UpwindOperator1d for Upwind9h2 {
     fn as_sbp(&self) -> &dyn SbpOperator1d {
         self
     }
+
+    #[cfg(feature = "sparse")]
+    fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
+        super::sparse_from_block(
+            Self::DISS_BLOCK,
+            Self::DISS_DIAG,
+            super::Symmetry::Symmetric,
+            super::OperatorType::H2,
+            n,
+        )
+    }
 }
 
 impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind9h2) {