From 0fddbe18ffc147aa8d9ab67dc3501ef3f38a5b09 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (l)" <j.wuttke@fz-juelich.de>
Date: Wed, 30 Dec 2020 21:12:35 +0100
Subject: [PATCH] + Polygon test
---
Sample/HardParticle/PolyhedralComponents.cpp | 15 +++++--
Sample/HardParticle/PolyhedralComponents.h | 2 +
.../Numeric/FormFactorTriangleTest.cpp | 41 +++++++++++++++++++
auto/Wrap/doxygenSample.i | 6 +++
4 files changed, 61 insertions(+), 3 deletions(-)
create mode 100644 Tests/UnitTests/Numeric/FormFactorTriangleTest.cpp
diff --git a/Sample/HardParticle/PolyhedralComponents.cpp b/Sample/HardParticle/PolyhedralComponents.cpp
index 6c19b2b83df..78319c48525 100644
--- a/Sample/HardParticle/PolyhedralComponents.cpp
+++ b/Sample/HardParticle/PolyhedralComponents.cpp
@@ -371,11 +371,10 @@ complex_t PolyhedralFace::ff_2D(cvector_t qpa) const
return m_area;
} else if (qpa_red < qpa_limit_series && !sym_S2) {
// summation of power series
- return m_area + expansion(1., 1., qpa, std::abs(m_area));
+ return ff_2D_expanded(qpa);
} else {
// direct evaluation of analytic formula
- complex_t ff = edge_sum_ff(qpa, qpa, false);
- complex_t ret = (sym_S2 ? 4. : 2. / I) * ff / qpa.mag2();
+ complex_t ret = ff_2D_direct(qpa);
#ifdef POLYHEDRAL_DIAGNOSTIC
if (diagnosis.debmsg >= 2)
std::cout << std::setprecision(16) << " ret=" << ret << " ff=" << ff << "\n";
@@ -384,6 +383,16 @@ complex_t PolyhedralFace::ff_2D(cvector_t qpa) const
}
}
+complex_t PolyhedralFace::ff_2D_direct(cvector_t qpa) const
+{
+ return (sym_S2 ? 4. : 2. / I) * edge_sum_ff(qpa, qpa, false) / qpa.mag2();
+}
+
+complex_t PolyhedralFace::ff_2D_expanded(cvector_t qpa) const
+{
+ return m_area + expansion(1., 1., qpa, std::abs(m_area));
+}
+
//! Throws if deviation from inversion symmetry is detected. Does not check vertices.
void PolyhedralFace::assert_Ci(const PolyhedralFace& other) const
diff --git a/Sample/HardParticle/PolyhedralComponents.h b/Sample/HardParticle/PolyhedralComponents.h
index 2010db55984..5ea62d689b9 100644
--- a/Sample/HardParticle/PolyhedralComponents.h
+++ b/Sample/HardParticle/PolyhedralComponents.h
@@ -61,6 +61,8 @@ public:
complex_t ff_n(int m, cvector_t q) const;
complex_t ff(cvector_t q, bool sym_Ci) const;
complex_t ff_2D(cvector_t qpa) const;
+ complex_t ff_2D_direct(cvector_t qpa) const; // for TestTriangle
+ complex_t ff_2D_expanded(cvector_t qpa) const; // for TestTriangle
void assert_Ci(const PolyhedralFace& other) const;
private:
diff --git a/Tests/UnitTests/Numeric/FormFactorTriangleTest.cpp b/Tests/UnitTests/Numeric/FormFactorTriangleTest.cpp
new file mode 100644
index 00000000000..351ba370d77
--- /dev/null
+++ b/Tests/UnitTests/Numeric/FormFactorTriangleTest.cpp
@@ -0,0 +1,41 @@
+#include "Base/Math/Constants.h"
+#include "Sample/HardParticle/PolyhedralComponents.h"
+#include "Tests/GTestWrapper/google_test.h"
+//#include "Tests/UnitTests/Numeric/FormFactorTest.h"
+#include <cstdio>
+
+//! Ad-hoc code for manuscript ffp - JWu dec20
+
+class FFTriangleTest : public ::testing::Test {};
+
+TEST_F(FFTriangleTest, Triangle)
+{
+ const double a = 1.;
+ const double as = a / 2;
+ const double ac = a / sqrt(3) / 2;
+ const double ah = a / sqrt(3);
+ const std::vector<kvector_t> V{{-ac, as, 0.}, {-ac, -as, 0.}, {ah, 0., 0.}};
+
+ const PolyhedralFace T(V, false);
+ EXPECT_NEAR(T.area(), sqrt(3)/4, 1e-15);
+
+ int failures = 0;
+ const int M=43;
+ for (int j=0; j < M; ++j) {
+ const double phi = M_PI_2*j/(M-1);
+ const cvector_t uQ{ sin(phi), cos(phi), 0. };
+ const int N=4000+j;
+ for (int i = 0; i < N; ++i) {
+ const double q = 1e-17*pow(1.7e17,i/(N-1.));
+ const cvector_t Q = uQ * q;
+ const double f1 = std::abs(T.ff_2D_direct(Q));
+ const double f2 = std::abs(T.ff_2D_expanded(Q));
+ const double relerr = std::abs(f1-f2)/f2;
+ if (relerr>7e-16) {
+ printf("%9.6f %16.11e %21.16e %21.16e %10.4e\n", phi, q, f1, f2, relerr);
+ ++failures;
+ }
+ }
+ }
+ EXPECT_EQ(failures, 0);
+}
diff --git a/auto/Wrap/doxygenSample.i b/auto/Wrap/doxygenSample.i
index 0b64464531d..c2e7db6d11c 100644
--- a/auto/Wrap/doxygenSample.i
+++ b/auto/Wrap/doxygenSample.i
@@ -6262,6 +6262,12 @@ Returns the contribution ff(q) of this face to the polyhedral form factor.
Returns the two-dimensional form factor of this face, for use in a prism.
";
+%feature("docstring") PolyhedralFace::ff_2D_direct "complex_t PolyhedralFace::ff_2D_direct(cvector_t qpa) const
+";
+
+%feature("docstring") PolyhedralFace::ff_2D_expanded "complex_t PolyhedralFace::ff_2D_expanded(cvector_t qpa) const
+";
+
%feature("docstring") PolyhedralFace::assert_Ci "void PolyhedralFace::assert_Ci(const PolyhedralFace &other) const
Throws if deviation from inversion symmetry is detected. Does not check vertices.
--
GitLab