VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
Namespaces | Classes | Typedefs | Enumerations | Functions | Variables
fmm Namespace Reference

Namespaces

namespace  detail
 
namespace  detail_interp
 

Classes

class  BinomNewton
 
class  FastMultipole
 
class  FastMultipole3d
 
class  MortonTree
 
struct  particle
 
struct  small_particle
 
struct  TreeCell
 
struct  treevector
 
class  Vector3
 

Typedefs

using Vector3d = Vector3< double >
 
using Matrix3d = std::array< Vector3d, 3 >
 
using Vector3cd = Vector3< std::complex< double > >
 
using point2d = std::complex< double >
 
using point3d = Vector3d
 
using particle2d = particle< point2d, double >
 
using particle3d = particle< point3d, double >
 
using particle3d3 = particle< point3d, point3d >
 
using TreeCell2d = TreeCell< point2d, double >
 
using TreeCell3d = TreeCell< point3d, double >
 
using TreeCell3d3 = TreeCell< point3d, point3d >
 
using MortonTree2d = MortonTree< 2, point2d, double >
 
using MortonTree3d = MortonTree< 3, point3d, double >
 
using MortonTree3d3 = MortonTree< 3, point3d, point3d >
 

Enumerations

enum class  InteractionType {
  Potential2d , Force2d , Potential3d , Force3d ,
  VelDipole3d , MomDipole3d
}
 

Functions

template<typename vec_type >
vec_type avx_zero_vec ()
 
template<typename vec_type >
void avx_store (double *data, vec_type &v)
 
template<typename vec_type >
vec_type avx_add (const vec_type &v1, const vec_type &v2)
 
template<typename vec_type >
vec_type avx_sub (const vec_type &v1, const vec_type &v2)
 
template<typename vec_type >
vec_type avx_mul (const vec_type &v1, const vec_type &v2)
 
double avx_hsum (__m128d v)
 
double avx_hsum (__m256d v)
 
double avx_hsum (__m512d v)
 
template<typename vec_type >
void avx_set_range (vec_type &vec, double *val)
 
template<typename vec_type >
void avx_set_constant (vec_type &vec, double val)
 
void avx_inv_dr (const __m512d &oneVec, const __m512d &epsVec, const __m512d &x_target, const __m512d &y_target, const __m512d &z_target, __m512d &temp1, __m512d &temp2, __m512d &temp3, __m512d &dx, __m512d &dy, __m512d &dz, __m512d &invdr, __m512d &invdr2)
 
void avx_inv_dr (const __m256d &oneVec, const __m256d &epsVec, const __m256d &x_target, const __m256d &y_target, const __m256d &z_target, __m256d &temp1, __m256d &temp2, __m256d &temp3, __m256d &dx, __m256d &dy, __m256d &dz, __m256d &invdr, __m256d &invdr2)
 
void avx_inv_dr (const __m128d &oneVec, const __m128d &epsVec, const __m128d &x_target, const __m128d &y_target, const __m128d &z_target, __m128d &temp1, __m128d &temp2, __m128d &temp3, __m128d &dx, __m128d &dy, __m128d &dz, __m128d &invdr, __m128d &invdr2)
 
bool IAmRoot ()
 
template<typename T >
constexpr T MyPow (T base, unsigned int exp)
 
constexpr double ni (int n)
 
Vector3d DecToSph (const Vector3d &vec)
 
double Potential2d (const particle2d &p1, const particle2d &p2)
 
std::complex< double > Force2d (const particle2d &p1, const particle2d &p2)
 
double Potential3d (const particle3d &p1, const particle3d &p2)
 
Vector3d Force3d (const particle3d &p1, const particle3d &p2)
 
Vector3d Force3d (const particle3d3 &p1, const particle3d3 &p2)
 
template<typename T >
std::pair< T, T > lambdaEta (T xi)
 
template<>
std::pair< double, double > lambdaEta< double > (double xi)
 
template<typename T >
std::pair< T, T > lambdaEtaA (T xi)
 
template<>
std::pair< double, double > lambdaEtaA< double > (double xi)
 
Vector3d velDipole3d (const particle3d3 &p1, const particle3d3 &p2)
 
Vector3d momDipole3d (const particle3d3 &p1, const particle3d3 &p2)
 
void Potential3dMutual (const particle3d &p1, const particle3d &p2, double &potential1, double &potential2)
 
void Force3dMutual (const particle3d &p1, const particle3d &p2, Vector3d &force1, Vector3d &force2)
 
void Force3dMutual (const particle3d3 &p1, const particle3d3 &p2, Vector3d &force1, Vector3d &force2)
 
void Potential2dMutual (const particle2d &p1, const particle2d &p2, double &potential1, double &potential2)
 
void Force2dMutual (const particle2d &p1, const particle2d &p2, std::complex< double > &force1, std::complex< double > &force2)
 
template<typename point_type , typename interaction_type , typename value_type >
void ComputeExact (const std::vector< particle< point_type, value_type > > &source_particles, const std::vector< particle< point_type, value_type > > &target_particles, std::function< interaction_type(const particle< point_type, value_type > &, const particle< point_type, value_type > &)> func, std::string filename)
 
template<InteractionType it, typename point_type , typename value_type >
void ComputeExact (const std::vector< particle< point_type, value_type > > &source_particles, const std::vector< particle< point_type, value_type > > &target_particles)
 
template<InteractionType it, typename point_type , typename value_type >
void ComputeExact (const std::vector< particle< point_type, value_type > > &particles)
 
template<InteractionType it, typename T >
void ReadError (const std::vector< T > &res)
 
void ReadErrorHyb (const std::vector< std::array< double, 2 > > &res)
 
treevector< double > ComputeKnm (int N)
 
treevector< double > ComputeAnm (int N)
 
std::vector< double > compute_wignerd_coef (int N)
 
std::vector< double > ComputeWignerD (int N, double theta)
 
std::vector< std::complex< double > > ComputeExp (int N, double phi)
 
void ComputePnm (int N, double theta, double *Pnm)
 
template<typename T >
constexpr Vector3< T > operator* (const double x, const Vector3< T > &rhs)
 
constexpr Vector3cd operator* (const std::complex< double > &x, const Vector3d &rhs)
 
constexpr Vector3cd operator* (const Vector3d &rhs, const std::complex< double > &x)
 
constexpr Vector3cd operator* (const std::complex< double > &x, const Vector3cd &rhs)
 
constexpr Vector3cd operator* (const Vector3cd &rhs, const std::complex< double > &x)
 
void operator*= (Vector3cd &lhs, const std::complex< double > &x)
 
void operator+= (Vector3cd &lhs, const Vector3d &rhs)
 
template<typename complex >
constexpr Vector3d real (const Vector3cd &rhs)
 
template<typename complex >
constexpr Vector3d imag (const Vector3cd &rhs)
 
template<typename T >
constexpr Vector3< T > cross (const Vector3< T > &lhs, const Vector3< T > &rhs)
 
constexpr double dot (const Vector3d &lhs, const Vector3d &rhs)
 
template<typename T >
std::ostream & operator<< (std::ostream &os, const Vector3< T > &rhs)
 
template<typename T >
std::istream & operator>> (std::istream &is, Vector3< T > &rhs)
 
template<typename T >
double norm (const Vector3< T > &vec)
 
template<typename T >
double abs (const Vector3< T > &vec)
 
int doublehash (double val)
 
int doublehash (double val1, double val2)
 
template<typename point , typename value >
std::pair< point, double > box_extent (const std::vector< particle< point, value > > &particles)
 
template<int dim>
size_t ExpandBits (double x, int tree_depth)
 
template<int dim, typename point >
size_t MortonCode (point val, int tree_depth, point shift, double scale)
 
template<typename point >
point PointFromMortonCode (size_t val, int depth, point shift, double scale)
 
template<>
size_t ExpandBits< 2 > (double x, int tree_depth)
 
template<>
size_t ExpandBits< 3 > (double x, int tree_depth)
 
template<>
size_t MortonCode< 2, point2d > (point2d val, int tree_depth, point2d shift, double scale)
 
template<>
size_t MortonCode< 3, point3d > (point3d val, int tree_depth, point3d shift, double scale)
 
template<>
point2d PointFromMortonCode (size_t val, int depth, point2d shift, double scale)
 
template<>
point3d PointFromMortonCode (size_t val, int depth, point3d shift, double scale)
 
template<typename value >
std::pair< point2d, double > box_extent (const std::vector< particle< point2d, value > > &particles)
 
template<typename value >
std::pair< point3d, double > box_extent (const std::vector< particle< point3d, value > > &particles)
 
std::complex< double > divComp (const std::complex< double > &a, const std::complex< double > &b)
 
void mv (const Matrix3d &A, const Vector3d &x, Vector3d &y)
 
void mm (const Matrix3d &A, const Matrix3d &B, Matrix3d &C)
 
void mmt (const Matrix3d &A, const Matrix3d &B, Matrix3d &C)
 
void vt (const Vector3d &x, const Matrix3d *T, Matrix3d &M)
 

Variables

const int FMM_AUTO = std::numeric_limits<int>::max()
 
constexpr double FORCE_EPS = 0.01
 
constexpr double FORCE_EPS2 = FORCE_EPS * FORCE_EPS
 
constexpr double FORCE_EPS3 = FORCE_EPS2 * FORCE_EPS
 
constexpr double FORCE_EPS4 = FORCE_EPS2 * FORCE_EPS2
 
