diff --git a/src/Solvers/EllipticDomain.cpp b/src/Solvers/EllipticDomain.cpp
index 5a4eb141cc55b032dfe8e961c7cc4447ae643926..452cab2dab896b0c3db79601278270c5ec646f30 100644
--- a/src/Solvers/EllipticDomain.cpp
+++ b/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
diff --git a/src/Solvers/EllipticDomain.h b/src/Solvers/EllipticDomain.h
index 9a1ec6da308d292a893f6394657526f4dff1872c..3562abf5c489b8fa43739f79740c31a6e3b390f0 100644
--- a/src/Solvers/EllipticDomain.h
+++ b/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
diff --git a/src/Solvers/RectangularDomain.cpp b/src/Solvers/RectangularDomain.cpp
index e8c7b7462345a5715407ae3af10a58d39ca73836..139fadfc21ae780e01dc8209981b8f5bec086bac 100644
--- a/src/Solvers/RectangularDomain.cpp
+++ b/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
diff --git a/src/Solvers/RectangularDomain.h b/src/Solvers/RectangularDomain.h
index de1888c48c38fb2383cfe2ec94f691692507bab5..1382d705db168a663dfd77045dc38c10693f5f51 100644
--- a/src/Solvers/RectangularDomain.h
+++ b/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
diff --git a/src/Solvers/RegularDomain.cpp b/src/Solvers/RegularDomain.cpp
index bed11b310c18e39a7afa3844108683f18bf07be6..38ad6464277626851be4c75bb745a668c5e1eea6 100644
--- a/src/Solvers/RegularDomain.cpp
+++ b/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]);
+ }
+}
diff --git a/src/Solvers/RegularDomain.h b/src/Solvers/RegularDomain.h
index 6bc0849dcbae700a00391e445b98193c9b1647bb..00cd69dd0f5739cc45b578fd03e82a6d271798d3 100644
--- a/src/Solvers/RegularDomain.h
+++ b/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;
};