33#ifndef MADNESS_MRA_DISPLACEMENTS_H__INCLUDED
34#define MADNESS_MRA_DISPLACEMENTS_H__INCLUDED
50 template <std::
size_t NDIM>
53 static std::vector< Key<NDIM> >
disp;
60 if (
NDIM == 1) bmax = 7;
61 else if (
NDIM == 2) bmax = 5;
62 else if (
NDIM == 3) bmax = 3;
63 else if (
NDIM == 4) bmax = 3;
64 else if (
NDIM == 5) bmax = 3;
65 else if (
NDIM == 6) bmax = 3;
72 return a.distsq() <
b.distsq();
84 for (std::size_t i=0; i<
NDIM; ++i) num *= (2*bmax + 1);
89 for (
d[0]=-bmax;
d[0]<=bmax; ++
d[0])
93 for (
d[0]=-bmax;
d[0]<=bmax; ++
d[0])
94 for (
d[1]=-bmax;
d[1]<=bmax; ++
d[1])
98 for (
d[0]=-bmax;
d[0]<=bmax; ++
d[0])
99 for (
d[1]=-bmax;
d[1]<=bmax; ++
d[1])
100 for (
d[2]=-bmax;
d[2]<=bmax; ++
d[2])
103 else if (
NDIM == 4) {
104 for (
d[0]=-bmax;
d[0]<=bmax; ++
d[0])
105 for (
d[1]=-bmax;
d[1]<=bmax; ++
d[1])
106 for (
d[2]=-bmax;
d[2]<=bmax; ++
d[2])
107 for (
d[3]=-bmax;
d[3]<=bmax; ++
d[3])
110 else if (
NDIM == 5) {
111 for (
d[0]=-bmax;
d[0]<=bmax; ++
d[0])
112 for (
d[1]=-bmax;
d[1]<=bmax; ++
d[1])
113 for (
d[2]=-bmax;
d[2]<=bmax; ++
d[2])
114 for (
d[3]=-bmax;
d[3]<=bmax; ++
d[3])
115 for (
d[4]=-bmax;
d[4]<=bmax; ++
d[4])
119 else if (
NDIM == 6) {
120 for (
d[0]=-bmax;
d[0]<=bmax; ++
d[0])
121 for (
d[1]=-bmax;
d[1]<=bmax; ++
d[1])
122 for (
d[2]=-bmax;
d[2]<=bmax; ++
d[2])
123 for (
d[3]=-bmax;
d[3]<=bmax; ++
d[3])
124 for (
d[4]=-bmax;
d[4]<=bmax; ++
d[4])
125 for (
d[5]=-bmax;
d[5]<=bmax; ++
d[5])
139 if (bmax > (twon-1)) bmax=twon-1;
148 if ((lx < 0) && (lx+twon > bmax)) bp[ip++] = lx + twon;
149 if ((lx > 0) && (lx-twon <-bmax)) bp[ip++] = lx - twon;
155 const int nbnp = inp;
162 for(
int i=0; i!=
NDIM; ++i) {
167 for (std::size_t i=0; i<
NDIM; ++i) {
198 if constexpr (
NDIM <= 3) {
216 if (kernel_lattice_sum_axes.
any()) {
219 if ((kernel_lattice_sum_axes &&
periodic_axes) != kernel_lattice_sum_axes) {
221 "Displacements<" + std::to_string(
NDIM) +
222 ">::get_disp(level, kernel_lattice_sum_axes): kernel_lattice_sum_axes is set for some axes that were not periodic in the FunctionDefault's boundary conditions active at the time when Displacements were initialized; invoke Displacements<NDIM>::reset_periodic_axes(kernel_lattice_sum_axes) to rebuild the periodic displacements";
258 for (
Level n = 0; n < nmax; ++n)
265 template <std::
size_t N, std::
size_t M>
266 constexpr std::enable_if_t<N>=M, std::array<std::size_t,
N-M>>
iota_array(std::array<std::size_t, M> values_to_skip_sorted) {
267 std::array<std::size_t,
N - M> result;
268 if constexpr (
N != M) {
269 std::size_t nadded = 0;
270 auto value_to_skip_it = values_to_skip_sorted.begin();
271 assert(*value_to_skip_it <
N);
272 auto value_to_skip = *value_to_skip_it;
273 for (std::size_t i = 0; i <
N; ++i) {
274 if (i < value_to_skip) {
275 result[nadded++] = i;
276 }
else if (value_to_skip_it != values_to_skip_sorted.end()) {
278 if (value_to_skip_it != values_to_skip_sorted.end()) {
279 value_to_skip = *value_to_skip_it;
293 template<std::
size_t NDIM>
305 using Box = std::array<std::pair<Translation, Translation>,
NDIM>;
331 mutable std::optional<Displacement>
disp;
369 auto increment_along_dim = [
this](
size_t dim) {
375 for (int64_t i =
NDIM - 1; i >= 0; --i) {
379 increment_along_dim(i);
387 const auto filtered_out = [&,
this]() {
389 const auto& filter = this->parent->
filter_;
393 std::optional<Displacement> nulldisp;
394 result = !filter(
point.level(), point_pattern, nulldisp);
430 for (
size_t i = 0; i <
NDIM; ++i) {
437 const auto filtered_out = [&]() ->
bool {
445 while (!done && filtered_out()) {
454 const auto is_fixed_dim = dim ==
fixed_dim;
469 l_dim_min = std::max(l_dim_min, l_dim_max - period + 1);
475 l_dim_min = std::max(
488 const auto filtered_out = [&,
this]() {
490 const auto& filter = this->parent->
filter_;
494 std::optional<Displacement> nulldisp;
495 result = !filter(
point.level(), point_pattern, nulldisp);
500 if (filtered_out()) {
501 bool have_another_surface_layer;
544 for (
size_t d = 0;
d !=
NDIM; ++
d) {
594 if (
a.done &&
b.done)
return true;
595 if (
a.done ||
b.done)
return false;
596 return a.fixed_dim ==
b.fixed_dim &&
633 bool has_finite_dimensions =
false;
634 const auto n =
center_.level();
639 r = (n == 0) ? (r+1)/2 : (r *
Translation(1) << (n-1));
642 if (!has_finite_dimensions)
643 probing_displacement_vec[
d] = r;
644 has_finite_dimensions =
true;
716 template <
size_t NDIM>
733 const std::array<KernelRange, NDIM>& range,
761 const auto twon = (
static_cast<Translation>(1) << level);
764 const auto map_to_range_twon = [&,
mask = ((~(
static_cast<std::uint64_t
>(0)) << (64-level)) >> (64-level))](std::int64_t x) -> std::int64_t {
765 const std::int64_t x_mapped = x &
mask;
770 const auto out_of_domain = [&](
const Translation& t) ->
bool {
771 return t < 0 || t >= twon;
775 const bool dest_is_in_domain = [&]() {
782 if (dest[
d].has_value() && out_of_domain(*dest[
d]))
return false;
787 if (dest_is_in_domain) {
795 bool among_standard_displacements =
true;
797 const auto disp_d = (*displacement)[
d];
802 auto disp_d_eff_abs =
std::abs(disp_d);
807 const std::int64_t disp_d_eff = map_to_range_twon(disp_d);
808 disp_d_eff_abs = std::min(disp_d_eff,
std::abs(disp_d_eff-twon));
812 if (dest[
d].has_value()) {
814 const auto dest_d_in_cell = map_to_range_twon(dest_d);
818 auto t = (*displacement).translation();
819 t[
d] += (dest_d_in_cell - dest_d);
827 if (disp_d_eff_abs > bmax_standard) {
828 among_standard_displacements =
false;
832 if (among_standard_displacements) {
Definition displacements.h:717
bool operator()(const Level level, const PointPattern &dest, std::optional< Displacement > &displacement) const
Apply filter to a displacement ending up at a point or a group of points (point pattern)
Definition displacements.h:755
array_of_bools< NDIM > domain_is_periodic_
Definition displacements.h:852
array_of_bools< NDIM > domain_is_infinite_
Definition displacements.h:851
std::array< KernelRange, NDIM > range_
Definition displacements.h:853
std::function< Translation(const Displacement &)> DistanceSquaredFunc
Definition displacements.h:723
Translation max_distsq_reached_
Definition displacements.h:855
Key< NDIM > Point
Definition displacements.h:719
DistanceSquaredFunc default_distance_squared_
Definition displacements.h:854
Vector< std::optional< Translation >, NDIM > PointPattern
Definition displacements.h:720
Key< NDIM > Displacement
Definition displacements.h:721
BoxSurfaceDisplacementFilter(const array_of_bools< NDIM > &domain_is_infinite, const array_of_bools< NDIM > &domain_is_periodic, const std::array< KernelRange, NDIM > &range, DistanceSquaredFunc default_distance_squared, Translation max_distsq_reached)
Definition displacements.h:730
Iterator class for lazy generation of surface points.
Definition displacements.h:325
std::optional< Displacement > disp
Memoized displacement from parent->center_ to point, computed by displacement(), reset by advance()
Definition displacements.h:331
Iterator(const BoxSurfaceDisplacementRange *p, Type type)
Constructs an iterator.
Definition displacements.h:536
Type
Definition displacements.h:327
@ End
Definition displacements.h:327
@ Begin
Definition displacements.h:327
size_t fixed_dim
Current fixed dimension (i.e. faces perpendicular to this axis are being iterated over)
Definition displacements.h:332
const std::optional< Displacement > & displacement() const
Definition displacements.h:515
friend bool operator!=(const Iterator &a, const Iterator &b)
Inequality comparison operator.
Definition displacements.h:606
bool done
Flag indicating iteration completion.
Definition displacements.h:334
const Point * pointer
Definition displacements.h:527
pointer operator->() const
Arrow operator for member access.
Definition displacements.h:566
Iterator operator++(int)
Post-increment operator.
Definition displacements.h:581
std::ptrdiff_t difference_type
Definition displacements.h:526
std::input_iterator_tag iterator_category
Definition displacements.h:524
Box box
updated box bounds in each dimension, used to avoid duplicate displacements by excluding the surface ...
Definition displacements.h:333
const BoxSurfaceDisplacementRange * parent
Pointer to parent box.
Definition displacements.h:329
bool next_surface_layer()
Definition displacements.h:337
reference operator*() const
Dereferences the iterator.
Definition displacements.h:560
Point value_type
Definition displacements.h:525
void reset_along_dim(size_t dim)
Definition displacements.h:453
void advance_till_valid()
Definition displacements.h:435
Point point
Current point.
Definition displacements.h:330
friend bool operator==(const Iterator &a, const Iterator &b)
Equality comparison operator.
Definition displacements.h:593
Iterator & operator++()
Pre-increment operator.
Definition displacements.h:572
const Point & reference
Definition displacements.h:528
void advance()
Advances the iterator to the next surface point.
Definition displacements.h:366
Definition displacements.h:294
Displacement probing_displacement_
displacement to a nearby point on the surface; it may not be able to pass the filter,...
Definition displacements.h:317
Key< NDIM > Point
Definition displacements.h:296
const Key< NDIM > & center() const
Definition displacements.h:684
std::array< std::optional< Translation >, NDIM > SurfaceThickness
Definition displacements.h:304
const std::array< std::optional< int64_t >, NDIM > & box_radius() const
Definition displacements.h:689
SurfaceThickness surface_thickness_
surface thickness in each dimension
Definition displacements.h:312
std::array< std::optional< Translation >, NDIM > BoxRadius
Definition displacements.h:303
const std::array< std::optional< int64_t >, NDIM > & surface_thickness() const
Definition displacements.h:694
BoxSurfaceDisplacementRange(const Key< NDIM > ¢er, const std::array< std::optional< std::int64_t >, NDIM > &box_radius, const std::array< std::optional< std::int64_t >, NDIM > &surface_thickness, const array_of_bools< NDIM > &is_periodic, Filter filter={})
Constructs a box with different sizes for each dimension.
Definition displacements.h:625
Hollowness hollowness_
does box contain non-surface points along each dimension?
Definition displacements.h:314
auto begin() const
Returns an iterator to the beginning of the surface points.
Definition displacements.h:663
Box box_
box bounds in each dimension
Definition displacements.h:313
Filter filter_
optional filter function
Definition displacements.h:316
const Displacement & probing_displacement() const
Definition displacements.h:704
Point center_
Center point of the box.
Definition displacements.h:309
Periodicity is_periodic_
which dimensions are periodic?
Definition displacements.h:315
std::function< bool(Level, const PointPattern &, std::optional< Displacement > &)> Filter
this callable filters out points and/or displacements; note that the displacement is optional (this u...
Definition displacements.h:300
Key< NDIM > Displacement
Definition displacements.h:298
std::array< bool, NDIM > Hollowness
Definition displacements.h:306
const array_of_bools< NDIM > & is_periodic() const
Definition displacements.h:699
auto end() const
Returns an iterator to the end of the surface points.
Definition displacements.h:669
Vector< std::optional< Translation >, NDIM > PointPattern
Definition displacements.h:297
std::array< std::pair< Translation, Translation >, NDIM > Box
Definition displacements.h:305
BoxRadius box_radius_
halved size of the box in each dimension
Definition displacements.h:310
Holds displacements for applying operators to avoid replicating for all operators.
Definition displacements.h:51
static std::vector< Key< NDIM > > disp
standard displacements to be used with standard kernels (range-unrestricted, no lattice sum)
Definition displacements.h:53
static bool cmp_keys(const Key< NDIM > &a, const Key< NDIM > &b)
Definition displacements.h:71
const std::vector< Key< NDIM > > & get_disp()
return the standard displacements appropriate for operators w/o lattice summation
Definition displacements.h:235
const std::vector< Key< NDIM > > & get_disp(Level n, const array_of_bools< NDIM > &kernel_lattice_sum_axes)
Definition displacements.h:211
static array_of_bools< NDIM > periodic_axes
along which axes lattice summation is performed?
Definition displacements.h:54
static void reset_periodic_axes(const array_of_bools< NDIM > &new_periodic_axes)
rebuilds periodic displacements so that they are optimal for the given set of periodic axes
Definition displacements.h:249
static void make_disp(int bmax)
Definition displacements.h:79
static std::vector< Key< NDIM > > disp_periodic[64]
displacements to be used with lattice-summed kernels
Definition displacements.h:55
static void make_disp_periodic(int bmax, Level n)
Definition displacements.h:135
static int bmax_default()
Definition displacements.h:58
static bool cmp_keys_periodic(const Key< NDIM > &a, const Key< NDIM > &b)
Definition displacements.h:75
Displacements()
Definition displacements.h:190
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:100
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:69
A simple, fixed dimension vector.
Definition vector.h:64
syntactic sugar for std::array<bool, N>
Definition array_of_bools.h:19
bool any() const
Definition array_of_bools.h:38
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
real_function_3d mask
Definition dirac-hatom.cc:27
Provides FunctionDefaults and utilities for coordinate transformation.
#define MADNESS_PRAGMA_CLANG(x)
Definition madness_config.h:200
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
int64_t Translation
Definition key.h:57
Key< NDIM > displacement(const Key< NDIM > &source, const Key< NDIM > &target)
given a source and a target, return the displacement in translation
Definition key.h:451
int Level
Definition key.h:58
static double pop(std::vector< double > &v)
Definition SCF.cc:113
constexpr std::array< std::size_t, N-M > iota_array(std::array< std::size_t, M > values_to_skip_sorted)
Definition displacements.h:266
std::string type(const PairType &n)
Definition PNOParameters.h:18
static long abs(long a)
Definition tensor.h:218
static const double b
Definition nonlinschro.cc:119
static const double d
Definition nonlinschro.cc:121
static const double a
Definition nonlinschro.cc:118
#define N
Definition testconv.cc:37
constexpr std::size_t NDIM
Definition testgconv.cc:54