constexpr double FORCE_EPS5 = FORCE_EPS3 * FORCE_EPS2
 
constexpr double FORCE_EPS6 = FORCE_EPS3 * FORCE_EPS3
 
const auto wignerd_auxillary_coef = compute_wignerd_coef(100)
 

Typedef Documentation

◆ Matrix3d

using fmm::Matrix3d = typedef std::array<Vector3d, 3>

Definition at line 64 of file utils.h.

◆ MortonTree2d

using fmm::MortonTree2d = typedef MortonTree<2, point2d, double>

Definition at line 50 of file morton_tree.h.

◆ MortonTree3d

using fmm::MortonTree3d = typedef MortonTree<3, point3d, double>

Definition at line 51 of file morton_tree.h.

◆ MortonTree3d3

using fmm::MortonTree3d3 = typedef MortonTree<3, point3d, point3d>

Definition at line 52 of file morton_tree.h.

◆ particle2d

using fmm::particle2d = typedef particle<point2d, double>

Definition at line 155 of file utils.h.

◆ particle3d

using fmm::particle3d = typedef particle<point3d, double>

Definition at line 156 of file utils.h.

◆ particle3d3

using fmm::particle3d3 = typedef particle<point3d, point3d>

Definition at line 157 of file utils.h.

◆ point2d

using fmm::point2d = typedef std::complex<double>

Definition at line 153 of file utils.h.

◆ point3d

using fmm::point3d = typedef Vector3d

Definition at line 154 of file utils.h.

◆ TreeCell2d

using fmm::TreeCell2d = typedef TreeCell<point2d, double>

Definition at line 47 of file morton_tree.h.

◆ TreeCell3d

using fmm::TreeCell3d = typedef TreeCell<point3d, double>

Definition at line 48 of file morton_tree.h.

◆ TreeCell3d3

using fmm::TreeCell3d3 = typedef TreeCell<point3d, point3d>

Definition at line 49 of file morton_tree.h.

◆ Vector3cd

using fmm::Vector3cd = typedef Vector3<std::complex<double> >

Definition at line 65 of file utils.h.

◆ Vector3d

using fmm::Vector3d = typedef Vector3<double>

Definition at line 63 of file utils.h.

Enumeration Type Documentation

◆ InteractionType

enum class fmm::InteractionType
strong
Enumerator
Potential2d 
Force2d 
Potential3d 
Force3d 
VelDipole3d 
MomDipole3d 

Definition at line 20 of file utils.h.

20 {
21 Potential2d,
22 Force2d,
23 Potential3d,
24 Force3d,
25 VelDipole3d,
26 MomDipole3d
27};
fmm::Force2d
std::complex< double > Force2d(const particle2d &p1, const particle2d &p2)
Definition simple_math.h:127
fmm::Force3d
Vector3d Force3d(const particle3d &p1, const particle3d &p2)
Definition simple_math.h:138
fmm::Potential3d
double Potential3d(const particle3d &p1, const particle3d &p2)
Definition simple_math.h:133
fmm::Potential2d
double Potential2d(const particle2d &p1, const particle2d &p2)
Definition simple_math.h:122

Function Documentation

◆ abs()

template<typename T >
double fmm::abs ( const Vector3< T > &  vec)
inline

Definition at line 122 of file utils.h.

122{return sqrt(norm(vec));}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ avx_add()

template<typename vec_type >
vec_type fmm::avx_add ( const vec_type &  v1,
const vec_type &  v2 
)
inline

Definition at line 42 of file avx.h.

43 {
44 if constexpr (std::is_same_v<vec_type, __m128d>)
45 return _mm_add_pd(v1, v2);
46 if constexpr (std::is_same_v<vec_type, __m256d>)
47 return _mm256_add_pd(v1, v2);
48 if constexpr (std::is_same_v<vec_type, __m512d>)
49 return _mm512_add_pd(v1, v2);
50 }

◆ avx_hsum() [1/3]

double fmm::avx_hsum ( __m128d  v)
inline

Definition at line 74 of file avx.h.

74 {
75 return _mm_cvtsd_f64(_mm_add_sd(v, _mm_unpackhi_pd(v, v)));
76 }
Here is the caller graph for this function:

◆ avx_hsum() [2/3]

double fmm::avx_hsum ( __m256d  v)
inline

Definition at line 78 of file avx.h.

78 {
79 __m128d vlow = _mm256_castpd256_pd128(v);
80 vlow = _mm_add_pd(vlow, _mm256_extractf128_pd(v, 1));
81 return _mm_cvtsd_f64(_mm_add_sd(vlow, _mm_unpackhi_pd(vlow, vlow)));
82 }

◆ avx_hsum() [3/3]

double fmm::avx_hsum ( __m512d  v)
inline

Definition at line 84 of file avx.h.

84 {
85 __m256d vlow = _mm512_castpd512_pd256(v);
86 return avx_hsum(_mm256_add_pd(vlow, _mm512_extractf64x4_pd(v, 1)));
87 }
fmm::avx_hsum
double avx_hsum(__m128d v)
Definition avx.h:74
Here is the call graph for this function:

◆ avx_inv_dr() [1/3]

void fmm::avx_inv_dr ( const __m128d &  oneVec,
const __m128d &  epsVec,
const __m128d &  x_target,
const __m128d &  y_target,
const __m128d &  z_target,
__m128d &  temp1,
__m128d &  temp2,
__m128d &  temp3,
__m128d &  dx,
__m128d &  dy,
__m128d &  dz,
__m128d &  invdr,
__m128d &  invdr2 
)
inline

Definition at line 149 of file avx.h.

153 {
154 dx = _mm_sub_pd(x_target, temp1);
155 dy = _mm_sub_pd(y_target, temp2);
156 dz = _mm_sub_pd(z_target, temp3);
157 temp1 = _mm_mul_pd(dx, dx);
158 temp2 = _mm_mul_pd(dy, dy);
159 temp3 = _mm_mul_pd(dz, dz);
160
161 invdr = _mm_add_pd(temp1, temp2);
162 invdr = _mm_add_pd(invdr, temp3);
163 invdr = _mm_max_pd(invdr, epsVec);
164 invdr2 = _mm_div_pd(oneVec, invdr);
165 invdr = _mm_sqrt_pd(invdr2);
166 }

◆ avx_inv_dr() [2/3]

void fmm::avx_inv_dr ( const __m256d &  oneVec,
const __m256d &  epsVec,
const __m256d &  x_target,
const __m256d &  y_target,
const __m256d &  z_target,
__m256d &  temp1,
__m256d &  temp2,
__m256d &  temp3,
__m256d &  dx,
__m256d &  dy,
__m256d &  dz,
__m256d &  invdr,
__m256d &  invdr2 
)
inline

Definition at line 129 of file avx.h.

133 {
134 dx = _mm256_sub_pd(x_target, temp1);
135 dy = _mm256_sub_pd(y_target, temp2);
136 dz = _mm256_sub_pd(z_target, temp3);
137 temp1 = _mm256_mul_pd(dx, dx);
138 temp2 = _mm256_mul_pd(dy, dy);
139 temp3 = _mm256_mul_pd(dz, dz);
140
141 invdr = _mm256_add_pd(temp1, temp2);
142 invdr = _mm256_add_pd(invdr, temp3);
143
144 invdr = _mm256_max_pd(invdr, epsVec);
145 invdr2 = _mm256_div_pd(oneVec, invdr);
146 invdr = _mm256_sqrt_pd(invdr2);
147 }

◆ avx_inv_dr() [3/3]

void fmm::avx_inv_dr ( const __m512d &  oneVec,
const __m512d &  epsVec,
const __m512d &  x_target,
const __m512d &  y_target,
const __m512d &  z_target,
__m512d &  temp1,
__m512d &  temp2,
__m512d &  temp3,
__m512d &  dx,
__m512d &  dy,
__m512d &  dz,
__m512d &  invdr,
__m512d &  invdr2 
)
inline

Definition at line 111 of file avx.h.

115 {
116 dx = _mm512_sub_pd(x_target, temp1);
117 dy = _mm512_sub_pd(y_target, temp2);
118 dz = _mm512_sub_pd(z_target, temp3);
119 temp1 = _mm512_mul_pd(dx, dx);
120 temp2 = _mm512_mul_pd(dy, dy);
121 temp3 = _mm512_mul_pd(dz, dz);
122 invdr = _mm512_add_pd(temp1, temp2);
123 invdr = _mm512_add_pd(invdr, temp3);
124 invdr = _mm512_max_pd(invdr, epsVec);
125 invdr2 = _mm512_div_pd(oneVec, invdr);
126 invdr = _mm512_sqrt_pd(invdr2);
127 }

◆ avx_mul()

template<typename vec_type >
vec_type fmm::avx_mul ( const vec_type &  v1,
const vec_type &  v2 
)
inline

Definition at line 64 of file avx.h.

65 {
66 if constexpr (std::is_same_v<vec_type, __m128d>)
67 return _mm_mul_pd(v1, v2);
68 if constexpr (std::is_same_v<vec_type, __m256d>)
69 return _mm256_mul_pd(v1, v2);
70 if constexpr (std::is_same_v<vec_type, __m512d>)
71 return _mm512_mul_pd(v1, v2);
72 }

◆ avx_set_constant()

template<typename vec_type >
void fmm::avx_set_constant ( vec_type &  vec,
double  val 
)
inline

Definition at line 101 of file avx.h.

