SAAMG: move common functions to RegularDomain

src/Solvers/EllipticDomain.cpp

@@ -138,42 +138,6 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
     }
 }
 
-void EllipticDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
-                                           double &scaleFactor) const
-{
-    scaleFactor = 1.0;
-
-    value.west  = -1.0 / (hr_m[0] * hr_m[0]);
-    value.east  = -1.0 / (hr_m[0] * hr_m[0]);
-    value.north = -1.0 / (hr_m[1] * hr_m[1]);
-    value.south = -1.0 / (hr_m[1] * hr_m[1]);
-    value.front = -1.0 / (hr_m[2] * hr_m[2]);
-    value.back  = -1.0 / (hr_m[2] * hr_m[2]);
-
-    value.center =  2.0 / (hr_m[0] * hr_m[0])
-                 +  2.0 / (hr_m[1] * hr_m[1])
-                 +  2.0 / (hr_m[2] * hr_m[2]);
-
-    // we are a right boundary point
-    if (!isInside(x + 1, y, z))
-        value.east = 0.0;
-
-    // we are a left boundary point
-    if (!isInside(x - 1, y, z))
-        value.west = 0.0;
-
-    // we are a upper boundary point
-    if (!isInside(x, y + 1, z))
-        value.north = 0.0;
-
-    // we are a lower boundary point
-    if (!isInside(x, y - 1, z))
-        value.south = 0.0;
-
-    // handle boundary condition in z direction
-    robinBoundaryStencil(z, value.front, value.back, value.center);
-}
-
 void EllipticDomain::linearInterpolation(int x, int y, int z, StencilValue_t& value,
                                          double &scaleFactor) const
 {
@@ -330,25 +294,6 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z,
     robinBoundaryStencil(z, value.front, value.back, value.center);
 }
 
-void EllipticDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) const {
-    if (z == 0 || z == nr_m[2] - 1) {
-
-        // case where we are on the Robin BC in Z-direction
-        // where we distinguish two cases
-        // IFF: this values should not matter because they
-        // never make it into the discretization matrix
-        if (z == 0)
-            F = 0.0;
-        else
-            B = 0.0;
-
-        // add contribution of Robin discretization to center point
-        // d the distance between the center of the bunch and the boundary
-        double d = 0.5 * hr_m[2] * (nr_m[2] - 1);
-        C += 2.0 / (d * hr_m[2]);
-    }
-}
-
 // vi: set et ts=4 sw=4 sts=4:
 // Local Variables:
 // mode:c

src/Solvers/EllipticDomain.h

@@ -84,17 +84,11 @@ private:
     }
 
     /// different interpolation methods for boundary points
-    void constantInterpolation(int x, int y, int z, StencilValue_t& value,
-                               double &scaleFactor) const override;
-
     void linearInterpolation(int x, int y, int z, StencilValue_t& value,
                              double &scaleFactor) const override;
 
     void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
                                 double &scaleFactor) const override;
-
-    /// function to handle the open boundary condition in longitudinal direction
-    void robinBoundaryStencil(int z, double &F, double &B, double &C) const;
 };
 
 #endif //#ifdef ELLIPTICAL_DOMAIN_H

src/Solvers/RectangularDomain.cpp

@@ -40,39 +40,4 @@ void RectangularDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
     setHr(hr);
 }
 
-void RectangularDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
-                                              double &scaleFactor) const
-{
-    scaleFactor = 1.0;
-    value.west  = -1.0 / (hr[0] * hr[0]);
-    value.east  = -1.0 / (hr[0] * hr[0]);
-    value.north = -1.0 / (hr[1] * hr[1]);
-    value.south = -1.0 / (hr[1] * hr[1]);
-    value.front = -1.0 / (hr[2] * hr[2]);
-    value.back  = -1.0 / (hr[2] * hr[2]);
-
-    value.center = 2.0 / (hr[0] * hr[0])
-                 + 2.0 / (hr[1] * hr[1])
-                 + 2.0 / (hr[2] * hr[2]);
-
-    if (!isInside(x + 1, y, z))
-        value.east = 0.0; //we are a right boundary point
-
-    if (!isInside(x - 1, y, z))
-        value.west = 0.0; //we are a left boundary point
-
-    if (!isInside(x, y + 1, z))
-        value.north = 0.0; //we are a upper boundary point
-
-    if (!isInside(x, y - 1, z))
-        value.south = 0.0; //we are a lower boundary point
-
-    //dirichlet
-    if (z == 0)
-        value.front = 0.0;
-
-    if (z == nr[2] - 1)
-        value.back = 0.0;
-}
-
 #endif //#ifdef HAVE_SAAMG_SOLVER
