From 64dac9f5e9bbf1e6a43e4e6365cad3b41f5209bf Mon Sep 17 00:00:00 2001
From: Matthias Frey <matthias.frey@psi.ch>
Date: Fri, 17 Jul 2020 16:56:27 +0200
Subject: [PATCH] SAAMG: move common functions to RegularDomain
---
src/Solvers/EllipticDomain.cpp | 55 ----------------------------
src/Solvers/EllipticDomain.h | 6 ----
src/Solvers/RectangularDomain.cpp | 35 ------------------
src/Solvers/RectangularDomain.h | 3 --
src/Solvers/RegularDomain.cpp | 59 ++++++++++++++++++++++++++++++-
src/Solvers/RegularDomain.h | 7 ++++
6 files changed, 65 insertions(+), 100 deletions(-)
diff --git a/src/Solvers/EllipticDomain.cpp b/src/Solvers/EllipticDomain.cpp
index 5a4eb141c..452cab2da 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 9a1ec6da3..3562abf5c 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 e8c7b7462..139fadfc2 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 de1888c48..1382d705d 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 bed11b310..38ad64642 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 6bc0849dc..00cd69dd0 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;
};
--
GitLab