102 {
103 if constexpr (std::is_same_v<vec_type, __m128d>)
104 vec = _mm_set_pd(val, val);
105 if constexpr (std::is_same_v<vec_type, __m256d>)
106 vec = _mm256_set_pd(val, val, val, val);
107 if constexpr (std::is_same_v<vec_type, __m512d>)
108 vec = _mm512_set_pd(val, val, val, val, val, val, val, val);
109 }
Here is the caller graph for this function:

◆ avx_set_range()

template<typename vec_type >
void fmm::avx_set_range ( vec_type &  vec,
double *  val 
)
inline

Definition at line 90 of file avx.h.

91 {
92 if constexpr (std::is_same_v<vec_type, __m128d>)
93 vec = _mm_set_pd(val[1], val[0]);
94 if constexpr (std::is_same_v<vec_type, __m256d>)
95 vec = _mm256_set_pd(val[3], val[2], val[1], val[0]);
96 if constexpr (std::is_same_v<vec_type, __m512d>)
97 vec = _mm512_set_pd(val[7], val[6], val[5], val[4], val[3], val[2], val[1], val[0]);
98 }

◆ avx_store()

template<typename vec_type >
void fmm::avx_store ( double *  data,
vec_type &  v 
)
inline

Definition at line 31 of file avx.h.

32 {
33 if constexpr (std::is_same_v<vec_type, __m128d>)
34 _mm_store_pd(data, v);
35 if constexpr (std::is_same_v<vec_type, __m256d>)
36 _mm256_store_pd(data, v);
37 if constexpr (std::is_same_v<vec_type, __m512d>)
38 _mm512_store_pd(data, v);
39 }

◆ avx_sub()

template<typename vec_type >
vec_type fmm::avx_sub ( const vec_type &  v1,
const vec_type &  v2 
)
inline

Definition at line 53 of file avx.h.

54 {
55 if constexpr (std::is_same_v<vec_type, __m128d>)
56 return _mm_sub_pd(v1, v2);
57 if constexpr (std::is_same_v<vec_type, __m256d>)
58 return _mm256_sub_pd(v1, v2);
59 if constexpr (std::is_same_v<vec_type, __m512d>)
60 return _mm512_sub_pd(v1, v2);
61 }

◆ avx_zero_vec()

template<typename vec_type >
vec_type fmm::avx_zero_vec ( )

Definition at line 20 of file avx.h.

21 {
22 if constexpr (std::is_same_v<vec_type, __m128d>)
23 return _mm_setzero_pd();
24 if constexpr (std::is_same_v<vec_type, __m256d>)
25 return _mm256_setzero_pd();
26 if constexpr (std::is_same_v<vec_type, __m512d>)
27 return _mm512_setzero_pd();
28 }

◆ box_extent() [1/3]

template<typename point , typename value >
std::pair< point, double > fmm::box_extent ( const std::vector< particle< point, value > > &  particles)
Here is the caller graph for this function:

◆ box_extent() [2/3]

template<typename value >
std::pair< point2d, double > fmm::box_extent ( const std::vector< particle< point2d, value > > &  particles)

Definition at line 219 of file morton_tree.cpp.

220{
221 tbb::enumerable_thread_specific<double> xmin(particles[0].center.real()), xmax(particles[0].center.real()),
222 ymin(particles[0].center.imag()), ymax(particles[0].center.imag());
223 tbb::parallel_for(tbb::blocked_range<size_t>(0, particles.size()), [&](const tbb::blocked_range<size_t>& r)
224 {
225 auto& Xmin = xmin.local();
226 auto& Xmax = xmax.local();
227 auto& Ymin = ymin.local();
228 auto& Ymax = ymax.local();
229 for (size_t i = r.begin(); i < r.end(); ++i)
230 {
231 double x = particles[i].center.real();
232 double y = particles[i].center.imag();
233 Xmin = std::min(Xmin, x);
234 Xmax = std::max(Xmax, x);
235 Ymin = std::min(Ymin, y);
236 Ymax = std::max(Ymax, y);
237 }
238 });
239 double Xmin = *std::min_element(xmin.begin(), xmin.end());
240 double Ymin = *std::min_element(ymin.begin(), ymin.end());
241 double Xmax = *std::max_element(xmax.begin(), xmax.end());
242 double Ymax = *std::max_element(ymax.begin(), ymax.end());
243
244 double dw = std::max(Xmax - Xmin, Ymax - Ymin);
245 return { {Xmin,Ymin}, dw };
246}
Here is the call graph for this function:

◆ box_extent() [3/3]

template<typename value >
std::pair< point3d, double > fmm::box_extent ( const std::vector< particle< point3d, value > > &  particles)

Definition at line 249 of file morton_tree.cpp.

250{
251 tbb::enumerable_thread_specific<double> xmin(particles[0].center[0]), xmax(particles[0].center[0]),
252 ymin(particles[0].center[1]), ymax(particles[0].center[1]),
253 zmin(particles[0].center[2]), zmax(particles[0].center[2]);
254 tbb::parallel_for(tbb::blocked_range<size_t>(0, particles.size()), [&](const tbb::blocked_range<size_t>& r)
255 {
256 auto& Xmin = xmin.local();
257 auto& Xmax = xmax.local();
258 auto& Ymin = ymin.local();
259 auto& Ymax = ymax.local();
260 auto& Zmin = zmin.local();
261 auto& Zmax = zmax.local();
262 for (size_t i = r.begin(); i < r.end(); ++i)
263 {
264 double x = particles[i].center[0];
265 double y = particles[i].center[1];
266 double z = particles[i].center[2];
267 Xmin = std::min(Xmin, x);
268 Xmax = std::max(Xmax, x);
269 Ymin = std::min(Ymin, y);
270 Ymax = std::max(Ymax, y);
271 Zmin = std::min(Zmin, z);
272 Zmax = std::max(Zmax, z);
273 }
274 });
275 double Xmin = *std::min_element(xmin.begin(), xmin.end());
276 double Ymin = *std::min_element(ymin.begin(), ymin.end());
277 double Zmin = *std::min_element(zmin.begin(), zmin.end());
278 double Xmax = *std::max_element(xmax.begin(), xmax.end());
279 double Ymax = *std::max_element(ymax.begin(), ymax.end());
280 double Zmax = *std::max_element(zmax.begin(), zmax.end());
281
282 double dw = std::max({Xmax - Xmin, Ymax - Ymin, Zmax - Zmin});
283 return { {Xmin, Ymin, Zmin}, dw };
284}
Here is the call graph for this function:

◆ compute_wignerd_coef()

std::vector< double > fmm::compute_wignerd_coef ( int  N)

Definition at line 194 of file special_functions.cpp.

195{
196 std::vector<double> g(N * (N + 1) / 2);
197 g[0] = 1;
198 for (int n = 1; n < N; ++n)
199 {
200 g[n * (n + 1) / 2] = sqrt((2 * n - 1.) / (2 * n)) * g[n * (n - 1) / 2];
201 for (int m = 1; m <= n; ++m)
202 {
203 g[n * (n + 1) / 2 + m] = sqrt((n - m + 1.) / (n + m)) * g[n * (n + 1) / 2 + m - 1];
204 }
205 }
206 return g;
207}

◆ ComputeAnm()

treevector< double > fmm::ComputeAnm ( int  N)

Definition at line 27 of file special_functions.cpp.

28{
29 treevector<double> Anm(2 * N);
30 Anm(0, 0) = 1.0;
31 for (int n = 1; n < 2 * N; ++n)
32 {
33 for (int m = 0; m <= n - 1; ++m)
34 Anm(n, m) = Anm(n - 1, m) * (n - m) * (n + m);
35 Anm(n, n) = Anm(n, n - 1) * 2 * n;
36 }
37 for (int n = 1; n < 2 * N; ++n)
38 {
39 for (int m = 0; m <= n; ++m)
40 Anm(n, m) = ni(n) / sqrt(Anm(n, m));
41 }
42 return Anm;
43}
fmm::detail::Anm
constexpr anm_wrapper< _3d_MAX_MULTIPOLE_NUM > Anm
Definition special_functions.h:116
Here is the call graph for this function:

◆ ComputeExact() [1/3]

template<InteractionType it, typename point_type , typename value_type >
void fmm::ComputeExact ( const std::vector< particle< point_type, value_type > > &  particles)

Definition at line 527 of file simple_math.h.

528{
529 ComputeExact<it, point_type, value_type>(particles, particles);
530}

◆ ComputeExact() [2/3]

template<InteractionType it, typename point_type , typename value_type >
void fmm::ComputeExact ( const std::vector< particle< point_type, value_type > > &  source_particles,
const std::vector< particle< point_type, value_type > > &  target_particles 
)

Definition at line 494 of file simple_math.h.

