1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
// ************************************************************************************************
//
// libformfactor: efficient and accurate computation of scattering form factors
//
//! @file ff/Prism.cpp
//! @brief Implements class Prism.
//!
//! @homepage https://jugit.fz-juelich.de/mlz/libformfactor
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2022
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
//! The mathematics implemented here is described in full detail in a paper
//! by Joachim Wuttke, entitled
//! "Form factor (Fourier shape transform) of polygon and polyhedron."
#include "ff/Prism.h"
#include "ff/Edge.h"
#include "ff/Face.h"
#include "ff/Polyhedron.h"
#include "ff/Topology.h"
#include <algorithm>
#include <numeric>
#include <stdexcept>
namespace {
complex_t sinc(const complex_t z) // cardinal sine function, sin(x)/x
{
if (z == complex_t(0., 0.))
return 1.0;
return std::sin(z) / z;
}
ff::Face precompute_base(bool symmetry_Ci, const std::vector<R3>& base_vertices)
{
try {
return ff::Face(base_vertices, symmetry_Ci);
} catch (std::invalid_argument& e) {
throw std::invalid_argument(std::string("Invalid parameterization of Prism: ") + e.what());
} catch (std::logic_error& e) {
throw std::logic_error(std::string("Bug in Prism: ") + e.what()
+ " [please report to the maintainers]");
} catch (std::exception& e) {
throw std::runtime_error(std::string("Unexpected exception in Prism: ") + e.what()
+ " [please report to the maintainers]");
}
}
} // namespace
ff::Prism::Prism(bool symmetry_Ci, double height, const std::vector<R3>& base_vertices,
const R3& location)
: IBody(location)
, m_height(height)
, m_base_vertices(base_vertices)
, m_base(std::make_unique<const ff::Face>(precompute_base(symmetry_Ci, base_vertices)))
, m_radius(std::hypot(height / 2, m_base->radius()))
, m_volume(height * m_base->area())
{
}
ff::Prism::~Prism() = default;
const std::vector<R3>& ff::Prism::vertices() const
{
if (my_vertices.empty()) {
const int N = m_base_vertices.size();
my_vertices.resize(2 * N);
for (int i = 0; i < N; i++) {
my_vertices[i] = m_base_vertices[i] + m_height / 2 * m_base->normal() + location();
my_vertices[i + N] = m_base_vertices[i] - m_height / 2 * m_base->normal() + location();
}
}
return my_vertices;
}
const ff::Topology& ff::Prism::topology() const
{
if (!my_topology) {
const int N = m_base_vertices.size();
std::vector<int> vertex_indices_floor(N);
std::vector<int> vertex_indices_ceil(N);
for (int i = 0; i < N; i++) {
vertex_indices_floor[i] = i;
vertex_indices_ceil[i] = (2 * N - 1) - i;
}
std::vector<FacialTopology> faces(N + 2);
faces[0] = {vertex_indices_floor, m_base->is_symmetric()};
faces[N + 1] = {vertex_indices_ceil, m_base->is_symmetric()};
for (int i = 1; i < N + 1; i++) {
int j = i % N;
faces[i] = {{i + N - 1, j + N, j, i - 1}, true};
}
if (m_base->is_symmetric())
std::reverse(faces.begin() + (N / 2) + 1, faces.end() - 1);
my_topology.reset(new Topology{faces, m_base->is_symmetric()});
}
return *my_topology;
}
double ff::Prism::z_bottom() const
{
return -m_height / 2;
}
complex_t ff::Prism::formfactor_at_center(const C3& q) const
{
try {
#ifdef ALGORITHM_DIAGNOSTIC
polyhedralDiagnosis.reset();
polyhedralDiagnosis.algo = 500;
#endif
C3 qxy(q.x(), q.y(), 0.);
return m_height * sinc(m_height / 2 * q.z()) * m_base->ff_2D(qxy);
} catch (std::logic_error& e) {
throw std::logic_error(std::string("Bug in Prism: ") + e.what()
+ " [please report to the maintainers]");
} catch (std::runtime_error& e) {
throw std::runtime_error(std::string("Numeric computation failed in Prism: ") + e.what()
+ " [please report to the maintainers]");
} catch (std::exception& e) {
throw std::runtime_error(std::string("Unexpected exception in Prism: ") + e.what()
+ " [please report to the maintainers]");
}
}
ff::Polyhedron* ff::Prism::asPolyhedron() const
{
std::vector<R3> vv = vertices();
for (R3& v : vv)
v -= location();
return new Polyhedron(topology(), vv, location());
}
|