Interpolate vectors instead of scalars

Same result, nicer code

2 files changed, 41 insertions, 92 deletions

diff --git a/src/refinement/coupler2D.h b/src/refinement/coupler2D.h
index db8ac15..d43a294 100644
--- a/src/refinement/coupler2D.h
+++ b/src/refinement/coupler2D.h
@@ -60,8 +60,8 @@ public:
 template <typename T, template<typename> class DESCRIPTOR>
 class FineCoupler2D : public Coupler2D<T,DESCRIPTOR> {
 private:
-  std::vector<T>           _c2f_rho;
-  std::vector<Vector<T,2>> _c2f_u;
+  std::vector<Vector<T,1>>                _c2f_rho;
+  std::vector<Vector<T,DESCRIPTOR<T>::d>> _c2f_u;
   std::vector<Vector<T,DESCRIPTOR<T>::q>> _c2f_fneq;
 
 public:
diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh
index 1726df4..2218c28 100644
--- a/src/refinement/coupler2D.hh
+++ b/src/refinement/coupler2D.hh
@@ -117,30 +117,36 @@ void FineCoupler2D<T,DESCRIPTOR>::store()
     lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, rho, u);
     lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u);
 
-    _c2f_rho[y]  = rho;
+    _c2f_rho[y]  = Vector<T,1>(rho);
     _c2f_u[y]    = Vector<T,2>(u);
     _c2f_fneq[y] = Vector<T,DESCRIPTOR<T>::q>(fNeq);
   }
 }
 
-template <typename T>
-T order4interpolation(T fm3, T fm1, T f1, T f3)
+template <typename T, unsigned N>
+Vector<T,N> order2interpolation(const Vector<T,N>& f0, const Vector<T,N>& f1)
 {
-  return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
+  return 0.5 * (f0 + f1);
 }
 
-template <typename T>
-T order3interpolation(T fm1, T f1, T f3)
+template <typename T, unsigned N>
+Vector<T,N> order3interpolation(const std::vector<Vector<T,N>>& data, int y, bool ascending)
 {
-  return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
+  if (ascending) {
+    return 3./8. * data[y] + 3./4. * data[y+1] - 1./8. * data[y+2];
+  }
+  else {
+    return 3./8. * data[y] + 3./4. * data[y-1] - 1./8. * data[y-2];
+  }
 }
 
-template <typename T>
-T order2interpolation(T f0, T f1)
+template <typename T, unsigned N>
+Vector<T,N> order4interpolation(const std::vector<Vector<T,N>>& data, int y)
 {
-  return 0.5 * (f0 + f1);
+  return 9./16. * (data[y] + data[y+1]) - 1./16. * (data[y-1] + data[y+2]);
 }
 
+
 template <typename T, template<typename> class DESCRIPTOR>
 void FineCoupler2D<T,DESCRIPTOR>::interpolate()
 {
@@ -153,15 +159,14 @@ void FineCoupler2D<T,DESCRIPTOR>::interpolate()
     T rho{};
     T u[2] {};
     lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, rho, u);
-    _c2f_rho[y]  = order2interpolation(rho, _c2f_rho[y]);
-    _c2f_u[y][0] = order2interpolation(u[0], _c2f_u[y][0]);
-    _c2f_u[y][1] = order2interpolation(u[1], _c2f_u[y][1]);
+
+    _c2f_rho[y] = order2interpolation(Vector<T,1>(rho), _c2f_rho[y]);
+    _c2f_u[y]   = order2interpolation(Vector<T,2>(u), _c2f_u[y]);
 
     T fNeq[DESCRIPTOR<T>::q] {};
     lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u);
-    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]);
-    }
+
+    _c2f_fneq[y] = order2interpolation(Vector<T,DESCRIPTOR<T>::q>(fNeq), _c2f_fneq[y]);
   }
 }
 