495{
496 if constexpr (it == fmm::InteractionType::Potential2d) {
497 auto foo = [](const particle2d& p1, const particle2d& p2) {return Potential2d(p1, p2); };
498 ComputeExact<point2d, double, value_type>(source_particles, target_particles, foo, "potential2d_exact.txt");
499 }
500 if constexpr (it == fmm::InteractionType::Force2d) {
501 auto foo = [](const particle2d& p1, const particle2d& p2) {return Force2d(p1, p2); };
502 ComputeExact<point2d, point2d, value_type>(source_particles, target_particles, foo, "force2d_exact.txt");
503 }
504 if constexpr (it == fmm::InteractionType::Potential3d)
505 {
506 auto foo = [](const particle3d& p1, const particle3d& p2) {return Potential3d(p1, p2); };
507 ComputeExact<point3d, double, value_type>(source_particles, target_particles, foo, "potential3d_exact.txt");
508 }
509 if constexpr (it == fmm::InteractionType::Force3d)
510 {
511 auto foo = [](const particle<point_type, value_type>& p1, const particle<point_type, value_type>& p2) {return Force3d(p1, p2); };
512 ComputeExact<point3d, Vector3d, value_type>(source_particles, target_particles, foo, "force3d_exact.txt");
513 }
514 if constexpr (it == fmm::InteractionType::VelDipole3d)
515 {
516 auto foo = [](const particle<point_type, value_type>& p1, const particle<point_type, value_type>& p2) {return velDipole3d(p1, p2); };
517 ComputeExact<point3d, Vector3d, value_type>(source_particles, target_particles, foo, "velDipole3d_exact.txt");
518 }
519 if constexpr (it == fmm::InteractionType::MomDipole3d)
520 {
521 auto foo = [](const particle<point_type, value_type>& p1, const particle<point_type, value_type>& p2) {return momDipole3d(p1, p2); };
522 ComputeExact<point3d, Vector3d, value_type>(source_particles, target_particles, foo, "momDipole3d_exact.txt");
523 }
524}
fmm::momDipole3d
Vector3d momDipole3d(const particle3d3 &p1, const particle3d3 &p2)
Definition simple_math.h:289
fmm::velDipole3d
Vector3d velDipole3d(const particle3d3 &p1, const particle3d3 &p2)
Definition simple_math.h:263
fmm::particle2d
particle< point2d, double > particle2d
Definition utils.h:155
fmm::particle3d
particle< point3d, double > particle3d
Definition utils.h:156
Here is the call graph for this function:

◆ ComputeExact() [3/3]

template<typename point_type , typename interaction_type , typename value_type >
void fmm::ComputeExact ( const std::vector< particle< point_type, value_type > > &  source_particles,
const std::vector< particle< point_type, value_type > > &  target_particles,
std::function< interaction_type(const particle< point_type, value_type > &, const particle< point_type, value_type > &)>  func,
std::string  filename 
)

Definition at line 398 of file simple_math.h.

401{
402 size_t num_particles = target_particles.size();
403#ifdef FMM_MPI
404 auto [shift, end_part] = LocalPart(0, num_particles);
405 size_t local_size = end_part - shift;
406#else
407 size_t shift = 0, local_size = num_particles;
408#endif
409 std::vector<interaction_type> exact(local_size);
410
411 auto particles_ptr = target_particles.data() + shift;
412 tbb::parallel_for(size_t(0), local_size, [&](size_t i)
413 {
414 const auto& p1 = particles_ptr[i];
415 for (const auto& p2 : source_particles)
416 {
417 exact[i] += func(p1, p2);
418 }
419 });
420
421#ifdef FMM_MPI
422 std::vector<int> sizes(NProc());
423 std::vector<int> displs(NProc());
424 std::vector<interaction_type> buf(num_particles);
425 for (int j = 0; j < NProc(); ++j)
426 {
427 sizes[j] = LocalPart(0, num_particles, j);
428 }
429 for (int j = 1; j < NProc(); ++j)
430 {
431 displs[j] = displs[j - 1] + sizes[j - 1];
432 }
433 if constexpr (std::is_same_v<interaction_type,double>)
434 MPI_Allgatherv(exact.data(), local_size, MPI_DOUBLE,
435 buf.data(), sizes.data(), displs.data(), MPI_DOUBLE, MPI_COMM_WORLD);
436 if constexpr (std::is_same_v<interaction_type,std::complex<double>>)
437 MPI_Allgatherv(exact.data(), local_size, MPI_COMPLEX16,
438 buf.data(), sizes.data(), displs.data(), MPI_COMPLEX16, MPI_COMM_WORLD);
439 if constexpr (std::is_same_v<interaction_type, Vector3d>)
440 {
441 MPI_Datatype MPI_VECTOR3;
442 MPI_Type_contiguous(3, MPI_DOUBLE, &MPI_VECTOR3);
443 MPI_Type_commit(&MPI_VECTOR3);
444 MPI_Allgatherv(exact.data(), local_size, MPI_VECTOR3,
445 buf.data(), sizes.data(), displs.data(), MPI_VECTOR3, MPI_COMM_WORLD);
446 }
447#endif
448
449 std::ofstream fout(filename);
450 fout.precision(12);
451#ifdef FMM_MPI
452 for (const auto& x : buf)
453 fout << x << "\n";
454#else
455 for (const auto& x : exact)
456 fout << x << "\n";
457#endif
458}

◆ ComputeExp()

std::vector< std::complex< double > > fmm::ComputeExp ( int  N,
double  phi 
)

Definition at line 353 of file special_functions.cpp.

354{
355 std::vector<std::complex<double>> res(N);
356 for (int i = 0; i < N; ++i)
357 res[i] = std::complex<double>(cos(i * phi), sin(i * phi));
358 return res;
359}
Here is the caller graph for this function:

◆ ComputeKnm()

treevector< double > fmm::ComputeKnm ( int  N)

Definition at line 9 of file special_functions.cpp.

10{
11 treevector<double> Knm(N);
12 Knm(0, 0) = 1.0;
13 for (int n = 1; n < N; ++n)
14 {
15 for (int m = 0; m <= n - 1; ++m)
16 Knm(n, m) = Knm(n - 1, m) * (n - m) / (n + m);
17 Knm(n, n) = Knm(n, n - 1) / (2 * n);
18 }
19 for (int n = 1; n < N; ++n)
20 {
21 for (int m = 0; m <= n; ++m)
22 Knm(n, m) = sqrt(Knm(n, m));
23 }
24 return Knm;
25}
fmm::detail::Knm
constexpr knm_wrapper< _3d_MAX_MULTIPOLE_NUM > Knm
Definition special_functions.h:85

◆ ComputePnm()

void fmm::ComputePnm ( int  N,
double  theta,
double *  Pnm 
)
inline

Definition at line 12 of file special_functions.h.

13{
14#ifdef __NVCC__
15 double x, y;
16 sincos(theta, &y, &x);
17#else
18 auto x = cos(theta);
19 auto y = sin(theta);
20#endif
21 Pnm[0] = 1;
22 for (int n = 1; n < N; ++n)
23 {
24 for (int m = 0; m <= n - 2; ++m)
25 {
26 Pnm[n * (n + 1) / 2 + m] = ((2 * n - 1) * x * Pnm[(n - 1) * n / 2 + m] - (n + m - 1) * Pnm[(n - 2) * (n - 1) / 2 + m]) / (n - m);
27 }
28 Pnm[n * (n + 1) / 2 + n - 1] = Pnm[(n - 1) * n / 2 + n - 1] * x * (2 * n - 1);
29 Pnm[n * (n + 1) / 2 + n] = Pnm[(n - 1) * n / 2 + n - 1] * y * (1 - 2 * n);
30 }
31}
Here is the caller graph for this function:

◆ ComputeWignerD()

std::vector< double > fmm::ComputeWignerD ( int  N,
double  theta 
)

Definition at line 211 of file special_functions.cpp.

