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
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
393 template<class Tvec>
394 class LatticeCubic {
395 public:
396 static constexpr u32 dim = 3;
397
398 static_assert(
399 dim == shambase::VectorProperties<Tvec>::dimension, "this lattice exist only in dim 3");
400
401 using Tscal = shambase::VecComponent<Tvec>;
402
410 static inline constexpr Tvec generator(Tscal dr, std::array<i32, dim> coord) noexcept {
411
412 i32 i = coord[0];
413 i32 j = coord[1];
414 i32 k = coord[2];
415
416 Tvec r_a = {i, j, k};
417
418 return dr * r_a;
419 }
420
421 static inline constexpr std::pair<std::array<i32, dim>, std::array<i32, dim>>
422 get_box_index_bounds(Tscal dr, Tvec box_min, Tvec box_max) {
423
424 Tvec coord_min;
425 Tvec coord_max;
426
427 coord_min[0] = box_min[0];
428 coord_max[0] = box_max[0];
429
430 coord_min[1] = box_min[1];
431 coord_max[1] = box_max[1];
432
433 coord_min[2] = box_min[2];
434 coord_max[2] = box_max[2];
435
436 coord_min /= dr;
437 coord_max /= dr;
438
439 std::array<i32, 3> ret_coord_min
440 = {i32(coord_min.x()) - 1, i32(coord_min.y()) - 1, i32(coord_min.z()) - 1};
441 std::array<i32, 3> ret_coord_max
442 = {i32(coord_max.x()) + 1, i32(coord_max.y()) + 1, i32(coord_max.z()) + 1};
443
444 return {ret_coord_min, ret_coord_max};
445 }
446
451 class IteratorDiscontinuous {
452 Tscal dr;
453
454 std::array<std::vector<size_t>, dim> remapped_indices;
455
456 std::array<size_t, dim> coord_delta;
457 size_t max_coord;
458
459 bool done = false;
460
461 public:
462 size_t current_idx;
463 IteratorDiscontinuous(
464 Tscal dr, std::array<i32, dim> coord_min, std::array<i32, dim> coord_max)
465 : dr(dr), current_idx(0), coord_delta({
466 size_t(coord_max[0] - coord_min[0]),
467 size_t(coord_max[1] - coord_min[1]),
468 size_t(coord_max[2] - coord_min[2]),
469 }) {
470
471 // must check for all axis otherwise we loop forever
472 for (int ax = 0; ax < dim; ax++) {
473 if (coord_min[ax] == coord_max[ax]) {
474 done = true;
475 }
476 }
477
478 max_coord = coord_delta[0] * coord_delta[1] * coord_delta[2];
479
480 for (int ax = 0; ax < dim; ax++) {
481 DiscontinuousIterator<i32> it(coord_min[ax], coord_max[ax]);
482 while (!it.is_done()) {
483 remapped_indices[ax].push_back(it.next());
484 }
485 }
486 }
487
488 inline bool is_done() { return done; }
489
490 inline Tvec next() {
491
492 // std::array<i32, 3> current = {
493 // coord_min[0] + i32(current_idx / (coord_delta[1] * coord_delta[2])),
494 // coord_min[1] + i32((current_idx / coord_delta[2]) % coord_delta[1]),
495 // coord_min[2] + i32(current_idx % coord_delta[2]),
496 // };
497
498 std::array<i32, 3> current = {
499 i32(current_idx % coord_delta[0]),
500 i32((current_idx / coord_delta[0]) % coord_delta[1]),
501 i32((current_idx / (coord_delta[0] * coord_delta[1]))),
502 };
503
504 current[0] = remapped_indices[0][current[0]];
505 current[1] = remapped_indices[1][current[1]];
506 current[2] = remapped_indices[2][current[2]];
507
508 // logger::raw_ln(current, current_idx, max_coord);
509
510 Tvec ret = generator(dr, current);
511
512 if (!done) {
513 current_idx++;
514 }
515 if (current_idx >= max_coord) {
516 done = true;
517 }
518
519 return ret;
520 }
521
522 inline std::vector<Tvec> next_n(u64 nmax) {
523 std::vector<Tvec> ret{};
524 for (u64 i = 0; i < nmax; i++) {
525 if (done) {
526 break;
527 }
528
529 ret.push_back(next());
530 }
531 shamlog_debug_ln("Discontinuous iterator", "next_n final idx", current_idx);
532 return ret;
533 }
534
535 inline void skip(u64 n) {
536 if (!done) {
537 current_idx += n;
538 }
539 if (current_idx >= max_coord) {
540 done = true;
541 }
542 shamlog_debug_ln("Discontinuous iterator", "skip final idx", current_idx);
543 }
544 };
545 };
546
547} // 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::IteratorDiscontinuous
Iterator utility to generate the lattice.
Definition crystalLattice.hpp:451
shammath::LatticeCubic
utility for generating Cubic crystal lattices
Definition crystalLattice.hpp:394
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:410
shammath::LatticeHCP::IteratorDiscontinuous
Iterator utility to generate the lattice.
Definition crystalLattice.hpp:292
shammath::LatticeHCP::Iterator
Iterator utility to generate the lattice.
Definition crystalLattice.hpp:205
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 Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8