@@ -189,110 +194,54 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
   }
 
   for (int y=1; y < this->_coarseSize-2; ++y) {
-    const T rho = order4interpolation(
-                    _c2f_rho[y-1],
-                    _c2f_rho[y],
-                    _c2f_rho[y+1],
-                    _c2f_rho[y+2]
-                  );
-    T u[2] {};
-    u[0] = order4interpolation(
-             _c2f_u[y-1][0],
-             _c2f_u[y][0],
-             _c2f_u[y+1][0],
-             _c2f_u[y+2][0]
-           );
-    u[1] = order4interpolation(
-             _c2f_u[y-1][1],
-             _c2f_u[y][1],
-             _c2f_u[y+1][1],
-             _c2f_u[y+2][1]
-           );
-    const T uSqr = u[0]*u[0] + u[1]*u[1];
+    const auto rho  = order4interpolation(_c2f_rho,  y);
+    const auto u    = order4interpolation(_c2f_u,    y);
+    const auto fneq = order4interpolation(_c2f_fneq, y);
+
+    const T uSqr = u*u;
 
     const auto finePos = this->getFineLatticeR(1+2*y);
     Cell<T,DESCRIPTOR> fineCell;
     fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      const T fneq = order4interpolation(
-                       _c2f_fneq[y-1][iPop],
-                       _c2f_fneq[y][iPop],
-                       _c2f_fneq[y+1][iPop],
-                       _c2f_fneq[y+2][iPop]
-                     );
-
-      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq;
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop];
     }
 
     fineLattice.set(finePos, fineCell);
   }
 
   {
-    const T rho = order3interpolation(
-                    _c2f_rho[0],
-                    _c2f_rho[1],
-                    _c2f_rho[2]
-                  );
-    T u[2] {};
-    u[0] = order3interpolation(
-             _c2f_u[0][0],
-             _c2f_u[1][0],
-             _c2f_u[2][0]
-           );
-    u[1] = order3interpolation(
-             _c2f_u[0][1],
-             _c2f_u[1][1],
-             _c2f_u[2][1]
-           );
-    const T uSqr = u[0]*u[0] + u[1]*u[1];
+    const auto rho  = order3interpolation(_c2f_rho,  0, true);
+    const auto u    = order3interpolation(_c2f_u,    0, true);
+    const auto fneq = order3interpolation(_c2f_fneq, 0, true);
+
+    const T uSqr = u*u;
 
     const auto& finePos = this->getFineLatticeR(1);
     Cell<T,DESCRIPTOR> fineCell;
     fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      const T fneq = order3interpolation(
-                       _c2f_fneq[0][iPop],
-                       _c2f_fneq[1][iPop],
-                       _c2f_fneq[2][iPop]
-                     );
-      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq;
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop];
     }
 
     fineLattice.set(finePos, fineCell);
   }
 
   {
-    const T rho = order3interpolation(
-                    _c2f_rho[this->_coarseSize-1],
-                    _c2f_rho[this->_coarseSize-2],
-                    _c2f_rho[this->_coarseSize-3]
-                  );
-    T u[2] {};
-    u[0] = order3interpolation(
-             _c2f_u[this->_coarseSize-1][0],
-             _c2f_u[this->_coarseSize-2][0],
-             _c2f_u[this->_coarseSize-3][0]
-           );
-    u[1] = order3interpolation(
-             _c2f_u[this->_coarseSize-1][1],
-             _c2f_u[this->_coarseSize-2][1],
-             _c2f_u[this->_coarseSize-3][1]
-           );
-    const T uSqr = u[0]*u[0] + u[1]*u[1];
+    const auto rho  = order3interpolation(_c2f_rho,  this->_coarseSize-1, false);
+    const auto u    = order3interpolation(_c2f_u,    this->_coarseSize-1, false);
+    const auto fneq = order3interpolation(_c2f_fneq, this->_coarseSize-1, false);
+
+    const T uSqr = u*u;
 
     const auto& finePos = this->getFineLatticeR(this->_fineSize-2);
     Cell<T,DESCRIPTOR> fineCell;
     fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      const T fneq = order3interpolation(
-                       _c2f_fneq[this->_coarseSize-1][iPop],
-                       _c2f_fneq[this->_coarseSize-2][iPop],
-                       _c2f_fneq[this->_coarseSize-3][iPop]
-                     );
-      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq;
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop];
     }
 
     fineLattice.set(finePos, fineCell);