212{
213 assert(std::abs(theta) <= std::numbers::pi);
214 // size can be computed as
215 // \sum_{n=0}^{N-1} \sum_{m=0}^n \sum_{k=-n}^n (1) = N * (1 + N) * (4 * N - 1) / 6
216 // linear index (n,m,k) can be found as
217 // \sum_{n'=0}^{n-1} \sum_{m=0}^{n'} \sum_{k=-n'}^{n'} 1 +
218 // + \sum_{m'=0}^{m-1} \sum_{k=-n}^n 1 +
219 // + \sum_{k'=-n}^{k-1} 1 =
220 // = 1/6 n (1 + n) (-1 + 4 n) + m (1 + 2 n) + (k + n)
221 std::vector<double> D(N * (1 + N) * (4 * N - 1) / 6);
222
223 // theta must be in [0;pi/2], other cases are considered using properties of d-matrix
224 // see properties here: https://en.wikipedia.org/wiki/Wigner_D-matrix
225 bool isnegative = theta < 0;
226 theta = std::abs(theta);
227
228 bool closetopi = false;
229 if (theta > std::numbers::pi / 2)
230 {
231 theta = std::numbers::pi - theta;
232 closetopi = true;
233 }
234
235 if (fabs(theta) < 1e-3) // theta = 0 => d^n_mm = 1, other zero
236 {
237 int idx0, idx1, idx2;
238 for (int n = 0; n < N; ++n)
239 {
240 idx0 = n * (1 + n) * (4 * n - 1) / 6 + n;
241 idx1 = 2 * n + 1;
242 for (int m = 0; m <= n; ++m)
243 {
244 idx2 = idx0 + idx1 * m;
245 D[idx2 + m] = 1.0;
246 }
247 }
248 }
249 else
250 {
251 const auto& g = wignerd_auxillary_coef;
252 int idx0, idx1, idx2, idx3;
253 double s = sin(theta);
254 double c = 1.0 + cos(theta);
255 double f = s / c;
256 std::vector<double> sn(2 * N), cn(2 * N);
257 sn[N] = cn[N] = 1.0;
258 for (int i = 1; i < N; ++i)
259 {
260 sn[N + i] = sn[N + i - 1] * s;
261 sn[N - i] = sn[N - i + 1] / s;
262 cn[N + i] = cn[N + i - 1] * c;
263 cn[N - i] = cn[N - i + 1] / c;
264 }
265
266
267 // 1) compute d^n_mn
268 for (int n = 0; n < N; ++n)
269 {
270 idx0 = n * (1 + n) * (4 * n - 1) / 6 + 2 * n;
271 idx1 = 1 + 2 * n;
272 for (int m = 0; m <= n; ++m)
273 {
274 D[idx0 + m * idx1] = ni(n+m) * g[n * (n + 1) / 2 + m] * cn[N + m] * sn[N + n - m];
275 }
276 }
277
278 // 2) compute d^n_n,m-1
279 for (int n = 0; n < N; ++n)
280 {
281 idx0 = n * (1 + n) * (4 * n - 1) / 6 + n * (2 * n + 1) + n;
282 for (int m = n; m > -n; --m)
283 {
284 D[idx0 + m - 1] = (n + m) / sqrt(n * (n + 1.) - m * (m - 1.)) * f * D[idx0 + m];
285 }
286 }
287
288 // 3) compute d^n_m,k-1
289 for (int n = 0; n < N; ++n)
290 {
291 idx0 = n * (1 + n) * (4 * n - 1) / 6 + n;
292 idx1 = 2 * n + 1;
293 for (int m = n - 1; m >= 0; --m)
294 {
295 idx2 = m * idx1;
296 idx3 = (m + 1) * idx1;
297 for (int k = n; k > -n; --k)
298 {
299 D[idx0 + idx2 + k - 1] = sqrt((n * (n + 1.) - m * (m + 1.)) / (n * (n + 1.) - k * (k - 1.))) * D[idx0 + idx3 + k] +
300 (m + k) / sqrt(n * (n + 1.) - k * (k - 1.)) * f * D[idx0 + idx2 + k];
301 }
302 }
303 }
304
305 }
306
307 int idx0, idx1, idx2;
308 auto Dcopy = D;
309 if (closetopi) // D^n_{m,k}(pi - theta) = (-1)^(n+m) * D^n_{m,-k}(theta)
310 {
311 for (int n = 1; n < N; ++n)
312 {
313 idx0 = n * (1 + n) * (4 * n - 1) / 6;
314 idx1 = 2 * n + 1;
315 for (int m = 0; m <= n; ++m)
316 {
317 idx2 = idx0 + idx1 * m + n;
318 for (int k = -n; k <= n; ++k)
319 {
320 D[idx2 + k] = ni(n + m) * Dcopy[idx2 - k];
321 }
322 }
323 }
324 }
325
326 for (int n = 1; n < N; ++n)
327 {
328 idx0 = n * (1 + n) * (4 * n - 1) / 6;
329 idx1 = 2 * n + 1;
330 for (int m = 0; m <= n; ++m)
331 {
332 idx2 = idx0 + idx1 * m + n;
333 for (int k = -n; k <= n; ++k)
334 {
335 // define eps(m) = (-1)^(m-|m|)/2, 1/eps(m) = eps(m),
336 // because Ynm in Laplace FMM is different from Ynm in physics (denote it Ynm_p)
337 // and they are connected as Ynm_p = eps(m) * sqrt((2n+1)/4pi) * Ynm,
338 // Wigner D-matrix are introduced for Ynm_p such that Ynm_p = sum_k Ynk_p*Dnkm,
339 // or equivalently eps(m)*Ynm = sum_k eps(k)*Ynk*Dnkm <=> Ynm = sum_k Ynk*[eps(m)*eps(k)*Dnkm]
340 // so Wigner D-matrix in fmm should have multiplier eps(m)*eps(k)
341 // since we consider only m>=0, we have only eps(k) left (eps(m) = 1 for m >= 0)
342 D[idx2 + k] *= ni((k - std::abs(k)) / 2);
343
344 if (isnegative) // D^n_{m,k}(-theta) = (-1)^(m - k) * D^n_{m,k}(theta)
345 D[idx2 + k] *= ni(m - k);
346 }
347 }
348 }
349
350 return D;
351}
fmm::ni
constexpr double ni(int n)
Definition simple_math.h:99
fmm::wignerd_auxillary_coef
const auto wignerd_auxillary_coef
Definition special_functions.cpp:209
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cross()

template<typename T >
constexpr Vector3< T > fmm::cross ( const Vector3< T > &  lhs,
const Vector3< T > &  rhs 
)
inlineconstexpr

Definition at line 108 of file utils.h.

108{ return Vector3<T>(lhs[1] * rhs[2] - lhs[2] * rhs[1], lhs[2] * rhs[0] - lhs[0] * rhs[2], lhs[0] * rhs[1] - lhs[1] * rhs[0]); }
Here is the caller graph for this function:

◆ DecToSph()

Vector3d fmm::DecToSph ( const Vector3d &  vec)
inline

Definition at line 104 of file simple_math.h.

105{
106 const double eps = 1e-12;
107 Vector3d c;
108 c[0] = abs(vec) + eps;
109 c[1] = acos(vec[2] / c[0]);
110 if (fabs(vec[0]) + fabs(vec[1]) < eps) {
111 c[2] = 0;
112 }
113 else if (fabs(vec[0]) < eps) {
114 c[2] = vec[1] / fabs(vec[1]) * /*std::numbers::pi*/ 3.14159265358979323846 * 0.5;
115 }
116 else {
117 c[2] = atan2(vec[1], vec[0]);
118 }
119 return c;
120}
fmm::abs
double abs(const Vector3< T > &vec)
Definition utils.h:122
Here is the call graph for this function:
Here is the caller graph for this function:

◆ divComp()

std::complex< double > fmm::divComp ( const std::complex< double > &  a,
const std::complex< double > &  b 
)
inline

Definition at line 58 of file multipole.cpp.

59{
60 const double& ax = a.real();
61 const double& bx = b.real();
62 const double& ay = a.imag();
63 const double& by = b.imag();
64
65 double zn = bx * bx + by * by;
66 return { (ax * bx + ay * by) / zn, (ay * bx - ax * by) / zn };
67}
Here is the caller graph for this function:

◆ dot()

constexpr double fmm::dot ( const Vector3d &  lhs,
const Vector3d &  rhs 
)
inlineconstexpr

Definition at line 110 of file utils.h.

110{ return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2]; }
Here is the caller graph for this function:

◆ doublehash() [1/2]

int fmm::doublehash ( double  val)
inline

Definition at line 124 of file utils.h.

125{
126 return int(val * 100);
127}
Here is the caller graph for this function:

◆ doublehash() [2/2]

int fmm::doublehash ( double  val1,
double  val2 
)
inline

Definition at line 129 of file utils.h.

130{
131 return 100 * int(val1 * 100) + int(val2 * 100);
132}

◆ ExpandBits()

template<int dim>
size_t fmm::ExpandBits ( double  x,
int  tree_depth 
)
inline

◆ ExpandBits< 2 >()

template<>
size_t fmm::ExpandBits< 2 > ( double  x,
int  tree_depth 
)
inline

Definition at line 49 of file morton_tree.cpp.

50{
51 assert(x > 0.0 && x < 1.0);
52#ifdef deterministic_behaviour
53 size_t bit_shift = (1ull << 30);
54#else
55 size_t bit_shift = (1ull << tree_depth);
56#endif
57 size_t val = (size_t)(bit_shift * x);
58 val = ((val << 16) | val) & 0xffff0000ffff;
59 val = ((val << 8) | val) & 0xff00ff00ff00ff;
60 val = ((val << 4) | val) & 0xf0f0f0f0f0f0f0f;
61 val = ((val << 2) | val) & 0x3333333333333333;
62 val = ((val << 1) | val) & 0x5555555555555555;
63 return val;
64}
Here is the caller graph for this function:

◆ ExpandBits< 3 >()

template<>
size_t fmm::ExpandBits< 3 > ( double  x,
int  tree_depth 
)
inline

Definition at line 67 of file morton_tree.cpp.

68{
69 assert(x > 0.0 && x < 1.0);
70#ifdef deterministic_behaviour
71 size_t bit_shift = (1ull << 20);
72#else
73 size_t bit_shift = (1ull << tree_depth);
74#endif
75 size_t val = (size_t)(bit_shift * x);
76 val = ((val << 32) | val) & 0xffff00000000ffff;
77 val = ((val << 16) | val) & 0xff0000ff0000ff;
78 val = ((val << 8) | val) & 0xf00f00f00f00f00f;
79 val = ((val << 4) | val) & 0x30c30c30c30c30c3;
80 val = ((val << 2) | val) & 0x9249249249249249;
81 return val;
82}
Here is the caller graph for this function:

◆ Force2d()

std::complex< double > fmm::Force2d ( const particle2d &  p1,
const particle2d &  p2 
)
inline

Definition at line 127 of file simple_math.h.

128{
129 auto dz = p1.center - p2.center;
130 return p2.q * std::conj(dz) / std::max(norm(dz), FORCE_EPS2);
131}
fmm::particle::q
value q
Definition utils.h:33
fmm::particle::center
point center
Definition utils.h:32
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Force2dMutual()

void fmm::Force2dMutual ( const particle2d &  p1,
const particle2d &  p2,
std::complex< double > &  force1,
std::complex< double > &  force2 
)
inline

Definition at line 355 of file simple_math.h.

