From 387ddfac098d47200fa9a3a0d6e071e78ee5997a Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (h)" <j.wuttke@fz-juelich.de>
Date: Wed, 1 Dec 2021 17:39:55 +0100
Subject: [PATCH] temporarily reintegrate libformfactor in BornAgain
---
Base/CMakeLists.txt | 2 +-
CMakeLists.txt | 2 +
cmake/BornAgain/Dependences.cmake | 4 +-
ff/CMakeLists.txt | 33 +++
ff/Factorial.h | 59 +++++
ff/PolyhedralComponents.cpp | 397 ++++++++++++++++++++++++++++++
ff/PolyhedralComponents.h | 103 ++++++++
ff/PolyhedralTopology.h | 44 ++++
ff/Polyhedron.cpp | 192 +++++++++++++++
ff/Polyhedron.h | 64 +++++
10 files changed, 897 insertions(+), 3 deletions(-)
create mode 100644 ff/CMakeLists.txt
create mode 100644 ff/Factorial.h
create mode 100644 ff/PolyhedralComponents.cpp
create mode 100644 ff/PolyhedralComponents.h
create mode 100644 ff/PolyhedralTopology.h
create mode 100644 ff/Polyhedron.cpp
create mode 100644 ff/Polyhedron.h
diff --git a/Base/CMakeLists.txt b/Base/CMakeLists.txt
index 9cded5346ed..e9c83e0b799 100644
--- a/Base/CMakeLists.txt
+++ b/Base/CMakeLists.txt
@@ -22,7 +22,7 @@ target_link_libraries(${lib}
PUBLIC
${GSL_LIBRARIES}
${FFTW3_LIBRARIES}
- formfactor
+ formfactor0
)
target_include_directories(${lib}
PUBLIC
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 10289ce2af6..89606bb2f13 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -150,6 +150,8 @@ add_subdirectory(Fit)
add_subdirectory(Fit/Test/Unit)
add_subdirectory(Fit/Test/Functional)
+# TEMPORARY formfactor
+add_subdirectory(ff)
## recurse into the given subdirectories
diff --git a/cmake/BornAgain/Dependences.cmake b/cmake/BornAgain/Dependences.cmake
index 583d35375c7..02f7790b4f4 100644
--- a/cmake/BornAgain/Dependences.cmake
+++ b/cmake/BornAgain/Dependences.cmake
@@ -8,8 +8,8 @@ find_package(LibHeinz REQUIRED)
message(STATUS "LibHeinz: found=${LibHeinz_FOUND}, include_dirs=${LibHeinz_INCLUDE_DIR}, "
"version=${LibHeinz_VERSION}")
-find_package(formfactor REQUIRED)
-message(STATUS "formfactor: found=${formfactor_FOUND}, version=${formfactor_VERSION}")
+#find_package(formfactor REQUIRED)
+#message(STATUS "formfactor: found=${formfactor_FOUND}, version=${formfactor_VERSION}")
find_package(Threads REQUIRED)
find_package(FFTW3 REQUIRED)
diff --git a/ff/CMakeLists.txt b/ff/CMakeLists.txt
new file mode 100644
index 00000000000..8ed490a99b5
--- /dev/null
+++ b/ff/CMakeLists.txt
@@ -0,0 +1,33 @@
+set(lib formfactor0)
+set(${lib}_LIBRARY ${lib} PARENT_SCOPE)
+
+file(GLOB src_files *.cpp)
+set(api_files Polyhedron.h PolyhedralTopology.h PolyhedralComponents.h)
+
+add_library(${lib} ${src_files})
+
+target_include_directories(${lib}
+ PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>
+ $<INSTALL_INTERFACE:include>
+ )
+target_include_directories(${lib} PRIVATE ${LibHeinz_INCLUDE_DIR})
+
+set_target_properties(
+ ${lib} PROPERTIES
+ OUTPUT_NAME ${lib}
+ VERSION ${formfactor_VERSION}
+ SOVERSION ${formfactor_VERSION_MAJOR})
+
+install(
+ TARGETS ${lib}
+ EXPORT formfactorTargets
+ LIBRARY DESTINATION lib
+ RUNTIME DESTINATION lib
+ ARCHIVE DESTINATION lib
+ COMPONENT Libraries)
+
+install(
+ FILES ${api_files}
+ DESTINATION include/ff
+ COMPONENT Headers)
diff --git a/ff/Factorial.h b/ff/Factorial.h
new file mode 100644
index 00000000000..19d41493f18
--- /dev/null
+++ b/ff/Factorial.h
@@ -0,0 +1,59 @@
+// ************************************************************************************************
+//
+// libformfactor: efficient and accurate computation of scattering form factors
+//
+//! @file lib/Factorial.h
+//! @brief Precomputes 1/n!
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see LICENSE)
+//! @copyright Forschungszentrum Jülich GmbH 2021
+//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
+//
+// ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif
+
+#ifndef USER_API
+#ifndef LIBFORMFACTOR_LIB_FACTORIAL_H
+#define LIBFORMFACTOR_LIB_FACTORIAL_H
+
+#include <array>
+#include <cstddef>
+#include <utility>
+
+namespace {
+
+template <size_t N> struct ReciprocalFactorial {
+ static constexpr double value = ReciprocalFactorial<N - 1>::value / N;
+};
+
+template <> struct ReciprocalFactorial<0> {
+ static constexpr double value = 1.0;
+};
+
+template <template <size_t> class F, size_t... I>
+constexpr std::array<double, sizeof...(I)> generateArrayHelper(std::index_sequence<I...>)
+{
+ return {F<I>::value...};
+}
+
+} // namespace
+
+
+namespace ff_aux {
+
+//! Returns a compile-time generated std::array of reciprocal factorials.
+
+template <size_t N, typename Indices = std::make_index_sequence<N>>
+constexpr std::array<double, N> generateReciprocalFactorialArray()
+{
+ return generateArrayHelper<ReciprocalFactorial>(Indices{});
+}
+
+} // namespace ff_aux
+
+#endif // LIBFORMFACTOR_LIB_FACTORIAL_H
+#endif // USER_API
diff --git a/ff/PolyhedralComponents.cpp b/ff/PolyhedralComponents.cpp
new file mode 100644
index 00000000000..11a4eae2d7f
--- /dev/null
+++ b/ff/PolyhedralComponents.cpp
@@ -0,0 +1,397 @@
+// ************************************************************************************************
+//
+// libformfactor: efficient and accurate computation of scattering form factors
+//
+//! @file lib/PolyhedralComponents.cpp
+//! @brief Implements classes PolyhedralEdge, PolyhedralFace
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see LICENSE)
+//! @copyright Forschungszentrum Jülich GmbH 2021
+//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
+//
+// ************************************************************************************************
+
+#include "ff/PolyhedralComponents.h"
+#include "ff/Factorial.h"
+
+#include <iomanip>
+#include <stdexcept>
+
+namespace {
+
+const double eps = 2e-16;
+constexpr auto ReciprocalFactorialArray = ff_aux::generateReciprocalFactorialArray<171>();
+
+complex_t sinc(const complex_t z) // cardinal sine function, sin(x)/x
+{
+ // This is an exception from the rule that we must not test floating-point numbers for equality.
+ // For small non-zero arguments, sin(z) returns quite accurately z or z-z^3/6.
+ // There is no loss of precision in computing sin(z)/z.
+ // Therefore there is no need for an expensive test like abs(z)<eps.
+ if (z == complex_t(0., 0.))
+ return 1.0;
+ return std::sin(z) / z;
+}
+
+} // namespace
+
+
+#ifdef ALGORITHM_DIAGNOSTIC
+void ff::PolyhedralDiagnosis::reset()
+{
+ order = 0;
+ algo = 0;
+ msg.clear();
+};
+std::string ff::PolyhedralDiagnosis::message() const
+{
+ std::string ret = "algo=" + std::to_string(algo) + ", order=" + std::to_string(order);
+ if (!msg.empty())
+ ret += ", msg:\n" + msg;
+ return ret;
+}
+bool ff::PolyhedralDiagnosis::operator==(const PolyhedralDiagnosis& other) const
+{
+ return order == other.order && algo == other.algo;
+}
+bool ff::PolyhedralDiagnosis::operator!=(const PolyhedralDiagnosis& other) const
+{
+ return !(*this == other);
+}
+#endif
+
+// ************************************************************************************************
+// PolyhedralEdge implementation
+// ************************************************************************************************
+
+ff::PolyhedralEdge::PolyhedralEdge(R3 Vlow, R3 Vhig) : m_E((Vhig - Vlow) / 2), m_R((Vhig + Vlow) / 2)
+{
+ if (m_E.mag2() == 0)
+ throw std::invalid_argument("At least one edge has zero length");
+}
+
+//! Returns sum_l=0^M/2 u^2l v^(M-2l) / (2l+1)!(M-2l)! - vperp^M/M!
+
+complex_t ff::PolyhedralEdge::contrib(int M, C3 qpa, complex_t qrperp) const
+{
+ complex_t u = qE(qpa);
+ complex_t v2 = m_R.dot(qpa);
+ complex_t v1 = qrperp;
+ complex_t v = v2 + v1;
+ // std::cout << std::scientific << std::showpos << std::setprecision(16) << "contrib: u=" << u
+ // << " v1=" << v1 << " v2=" << v2 << "\n";
+ if (v == 0.) { // only 2l=M contributes
+ if (M & 1) // M is odd
+ return 0.;
+ return ReciprocalFactorialArray[M] * (pow(u, M) / (M + 1.) - pow(v1, M));
+ }
+ complex_t ret = 0;
+ // the l=0 term, minus (qperp.R)^M, which cancels under the sum over E*contrib()
+ if (v1 == 0.) {
+ ret = ReciprocalFactorialArray[M] * pow(v2, M);
+ } else if (v2 == 0.) {
+ ; // leave ret=0
+ } else {
+ // binomial expansion
+ for (int mm = 1; mm <= M; ++mm) {
+ complex_t term = ReciprocalFactorialArray[mm] * ReciprocalFactorialArray[M - mm]
+ * pow(v2, mm) * pow(v1, M - mm);
+ ret += term;
+ // std::cout << "contrib mm=" << mm << " t=" << term << " s=" << ret << "\n";
+ }
+ }
+ if (u == 0.)
+ return ret;
+ for (int l = 1; l <= M / 2; ++l) {
+ complex_t term = ReciprocalFactorialArray[M - 2 * l] * ReciprocalFactorialArray[2 * l + 1]
+ * pow(u, 2 * l) * pow(v, M - 2 * l);
+ ret += term;
+ // std::cout << "contrib l=" << l << " t=" << term << " s=" << ret << "\n";
+ }
+ return ret;
+}
+
+// ************************************************************************************************
+// PolyhedralFace implementation
+// ************************************************************************************************
+
+double ff::PolyhedralFace::qpa_limit_series = 1e-2;
+int ff::PolyhedralFace::n_limit_series = 20;
+
+//! Static method, returns diameter of circle that contains all vertices.
+
+double ff::PolyhedralFace::diameter(const std::vector<R3>& V)
+{
+ double diameterFace = 0;
+ for (size_t j = 0; j < V.size(); ++j)
+ for (size_t jj = j + 1; jj < V.size(); ++jj)
+ diameterFace = std::max(diameterFace, (V[j] - V[jj]).mag());
+ return diameterFace;
+}
+
+//! Sets internal variables for given vertex chain.
+
+//! @param V oriented vertex list
+//! @param _sym_S2 true if face has a perpedicular two-fold symmetry axis
+
+ff::PolyhedralFace::PolyhedralFace(const std::vector<R3>& V, bool _sym_S2) : sym_S2(_sym_S2)
+{
+ size_t NV = V.size();
+ if (!NV)
+ throw std::runtime_error("Invalid polyhedral face: no edges given");
+ if (NV < 3)
+ throw std::runtime_error("Invalid polyhedral face: less than three edges");
+
+ // compute radius in 2d and 3d
+ m_radius_2d = diameter(V) / 2;
+ m_radius_3d = 0;
+ for (const R3& v : V)
+ m_radius_3d = std::max(m_radius_3d, v.mag());
+
+ // Initialize list of 'edges'.
+ // Do not create an edge if two vertices are too close to each other.
+ // TODO This is implemented in a somewhat sloppy way: we just skip an edge if it would
+ // be too short. This leaves tiny open edges. In a clean implementation, we
+ // rather should merge adjacent vertices before generating edges.
+ for (size_t j = 0; j < NV; ++j) {
+ size_t jj = (j + 1) % NV;
+ if ((V[j] - V[jj]).mag() < 1e-14 * m_radius_2d)
+ continue; // distance too short -> skip this edge
+ edges.emplace_back(V[j], V[jj]);
+ }
+ size_t NE = edges.size();
+ if (NE < 3)
+ throw std::invalid_argument("Face has less than three non-vanishing edges");
+
+ // compute n_k, rperp
+ m_normal = R3();
+ for (size_t j = 0; j < NE; ++j) {
+ size_t jj = (j + 1) % NE;
+ R3 ee = edges[j].E().cross(edges[jj].E());
+ if (ee.mag2() == 0) {
+ throw std::runtime_error("Invalid polyhedral face: two adjacent edges are parallel");
+ }
+ m_normal += ee.unit();
+ }
+ m_normal /= NE;
+ m_rperp = 0;
+ for (size_t j = 0; j < NV; ++j)
+ m_rperp += V[j].dot(m_normal);
+ m_rperp /= NV;
+ // assert that the vertices lay in a plane
+ for (size_t j = 1; j < NV; ++j)
+ if (std::abs(V[j].dot(m_normal) - m_rperp) > 1e-14 * m_radius_3d)
+ throw std::runtime_error("Invalid polyhedral face: not planar");
+ // compute m_area
+ m_area = 0;
+ for (size_t j = 0; j < NV; ++j) {
+ size_t jj = (j + 1) % NV;
+ m_area += m_normal.dot(V[j].cross(V[jj])) / 2;
+ }
+ // only now deal with inversion symmetry
+ if (sym_S2) {
+ if (NE & 1)
+ throw std::runtime_error("Invalid polyhedral face: odd #edges violates symmetry S2");
+ NE /= 2;
+ for (size_t j = 0; j < NE; ++j) {
+ if (((edges[j].R() - m_rperp * m_normal) + (edges[j + NE].R() - m_rperp * m_normal))
+ .mag()
+ > 1e-12 * m_radius_2d)
+ throw std::runtime_error(
+ "Invalid polyhedral face: edge centers violate symmetry S2");
+ if ((edges[j].E() + edges[j + NE].E()).mag() > 1e-12 * m_radius_2d)
+ throw std::runtime_error(
+ "Invalid polyhedral face: edge vectors violate symmetry S2");
+ }
+ // keep only half of the egdes
+ edges.erase(edges.begin() + NE, edges.end());
+ }
+}
+
+//! Sets qperp and qpa according to argument q and to this polygon's normal.
+
+void ff::PolyhedralFace::decompose_q(C3 q, complex_t& qperp, C3& qpa) const
+{
+ qperp = m_normal.dot(q);
+ qpa = q - qperp * m_normal;
+ // improve numeric accuracy:
+ qpa -= m_normal.dot(qpa) * m_normal;
+ if (qpa.mag() < eps * std::abs(qperp))
+ qpa = C3(0., 0., 0.);
+}
+
+//! Returns core contribution to f_n
+
+complex_t ff::PolyhedralFace::ff_n_core(int m, C3 qpa, complex_t qperp) const
+{
+ const C3 prevec = 2. * m_normal.cross(qpa); // complex conjugation not here but in .dot
+ complex_t ret = 0;
+ const complex_t qrperp = qperp * m_rperp;
+ for (size_t i = 0; i < edges.size(); ++i) {
+ const PolyhedralEdge& e = edges[i];
+ const complex_t vfac = prevec.dot(e.E());
+ const complex_t tmp = e.contrib(m + 1, qpa, qrperp);
+ ret += vfac * tmp;
+ // std::cout << std::scientific << std::showpos << std::setprecision(16)
+ // << "DBX ff_n_core " << m << " " << vfac << " " << tmp
+ // << " term=" << vfac * tmp << " sum=" << ret << "\n";
+ }
+ return ret;
+}
+
+//! Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor.
+
+complex_t ff::PolyhedralFace::ff_n(int n, C3 q) const
+{
+ complex_t qn = q.dot(m_normal); // conj(q)*normal (dot is antilinear in 'this' argument)
+ if (std::abs(qn) < eps * q.mag())
+ return 0.;
+ complex_t qperp;
+ C3 qpa;
+ decompose_q(q, qperp, qpa);
+ double qpa_mag2 = qpa.mag2();
+ if (qpa_mag2 == 0.)
+ return qn * pow(qperp * m_rperp, n) * m_area * ReciprocalFactorialArray[n];
+ if (sym_S2)
+ return qn * (ff_n_core(n, qpa, qperp) + ff_n_core(n, -qpa, qperp)) / qpa_mag2;
+ complex_t tmp = ff_n_core(n, qpa, qperp);
+ // std::cout << "DBX ff_n " << n << " " << qn << " " << tmp << " " << qpa_mag2 << "\n";
+ return qn * tmp / qpa_mag2;
+}
+
+//! Returns sum of n>=1 terms of qpa expansion of 2d form factor
+
+complex_t ff::PolyhedralFace::expansion(complex_t fac_even, complex_t fac_odd, C3 qpa,
+ double abslevel) const
+{
+#ifdef ALGORITHM_DIAGNOSTIC
+ polyhedralDiagnosis.algo += 1;
+#endif
+ complex_t sum = 0;
+ complex_t n_fac = I;
+ int count_return_condition = 0;
+ for (int n = 1; n < n_limit_series; ++n) {
+#ifdef ALGORITHM_DIAGNOSTIC
+ polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
+#endif
+ complex_t term = n_fac * (n & 1 ? fac_odd : fac_even) * ff_n_core(n, qpa, 0) / qpa.mag2();
+ // std::cout << std::setprecision(16) << " sum=" << sum << " +term=" << term << "\n";
+ sum += term;
+ if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * abslevel)
+ ++count_return_condition;
+ else
+ count_return_condition = 0;
+ if (count_return_condition > 2)
+ return sum; // regular exit
+ n_fac = mul_I(n_fac);
+ }
+ throw std::runtime_error("Numeric error in polyhedral face: series f(q_pa) not converged");
+}
+
+//! Returns core contribution to analytic 2d form factor.
+
+complex_t ff::PolyhedralFace::edge_sum_ff(C3 q, C3 qpa, bool sym_Ci) const
+{
+ C3 prevec = m_normal.cross(qpa); // complex conjugation will take place in .dot
+ complex_t sum = 0;
+ complex_t vfacsum = 0;
+ for (size_t i = 0; i < edges.size(); ++i) {
+ const ff::PolyhedralEdge& e = edges[i];
+ complex_t qE = e.qE(qpa);
+ complex_t qR = e.qR(qpa);
+ complex_t Rfac = sym_S2 ? sin(qR) : (sym_Ci ? cos(e.qR(q)) : exp_I(qR));
+ complex_t vfac;
+ if (sym_S2 || i < edges.size() - 1) {
+ vfac = prevec.dot(e.E());
+ vfacsum += vfac;
+ } else {
+ vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
+ }
+ complex_t term = vfac * sinc(qE) * Rfac;
+ sum += term;
+ // std::cout << std::scientific << std::showpos << std::setprecision(16)
+ // << " sum=" << sum << " term=" << term << " vf=" << vfac << " qE=" << qE
+ // << " qR=" << qR << " sinc=" << sinc(qE) << " Rfac=" << Rfac << "\n";
+ }
+ return sum;
+}
+
+//! Returns the contribution ff(q) of this face to the polyhedral form factor.
+
+complex_t ff::PolyhedralFace::ff(C3 q, bool sym_Ci) const
+{
+ complex_t qperp;
+ C3 qpa;
+ decompose_q(q, qperp, qpa);
+ double qpa_red = m_radius_2d * qpa.mag();
+ complex_t qr_perp = qperp * m_rperp;
+ complex_t ff0 = (sym_Ci ? 2. * I * sin(qr_perp) : exp_I(qr_perp)) * m_area;
+ if (qpa_red == 0)
+ return ff0;
+ if (qpa_red < qpa_limit_series && !sym_S2) {
+ // summation of power series
+ complex_t fac_even;
+ complex_t fac_odd;
+ if (sym_Ci) {
+ fac_even = 2. * mul_I(sin(qr_perp));
+ fac_odd = 2. * cos(qr_perp);
+ } else {
+ fac_even = exp_I(qr_perp);
+ fac_odd = fac_even;
+ }
+ return ff0 + expansion(fac_even, fac_odd, qpa, std::abs(ff0));
+ }
+ // direct evaluation of analytic formula
+ complex_t prefac;
+ if (sym_S2)
+ prefac = sym_Ci ? -8. * sin(qr_perp) : 4. * mul_I(exp_I(qr_perp));
+ else
+ prefac = sym_Ci ? 4. : 2. * exp_I(qr_perp);
+ // std::cout << " qrperp=" << qr_perp << " => prefac=" << prefac << "\n";
+ return prefac * edge_sum_ff(q, qpa, sym_Ci) / mul_I(qpa.mag2());
+}
+
+//! Two-dimensional form factor, for use in prism, from power series.
+
+complex_t ff::PolyhedralFace::ff_2D_expanded(C3 qpa) const
+{
+ return m_area + expansion(1., 1., qpa, std::abs(m_area));
+}
+
+//! Two-dimensional form factor, for use in prism, from sum over edge form factors.
+
+complex_t ff::PolyhedralFace::ff_2D_direct(C3 qpa) const
+{
+ return (sym_S2 ? 4. : 2. / I) * edge_sum_ff(qpa, qpa, false) / qpa.mag2();
+}
+
+//! Returns the two-dimensional form factor of this face, for use in a prism.
+
+complex_t ff::PolyhedralFace::ff_2D(C3 qpa) const
+{
+ if (std::abs(qpa.dot(m_normal)) > eps * qpa.mag())
+ throw std::runtime_error(
+ "Numeric error in polyhedral formfactor: ff_2D called with perpendicular q component");
+ double qpa_red = m_radius_2d * qpa.mag();
+ if (qpa_red == 0)
+ return m_area;
+ if (qpa_red < qpa_limit_series && !sym_S2)
+ return ff_2D_expanded(qpa);
+ return ff_2D_direct(qpa);
+}
+
+//! Throws if deviation from inversion symmetry is detected. Does not check vertices.
+
+void ff::PolyhedralFace::assert_Ci(const PolyhedralFace& other) const
+{
+ if (std::abs(m_rperp - other.m_rperp) > 1e-15 * (m_rperp + other.m_rperp))
+ throw std::runtime_error(
+ "Invalid polyhedron: faces with different distance from origin violate symmetry Ci");
+ if (std::abs(m_area - other.m_area) > 1e-15 * (m_area + other.m_area))
+ throw std::runtime_error(
+ "Invalid polyhedron: faces with different areas violate symmetry Ci");
+ if ((m_normal + other.m_normal).mag() > 1e-14)
+ throw std::runtime_error(
+ "Invalid polyhedron: faces do not have opposite orientation, violating symmetry Ci");
+}
diff --git a/ff/PolyhedralComponents.h b/ff/PolyhedralComponents.h
new file mode 100644
index 00000000000..41405221ce0
--- /dev/null
+++ b/ff/PolyhedralComponents.h
@@ -0,0 +1,103 @@
+// ************************************************************************************************
+//
+// libformfactor: efficient and accurate computation of scattering form factors
+//
+//! @file lib/PolyhedralComponents.h
+//! @brief Defines classes PolyhedralEdge, PolyhedralFace
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see LICENSE)
+//! @copyright Forschungszentrum Jülich GmbH 2021
+//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
+//
+// ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif
+
+#ifndef USER_API
+#ifndef LIBFORMFACTOR_LIB_POLYHEDRALCOMPONENTS_H
+#define LIBFORMFACTOR_LIB_POLYHEDRALCOMPONENTS_H
+
+#include <heinz/Complex.h>
+#include <heinz/Vectors3D.h>
+#include <vector>
+
+namespace ff {
+
+#ifdef ALGORITHM_DIAGNOSTIC
+#include <string>
+
+struct PolyhedralDiagnosis {
+ int algo;
+ int order;
+ std::string msg;
+ void reset();
+ std::string message() const;
+ bool operator==(const PolyhedralDiagnosis&) const;
+ bool operator!=(const PolyhedralDiagnosis&) const;
+};
+inline PolyhedralDiagnosis polyhedralDiagnosis;
+#endif
+
+//! One edge of a polygon, for form factor computation.
+
+class PolyhedralEdge {
+public:
+ PolyhedralEdge(R3 Vlow, R3 Vhig);
+
+ R3 E() const { return m_E; }
+ R3 R() const { return m_R; }
+ complex_t qE(C3 q) const { return m_E.dot(q); }
+ complex_t qR(C3 q) const { return m_R.dot(q); }
+
+ complex_t contrib(int m, C3 qpa, complex_t qrperp) const;
+
+private:
+ R3 m_E; //!< vector pointing from mid of edge to upper vertex
+ R3 m_R; //!< position vector of edge midpoint
+};
+
+//! A polygon, for form factor computation.
+
+class PolyhedralFace {
+public:
+ static double diameter(const std::vector<R3>& V);
+
+ PolyhedralFace(const std::vector<R3>& _V = std::vector<R3>(), bool _sym_S2 = false);
+
+ double area() const { return m_area; }
+ double pyramidalVolume() const { return m_rperp * m_area / 3; }
+ double radius3d() const { return m_radius_3d; }
+ //! Returns conj(q)*normal [BasicVector3D::dot is antilinear in 'this' argument]
+ complex_t normalProjectionConj(C3 q) const { return q.dot(m_normal); }
+ complex_t ff_n(int n, C3 q) const;
+ complex_t ff(C3 q, bool sym_Ci) const;
+ complex_t ff_2D(C3 qpa) const;
+ complex_t ff_2D_direct(C3 qpa) const; // for TestTriangle
+ complex_t ff_2D_expanded(C3 qpa) const; // for TestTriangle
+ void assert_Ci(const PolyhedralFace& other) const;
+
+private:
+ static double qpa_limit_series; //!< determines when use power series
+ static int n_limit_series;
+
+ bool sym_S2; //!< if true, then edges obtainable by inversion are not provided
+ std::vector<PolyhedralEdge> edges;
+ double m_area;
+ R3 m_normal; //!< normal vector of this polygon's plane
+ double m_rperp; //!< distance of this polygon's plane from the origin, along 'm_normal'
+ double m_radius_2d; //!< radius of enclosing cylinder
+ double m_radius_3d; //!< radius of enclosing sphere
+
+ void decompose_q(C3 q, complex_t& qperp, C3& qpa) const;
+ complex_t ff_n_core(int m, C3 qpa, complex_t qperp) const;
+ complex_t edge_sum_ff(C3 q, C3 qpa, bool sym_Ci) const;
+ complex_t expansion(complex_t fac_even, complex_t fac_odd, C3 qpa, double abslevel) const;
+};
+
+} // namespace ff
+
+#endif // LIBFORMFACTOR_LIB_POLYHEDRALCOMPONENTS_H
+#endif // USER_API
diff --git a/ff/PolyhedralTopology.h b/ff/PolyhedralTopology.h
new file mode 100644
index 00000000000..b3e4b99ba2d
--- /dev/null
+++ b/ff/PolyhedralTopology.h
@@ -0,0 +1,44 @@
+// ************************************************************************************************
+//
+// libformfactor: efficient and accurate computation of scattering form factors
+//
+//! @file lib/PolyhedralTopology.h
+//! @brief Defines classes PolygonalTopology, PolyhedralTopology
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see LICENSE)
+//! @copyright Forschungszentrum Jülich GmbH 2021
+//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
+//
+// ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif
+
+#ifndef USER_API
+#ifndef LIBFORMFACTOR_LIB_POLYHEDRALTOPOLOGY_H
+#define LIBFORMFACTOR_LIB_POLYHEDRALTOPOLOGY_H
+
+#include <vector>
+
+namespace ff {
+
+//! For internal use in PolyhedralFace.
+class PolygonalTopology {
+public:
+ std::vector<int> vertexIndices;
+ bool symmetry_S2;
+};
+
+//! For internal use in IFormFactorPolyhedron.
+class PolyhedralTopology {
+public:
+ std::vector<PolygonalTopology> faces;
+ bool symmetry_Ci;
+};
+
+} // namespace ff
+
+#endif // LIBFORMFACTOR_LIB_POLYHEDRALTOPOLOGY_H
+#endif // USER_API
diff --git a/ff/Polyhedron.cpp b/ff/Polyhedron.cpp
new file mode 100644
index 00000000000..87584be8bc5
--- /dev/null
+++ b/ff/Polyhedron.cpp
@@ -0,0 +1,192 @@
+// ************************************************************************************************
+//
+// libformfactor: efficient and accurate computation of scattering form factors
+//
+//! @file lib/Polyhedron.cpp
+//! @brief Implements class Polyhedron.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see LICENSE)
+//! @copyright Forschungszentrum Jülich GmbH 2021
+//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
+//
+// ************************************************************************************************
+
+//! The mathematics implemented here is described in full detail in a paper
+//! by Joachim Wuttke, entitled
+//! "Numerically stable form factor of any polygon and polyhedron"
+
+#include "ff/Polyhedron.h"
+#include "ff/PolyhedralComponents.h"
+#include "ff/PolyhedralTopology.h"
+#include <stdexcept>
+
+#ifdef ALGORITHM_DIAGNOSTIC_LEVEL2
+#include <boost/format.hpp>
+#endif
+
+namespace {
+
+const double eps = 2e-16;
+const double q_limit_series = 1e-2;
+const int n_limit_series = 20;
+
+} // namespace
+
+
+ff::Polyhedron::Polyhedron(const PolyhedralTopology& topology, double z_bottom,
+ const std::vector<R3>& vertices)
+{
+ m_vertices.clear();
+ for (const R3& vertex : vertices)
+ m_vertices.push_back(vertex - R3{0, 0, z_bottom});
+
+ m_z_bottom = z_bottom;
+ m_sym_Ci = topology.symmetry_Ci;
+
+ double diameter = 0;
+ for (size_t j = 0; j < vertices.size(); ++j)
+ for (size_t jj = j + 1; jj < vertices.size(); ++jj)
+ diameter = std::max(diameter, (vertices[j] - vertices[jj]).mag());
+
+ m_faces.clear();
+ for (const PolygonalTopology& tf : topology.faces) {
+ std::vector<R3> corners; // of one face
+ for (int i : tf.vertexIndices)
+ corners.push_back(vertices[i]);
+ if (PolyhedralFace::diameter(corners) <= 1e-14 * diameter)
+ continue; // skip ridiculously small face
+ m_faces.emplace_back(corners, tf.symmetry_S2);
+ }
+ if (m_faces.size() < 4)
+ throw std::runtime_error("Invalid polyhedron: less than four non-vanishing faces");
+
+ m_radius = 0;
+ m_volume = 0;
+ for (const PolyhedralFace& Gk : m_faces) {
+ m_radius = std::max(m_radius, Gk.radius3d());
+ m_volume += Gk.pyramidalVolume();
+ }
+ if (m_sym_Ci) {
+ if (m_faces.size() & 1)
+ throw std::runtime_error("Invalid polyhedron: odd #faces violates symmetry Ci");
+ size_t N = m_faces.size() / 2;
+ // for this tests, m_faces must be in a specific order
+ for (size_t k = 0; k < N; ++k)
+ m_faces[k].assert_Ci(m_faces[2 * N - 1 - k]);
+ // keep only half of the faces
+ m_faces.erase(m_faces.begin() + N, m_faces.end());
+ }
+}
+
+void ff::Polyhedron::assert_platonic() const
+{
+ // just one test; one could do much more ...
+ double pyramidal_volume = 0;
+ for (const auto& Gk : m_faces)
+ pyramidal_volume += Gk.pyramidalVolume();
+ pyramidal_volume /= m_faces.size();
+ for (const auto& Gk : m_faces)
+ if (std::abs(Gk.pyramidalVolume() - pyramidal_volume) > 160 * eps * pyramidal_volume)
+ throw std::runtime_error(
+ "Invalid Polyhedron: declared platonic but not sufficiently uniform");
+}
+
+double ff::Polyhedron::volume() const
+{
+ return m_volume;
+}
+
+double ff::Polyhedron::radius() const
+{
+ return m_radius;
+}
+
+const std::vector<R3> ff::Polyhedron::vertices() const
+{
+ std::vector<R3> ret;
+ ret.reserve(m_vertices.size());
+ for (const auto& vertex : m_vertices)
+ ret.emplace_back(R3{vertex});
+ return ret;
+}
+
+//! Returns the form factor F(q) of this polyhedron, respecting the offset z_bottom.
+
+complex_t ff::Polyhedron::evaluate_for_q(const C3& _q) const
+{
+ C3 q{_q};
+ return exp_I(-m_z_bottom * q.z()) * evaluate_centered(q);
+}
+
+//! Returns the form factor F(q) of this polyhedron, with origin at z=0.
+
+complex_t ff::Polyhedron::evaluate_centered(const C3& _q) const
+{
+ C3 q{_q};
+ double q_red = m_radius * q.mag();
+#ifdef ALGORITHM_DIAGNOSTIC
+ polyhedralDiagnosis.reset();
+#endif
+ if (q_red == 0)
+ return m_volume;
+ if (q_red < q_limit_series) {
+ // summation of power series
+#ifdef ALGORITHM_DIAGNOSTIC
+ polyhedralDiagnosis.algo = 100;
+#endif
+ complex_t sum = 0;
+ complex_t n_fac = (m_sym_Ci ? -2 : -1) / q.mag2();
+ int count_return_condition = 0;
+ for (int n = 2; n < n_limit_series; ++n) {
+ if (m_sym_Ci && n & 1)
+ continue;
+#ifdef ALGORITHM_DIAGNOSTIC
+ polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
+#endif
+ complex_t term = 0;
+ for (const PolyhedralFace& Gk : m_faces) {
+ complex_t tmp = Gk.ff_n(n + 1, q);
+ term += tmp;
+ }
+ term *= n_fac;
+#ifdef ALGORITHM_DIAGNOSTIC_LEVEL2
+ polyhedralDiagnosis.msg +=
+ boost::str(boost::format(" + term(n=%2i) = %23.17e+i*%23.17e\n") % n % term.real()
+ % term.imag());
+#endif
+ sum += term;
+ if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * m_volume)
+ ++count_return_condition;
+ else
+ count_return_condition = 0;
+ if (count_return_condition > 2)
+ return m_volume + sum; // regular exit
+ n_fac = m_sym_Ci ? -n_fac : mul_I(n_fac);
+ }
+ throw std::runtime_error("Numeric failure in polyhedron: series F(q) not converged");
+ }
+
+ // direct evaluation of analytic formula (coefficients may involve series)
+#ifdef ALGORITHM_DIAGNOSTIC
+ polyhedralDiagnosis.algo = 200;
+#endif
+ complex_t sum = 0;
+ for (const PolyhedralFace& Gk : m_faces) {
+ complex_t qn = Gk.normalProjectionConj(q); // conj(q)*normal
+ if (std::abs(qn) < eps * q.mag())
+ continue;
+ complex_t term = qn * Gk.ff(q, m_sym_Ci);
+#ifdef ALGORITHM_DIAGNOSTIC //_LEVEL2
+ polyhedralDiagnosis.msg += boost::str(boost::format(" + face_ff = %23.17e+i*%23.17e\n")
+ % term.real() % term.imag());
+#endif
+ sum += term;
+ }
+#ifdef ALGORITHM_DIAGNOSTIC //_LEVEL2
+ polyhedralDiagnosis.msg +=
+ boost::str(boost::format(" -> raw sum = %23.17e+i*%23.17e; divisor = %23.17e\n")
+ % sum.real() % sum.imag() % q.mag2());
+#endif
+ return sum / I / q.mag2();
+}
diff --git a/ff/Polyhedron.h b/ff/Polyhedron.h
new file mode 100644
index 00000000000..c8cc56a2438
--- /dev/null
+++ b/ff/Polyhedron.h
@@ -0,0 +1,64 @@
+// ************************************************************************************************
+//
+// libformfactor: efficient and accurate computation of scattering form factors
+//
+//! @file lib/Polyhedron.h
+//! @brief Defines class Polyhedron.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see LICENSE)
+//! @copyright Forschungszentrum Jülich GmbH 2021
+//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
+//
+// ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif
+
+#ifndef USER_API
+#ifndef LIBFORMFACTOR_LIB_POLYHEDRON_H
+#define LIBFORMFACTOR_LIB_POLYHEDRON_H
+
+#include <array>
+#include <memory>
+#include <vector>
+
+#include <heinz/Complex.h>
+#include <heinz/Vectors3D.h>
+
+namespace ff {
+
+class PolyhedralFace;
+class PolyhedralTopology;
+
+//! A polyhedron, implementation class for use in IFormFactorPolyhedron
+
+class Polyhedron {
+public:
+ Polyhedron(const PolyhedralTopology& topology, double z_bottom,
+ const std::vector<R3>& vertices);
+ Polyhedron(const Polyhedron&) = delete;
+
+ void assert_platonic() const;
+ double volume() const;
+ double radius() const;
+
+ const std::vector<R3> vertices() const; //! needed for topZ, bottomZ computation
+ complex_t evaluate_for_q(const C3& q) const;
+ complex_t evaluate_centered(const C3& q) const;
+
+private:
+ double m_z_bottom;
+ bool m_sym_Ci; //!< if true, then faces obtainable by inversion are not provided
+
+ std::vector<PolyhedralFace> m_faces;
+ double m_radius;
+ double m_volume;
+ std::vector<R3> m_vertices; //! for topZ, bottomZ computation only
+};
+
+} // namespace ff
+
+#endif // LIBFORMFACTOR_LIB_POLYHEDRON_H
+#endif // USER_API
--
GitLab