Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
crystalLattice.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
18
19#include "shambase/aliases_int.hpp"
20#include "shambackends/sycl_utils.hpp"
21#include "shambackends/vec.hpp"
22#include "shamcomm/logs.hpp"
23#include "shammath/CoordRange.hpp"
24#include "shammath/DiscontinuousIterator.hpp"
25#include <array>
26#include <functional>
27#include <utility>
28#include <vector>
29
30namespace shammath {
31
32 class LatticeError : public std::exception {
33 public:
34 explicit LatticeError(const char *message) : msg_(message) {}
35
36 explicit LatticeError(const std::string &message) : msg_(message) {}
37
38 ~LatticeError() noexcept override = default;
39
40 [[nodiscard]] const char *what() const noexcept override { return msg_.c_str(); }
41
42 protected:
43 std::string msg_;
44 };
45
51 template<class Tvec>
52 class LatticeHCP {
53 public:
54 static constexpr u32 dim = 3;
55
56 static_assert(
57 dim == shambase::VectorProperties<Tvec>::dimension, "this lattice exist only in dim 3");
58
59 using Tscal = shambase::VecComponent<Tvec>;
60
68 static inline constexpr Tvec generator(Tscal dr, std::array<i32, dim> coord) noexcept {
69
70 i32 i = coord[0];
71 i32 j = coord[1];
72 i32 k = coord[2];
73
74 Tvec r_a
75 = {2 * i + (sycl::abs(j + k) % 2),
76 sycl::sqrt(3.) * (j + (1. / 3.) * (sycl::abs(k) % 2)),
77 2 * sycl::sqrt(6.) * k / 3};
78
79 return dr * r_a;
80 }
81
90 constexpr static bool can_make_periodic_box(
91 std::array<i32, dim> coord_min, std::array<i32, dim> coord_max) {
92 if (coord_max[0] - coord_min[0] < 2) {
93 return false;
94 }
95
96 if ((coord_max[1] + coord_min[1]) % 2 != 0) {
97 return false;
98 }
99
100 if ((coord_max[2] + coord_min[2]) % 2 != 0) {
101 return false;
102 }
103
104 return true;
105 }
106
116 static inline constexpr CoordRange<Tvec> get_periodic_box(
117 Tscal dr, std::array<i32, dim> coord_min, std::array<i32, dim> coord_max) {
118 Tscal xmin, xmax, ymin, ymax, zmin, zmax;
119
120 xmin = 2 * coord_min[0];
121 xmax = 2 * coord_max[0];
122
123 ymin = sycl::sqrt(3.) * coord_min[1];
124 ymax = sycl::sqrt(3.) * coord_max[1];
125
126 zmin = 2 * sycl::sqrt(6.) * coord_min[2] / 3;
127 zmax = 2 * sycl::sqrt(6.) * coord_max[2] / 3;
128
129 if (!can_make_periodic_box(coord_min, coord_max)) {
130 throw LatticeError(
131 "x axis count should be greater than 1\n"
132 "y axis count should be even\n"
133 "z axis count should be even");
134 }
135
136 return {Tvec{xmin, ymin, zmin} * dr, Tvec{xmax, ymax, zmax} * dr};
137 }
138
139 static inline constexpr std::pair<std::array<i32, dim>, std::array<i32, dim>>
140 get_box_index_bounds(Tscal dr, Tvec box_min, Tvec box_max) {
141
142 Tvec coord_min;
143 Tvec coord_max;
144
145 coord_min[0] = box_min[0] / 2.;
146 coord_max[0] = box_max[0] / 2.;
147
148 coord_min[1] = box_min[1] / sycl::sqrt(3.);
149 coord_max[1] = box_max[1] / sycl::sqrt(3.);
150
151 coord_min[2] = box_min[2] / (2 * sycl::sqrt(6.) / 3);
152 coord_max[2] = box_max[2] / (2 * sycl::sqrt(6.) / 3);
153
154 coord_min /= dr;
155 coord_max /= dr;
156
157 std::array<i32, 3> ret_coord_min
158 = {i32(coord_min.x()) - 1, i32(coord_min.y()) - 1, i32(coord_min.z()) - 1};
159 std::array<i32, 3> ret_coord_max
160 = {i32(coord_max.x()) + 1, i32(coord_max.y()) + 1, i32(coord_max.z()) + 1};
161
162 return {ret_coord_min, ret_coord_max};
163 }
164
172 static inline constexpr std::pair<std::array<i32, dim>, std::array<i32, dim>>
173 nearest_periodic_box_indices(
174 std::array<i32, dim> coord_min, std::array<i32, dim> coord_max) {
175 std::array<i32, dim> ret_coord_min;
176 std::array<i32, dim> ret_coord_max;
177
178 ret_coord_min[0] = coord_min[0];
179 ret_coord_min[1] = coord_min[1];
180 ret_coord_min[2] = coord_min[2];
181
182 ret_coord_max[0] = coord_max[0];
183 ret_coord_max[1] = coord_max[1];
184 ret_coord_max[2] = coord_max[2];
185
186 if (coord_max[0] - coord_min[0] < 2) {
187 ret_coord_max[0]++;
188 }
189
190 if ((coord_max[1] + coord_min[1]) % 2 != 0) {
191 ret_coord_max[1]++;
192 }
193
194 if ((coord_max[2] + coord_min[2]) % 2 != 0) {
195 ret_coord_max[2]++;
196 }
197
198 return {ret_coord_min, ret_coord_max};
199 }
200
205 class Iterator {
206 Tscal dr;
207 std::array<i32, dim> coord_min;
208
209 std::array<size_t, dim> coord_delta;
210 size_t current_idx;
211 size_t max_coord;
212
213 bool done = false;
214
215 public:
216 Iterator(Tscal dr, std::array<i32, dim> coord_min, std::array<i32, dim> coord_max)
217 : dr(dr), coord_min(coord_min), current_idx(0),
218 coord_delta({
219 size_t(coord_max[0] - coord_min[0]),
220 size_t(coord_max[1] - coord_min[1]),
221 size_t(coord_max[2] - coord_min[2]),
222 }) {
223
224 // must check for all axis otherwise we loop forever
225 for (int ax = 0; ax < dim; ax++) {
226 if (coord_min[ax] == coord_max[ax]) {
227 done = true;
228 }
229 }
230
231 max_coord = coord_delta[0] * coord_delta[1] * coord_delta[2];
232 }
233
234 inline bool is_done() { return done; }
235
236 inline Tvec next() {
237
238 // std::array<i32, 3> current = {
239 // coord_min[0] + i32(current_idx / (coord_delta[1] * coord_delta[2])),
240 // coord_min[1] + i32((current_idx / coord_delta[2]) % coord_delta[1]),
241 // coord_min[2] + i32(current_idx % coord_delta[2]),
242 // };
243
244 std::array<i32, 3> current = {
245 coord_min[0] + i32(current_idx % coord_delta[0]),
246 coord_min[1] + i32((current_idx / coord_delta[0]) % coord_delta[1]),
247 coord_min[2] + i32((current_idx / (coord_delta[0] * coord_delta[1]))),
248 };
249
250 // logger::raw_ln(current, current_idx, max_coord);
251
252 Tvec ret = generator(dr, current);
253
254 if (!done) {
255 current_idx++;
256 }
257 if (current_idx >= max_coord) {
258 done = true;
259 }
260
261 return ret;
262 }
263
264 inline std::vector<Tvec> next_n(u64 nmax) {
265 std::vector<Tvec> ret{};
266 for (u64 i = 0; i < nmax; i++) {
267 if (done) {
268 break;
269 }
270
271 ret.push_back(next());
272 }
273 shamlog_debug_ln("Discontinuous iterator", "next_n final idx", current_idx);
274 return ret;
275 }
276
277 inline void skip(u64 n) {
278 if (!done) {
279 current_idx += n;
280 }
281 if (current_idx >= max_coord) {
282 done = true;
283 }
284 shamlog_debug_ln("Discontinuous iterator", "skip final idx", current_idx);
285 }
286 };
287
292 class IteratorDiscontinuous {
293 Tscal dr;
294
295 std::array<std::vector<size_t>, dim> remapped_indices;
296
297 std::array<size_t, dim> coord_delta;
298 size_t max_coord;
299
300 bool done = false;
301
302 public:
303 size_t current_idx;
304 IteratorDiscontinuous(
305 Tscal dr, std::array<i32, dim> coord_min, std::array<i32, dim> coord_max)
306 : dr(dr), current_idx(0), coord_delta({
307 size_t(coord_max[0] - coord_min[0]),
308 size_t(coord_max[1] - coord_min[1]),
309 size_t(coord_max[2] - coord_min[2]),
310 }) {
311
312 // must check for all axis otherwise we loop forever
313 for (int ax = 0; ax < dim; ax++) {
314 if (coord_min[ax] == coord_max[ax]) {
315 done = true;
316 }
317 }
318
319 max_coord = coord_delta[0] * coord_delta[1] * coord_delta[2];
320
321 for (int ax = 0; ax < dim; ax++) {
322 DiscontinuousIterator<i32> it(coord_min[ax], coord_max[ax]);
323 while (!it.is_done()) {
324 remapped_indices[ax].push_back(it.next());
325 }
326 }
327 }
328
329 inline bool is_done() { return done; }
330
331 inline Tvec next() {
332
333 // std::array<i32, 3> current = {
334 // coord_min[0] + i32(current_idx / (coord_delta[1] * coord_delta[2])),
335 // coord_min[1] + i32((current_idx / coord_delta[2]) % coord_delta[1]),
336 // coord_min[2] + i32(current_idx % coord_delta[2]),
337 // };
338
339 std::array<i32, 3> current = {
340 i32(current_idx % coord_delta[0]),
341 i32((current_idx / coord_delta[0]) % coord_delta[1]),
342 i32((current_idx / (coord_delta[0] * coord_delta[1]))),
343 };
344
345 current[0] = remapped_indices[0][current[0]];
346 current[1] = remapped_indices[1][current[1]];
347 current[2] = remapped_indices[2][current[2]];
348
349 // logger::raw_ln(current, current_idx, max_coord);
350
351 Tvec ret = generator(dr, current);
352
353 if (!done) {
354 current_idx++;
355 }
356 if (current_idx >= max_coord) {
357 done = true;
358 }
359
360 return ret;
361 }
362
363 inline std::vector<Tvec> next_n(u64 nmax) {
364 std::vector<Tvec> ret{};
365 for (u64 i = 0; i < nmax; i++) {
366 if (done) {
367 break;
368 }
369
370 ret.push_back(next());
371 }
372 shamlog_debug_ln("Discontinuous iterator", "next_n final idx", current_idx);
373 return ret;
374 }
375
376 inline void skip(u64 n) {
377 if (!done) {
378 current_idx += n;
379 }
380 if (current_idx >= max_coord) {
381 done = true;
382 }
383 shamlog_debug_ln("Discontinuous iterator", "skip final idx", current_idx);
384 }
385 };
386
387 inline static std::tuple<Tvec, Tvec> get_ideal_hcp_box(
388 Tscal dr, std::tuple<Tvec, Tvec> box) {
389
390 auto [bmin, bmax] = box;
391
392 auto [idxs_min, idxs_max] = LatticeHCP::get_box_index_bounds(dr, bmin, bmax);
393
394 auto [idxs_min_per, idxs_max_per]
395 = LatticeHCP::nearest_periodic_box_indices(idxs_min, idxs_max);
396
397 shammath::CoordRange<Tvec> ret
398 = LatticeHCP::get_periodic_box(dr, idxs_min_per, idxs_max_per);
399
400 return {ret.lower, ret.upper};
401 }
402 };
403
409 template<class Tvec>
410 class LatticeCubic {
411 public:
412 static constexpr u32 dim = 3;
413
414 static_assert(
415 dim == shambase::VectorProperties<Tvec>::dimension, "this lattice exist only in dim 3");
416
417 using Tscal = shambase::VecComponent<Tvec>;
418
426 static inline constexpr Tvec generator(Tscal dr, std::array<i32, dim> coord) noexcept {
427
428 i32 i = coord[0];
429 i32 j = coord[1];
430 i32 k = coord[2];
431
432 Tvec r_a = {i, j, k};
433
434 return dr * r_a;
435 }
436
437 static inline constexpr std::pair<std::array<i32, dim>, std::array<i32, dim>>
438 get_box_index_bounds(Tscal dr, Tvec box_min, Tvec box_max) {
439
440 Tvec coord_min;
441 Tvec coord_max;
442
443 coord_min[0] = box_min[0];
444 coord_max[0] = box_max[0];
445
446 coord_min[1] = box_min[1];
447 coord_max[1] = box_max[1];
448
449 coord_min[2] = box_min[2];
450 coord_max[2] = box_max[2];
451
452 coord_min /= dr;
453 coord_max /= dr;
454
455 std::array<i32, 3> ret_coord_min
456 = {i32(coord_min.x()) - 1, i32(coord_min.y()) - 1, i32(coord_min.z()) - 1};
457 std::array<i32, 3> ret_coord_max
458 = {i32(coord_max.x()) + 1, i32(coord_max.y()) + 1, i32(coord_max.z()) + 1};
459
460 return {ret_coord_min, ret_coord_max};
461 }
462
467 class IteratorDiscontinuous {
468 Tscal dr;
469
470 std::array<std::vector<size_t>, dim> remapped_indices;
471
472 std::array<size_t, dim> coord_delta;
473 size_t max_coord;
474
475 bool done = false;
476
477 public:
478 size_t current_idx;
479 IteratorDiscontinuous(
480 Tscal dr, std::array<i32, dim> coord_min, std::array<i32, dim> coord_max)
481 : dr(dr), current_idx(0), coord_delta({
482 size_t(coord_max[0] - coord_min[0]),
483 size_t(coord_max[1] - coord_min[1]),
484 size_t(coord_max[2] - coord_min[2]),
485 }) {
486
487 // must check for all axis otherwise we loop forever
488 for (int ax = 0; ax < dim; ax++) {
489 if (coord_min[ax] == coord_max[ax]) {
490 done = true;
491 }
492 }
493
494 max_coord = coord_delta[0] * coord_delta[1] * coord_delta[2];
495
496 for (int ax = 0; ax < dim; ax++) {
497 DiscontinuousIterator<i32> it(coord_min[ax], coord_max[ax]);
498 while (!it.is_done()) {
499 remapped_indices[ax].push_back(it.next());
500 }
501 }
502 }
503
504 inline bool is_done() { return done; }
505
506 inline Tvec next() {
507
508 // std::array<i32, 3> current = {
509 // coord_min[0] + i32(current_idx / (coord_delta[1] * coord_delta[2])),
510 // coord_min[1] + i32((current_idx / coord_delta[2]) % coord_delta[1]),
511 // coord_min[2] + i32(current_idx % coord_delta[2]),
512 // };
513
514 std::array<i32, 3> current = {
515 i32(current_idx % coord_delta[0]),
516 i32((current_idx / coord_delta[0]) % coord_delta[1]),
517 i32((current_idx / (coord_delta[0] * coord_delta[1]))),
518 };
519
520 current[0] = remapped_indices[0][current[0]];
521 current[1] = remapped_indices[1][current[1]];
522 current[2] = remapped_indices[2][current[2]];
523
524 // logger::raw_ln(current, current_idx, max_coord);
525
526 Tvec ret = generator(dr, current);
527
528 if (!done) {
529 current_idx++;
530 }
531 if (current_idx >= max_coord) {
532 done = true;
533 }
534
535 return ret;
536 }
537
538 inline std::vector<Tvec> next_n(u64 nmax) {
539 std::vector<Tvec> ret{};
540 for (u64 i = 0; i < nmax; i++) {
541 if (done) {
542 break;
543 }
544
545 ret.push_back(next());
546 }
547 shamlog_debug_ln("Discontinuous iterator", "next_n final idx", current_idx);
548 return ret;
549 }
550
551 inline void skip(u64 n) {
552 if (!done) {
553 current_idx += n;
554 }
555 if (current_idx >= max_coord) {
556 done = true;
557 }
558 shamlog_debug_ln("Discontinuous iterator", "skip final idx", current_idx);
559 }
560 };
561 };
562
563} // namespace shammath
DiscontinuousIterator.hpp
header for PatchData related function and declaration
aliases_int.hpp
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
u64
std::uint64_t u64
64 bit unsigned integer
Definition aliases_int.hpp:26
i32
std::int32_t i32
32 bit integer
Definition aliases_int.hpp:23
shammath::DiscontinuousIterator
Discontinuous Iterator will iterate over every values in an integer set in the most discontinuous way...
Definition DiscontinuousIterator.hpp:45
shammath::LatticeCubic
utility for generating Cubic crystal lattices
Definition crystalLattice.hpp:410
shammath::LatticeCubic::generator
static constexpr Tvec generator(Tscal dr, std::array< i32, dim > coord) noexcept
generate a Cubic lattice centered on (0,0,0)
Definition crystalLattice.hpp:426
shammath::LatticeHCP
utility for generating HCP crystal lattices
Definition crystalLattice.hpp:52
shammath::LatticeHCP::can_make_periodic_box
static constexpr bool can_make_periodic_box(std::array< i32, dim > coord_min, std::array< i32, dim > coord_max)
check if the given lattice coordinates bounds can make a periodic box
Definition crystalLattice.hpp:90
shammath::LatticeHCP::generator
static constexpr Tvec generator(Tscal dr, std::array< i32, dim > coord) noexcept
generate a HCP lattice centered on (0,0,0)
Definition crystalLattice.hpp:68
shammath::LatticeHCP::nearest_periodic_box_indices
static constexpr std::pair< std::array< i32, dim >, std::array< i32, dim > > nearest_periodic_box_indices(std::array< i32, dim > coord_min, std::array< i32, dim > coord_max)
get the nearest integer triplets bound that gives a periodic box
Definition crystalLattice.hpp:173
shammath::LatticeHCP::get_periodic_box
static constexpr CoordRange< Tvec > get_periodic_box(Tscal dr, std::array< i32, dim > coord_min, std::array< i32, dim > coord_max)
Get the periodic box corresponding to integer lattice coordinates this function will throw if the coo...
Definition crystalLattice.hpp:116
shammath
namespace for math utility
Definition AABB.hpp:26
sycl_utils.hpp
Generated on for Shamrock by doxygen 1.15.0