356{
357 auto dz = p1.center - p2.center;
358 auto invdz = std::conj(dz) / std::max(norm(dz), FORCE_EPS2);
359 force1 += p2.q * invdz;
360 force2 -= p1.q * invdz;
361}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Force3d() [1/2]

Vector3d fmm::Force3d ( const particle3d &  p1,
const particle3d &  p2 
)
inline

Definition at line 138 of file simple_math.h.

139{
140 auto dr = p1.center - p2.center;
141 return p2.q * dr / MyPow(std::max(abs(dr), FORCE_EPS), 3);
142}
fmm::MyPow
constexpr T MyPow(T base, unsigned int exp)
Definition simple_math.h:85
Here is the call graph for this function:

◆ Force3d() [2/2]

Vector3d fmm::Force3d ( const particle3d3 &  p1,
const particle3d3 &  p2 
)
inline

Definition at line 144 of file simple_math.h.

145{
146 auto dr = p1.center - p2.center;
147 return cross(p2.q, dr) / MyPow(std::max(abs(dr), FORCE_EPS), 3);
148}
fmm::cross
constexpr Vector3< T > cross(const Vector3< T > &lhs, const Vector3< T > &rhs)
Definition utils.h:108
Here is the call graph for this function:

◆ Force3dMutual() [1/2]

void fmm::Force3dMutual ( const particle3d &  p1,
const particle3d &  p2,
Vector3d &  force1,
Vector3d &  force2 
)
inline

Definition at line 332 of file simple_math.h.

333{
334 auto dr = p1.center - p2.center;
335 auto invdr = dr / MyPow(std::max(abs(dr), FORCE_EPS), 3);
336 force1 += p2.q * invdr;
337 force2 -= p1.q * invdr;
338}
Here is the call graph for this function:

◆ Force3dMutual() [2/2]

void fmm::Force3dMutual ( const particle3d3 &  p1,
const particle3d3 &  p2,
Vector3d &  force1,
Vector3d &  force2 
)
inline

Definition at line 340 of file simple_math.h.

341{
342 auto dr = p1.center - p2.center;
343 auto invdr = dr / MyPow(std::max(abs(dr), FORCE_EPS), 3);
344 force1 += cross(p2.q, invdr);
345 force2 -= cross(p1.q, invdr);
346}
Here is the call graph for this function:

◆ IAmRoot()

bool fmm::IAmRoot ( )
inline

Definition at line 101 of file mpi_utils.h.

102{
103 return true;
104}
Here is the caller graph for this function:

◆ imag()

template<typename complex >
constexpr Vector3d fmm::imag ( const Vector3cd &  rhs)
inlineconstexpr

Definition at line 103 of file utils.h.

103{ return Vector3d(rhs[0].imag(), rhs[1].imag(), rhs[2].imag()); }
fmm::imag
constexpr Vector3d imag(const Vector3cd &rhs)
Definition utils.h:103
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lambdaEta()

template<typename T >
std::pair< T, T > fmm::lambdaEta ( T  xi)

Definition at line 211 of file simple_math.h.

212{
213 T sqrtPi = sqrt((T)3.141592653589793);
214 if (fabs(xi) < 1e-10)
215 return { 8.0 / (3.0 * sqrtPi), 8.0 / (5.0 * sqrtPi) };
216 T xi2 = xi * xi, xi3 = xi2 * xi, xi5 = xi3 * xi2;
217 T e = std::erf(xi), ex = std::exp(-xi2);
218 return {
219 -e/xi3 + (2/sqrtPi) *(2 + 1/xi2)*ex,
220 3 * e/xi5 - (2/sqrtPi/xi2)*(2 + 3/xi2)*ex
221 };
222};
Here is the caller graph for this function:

◆ lambdaEta< double >()

template<>
std::pair< double, double > fmm::lambdaEta< double > ( double  xi)
inline

Definition at line 211 of file simple_math.h.

227{
228 if (xi >= detail_interp::XI_MAX_INTERP) {
229 // far-field fallback (power-law, no erf/exp needed)
230 double xi3 = xi*xi*xi, xi5 = xi3*xi*xi;
231 return { -1.0/xi3, 3.0/xi5 };
232 }
233 const auto& t = detail_interp::tables();
234 return detail_interp::lookup2(t.lambda, t.eta, xi);
235}

◆ lambdaEtaA()

template<typename T >
std::pair< T, T > fmm::lambdaEtaA ( T  xi)

Definition at line 238 of file simple_math.h.

239{
240 T sqrtPi = sqrt((T)3.141592653589793);
241 if (fabs(xi) < 1e-10) return { 0.0, 0.0 };
242 T xi2 = xi*xi, xi3 = xi2*xi, xi4 = xi2*xi2, xi5 = xi3*xi2, xi6 = xi4*xi2;
243 T e = std::erf(xi), ex = std::exp(-xi2);
244 return {
245 3*e/xi4 - (8/sqrtPi)*(3/(4*xi3) + 1/(2*xi) + xi)*ex,
246 -15*e/xi6 + (1/sqrtPi)*(8/xi + 20/xi3 + 30/xi5)*ex
247 };
248};
Here is the caller graph for this function:

◆ lambdaEtaA< double >()

template<>
std::pair< double, double > fmm::lambdaEtaA< double > ( double  xi)
inline

Definition at line 238 of file simple_math.h.

253{
254 if (xi >= detail_interp::XI_MAX_INTERP) {
255 double xi4 = xi*xi*xi*xi, xi6 = xi4*xi*xi;
256 return { 3.0/xi4, -15.0/xi6 };
257 }
258 const auto& t = detail_interp::tables();
259 return detail_interp::lookup2(t.lambdaA, t.etaA, xi);
260}

◆ mm()

void fmm::mm ( const Matrix3d &  A,
const Matrix3d &  B,
Matrix3d &  C 
)
inline

Definition at line 424 of file multipole3d.cpp.

425{
426 C[0][0] = A[0][0] * B[0][0] + A[0][1] * B[1][0] + A[0][2] * B[2][0];
427 C[0][1] = A[0][0] * B[0][1] + A[0][1] * B[1][1] + A[0][2] * B[2][1];
428 C[0][2] = A[0][0] * B[0][2] + A[0][1] * B[1][2] + A[0][2] * B[2][2];
429
430 C[1][0] = A[1][0] * B[0][0] + A[1][1] * B[1][0] + A[1][2] * B[2][0];
431 C[1][1] = A[1][0] * B[0][1] + A[1][1] * B[1][1] + A[1][2] * B[2][1];
432 C[1][2] = A[1][0] * B[0][2] + A[1][1] * B[1][2] + A[1][2] * B[2][2];
433
434 C[2][0] = A[2][0] * B[0][0] + A[2][1] * B[1][0] + A[2][2] * B[2][0];
435 C[2][1] = A[2][0] * B[0][1] + A[2][1] * B[1][1] + A[2][2] * B[2][1];
436 C[2][2] = A[2][0] * B[0][2] + A[2][1] * B[1][2] + A[2][2] * B[2][2];
437}
Here is the caller graph for this function:

◆ mmt()

void fmm::mmt ( const Matrix3d &  A,
const Matrix3d &  B,
Matrix3d &  C 
)
inline

Definition at line 439 of file multipole3d.cpp.

440{
441 C[0][0] = A[0][0] * B[0][0] + A[0][1] * B[0][1] + A[0][2] * B[0][2];
442 C[0][1] = A[0][0] * B[1][0] + A[0][1] * B[1][1] + A[0][2] * B[1][2];
443 C[0][2] = A[0][0] * B[2][0] + A[0][1] * B[2][1] + A[0][2] * B[2][2];
444
445 C[1][0] = A[1][0] * B[0][0] + A[1][1] * B[0][1] + A[1][2] * B[0][2];
446 C[1][1] = A[1][0] * B[1][0] + A[1][1] * B[1][1] + A[1][2] * B[1][2];
447 C[1][2] = A[1][0] * B[2][0] + A[1][1] * B[2][1] + A[1][2] * B[2][2];
448
449 C[2][0] = A[2][0] * B[0][0] + A[2][1] * B[0][1] + A[2][2] * B[0][2];
450 C[2][1] = A[2][0] * B[1][0] + A[2][1] * B[1][1] + A[2][2] * B[1][2];
451 C[2][2] = A[2][0] * B[2][0] + A[2][1] * B[2][1] + A[2][2] * B[2][2];
452}
Here is the caller graph for this function:

◆ momDipole3d()

Vector3d fmm::momDipole3d ( const particle3d3 &  p1,
const particle3d3 &  p2 
)
inline

Definition at line 289 of file simple_math.h.