\ No newline at end of file

src/Solvers/RectangularDomain.h

@@ -60,9 +60,6 @@ private:
     int coordAccess(int idx) const {
         return idx % getNumXY();
     }
-
-    void constantInterpolation(int x, int y, int z, StencilValue_t& value,
-                               double &scaleFactor) const override;
 };
 
 #endif

src/Solvers/RegularDomain.cpp

@@ -25,6 +25,7 @@ RegularDomain::RegularDomain(const IntVector_t& nr,
     , nxy_m(nr[0] * nr[1])
 { }
 
+
 void RegularDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
                                 const Vector_t& rmax, double dh)
 {
@@ -41,4 +42,60 @@ void RegularDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& r
 
     for (int i = 0; i < 3; ++i)
         hr[i] = (mymax[i] - origin[i]) / nr_m[i];
-}
+}
\ No newline at end of file
+
+
+void RegularDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
+                                          double &scaleFactor) const
+{
+    scaleFactor = 1.0;
+    value.west  = -1.0 / (hr_m[0] * hr_m[0]);
+    value.east  = -1.0 / (hr_m[0] * hr_m[0]);
+    value.north = -1.0 / (hr_m[1] * hr_m[1]);
+    value.south = -1.0 / (hr_m[1] * hr_m[1]);
+    value.front = -1.0 / (hr_m[2] * hr_m[2]);
+    value.back  = -1.0 / (hr_m[2] * hr_m[2]);
+
+    value.center = 2.0 / (hr_m[0] * hr_m[0])
+                 + 2.0 / (hr_m[1] * hr_m[1])
+                 + 2.0 / (hr_m[2] * hr_m[2]);
+
+    // we are a right boundary point
+    if (!isInside(x + 1, y, z))
+        value.east = 0.0;
+
+    // we are a left boundary point
+    if (!isInside(x - 1, y, z))
+        value.west = 0.0;
+
+    // we are a upper boundary point
+    if (!isInside(x, y + 1, z))
+        value.north = 0.0;
+
+    // we are a lower boundary point
+    if (!isInside(x, y - 1, z))
+        value.south = 0.0;
+
+    // handle boundary condition in z direction
+    robinBoundaryStencil(z, value.front, value.back, value.center);
+}
+
+
+void RegularDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) const {
+    if (z == 0 || z == nr_m[2] - 1) {
+
+        // case where we are on the Robin BC in Z-direction
+        // where we distinguish two cases
+        // IFF: this values should not matter because they
+        // never make it into the discretization matrix
+        if (z == 0)
+            F = 0.0;
+        else
+            B = 0.0;
+
+        // add contribution of Robin discretization to center point
+        // d the distance between the center of the bunch and the boundary
+        double d = 0.5 * hr_m[2] * (nr_m[2] - 1);
+        C += 2.0 / (d * hr_m[2]);
+    }
+}

src/Solvers/RegularDomain.h

@@ -38,7 +38,14 @@ public:
     void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
                     const Vector_t& rmax, double dh) override;
 
+protected:
+    /// function to handle the open boundary condition in longitudinal direction
+    void robinBoundaryStencil(int z, double &F, double &B, double &C) const;
+
 private:
+    void constantInterpolation(int x, int y, int z, StencilValue_t& value,
+                               double &scaleFactor) const override;
+
     /// number of nodes in the xy plane (for this case: independent of the z coordinate)
     int nxy_m;
 };