diff --git a/sbp/benches/euler.rs b/sbp/benches/euler.rs
index 2ed59b8..d912573 100644
--- a/sbp/benches/euler.rs
+++ b/sbp/benches/euler.rs
@@ -1,14 +1,14 @@
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use sbp::euler::System;
-use sbp::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4};
+use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4};
 
-fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) {
+fn advance_system<SBP: SbpOperator2d>(universe: &mut System<SBP>, n: usize) {
     for _ in 0..n {
         universe.advance(1.0 / 40.0 * 0.2);
     }
 }
 
-fn advance_system_upwind<UO: UpwindOperator>(universe: &mut System<UO>, n: usize) {
+fn advance_system_upwind<UO: UpwindOperator2d>(universe: &mut System<UO>, n: usize) {
     for _ in 0..n {
         universe.advance_upwind(1.0 / 40.0 * 0.2);
     }
diff --git a/sbp/benches/maxwell.rs b/sbp/benches/maxwell.rs
index 70005ca..8dbbf99 100644
--- a/sbp/benches/maxwell.rs
+++ b/sbp/benches/maxwell.rs
@@ -1,15 +1,15 @@
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use sbp::maxwell::System;
-use sbp::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4};
+use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4};
 use sbp::Float;
 
-fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) {
+fn advance_system<SBP: SbpOperator2d>(universe: &mut System<SBP>, n: usize) {
     for _ in 0..n {
         universe.advance(0.01);
     }
 }
 
-fn advance_system_upwind<UO: UpwindOperator>(universe: &mut System<UO>, n: usize) {
+fn advance_system_upwind<UO: UpwindOperator2d>(universe: &mut System<UO>, n: usize) {
     for _ in 0..n {
         universe.advance_upwind(0.01);
     }
diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs
index a988159..a1be52c 100644
--- a/sbp/src/euler.rs
+++ b/sbp/src/euler.rs
@@ -1,6 +1,6 @@
 use super::grid::{Grid, Metrics};
 use super::integrate;
-use super::operators::{InterpolationOperator, SbpOperator, UpwindOperator};
+use super::operators::{InterpolationOperator, SbpOperator2d, UpwindOperator2d};
 use super::Float;
 use ndarray::azip;
 use ndarray::prelude::*;
@@ -10,7 +10,7 @@ pub const GAMMA: Float = 1.4;
 // A collection of buffers that allows one to efficiently
 // move to the next state
 #[derive(Debug)]
-pub struct System<SBP: SbpOperator> {
+pub struct System<SBP: SbpOperator2d> {
     sys: (Field, Field),
     k: [Field; 4],
     wb: WorkBuffers,
@@ -18,12 +18,12 @@ pub struct System<SBP: SbpOperator> {
     op: SBP,
 }
 
-impl<SBP: SbpOperator> System<SBP> {
+impl<SBP: SbpOperator2d> System<SBP> {
     pub fn new(x: ndarray::Array2<Float>, y: ndarray::Array2<Float>, op: SBP) -> Self {
         let grid = Grid::new(x, y).expect(
             "Could not create grid. Different number of elements compared to width*height?",
         );
-        let metrics = grid.metrics(op, op).unwrap();
+        let metrics = grid.metrics(op).unwrap();
         let nx = grid.nx();
         let ny = grid.ny();
         Self {
@@ -51,7 +51,7 @@ impl<SBP: SbpOperator> System<SBP> {
         let rhs_trad = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| {
             let (grid, metrics) = gm;
             let boundaries = boundary_extractor(y, grid, &bc);
-            RHS_trad((op, op), k, y, metrics, &boundaries, wb)
+            RHS_trad(op, k, y, metrics, &boundaries, wb)
         };
         integrate::integrate::<integrate::Rk4, _, _, _, _>(
             rhs_trad,
@@ -100,7 +100,7 @@ impl<SBP: SbpOperator> System<SBP> {
     }
 }
 
-impl<UO: UpwindOperator> System<UO> {
+impl<UO: UpwindOperator2d> System<UO> {
     pub fn advance_upwind(&mut self, dt: Float) {
         let bc = BoundaryCharacteristics {
             north: BoundaryCharacteristic::This,
@@ -112,7 +112,7 @@ impl<UO: UpwindOperator> System<UO> {
         let rhs_upwind = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| {
             let (grid, metrics) = gm;
             let boundaries = boundary_extractor(y, grid, &bc);
-            RHS_upwind((op, op), k, y, metrics, &boundaries, wb)
+            RHS_upwind(op, k, y, metrics, &boundaries, wb)
         };
         integrate::integrate::<integrate::Rk4, _, _, _, _>(
             rhs_upwind,
@@ -272,11 +272,7 @@ impl Field {
 
 impl Field {
     /// sqrt((self-other)^T*H*(self-other))
-    pub fn h2_err<SBPeta: SbpOperator, SBPxi: SbpOperator>(
-        &self,
-        (opeta, opxi): (SBPeta, SBPxi),
-        other: &Self,
-    ) -> Float {
+    pub fn h2_err<SBP: SbpOperator2d>(&self, op: SBP, other: &Self) -> Float {
         assert_eq!(self.nx(), other.nx());
         assert_eq!(self.ny(), other.ny());
 
@@ -297,9 +293,9 @@ impl Field {
         };
 
         let hxiterator = itermaker(
-            opxi.h(),
+            op.hxi(),
             self.nx(),
-            if opxi.is_h2() {
+            if op.is_h2xi() {
                 1.0 / (self.nx() - 2) as Float
             } else {
                 1.0 / (self.nx() - 1) as Float
@@ -310,9 +306,9 @@ impl Field {
         let hxiterator = hxiterator.cycle().take(self.nx() * self.ny());
 
         let hyiterator = itermaker(
-            opeta.h(),
+            op.heta(),
             self.ny(),
-            1.0 / if opeta.is_h2() {
+            1.0 / if op.is_h2eta() {
                 (self.ny() - 2) as Float
             } else {
                 (self.ny() - 1) as Float
@@ -410,8 +406,8 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
 }
 
 #[allow(non_snake_case)]
-pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
-    (opeta, opxi): (SBPeta, SBPxi),
+pub fn RHS_trad<SBP: SbpOperator2d>(
+    op: SBP,
     k: &mut Field,
     y: &Field,
     metrics: &Metrics,
@@ -424,15 +420,15 @@ pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     let dE = &mut tmp.2;
     let dF = &mut tmp.3;
 
-    opxi.diffxi(ehat.rho(), dE.rho_mut());
-    opxi.diffxi(ehat.rhou(), dE.rhou_mut());
-    opxi.diffxi(ehat.rhov(), dE.rhov_mut());
-    opxi.diffxi(ehat.e(), dE.e_mut());
+    op.diffxi(ehat.rho(), dE.rho_mut());
+    op.diffxi(ehat.rhou(), dE.rhou_mut());
+    op.diffxi(ehat.rhov(), dE.rhov_mut());
+    op.diffxi(ehat.e(), dE.e_mut());
 
-    opeta.diffeta(fhat.rho(), dF.rho_mut());
-    opeta.diffeta(fhat.rhou(), dF.rhou_mut());
-    opeta.diffeta(fhat.rhov(), dF.rhov_mut());
-    opeta.diffeta(fhat.e(), dF.e_mut());
+    op.diffeta(fhat.rho(), dF.rho_mut());
+    op.diffeta(fhat.rhou(), dF.rhou_mut());
+    op.diffeta(fhat.rhov(), dF.rhov_mut());
+    op.diffeta(fhat.e(), dF.e_mut());
 
     azip!((out in &mut k.0,
                     eflux in &dE.0,
@@ -441,12 +437,12 @@ pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
         *out = (-eflux - fflux)/detj
     });
 
-    SAT_characteristics((opeta, opxi), k, y, metrics, boundaries);
+    SAT_characteristics(op, k, y, metrics, boundaries);
 }
 
 #[allow(non_snake_case)]
-pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
-    (opeta, opxi): (UOeta, UOxi),
+pub fn RHS_upwind<UO: UpwindOperator2d>(
+    op: UO,
     k: &mut Field,
     y: &Field,
     metrics: &Metrics,
@@ -459,25 +455,19 @@ pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
     let dE = &mut tmp.2;
     let dF = &mut tmp.3;
 
-    opxi.diffxi(ehat.rho(), dE.rho_mut());
-    opxi.diffxi(ehat.rhou(), dE.rhou_mut());
-    opxi.diffxi(ehat.rhov(), dE.rhov_mut());
-    opxi.diffxi(ehat.e(), dE.e_mut());
+    op.diffxi(ehat.rho(), dE.rho_mut());
+    op.diffxi(ehat.rhou(), dE.rhou_mut());
+    op.diffxi(ehat.rhov(), dE.rhov_mut());
+    op.diffxi(ehat.e(), dE.e_mut());
 
-    opeta.diffeta(fhat.rho(), dF.rho_mut());
-    opeta.diffeta(fhat.rhou(), dF.rhou_mut());
-    opeta.diffeta(fhat.rhov(), dF.rhov_mut());
-    opeta.diffeta(fhat.e(), dF.e_mut());
+    op.diffeta(fhat.rho(), dF.rho_mut());
+    op.diffeta(fhat.rhou(), dF.rhou_mut());
+    op.diffeta(fhat.rhov(), dF.rhov_mut());
+    op.diffeta(fhat.e(), dF.e_mut());
 
     let ad_xi = &mut tmp.4;
     let ad_eta = &mut tmp.5;
-    upwind_dissipation(
-        (opeta, opxi),
-        (ad_xi, ad_eta),
-        y,
-        metrics,
-        (&mut tmp.0, &mut tmp.1),
-    );
+    upwind_dissipation(op, (ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1));
 
     azip!((out in &mut k.0,
                     eflux in &dE.0,
@@ -488,12 +478,12 @@ pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
         *out = (-eflux - fflux + ad_xi + ad_eta)/detj
     });
 
-    SAT_characteristics((opeta, opxi), k, y, metrics, boundaries);
+    SAT_characteristics(op, k, y, metrics, boundaries);
 }
 
 #[allow(clippy::many_single_char_names)]
-fn upwind_dissipation<UOeta: UpwindOperator, UOxi: UpwindOperator>(
-    (opeta, opxi): (UOeta, UOxi),
+fn upwind_dissipation<UO: UpwindOperator2d>(
+    op: UO,
     k: (&mut Field, &mut Field),
     y: &Field,
     metrics: &Metrics,
@@ -549,15 +539,15 @@ fn upwind_dissipation<UOeta: UpwindOperator, UOxi: UpwindOperator>(
         tmp1[3] = alpha_v * e * detj;
     }
 
-    opxi.dissxi(tmp.0.rho(), k.0.rho_mut());
-    opxi.dissxi(tmp.0.rhou(), k.0.rhou_mut());
-    opxi.dissxi(tmp.0.rhov(), k.0.rhov_mut());
-    opxi.dissxi(tmp.0.e(), k.0.e_mut());
+    op.dissxi(tmp.0.rho(), k.0.rho_mut());
+    op.dissxi(tmp.0.rhou(), k.0.rhou_mut());
+    op.dissxi(tmp.0.rhov(), k.0.rhov_mut());
+    op.dissxi(tmp.0.e(), k.0.e_mut());
 
-    opeta.disseta(tmp.1.rho(), k.1.rho_mut());
-    opeta.disseta(tmp.1.rhou(), k.1.rhou_mut());
-    opeta.disseta(tmp.1.rhov(), k.1.rhov_mut());
-    opeta.disseta(tmp.1.e(), k.1.e_mut());
+    op.disseta(tmp.1.rho(), k.1.rho_mut());
+    op.disseta(tmp.1.rhou(), k.1.rhou_mut());
+    op.disseta(tmp.1.rhov(), k.1.rhov_mut());
+    op.disseta(tmp.1.e(), k.1.e_mut());
 }
 
 fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics) {
@@ -853,8 +843,8 @@ fn vortexify(
 
 #[allow(non_snake_case)]
 /// Boundary conditions (SAT)
-fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
-    (opeta, opxi): (SBPeta, SBPxi),
+fn SAT_characteristics<SBP: SbpOperator2d>(
+    op: SBP,
     k: &mut Field,
     y: &Field,
     metrics: &Metrics,
@@ -862,10 +852,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
 ) {
     // North boundary
     {
-        let hi = if opeta.is_h2() {
-            (k.ny() - 2) as Float / opeta.h()[0]
+        let hi = if op.is_h2eta() {
+            (k.ny() - 2) as Float / op.heta()[0]
         } else {
-            (k.ny() - 1) as Float / opeta.h()[0]
+            (k.ny() - 1) as Float / op.heta()[0]
         };
         let sign = -1.0;
         let tau = 1.0;
@@ -884,10 +874,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     }
     // South boundary
     {
-        let hi = if opeta.is_h2() {
-            (k.ny() - 2) as Float / opeta.h()[0]
+        let hi = if op.is_h2eta() {
+            (k.ny() - 2) as Float / op.heta()[0]
         } else {
-            (k.ny() - 1) as Float / opeta.h()[0]
+            (k.ny() - 1) as Float / op.heta()[0]
         };
         let sign = 1.0;
         let tau = -1.0;
@@ -906,10 +896,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     }
     // West Boundary
     {
-        let hi = if opxi.is_h2() {
-            (k.nx() - 2) as Float / opxi.h()[0]
+        let hi = if op.is_h2xi() {
+            (k.nx() - 2) as Float / op.hxi()[0]
         } else {
-            (k.nx() - 1) as Float / opxi.h()[0]
+            (k.nx() - 1) as Float / op.hxi()[0]
         };
         let sign = 1.0;
         let tau = -1.0;
@@ -928,10 +918,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     }
     // East Boundary
     {
-        let hi = if opxi.is_h2() {
-            (k.nx() - 2) as Float / opxi.h()[0]
+        let hi = if op.is_h2xi() {
+            (k.nx() - 2) as Float / op.hxi()[0]
         } else {
-            (k.nx() - 1) as Float / opxi.h()[0]
+            (k.nx() - 1) as Float / op.hxi()[0]
         };
         let sign = -1.0;
         let tau = 1.0;
diff --git a/sbp/src/grid.rs b/sbp/src/grid.rs
index c99e6a3..76fb5f8 100644
--- a/sbp/src/grid.rs
+++ b/sbp/src/grid.rs
@@ -36,12 +36,11 @@ impl Grid {
         self.y.view()
     }
 
-    pub fn metrics<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>(
+    pub fn metrics<SBP: super::operators::SbpOperator2d>(
         &self,
-        opeta: SBPeta,
-        opxi: SBPxi,
+        op: SBP,
     ) -> Result<Metrics, ndarray::ShapeError> {
-        Metrics::new(self, opeta, opxi)
+        Metrics::new(self, op)
     }
 
     pub fn north(&self) -> (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>) {
@@ -71,10 +70,9 @@ impl Grid {
 }
 
 impl Metrics {
-    fn new<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>(
+    fn new<SBP: super::operators::SbpOperator2d>(
         grid: &Grid,
-        opeta: SBPeta,
-        opxi: SBPxi,
+        op: SBP,
     ) -> Result<Self, ndarray::ShapeError> {
         let ny = grid.ny();
         let nx = grid.nx();
@@ -82,13 +80,13 @@ impl Metrics {
         let y = &grid.y;
 
         let mut dx_dxi = Array2::zeros((ny, nx));
-        opxi.diffxi(x.view(), dx_dxi.view_mut());
+        op.diffxi(x.view(), dx_dxi.view_mut());
         let mut dx_deta = Array2::zeros((ny, nx));
-        opeta.diffeta(x.view(), dx_deta.view_mut());
+        op.diffeta(x.view(), dx_deta.view_mut());
         let mut dy_dxi = Array2::zeros((ny, nx));
-        opxi.diffxi(y.view(), dy_dxi.view_mut());
+        op.diffxi(y.view(), dy_dxi.view_mut());
         let mut dy_deta = Array2::zeros((ny, nx));
-        opeta.diffeta(y.view(), dy_deta.view_mut());
+        op.diffeta(y.view(), dy_deta.view_mut());
 
         let mut detj = Array2::zeros((ny, nx));
         ndarray::azip!((detj in &mut detj,
diff --git a/sbp/src/maxwell.rs b/sbp/src/maxwell.rs
index de2519c..1bcab95 100644
--- a/sbp/src/maxwell.rs
+++ b/sbp/src/maxwell.rs
@@ -1,6 +1,6 @@
 use super::grid::{Grid, Metrics};
 use super::integrate;
-use super::operators::{SbpOperator, UpwindOperator};
+use super::operators::{SbpOperator2d, UpwindOperator2d};
 use crate::Float;
 use ndarray::azip;
 use ndarray::prelude::*;
@@ -74,7 +74,7 @@ impl Field {
 }
 
 #[derive(Debug, Clone)]
-pub struct System<SBP: SbpOperator> {
+pub struct System<SBP: SbpOperator2d> {
     sys: (Field, Field),
     wb: WorkBuffers,
     grid: Grid,
@@ -82,14 +82,14 @@ pub struct System<SBP: SbpOperator> {
     op: SBP,
 }
 
-impl<SBP: SbpOperator> System<SBP> {
+impl<SBP: SbpOperator2d> System<SBP> {
     pub fn new(x: Array2<Float>, y: Array2<Float>, op: SBP) -> Self {
         assert_eq!(x.shape(), y.shape());
         let ny = x.shape()[0];
         let nx = x.shape()[1];
 
         let grid = Grid::new(x, y).unwrap();
-        let metrics = grid.metrics(op, op).unwrap();
+        let metrics = grid.metrics(op).unwrap();
 
         Self {
             op,
@@ -146,7 +146,7 @@ impl<SBP: SbpOperator> System<SBP> {
     }
 }
 
-impl<UO: UpwindOperator> System<UO> {
+impl<UO: UpwindOperator2d> System<UO> {
     /// Using artificial dissipation with the upwind operator
     pub fn advance_upwind(&mut self, dt: Float) {
         let op = self.op;
@@ -205,7 +205,7 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
 /// where J is the grid determinant
 ///
 /// This is used both in fluxes and SAT terms
-fn RHS<SBP: SbpOperator>(
+fn RHS<SBP: SbpOperator2d>(
     op: SBP,
     k: &mut Field,
     y: &Field,
@@ -230,7 +230,7 @@ fn RHS<SBP: SbpOperator>(
 }
 
 #[allow(non_snake_case)]
-fn RHS_upwind<UO: UpwindOperator>(
+fn RHS_upwind<UO: UpwindOperator2d>(
     op: UO,
     k: &mut Field,
     y: &Field,
@@ -255,7 +255,7 @@ fn RHS_upwind<UO: UpwindOperator>(
     });
 }
 
-fn fluxes<SBP: SbpOperator>(
+fn fluxes<SBP: super::operators::SbpOperator2d>(
     op: SBP,
     k: &mut Field,
     y: &Field,
@@ -330,7 +330,7 @@ fn fluxes<SBP: SbpOperator>(
     }
 }
 
-fn dissipation<UO: UpwindOperator>(
+fn dissipation<UO: UpwindOperator2d>(
     op: UO,
     k: &mut Field,
     y: &Field,
@@ -432,7 +432,7 @@ pub struct BoundaryTerms {
 
 #[allow(non_snake_case)]
 /// Boundary conditions (SAT)
-fn SAT_characteristics<SBP: SbpOperator>(
+fn SAT_characteristics<SBP: SbpOperator2d>(
     op: SBP,
     k: &mut Field,
     y: &Field,
@@ -464,10 +464,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
             Boundary::This => y.slice(s![.., .., 0]),
         };
         // East boundary
-        let hinv = if op.is_h2() {
-            (nx - 2) as Float / op.h()[0]
+        let hinv = if op.is_h2xi() {
+            (nx - 2) as Float / op.hxi()[0]
         } else {
-            (nx - 1) as Float / op.h()[0]
+            (nx - 1) as Float / op.hxi()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., .., nx - 1])
@@ -502,10 +502,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let g = match boundaries.east {
             Boundary::This => y.slice(s![.., .., nx - 1]),
         };
-        let hinv = if op.is_h2() {
-            (nx - 2) as Float / op.h()[0]
+        let hinv = if op.is_h2xi() {
+            (nx - 2) as Float / op.hxi()[0]
         } else {
-            (nx - 1) as Float / op.h()[0]
+            (nx - 1) as Float / op.hxi()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., .., 0])
@@ -545,10 +545,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let g = match boundaries.north {
             Boundary::This => y.slice(s![.., 0, ..]),
         };
-        let hinv = if op.is_h2() {
-            (ny - 2) as Float / op.h()[0]
+        let hinv = if op.is_h2eta() {
+            (ny - 2) as Float / op.heta()[0]
         } else {
-            (ny - 1) as Float / op.h()[0]
+            (ny - 1) as Float / op.heta()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., ny - 1, ..])
@@ -582,10 +582,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let g = match boundaries.south {
             Boundary::This => y.slice(s![.., ny - 1, ..]),
         };
-        let hinv = if op.is_h2() {
-            (ny - 2) as Float / op.h()[0]
+        let hinv = if op.is_h2eta() {
+            (ny - 2) as Float / op.heta()[0]
         } else {
-            (ny - 1) as Float / op.h()[0]
+            (ny - 1) as Float / op.heta()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., 0, ..])
diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs
index 5ad2742..8b02820 100644
--- a/sbp/src/operators.rs
+++ b/sbp/src/operators.rs
@@ -2,32 +2,104 @@
 #![allow(clippy::unreadable_literal)]
 
 use crate::Float;
-
 use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
 
-pub trait SbpOperator: Send + Sync + Copy {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
-    fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
-        assert_eq!(prev.shape(), fut.shape());
-        for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
-            self.diff1d(r0, r1);
-        }
-    }
-    fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
-        self.diffxi(prev.reversed_axes(), fut.reversed_axes())
-    }
+pub trait SbpOperator1d: Copy + Clone + core::fmt::Debug {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
+
     fn h(&self) -> &'static [Float];
     fn is_h2(&self) -> bool {
         false
     }
 }
 
-pub trait UpwindOperator: SbpOperator {
-    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
+pub trait SbpOperator2d: Copy + Clone {
+    fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
+    fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
+
+    fn hxi(&self) -> &'static [Float];
+    fn heta(&self) -> &'static [Float];
+
+    fn is_h2xi(&self) -> bool;
+    fn is_h2eta(&self) -> bool;
+}
+
+impl<SBPeta: SbpOperator1d, SBPxi: SbpOperator1d> SbpOperator2d for (SBPeta, SBPxi) {
+    fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+        assert_eq!(prev.shape(), fut.shape());
+        for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
+            self.1.diff(r0, r1)
+        }
+    }
+    fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
+        let ba = (self.1, self.0);
+        ba.diffxi(prev.reversed_axes(), fut.reversed_axes())
+    }
+    fn hxi(&self) -> &'static [Float] {
+        self.1.h()
+    }
+    fn heta(&self) -> &'static [Float] {
+        self.0.h()
+    }
+    fn is_h2xi(&self) -> bool {
+        self.1.is_h2()
+    }
+    fn is_h2eta(&self) -> bool {
+        self.0.is_h2()
+    }
+}
+
+impl<SBP: SbpOperator1d> SbpOperator2d for SBP {
+    fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+        assert_eq!(prev.shape(), fut.shape());
+        for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
+            self.diff(r0, r1)
+        }
+    }
+    fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
+        self.diffxi(prev.reversed_axes(), fut.reversed_axes())
+    }
+    fn hxi(&self) -> &'static [Float] {
+        self.h()
+    }
+    fn heta(&self) -> &'static [Float] {
+        self.h()
+    }
+    fn is_h2xi(&self) -> bool {
+        self.is_h2()
+    }
+    fn is_h2eta(&self) -> bool {
+        self.is_h2()
+    }
+}
+
+pub trait UpwindOperator1d: SbpOperator1d + Copy + Clone {
+    fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
+}
+
+pub trait UpwindOperator2d: SbpOperator2d + Copy + Clone {
+    fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
+    fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
+}
+
+impl<UOeta: UpwindOperator1d, UOxi: UpwindOperator1d> UpwindOperator2d for (UOeta, UOxi) {
     fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
-            self.diss1d(r0, r1);
+            self.1.diss(r0, r1);
+        }
+    }
+    fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
+        let ba = (self.1, self.0);
+        ba.dissxi(prev.reversed_axes(), fut.reversed_axes())
+    }
+}
+
+impl<UO: UpwindOperator1d> UpwindOperator2d for UO {
+    fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+        assert_eq!(prev.shape(), fut.shape());
+        for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
+            self.diss(r0, r1);
         }
     }
     fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
@@ -140,7 +212,7 @@ pub(crate) mod testing {
         dfdy: FY,
         eps: Float,
     ) where
-        SBP: SbpOperator,
+        SBP: SbpOperator2d,
         F: Fn(Float, Float) -> Float,
         FX: Fn(Float, Float) -> Float,
         FY: Fn(Float, Float) -> Float,
diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs
index 4dd9917..100a0f6 100644
--- a/sbp/src/operators/traditional4.rs
+++ b/sbp/src/operators/traditional4.rs
@@ -1,4 +1,4 @@
-use super::SbpOperator;
+use super::*;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
@@ -23,8 +23,8 @@ impl SBP4 {
     ];
 }
 
-impl SbpOperator for SBP4 {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl SbpOperator1d for SBP4 {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs
index 8fe19b6..10ffe03 100644
--- a/sbp/src/operators/traditional8.rs
+++ b/sbp/src/operators/traditional8.rs
@@ -1,4 +1,4 @@
-use super::SbpOperator;
+use super::*;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
@@ -27,8 +27,8 @@ impl SBP8 {
     ];
 }
 
-impl SbpOperator for SBP8 {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl SbpOperator1d for SBP8 {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs
index 834fbeb..ba0f362 100644
--- a/sbp/src/operators/upwind4.rs
+++ b/sbp/src/operators/upwind4.rs
@@ -1,4 +1,4 @@
-use super::{SbpOperator, UpwindOperator};
+use super::*;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis};
 
@@ -275,8 +275,8 @@ impl Upwind4 {
     ];
 }
 
-impl SbpOperator for Upwind4 {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl SbpOperator1d for Upwind4 {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -286,6 +286,12 @@ impl SbpOperator for Upwind4 {
             fut,
         )
     }
+    fn h(&self) -> &'static [Float] {
+        Self::HBLOCK
+    }
+}
+
+/*
     fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
@@ -307,10 +313,8 @@ impl SbpOperator for Upwind4 {
         }
     }
 
-    fn h(&self) -> &'static [Float] {
-        Self::HBLOCK
-    }
 }
+*/
 
 #[test]
 fn upwind4_test() {
@@ -326,7 +330,7 @@ fn upwind4_test() {
         target[i] = 1.0;
     }
     res.fill(0.0);
-    Upwind4.diff1d(source.view(), res.view_mut());
+    Upwind4.diff(source.view(), res.view_mut());
     approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
     {
         let source = source.to_owned().insert_axis(ndarray::Axis(0));
@@ -352,7 +356,7 @@ fn upwind4_test() {
         target[i] = 2.0 * x;
     }
     res.fill(0.0);
-    Upwind4.diff1d(source.view(), res.view_mut());
+    Upwind4.diff(source.view(), res.view_mut());
     approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
     {
         let source = source.to_owned().insert_axis(ndarray::Axis(0));
@@ -378,7 +382,7 @@ fn upwind4_test() {
         target[i] = 3.0 * x * x;
     }
     res.fill(0.0);
-    Upwind4.diff1d(source.view(), res.view_mut());
+    Upwind4.diff(source.view(), res.view_mut());
     approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
 
     {
@@ -400,8 +404,8 @@ fn upwind4_test() {
     }
 }
 
-impl UpwindOperator for Upwind4 {
-    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl UpwindOperator1d for Upwind4 {
+    fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
@@ -411,6 +415,9 @@ impl UpwindOperator for Upwind4 {
             fut,
         )
     }
+}
+
+/*
     fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
@@ -432,6 +439,7 @@ impl UpwindOperator for Upwind4 {
         }
     }
 }
+*/
 
 #[test]
 fn upwind4_test2() {
diff --git a/sbp/src/operators/upwind4h2.rs b/sbp/src/operators/upwind4h2.rs
index d9304ba..837fce6 100644
--- a/sbp/src/operators/upwind4h2.rs
+++ b/sbp/src/operators/upwind4h2.rs
@@ -1,4 +1,4 @@
-use super::{SbpOperator, UpwindOperator};
+use super::*;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
@@ -36,8 +36,8 @@ impl Upwind4h2 {
     ];
 }
 
-impl SbpOperator for Upwind4h2 {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl SbpOperator1d for Upwind4h2 {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -64,25 +64,25 @@ fn upwind4h2_test() {
 
     let mut res = ndarray::Array1::zeros(nx);
 
-    Upwind4h2.diff1d(x.view(), res.view_mut());
+    Upwind4h2.diff(x.view(), res.view_mut());
     let ans = &x * 0.0 + 1.0;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x / 2.0;
-    Upwind4h2.diff1d(y.view(), res.view_mut());
+    Upwind4h2.diff(y.view(), res.view_mut());
     let ans = &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x * &x / 3.0;
-    Upwind4h2.diff1d(y.view(), res.view_mut());
+    Upwind4h2.diff(y.view(), res.view_mut());
     let ans = &x * &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2);
 }
 
-impl UpwindOperator for Upwind4h2 {
-    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl UpwindOperator1d for Upwind4h2 {
+    fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
diff --git a/sbp/src/operators/upwind9.rs b/sbp/src/operators/upwind9.rs
index 9c11143..c6d6dce 100644
--- a/sbp/src/operators/upwind9.rs
+++ b/sbp/src/operators/upwind9.rs
@@ -1,4 +1,4 @@
-use super::{SbpOperator, UpwindOperator};
+use super::*;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
@@ -44,8 +44,8 @@ impl Upwind9 {
     ];
 }
 
-impl SbpOperator for Upwind9 {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl SbpOperator1d for Upwind9 {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -61,8 +61,8 @@ impl SbpOperator for Upwind9 {
     }
 }
 
-impl UpwindOperator for Upwind9 {
-    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl UpwindOperator1d for Upwind9 {
+    fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
diff --git a/sbp/src/operators/upwind9h2.rs b/sbp/src/operators/upwind9h2.rs
index b65296d..811d53a 100644
--- a/sbp/src/operators/upwind9h2.rs
+++ b/sbp/src/operators/upwind9h2.rs
@@ -1,4 +1,4 @@
-use super::{SbpOperator, UpwindOperator};
+use super::*;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
@@ -44,8 +44,8 @@ impl Upwind9h2 {
     ];
 }
 
-impl SbpOperator for Upwind9h2 {
-    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl SbpOperator1d for Upwind9h2 {
+    fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -72,25 +72,25 @@ fn upwind9h2_test() {
 
     let mut res = ndarray::Array1::zeros(nx);
 
-    Upwind9h2.diff1d(x.view(), res.view_mut());
+    Upwind9h2.diff(x.view(), res.view_mut());
     let ans = &x * 0.0 + 1.0;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x / 2.0;
-    Upwind9h2.diff1d(y.view(), res.view_mut());
+    Upwind9h2.diff(y.view(), res.view_mut());
     let ans = &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x * &x / 3.0;
-    Upwind9h2.diff1d(y.view(), res.view_mut());
+    Upwind9h2.diff(y.view(), res.view_mut());
     let ans = &x * &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2);
 }
 
-impl UpwindOperator for Upwind9h2 {
-    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+impl UpwindOperator1d for Upwind9h2 {
+    fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),