290{
291 //p1 - �����������, p2 - ��������
292 auto dr = (p1.center - p2.center);
293 auto Ldr = abs(dr);
294
295 auto argLE = (Ldr / FORCE_EPS);
296 if (argLE < 3.0)
297 {
298 auto [lambda, eta] = lambdaEta(argLE);
299 auto [lambdaA, etaA] = lambdaEtaA(argLE);
300
301 auto result = -1 * ((dot(p1.q, p2.q) * lambdaA / FORCE_EPS4 + dot(p2.q, dr) * dot(p1.q, dr) * etaA / FORCE_EPS6) / std::max(Ldr, 1e-10) * dr + \
302 (p2.q * dot(p1.q, dr) + p1.q * dot(p2.q, dr)) * eta / FORCE_EPS5);
303
304 //if (std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2]))
305 // std::cout << "AAA" << std::endl;
306 return result;
307 }
308 else
309 {
310 auto Ldr2 = Ldr * Ldr;
311 auto Ldr3 = Ldr2 * Ldr;
312 auto Ldr4 = Ldr2 * Ldr2;
313 auto Ldr5 = Ldr3 * Ldr2;
314 auto Ldr6 = Ldr3 * Ldr3;
315 auto result = -1 * ((dot(p1.q, p2.q) * (3 / Ldr4) + dot(p2.q, dr) * dot(p1.q, dr) * (-15 / Ldr6)) / std::max(Ldr, 1e-10) * dr + \
316 (p2.q * dot(p1.q, dr) + p1.q * dot(p2.q, dr)) * 3 / Ldr5);
317
318 //if (std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2]))
319 // std::cout << "BBB" << std::endl;
320 return result;
321 }
322}
fmm::FORCE_EPS6
constexpr double FORCE_EPS6
Definition defs.h:33
fmm::FORCE_EPS4
constexpr double FORCE_EPS4
Definition defs.h:31
fmm::lambdaEtaA
std::pair< T, T > lambdaEtaA(T xi)
Definition simple_math.h:238
fmm::dot
constexpr double dot(const Vector3d &lhs, const Vector3d &rhs)
Definition utils.h:110
fmm::lambdaEta
std::pair< T, T > lambdaEta(T xi)
Definition simple_math.h:211
fmm::FORCE_EPS
constexpr double FORCE_EPS
Definition defs.h:28
Here is the call graph for this function:
Here is the caller graph for this function:

◆ MortonCode()

template<int dim, typename point >
size_t fmm::MortonCode ( point  val,
int  tree_depth,
point  shift,
double  scale 
)
inline

◆ MortonCode< 2, point2d >()

template<>
size_t fmm::MortonCode< 2, point2d > ( point2d  val,
int  tree_depth,
point2d  shift,
double  scale 
)
inline

Definition at line 85 of file morton_tree.cpp.

86{
87 double x = (val.real() - shift.real()) / scale;
88 double y = (val.imag() - shift.imag()) / scale;
89 return (ExpandBits<2>(x, tree_depth) << 1) | ExpandBits<2>(y, tree_depth);
90}
Here is the call graph for this function:

◆ MortonCode< 3, point3d >()

template<>
size_t fmm::MortonCode< 3, point3d > ( point3d  val,
int  tree_depth,
point3d  shift,
double  scale 
)
inline

Definition at line 93 of file morton_tree.cpp.

94{
95 double x = (val[0] - shift[0]) / scale;
96 double y = (val[1] - shift[1]) / scale;
97 double z = (val[2] - shift[2]) / scale;
98 return (ExpandBits<3>(x, tree_depth) << 2) | (ExpandBits<3>(y, tree_depth) << 1) | ExpandBits<3>(z, tree_depth);
99}
Here is the call graph for this function:

◆ mv()

void fmm::mv ( const Matrix3d &  A,
const Vector3d &  x,
Vector3d &  y 
)
inline

Definition at line 417 of file multipole3d.cpp.

418{
419 y[0] = A[0][0] * x[0] + A[0][1] * x[1] + A[0][2] * x[2];
420 y[1] = A[1][0] * x[0] + A[1][1] * x[1] + A[1][2] * x[2];
421 y[2] = A[2][0] * x[0] + A[2][1] * x[1] + A[2][2] * x[2];
422}
Here is the caller graph for this function:

◆ MyPow()

template<typename T >
constexpr T fmm::MyPow ( T  base,
unsigned int  exp 
)
constexpr

Definition at line 85 of file simple_math.h.

86{
87 T res = 1;
88 while (exp) {
89 if (exp & 1)
90 {
91 res *= base;
92 }
93 exp >>= 1;
94 base *= base;
95 }
96 return res;
97}
Here is the caller graph for this function:

◆ ni()

constexpr double fmm::ni ( int  n)
inlineconstexpr

Definition at line 99 of file simple_math.h.

100{
101 return ((n & 1) == 1) ? -1.0 : 1.0;
102}
Here is the caller graph for this function:

◆ norm()

template<typename T >
double fmm::norm ( const Vector3< T > &  vec)
inline

Definition at line 119 of file utils.h.

119{return vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2];}
Here is the caller graph for this function:

◆ operator*() [1/5]

template<typename T >
constexpr Vector3< T > fmm::operator* ( const double  x,
const Vector3< T > &  rhs 
)
inlineconstexpr

Definition at line 69 of file utils.h.

69{return Vector3<T>(rhs[0]*x, rhs[1]*x, rhs[2]*x);}

◆ operator*() [2/5]

constexpr Vector3cd fmm::operator* ( const std::complex< double > &  x,
const Vector3cd &  rhs 
)
inlineconstexpr

Definition at line 73 of file utils.h.

73{ return Vector3cd(rhs[0] * x, rhs[1] * x, rhs[2] * x); }
fmm::Vector3cd
Vector3< std::complex< double > > Vector3cd
Definition utils.h:65

◆ operator*() [3/5]

constexpr Vector3cd fmm::operator* ( const std::complex< double > &  x,
const Vector3d &  rhs 
)
inlineconstexpr

Definition at line 71 of file utils.h.

71{ return Vector3cd(rhs[0] * x, rhs[1] * x, rhs[2] * x); }

◆ operator*() [4/5]

constexpr Vector3cd fmm::operator* ( const Vector3cd &  rhs,
const std::complex< double > &  x 
)
inlineconstexpr

Definition at line 74 of file utils.h.

74{ return Vector3cd(rhs[0] * x, rhs[1] * x, rhs[2] * x); }

◆ operator*() [5/5]

constexpr Vector3cd fmm::operator* ( const Vector3d &  rhs,
const std::complex< double > &  x 
)
inlineconstexpr

Definition at line 72 of file utils.h.

72{ return Vector3cd(rhs[0] * x, rhs[1] * x, rhs[2] * x); }

◆ operator*=()

void fmm::operator*= ( Vector3cd &  lhs,
const std::complex< double > &  x 
)
inline

Definition at line 76 of file utils.h.

76{ lhs[0] *= x; lhs[1] *= x; lhs[2] *= x; }

◆ operator+=()

void fmm::operator+= ( Vector3cd &  lhs,
const Vector3d &  rhs 
)
inline

Definition at line 85 of file utils.h.

85 {
86 lhs[0].real(lhs[0].real() + rhs[0]);
87 lhs[1].real(lhs[1].real() + rhs[1]);
88 lhs[2].real(lhs[2].real() + rhs[2]);
89}
fmm::real
constexpr Vector3d real(const Vector3cd &rhs)
Definition utils.h:98
Here is the call graph for this function:

◆ operator<<()

template<typename T >
std::ostream & fmm::operator<< ( std::ostream &  os,
const Vector3< T > &  rhs 
)
inline

Definition at line 113 of file utils.h.

113{ os << rhs[0] << " " << rhs[1] << " " << rhs[2]; return os; }

◆ operator>>()

template<typename T >
std::istream & fmm::operator>> ( std::istream &  is,
Vector3< T > &  rhs 
)
inline

Definition at line 116 of file utils.h.

116{is >> rhs[0] >> rhs[1] >> rhs[2]; return is;}

◆ PointFromMortonCode() [1/3]

template<typename point >
point fmm::PointFromMortonCode ( size_t  val,
int  depth,
point  shift,
double  scale 
)
Here is the caller graph for this function:

◆ PointFromMortonCode() [2/3]

template<>
point2d fmm::PointFromMortonCode ( size_t  val,
int  depth,
point2d  shift,
double  scale 
)

Definition at line 102 of file morton_tree.cpp.

103{
104 double step = 1.0 / (1 << (depth + 1));
105 double x = step, y = step;
106 step *= 2;
107 for (int i = 0; i < depth; ++i)
108 {
109 y += step * (val & 1);
110 val >>= 1;
111 x += step * (val & 1);
112 val >>= 1;
113 step *= 2;
114 }
115 x *= scale;
116 x += shift.real();
117 y *= scale;
118 y += shift.imag();
119 return { x, y };
120}

◆ PointFromMortonCode() [3/3]

template<>
point3d fmm::PointFromMortonCode ( size_t  val,
int  depth,
point3d  shift,
double  scale 
)

Definition at line 123 of file morton_tree.cpp.

124{
125 double step = 1.0 / (1 << (depth + 1));
126 double x = step, y = step, z = step;
127 step *= 2;
128 for (int i = 0; i < depth; ++i)
129 {
130 z += step * (val & 1);
131 val >>= 1;
132 y += step * (val & 1);
133 val >>= 1;
134 x += step * (val & 1);
135 val >>= 1;
136 step *= 2;
137 }
138 x *= scale;
139 x += shift[0];
140 y *= scale;
141 y += shift[1];
142 z *= scale;
143 z += shift[2];
144 return { x, y, z };
145}

◆ Potential2d()

double fmm::Potential2d ( const particle2d &  p1,
const particle2d &  p2 
)
inline

Definition at line 122 of file simple_math.h.

123{
124 return p2.q * 0.5 * log(std::max(norm(p1.center - p2.center), FORCE_EPS2));
125}
Here is the call graph for this function:

◆ Potential2dMutual()

void fmm::Potential2dMutual ( const particle2d &  p1,
const particle2d &  p2,
double &  potential1,
double &  potential2 
)
inline

Definition at line 348 of file simple_math.h.

349{
350 double dz = 0.5 * log(std::max(norm(p1.center - p2.center), FORCE_EPS2));
351 potential1 += p2.q * dz;
352 potential2 += p1.q * dz;
353}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Potential3d()

double fmm::Potential3d ( const particle3d &  p1,
const particle3d &  p2 
)
inline

Definition at line 133 of file simple_math.h.

134{
135 return p2.q / std::max(abs(p1.center - p2.center), FORCE_EPS);
136}
Here is the call graph for this function:

◆ Potential3dMutual()

void fmm::Potential3dMutual ( const particle3d &  p1,
const particle3d &  p2,
double &  potential1,
double &  potential2 
)
inline

Definition at line 325 of file simple_math.h.

326{
327 auto invdr = 1.0 / std::max(abs(p1.center - p2.center), FORCE_EPS);
328 potential1 += p2.q * invdr;
329 potential2 += p1.q * invdr;
330}
Here is the call graph for this function:

◆ ReadError()

template<InteractionType it, typename T >
void fmm::ReadError ( const std::vector< T > &  res)

Definition at line 533 of file simple_math.h.

534{
535 size_t num_particles = res.size();
536 std::vector<T> exact(num_particles);
537
538 std::string filename;
539 switch (it)
540 {
541 case fmm::InteractionType::Potential2d:
542 filename = "potential2d_exact.txt";
543 break;
544 case fmm::InteractionType::Force2d:
545 filename = "force2d_exact.txt";
546 break;
547 case fmm::InteractionType::Potential3d:
548 filename = "potential3d_exact.txt";
549 break;
550 case fmm::InteractionType::Force3d:
551 filename = "force3d_exact.txt";
552 break;
553 case fmm::InteractionType::VelDipole3d:
554 filename = "velDipole3d_exact.txt";
555 break;
556 default:
557 break;
558 }
559
560 std::ifstream fin(filename);
561 for (auto& x : exact)
562 fin >> x;
563
564 double max_error = 0.0;
565 double err1 = 0.0, err2 = 0.0;
566 double l2error = 0.0;
567//#pragma omp parallel for reduction(+: err1, err2, l2error) reduction(max: max_error)
568 for (int i = 0; i < num_particles; ++i)
569 {
570 double err = 0;
571 if constexpr ((it != InteractionType::Force3d) && (it != InteractionType::VelDipole3d))
572 err = std::abs(exact[i] - res[i]);
573 else
574 err = abs(exact[i] - res[i]);
575 max_error = std::max(max_error, err);
576
577 err1 += err;
578
579 if constexpr ((it != InteractionType::Force3d) && (it != InteractionType::VelDipole3d))
580 {
581 err2 += std::abs(exact[i]);
582 l2error += (err / std::abs(exact[i])) * (err / std::abs(exact[i]));
583 }
584 else
585 {
586 err2 += abs(exact[i]);
587 l2error += (err / abs(exact[i])) * (err / abs(exact[i]));
588 }
589
590 }
591 l2error /= num_particles;
592 //std::cout << "--------------------------------" << std::endl;
593 //std::cout << "max error = " << max_error << std::endl;
594 //std::cout << "l2 error = " << sqrt(l2error) << std::endl;
595 std::cout << "FMM relative error = " << err1 / err2 << std::endl;
596 //std::cout << "--------------------------------" << std::endl;
597}
Here is the call graph for this function:

◆ ReadErrorHyb()

void fmm::ReadErrorHyb ( const std::vector< std::array< double, 2 > > &  res)
inline

Definition at line 600 of file simple_math.h.

601{
602 size_t num_particles = res.size();
603 std::vector<std::complex<double>> exact(num_particles);
604
605 std::string filename = "force2d_exact.txt";
606
607 std::ifstream fin(filename);
608 for (auto& x : exact)
609 fin >> x;
610
611 double max_error = 0.0;
612 double err1 = 0.0, err2 = 0.0;
613 double l2error = 0.0;
614 //#pragma omp parallel for reduction(+: err1, err2, l2error) reduction(max: max_error)
615
616 for (int i = 0; i < num_particles; ++i)
617 {
618 double err = 0;
619
620 std::array<double, 2> hyb = { res[i][1] * 6.283185307179586476925286766559, res[i][0] * 6.283185307179586476925286766559 };
621
622 err = sqrt((exact[i].real() - hyb[0]) * (exact[i].real() - hyb[0]) +
623 (exact[i].imag() - hyb[1]) * (exact[i].imag() - hyb[1]));
624
625 max_error = std::max(max_error, err);
626
627 err1 += err;
628
629
630 err2 += std::abs(exact[i]);
631 l2error += (err / std::abs(exact[i])) * (err / std::abs(exact[i]));
632
633
634 }
635 l2error /= num_particles;
636 //std::cout << "--------------------------------" << std::endl;
637 //std::cout << "max error = " << max_error << std::endl;
638 //std::cout << "l2 error = " << sqrt(l2error) << std::endl;
639 std::cout << "Hybrid relative error = " << err1 / err2 << std::endl;
640 std::cout << "--------------------------------" << std::endl;
641}
Here is the call graph for this function:

◆ real()

template<typename complex >
constexpr Vector3d fmm::real ( const Vector3cd &  rhs)
inlineconstexpr

Definition at line 98 of file utils.h.

98{ return Vector3d(rhs[0].real(), rhs[1].real(), rhs[2].real()); }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ velDipole3d()

Vector3d fmm::velDipole3d ( const particle3d3 &  p1,
const particle3d3 &  p2 
)
inline

Definition at line 263 of file simple_math.h.

264{
265 auto dr = (p1.center - p2.center);
266 auto Ldr = abs(dr);
267
268 auto argLE = (Ldr / FORCE_EPS);
269 if (argLE < 3.0)
270 {
271 auto [lambda, eta] = lambdaEta(argLE);
272 auto result = p2.q * lambda / FORCE_EPS3 + dot(p2.q, dr) * dr * eta / FORCE_EPS5;
273 //if (std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2]))
274 // std::cout << "AAA" << std::endl;
275 return result;
276 }
277 else
278 {
279 auto Ldr2 = Ldr * Ldr;
280 auto Ldr3 = Ldr2 * Ldr;
281 auto Ldr5 = Ldr3 * Ldr2;
282 auto result = p2.q / (-Ldr3) + 3 * dot(p2.q, dr) * dr / Ldr5;
283 //if (std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2]))
284 // std::cout << "BBB" << std::endl;
285 return result;
286 }
287}
fmm::FORCE_EPS3
constexpr double FORCE_EPS3
Definition defs.h:30
fmm::FORCE_EPS5
constexpr double FORCE_EPS5
Definition defs.h:32
Here is the call graph for this function:
Here is the caller graph for this function:

◆ vt()

void fmm::vt ( const Vector3d &  x,
const Matrix3d *  T,
Matrix3d &  M 
)
inline

Definition at line 454 of file multipole3d.cpp.

455{
456 M[0][0] = x[0] * T[0][0][0] + x[1] * T[1][0][0] + x[2] * T[2][0][0];
457 M[0][1] = x[0] * T[0][0][1] + x[1] * T[1][0][1] + x[2] * T[2][0][1];
458 M[0][2] = x[0] * T[0][0][2] + x[1] * T[1][0][2] + x[2] * T[2][0][2];
459
460 M[1][0] = x[0] * T[0][1][0] + x[1] * T[1][1][0] + x[2] * T[2][1][0];
461 M[1][1] = x[0] * T[0][1][1] + x[1] * T[1][1][1] + x[2] * T[2][1][1];
462 M[1][2] = x[0] * T[0][1][2] + x[1] * T[1][1][2] + x[2] * T[2][1][2];
463
464 M[2][0] = x[0] * T[0][2][0] + x[1] * T[1][2][0] + x[2] * T[2][2][0];
465 M[2][1] = x[0] * T[0][2][1] + x[1] * T[1][2][1] + x[2] * T[2][2][1];
466 M[2][2] = x[0] * T[0][2][2] + x[1] * T[1][2][2] + x[2] * T[2][2][2];
467}
Here is the caller graph for this function:

Variable Documentation

◆ FMM_AUTO

const int fmm::FMM_AUTO = std::numeric_limits<int>::max()

Definition at line 26 of file defs.h.

◆ FORCE_EPS

constexpr double fmm::FORCE_EPS = 0.01
constexpr

Definition at line 28 of file defs.h.

◆ FORCE_EPS2

constexpr double fmm::FORCE_EPS2 = FORCE_EPS * FORCE_EPS
constexpr

Definition at line 29 of file defs.h.

◆ FORCE_EPS3

constexpr double fmm::FORCE_EPS3 = FORCE_EPS2 * FORCE_EPS
constexpr

Definition at line 30 of file defs.h.

◆ FORCE_EPS4

constexpr double fmm::FORCE_EPS4 = FORCE_EPS2 * FORCE_EPS2
constexpr

Definition at line 31 of file defs.h.

◆ FORCE_EPS5

constexpr double fmm::FORCE_EPS5 = FORCE_EPS3 * FORCE_EPS2
constexpr

Definition at line 32 of file defs.h.

◆ FORCE_EPS6

constexpr double fmm::FORCE_EPS6 = FORCE_EPS3 * FORCE_EPS3
constexpr

Definition at line 33 of file defs.h.

◆ wignerd_auxillary_coef

const auto fmm::wignerd_auxillary_coef = compute_wignerd_coef(100)

Definition at line 209 of file special_functions.cpp.