Version: SMASH-3.1
Public Types | Public Member Functions | Protected Attributes | Private Member Functions | List of all members
smash::RectangularLattice< T > Class Template Reference

#include <lattice.h>

template<typename T>
class smash::RectangularLattice< T >

A container class to hold all the arrays on the lattice and access them.

Template Parameters
TThe type of the contained values.

Definition at line 47 of file lattice.h.

Public Types

using iterator = typename std::vector< T >::iterator
 Iterator of lattice. More...
 
using const_iterator = typename std::vector< T >::const_iterator
 Const interator of lattice. More...
 

Public Member Functions

 RectangularLattice (const std::array< double, 3 > &l, const std::array< int, 3 > &n, const std::array< double, 3 > &orig, bool per, const LatticeUpdate upd)
 Rectangular lattice constructor. More...
 
 RectangularLattice (RectangularLattice< T > const &rl)
 Copy-constructor. More...
 
void reset ()
 Sets all values on lattice to zeros. More...
 
bool out_of_bounds (int ix, int iy, int iz) const
 Checks if 3D index is out of lattice bounds. More...
 
ThreeVector cell_center (int ix, int iy, int iz) const
 Find the coordinates of a given cell. More...
 
ThreeVector cell_center (int index) const
 Find the coordinate of cell center given the 1d index of the cell. More...
 
const std::array< double, 3 > & lattice_sizes () const
 
const std::array< int, 3 > & n_cells () const
 
const std::array< double, 3 > & cell_sizes () const
 
const std::array< double, 3 > & origin () const
 
bool periodic () const
 
LatticeUpdate when_update () const
 
iterator begin ()
 
const_iterator begin () const
 
iterator end ()
 
const_iterator end () const
 
T & operator[] (std::size_t i)
 
const T & operator[] (std::size_t i) const
 
std::size_t size () const
 
void assign_value (int lattice_index, T value)
 Overwrite with a template value T at a given node. More...
 
int index1d (int ix, int iy, int iz)
 Return the index of a given cell. More...
 
int index_left (int ix, int iy, int iz)
 Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -x direction ("left"). More...
 
int index_right (int ix, int iy, int iz)
 Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +x direction ("right"). More...
 
int index_down (int ix, int iy, int iz)
 Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -y direction ("down"). More...
 
int index_up (int ix, int iy, int iz)
 Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +y direction ("up"). More...
 
int index_backward (int ix, int iy, int iz)
 Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -z direction ("backward"). More...
 
int index_forward (int ix, int iy, int iz)
 Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +z direction ("forward"). More...
 
void compute_gradient_lattice (RectangularLattice< ThreeVector > &grad_lat) const
 Compute a gradient on a lattice of doubles via the finite difference method. More...
 
void compute_four_gradient_lattice (RectangularLattice< FourVector > &old_lat, double time_step, RectangularLattice< std::array< FourVector, 4 >> &grad_lat) const
 Compute a fourgradient on a lattice of FourVectors jmu via the finite difference method. More...
 
T & node (int ix, int iy, int iz)
 Take the value of a cell given its 3-D indices. More...
 
bool value_at (const ThreeVector &r, T &value)
 Interpolates lattice quantity to coordinate r. More...
 
template<typename F >
void iterate_sublattice (const std::array< int, 3 > &lower_bounds, const std::array< int, 3 > &upper_bounds, F &&func)
 A sub-lattice iterator, which iterates in a 3D-structured manner and calls a function on every cell. More...
 
template<typename F >
void iterate_in_cube (const ThreeVector &point, const double r_cut, F &&func)
 Iterates only nodes whose cell centers lie not further than r_cut in x, y, z directions from the given point, that is iterates within a cube of side length 2*r_cut, and applies a function to each node. More...
 
template<typename F >
void integrate_volume (F &integral, F(*integrand)(ThreeVector, T &, ThreeVector), const double rcut, const ThreeVector &point)
 Calculate a volume integral with given integrand. More...
 
template<typename F >
void iterate_in_rectangle (const ThreeVector &point, const std::array< double, 3 > &rectangle, F &&func)
 Iterates only nodes whose cell centers lie not further than d_x in x-, d_y in y-, and d_z in z-direction from the given point, that is iterates within a rectangle of side lengths (2*dx, 2*dy, 2*dz), and applies a function to each node. More...
 
template<typename F >
void iterate_nearest_neighbors (const ThreeVector &point, F &&func)
 Iterates only over nodes corresponding to the center cell (the cell containing the given point) and its nearest neighbors in the -x, +x, -y, +y, -z, +z directions, and applies a function to each node. More...
 
template<typename L >
bool identical_to_lattice (const L *lat) const
 Checks if lattices of possibly different types have identical structure. More...
 

Protected Attributes

std::vector< T > lattice_
 The lattice itself, array containing physical quantities. More...
 
const std::array< double, 3 > lattice_sizes_
 Lattice sizes in x, y, z directions. More...
 
const std::array< int, 3 > n_cells_
 Number of cells in x,y,z directions. More...
 
const std::array< double, 3 > cell_sizes_
 Cell sizes in x, y, z directions. More...
 
const double cell_volume_
 Volume of a cell. More...
 
const std::array< double, 3 > origin_
 Coordinates of the left down nearer corner. More...
 
const bool periodic_
 Whether the lattice is periodic. More...
 
const LatticeUpdate when_update_
 When the lattice should be recalculated. More...
 

Private Member Functions

int positive_modulo (int i, int n) const
 Returns division modulo, which is always between 0 and n-1 in is not suitable, because it returns results from -(n-1) to n-1. More...
 

Member Typedef Documentation

◆ iterator

template<typename T >
using smash::RectangularLattice< T >::iterator = typename std::vector<T>::iterator

Iterator of lattice.

Definition at line 172 of file lattice.h.

◆ const_iterator

template<typename T >
using smash::RectangularLattice< T >::const_iterator = typename std::vector<T>::const_iterator

Const interator of lattice.

Definition at line 174 of file lattice.h.

Constructor & Destructor Documentation

◆ RectangularLattice() [1/2]

template<typename T >
smash::RectangularLattice< T >::RectangularLattice ( const std::array< double, 3 > &  l,
const std::array< int, 3 > &  n,
const std::array< double, 3 > &  orig,
bool  per,
const LatticeUpdate  upd 
)
inline

Rectangular lattice constructor.

Parameters
[in]l3-dimensional array (lx,ly,lz) indicates the size of the lattice in x, y, z directions respectively [fm].
[in]n3-dimensional array (nx,ny,nz) indicates the number of cells of the lattice in x, y, z directions respectively. Each cell in the lattice is labeled by three integers i, j, k where \(i\in[0, nx-1]\), \(j\in[0, ny-1]\), \(k\in[0, nz-1]\). The sizes of each cell are given by lx/nx, ly/ny, lz/nz in x,y,z directions respectively.
[in]origA 3-dimensional array indicating the coordinates of the origin [fm].
[in]perBoolean indicating whether a periodic boundary condition is applied.
[in]updEnum indicating how frequently the lattice is updated.

Definition at line 66 of file lattice.h.

70  : lattice_sizes_(l),
71  n_cells_(n),
72  cell_sizes_{l[0] / n[0], l[1] / n[1], l[2] / n[2]},
73  cell_volume_{cell_sizes_[0] * cell_sizes_[1] * cell_sizes_[2]},
74  origin_(orig),
75  periodic_(per),
76  when_update_(upd) {
77  lattice_.resize(n_cells_[0] * n_cells_[1] * n_cells_[2]);
78  logg[LLattice].debug(
79  "Rectangular lattice created: sizes[fm] = (", lattice_sizes_[0], ",",
80  lattice_sizes_[1], ",", lattice_sizes_[2], "), dims = (", n_cells_[0],
81  ",", n_cells_[1], ",", n_cells_[2], "), origin = (", origin_[0], ",",
82  origin_[1], ",", origin_[2], "), periodic: ", periodic_);
83  if (n_cells_[0] < 1 || n_cells_[1] < 1 || n_cells_[2] < 1 ||
84  lattice_sizes_[0] < 0.0 || lattice_sizes_[1] < 0.0 ||
85  lattice_sizes_[2] < 0.0) {
86  throw std::invalid_argument(
87  "Lattice sizes should be positive, "
88  "number of lattice cells should be > 0.");
89  }
90  }
smash::RectangularLattice::cell_volume_
const double cell_volume_
Volume of a cell.
Definition: lattice.h:796
smash::RectangularLattice::n_cells_
const std::array< int, 3 > n_cells_
Number of cells in x,y,z directions.
Definition: lattice.h:792
smash::RectangularLattice::periodic_
const bool periodic_
Whether the lattice is periodic.
Definition: lattice.h:800
smash::RectangularLattice::origin_
const std::array< double, 3 > origin_
Coordinates of the left down nearer corner.
Definition: lattice.h:798
smash::RectangularLattice::when_update_
const LatticeUpdate when_update_
When the lattice should be recalculated.
Definition: lattice.h:802
smash::RectangularLattice::lattice_sizes_
const std::array< double, 3 > lattice_sizes_
Lattice sizes in x, y, z directions.
Definition: lattice.h:790
smash::RectangularLattice::cell_sizes_
const std::array< double, 3 > cell_sizes_
Cell sizes in x, y, z directions.
Definition: lattice.h:794
smash::RectangularLattice::lattice_
std::vector< T > lattice_
The lattice itself, array containing physical quantities.
Definition: lattice.h:788
smash::logg
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
smash::pdg::n
constexpr int n
Neutron.
Definition: pdgcode_constants.h:30
smash::LLattice
static constexpr int LLattice
Definition: lattice.h:25

◆ RectangularLattice() [2/2]

template<typename T >
smash::RectangularLattice< T >::RectangularLattice ( RectangularLattice< T > const &  rl)
inline

Copy-constructor.

Definition at line 93 of file lattice.h.

94  : lattice_(rl.lattice_),
95  lattice_sizes_(rl.lattice_sizes_),
96  n_cells_(rl.n_cells_),
97  cell_sizes_(rl.cell_sizes_),
98  cell_volume_(rl.cell_volume_),
99  origin_(rl.origin_),
100  periodic_(rl.periodic_),
101  when_update_(rl.when_update_) {}

Member Function Documentation

◆ reset()

template<typename T >
void smash::RectangularLattice< T >::reset ( )
inline

Sets all values on lattice to zeros.

Definition at line 104 of file lattice.h.

104 { std::fill(lattice_.begin(), lattice_.end(), T()); }

◆ out_of_bounds()

template<typename T >
bool smash::RectangularLattice< T >::out_of_bounds ( int  ix,
int  iy,
int  iz 
) const
inline

Checks if 3D index is out of lattice bounds.

Parameters
[in]ixThe index of the cell in x direction.
[in]iyThe index of the cell in y direction.
[in]izThe index of the cell in z direction.
Returns
Whether the cell is out of the lattice.

Definition at line 114 of file lattice.h.

114  {
115  // clang-format off
116  return !periodic_ &&
117  (ix < 0 || ix >= n_cells_[0] ||
118  iy < 0 || iy >= n_cells_[1] ||
119  iz < 0 || iz >= n_cells_[2]);
120  // clang-format on
121  }

◆ cell_center() [1/2]

template<typename T >
ThreeVector smash::RectangularLattice< T >::cell_center ( int  ix,
int  iy,
int  iz 
) const
inline

Find the coordinates of a given cell.

Parameters
[in]ixThe index of the cell in x direction.
[in]iyThe index of the cell in y direction.
[in]izThe index of the cell in z direction.
Returns
Coordinates of the center of the given cell [fm].

Definition at line 131 of file lattice.h.

131  {
132  return ThreeVector(origin_[0] + cell_sizes_[0] * (ix + 0.5),
133  origin_[1] + cell_sizes_[1] * (iy + 0.5),
134  origin_[2] + cell_sizes_[2] * (iz + 0.5));
135  }

◆ cell_center() [2/2]

template<typename T >
ThreeVector smash::RectangularLattice< T >::cell_center ( int  index) const
inline

Find the coordinate of cell center given the 1d index of the cell.

Parameters
[in]index1-dimensional index of the given cell. It can be related to a 3-dimensional one by index = ix + nx (iy + iz * ny).
Returns
Coordinates of the center of the given cell [fm].

Definition at line 145 of file lattice.h.

145  {
146  const int ix = index % n_cells_[0];
147  index = index / n_cells_[0];
148  const int iy = index % n_cells_[1];
149  const int iz = index / n_cells_[1];
150  return cell_center(ix, iy, iz);
151  }
smash::RectangularLattice::cell_center
ThreeVector cell_center(int ix, int iy, int iz) const
Find the coordinates of a given cell.
Definition: lattice.h:131

◆ lattice_sizes()

template<typename T >
const std::array<double, 3>& smash::RectangularLattice< T >::lattice_sizes ( ) const
inline
Returns
Lengths of the lattice in x, y, z directions.

Definition at line 154 of file lattice.h.

154 { return lattice_sizes_; }

◆ n_cells()

template<typename T >
const std::array<int, 3>& smash::RectangularLattice< T >::n_cells ( ) const
inline
Returns
Number of cells in x, y, z directions.

Definition at line 157 of file lattice.h.

157 { return n_cells_; }

◆ cell_sizes()

template<typename T >
const std::array<double, 3>& smash::RectangularLattice< T >::cell_sizes ( ) const
inline
Returns
Lengths of one cell in x, y, z directions.

Definition at line 160 of file lattice.h.

160 { return cell_sizes_; }

◆ origin()

template<typename T >
const std::array<double, 3>& smash::RectangularLattice< T >::origin ( ) const
inline
Returns
Lattice origin: left, down, near corner coordinates.

Definition at line 163 of file lattice.h.

163 { return origin_; }

◆ periodic()

template<typename T >
bool smash::RectangularLattice< T >::periodic ( ) const
inline
Returns
If lattice is periodic or not.

Definition at line 166 of file lattice.h.

166 { return periodic_; }

◆ when_update()

template<typename T >
LatticeUpdate smash::RectangularLattice< T >::when_update ( ) const
inline
Returns
The enum, which tells at which time lattice needs to be updated.

Definition at line 169 of file lattice.h.

169 { return when_update_; }

◆ begin() [1/2]

template<typename T >
iterator smash::RectangularLattice< T >::begin ( )
inline
Returns
First element of lattice.

Definition at line 176 of file lattice.h.

176 { return lattice_.begin(); }

◆ begin() [2/2]

template<typename T >
const_iterator smash::RectangularLattice< T >::begin ( ) const
inline
Returns
First element of lattice (const).

Definition at line 178 of file lattice.h.

178 { return lattice_.begin(); }

◆ end() [1/2]

template<typename T >
iterator smash::RectangularLattice< T >::end ( )
inline
Returns
Last element of lattice.

Definition at line 180 of file lattice.h.

180 { return lattice_.end(); }

◆ end() [2/2]

template<typename T >
const_iterator smash::RectangularLattice< T >::end ( ) const
inline
Returns
Last element of lattice (const).

Definition at line 182 of file lattice.h.

182 { return lattice_.end(); }

◆ operator[]() [1/2]

template<typename T >
T& smash::RectangularLattice< T >::operator[] ( std::size_t  i)
inline
Returns
ith element of lattice.

Definition at line 184 of file lattice.h.

184 { return lattice_[i]; }

◆ operator[]() [2/2]

template<typename T >
const T& smash::RectangularLattice< T >::operator[] ( std::size_t  i) const
inline
Returns
ith element of lattice (const).

Definition at line 186 of file lattice.h.

186 { return lattice_[i]; }

◆ size()

template<typename T >
std::size_t smash::RectangularLattice< T >::size ( ) const
inline
Returns
Size of lattice.

Definition at line 188 of file lattice.h.

188 { return lattice_.size(); }

◆ assign_value()

template<typename T >
void smash::RectangularLattice< T >::assign_value ( int  lattice_index,
T  value 
)
inline

Overwrite with a template value T at a given node.

Definition at line 193 of file lattice.h.

193  {
194  lattice_[lattice_index] = value;
195  }

◆ index1d()

template<typename T >
int smash::RectangularLattice< T >::index1d ( int  ix,
int  iy,
int  iz 
)
inline

Return the index of a given cell.

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of cell

Definition at line 205 of file lattice.h.

205  {
206  assert(ix < n_cells_[0]);
207  assert(iy < n_cells_[1]);
208  assert(iz < n_cells_[2]);
209  return (ix + n_cells_[0] * (iy + iz * n_cells_[1]));
210  }

◆ index_left()

template<typename T >
int smash::RectangularLattice< T >::index_left ( int  ix,
int  iy,
int  iz 
)
inline

Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -x direction ("left").

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of the nearest cell in the "left" direction (ix-1) with respect to the current cell

Definition at line 221 of file lattice.h.

221  {
222  if (unlikely(ix == 0)) {
223  return periodic_ ? index1d(n_cells_[0] - 1, iy, iz) : index1d(ix, iy, iz);
224  } else {
225  return index1d(ix - 1, iy, iz);
226  }
227  }
smash::RectangularLattice::index1d
int index1d(int ix, int iy, int iz)
Return the index of a given cell.
Definition: lattice.h:205
unlikely
#define unlikely(x)
Tell the branch predictor that this expression is likely false.
Definition: macros.h:16

◆ index_right()

template<typename T >
int smash::RectangularLattice< T >::index_right ( int  ix,
int  iy,
int  iz 
)
inline

Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +x direction ("right").

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of the nearest cell in the "right" direction (ix+1) with respect to the current cell

Definition at line 238 of file lattice.h.

238  {
239  if (unlikely(ix == (n_cells_[0] - 1))) {
240  return periodic_ ? index1d(0, iy, iz) : index1d(ix, iy, iz);
241  } else {
242  return index1d(ix + 1, iy, iz);
243  }
244  }

◆ index_down()

template<typename T >
int smash::RectangularLattice< T >::index_down ( int  ix,
int  iy,
int  iz 
)
inline

Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -y direction ("down").

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of the nearest cell in the "down" direction (iy-1) with respect to the current cell

Definition at line 255 of file lattice.h.

255  {
256  if (unlikely(iy == 0)) {
257  return periodic_ ? index1d(ix, n_cells_[1] - 1, iz) : index1d(ix, iy, iz);
258  } else {
259  return index1d(ix, iy - 1, iz);
260  }
261  }

◆ index_up()

template<typename T >
int smash::RectangularLattice< T >::index_up ( int  ix,
int  iy,
int  iz 
)
inline

Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +y direction ("up").

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of the nearest cell in the "up" direction (iy+1) with respect to the current cell

Definition at line 272 of file lattice.h.

272  {
273  if (unlikely(iy == (n_cells_[1] - 1))) {
274  return periodic_ ? index1d(ix, 0, iz) : index1d(ix, iy, iz);
275  } else {
276  return index1d(ix, iy + 1, iz);
277  }
278  }

◆ index_backward()

template<typename T >
int smash::RectangularLattice< T >::index_backward ( int  ix,
int  iy,
int  iz 
)
inline

Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -z direction ("backward").

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of the nearest cell in the "backward" direction (iz-1) with respect to the current cell

Definition at line 289 of file lattice.h.

289  {
290  if (unlikely(iz == 0)) {
291  return periodic_ ? index1d(ix, iy, n_cells_[2] - 1) : index1d(ix, iy, iz);
292  } else {
293  return index1d(ix, iy, iz - 1);
294  }
295  }

◆ index_forward()

template<typename T >
int smash::RectangularLattice< T >::index_forward ( int  ix,
int  iy,
int  iz 
)
inline

Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +z direction ("forward").

Parameters
[in]ixindex of a cell in the x-direction
[in]iyindex of a cell in the y-direction
[in]izindex of a cell in the z-direction
Returns
index of the nearest cell in the "forward" direction (iz+1) with respect to the current cell

Definition at line 306 of file lattice.h.

306  {
307  if (unlikely(iz == (n_cells_[2] - 1))) {
308  return periodic_ ? index1d(ix, iy, 0) : index1d(ix, iy, iz);
309  } else {
310  return index1d(ix, iy, iz + 1);
311  }
312  }

◆ compute_gradient_lattice()

template<typename T >
void smash::RectangularLattice< T >::compute_gradient_lattice ( RectangularLattice< ThreeVector > &  grad_lat) const
inline

Compute a gradient on a lattice of doubles via the finite difference method.

return a lattice of ThreeVectors which are gradients of the values on the original lattice

Definition at line 320 of file lattice.h.

321  {
322  if (n_cells_[0] < 2 || n_cells_[1] < 2 || n_cells_[2] < 2) {
323  // Gradient calculation is impossible
324  throw std::runtime_error(
325  "Lattice is too small for gradient calculation"
326  " (should be at least 2x2x2)");
327  }
328  if (!identical_to_lattice(&grad_lat)) {
329  // Lattice for gradient should have identical origin/dims/periodicity
330  throw std::invalid_argument(
331  "Lattice for gradient should have the"
332  " same origin/dims/periodicity as the original one.");
333  }
334  const double inv_2dx = 0.5 / cell_sizes_[0];
335  const double inv_2dy = 0.5 / cell_sizes_[1];
336  const double inv_2dz = 0.5 / cell_sizes_[2];
337  const int dix = 1;
338  const int diy = n_cells_[0];
339  const int diz = n_cells_[0] * n_cells_[1];
340  const int d = diz * n_cells_[2];
341 
342  for (int iz = 0; iz < n_cells_[2]; iz++) {
343  const int z_offset = diz * iz;
344  for (int iy = 0; iy < n_cells_[1]; iy++) {
345  const int y_offset = diy * iy + z_offset;
346  for (int ix = 0; ix < n_cells_[0]; ix++) {
347  const int index = ix + y_offset;
348  if (unlikely(ix == 0)) {
349  (grad_lat)[index].set_x1(
350  periodic_
351  ? (lattice_[index + dix] - lattice_[index + diy - dix]) *
352  inv_2dx
353  : (lattice_[index + dix] - lattice_[index]) * 2.0 *
354  inv_2dx);
355  } else if (unlikely(ix == n_cells_[0] - 1)) {
356  (grad_lat)[index].set_x1(
357  periodic_
358  ? (lattice_[index - diy + dix] - lattice_[index - dix]) *
359  inv_2dx
360  : (lattice_[index] - lattice_[index - dix]) * 2.0 *
361  inv_2dx);
362  } else {
363  (grad_lat)[index].set_x1(
364  (lattice_[index + dix] - lattice_[index - dix]) * inv_2dx);
365  }
366 
367  if (unlikely(iy == 0)) {
368  (grad_lat)[index].set_x2(
369  periodic_
370  ? (lattice_[index + diy] - lattice_[index + diz - diy]) *
371  inv_2dy
372  : (lattice_[index + diy] - lattice_[index]) * 2.0 *
373  inv_2dy);
374  } else if (unlikely(iy == n_cells_[1] - 1)) {
375  (grad_lat)[index].set_x2(
376  periodic_
377  ? (lattice_[index - diz + diy] - lattice_[index - diy]) *
378  inv_2dy
379  : (lattice_[index] - lattice_[index - diy]) * 2.0 *
380  inv_2dy);
381  } else {
382  (grad_lat)[index].set_x2(
383  (lattice_[index + diy] - lattice_[index - diy]) * inv_2dy);
384  }
385 
386  if (unlikely(iz == 0)) {
387  (grad_lat)[index].set_x3(
388  periodic_
389  ? (lattice_[index + diz] - lattice_[index + d - diz]) *
390  inv_2dz
391  : (lattice_[index + diz] - lattice_[index]) * 2.0 *
392  inv_2dz);
393  } else if (unlikely(iz == n_cells_[2] - 1)) {
394  (grad_lat)[index].set_x3(
395  periodic_
396  ? (lattice_[index - d + diz] - lattice_[index - diz]) *
397  inv_2dz
398  : (lattice_[index] - lattice_[index - diz]) * 2.0 *
399  inv_2dz);
400  } else {
401  (grad_lat)[index].set_x3(
402  (lattice_[index + diz] - lattice_[index - diz]) * inv_2dz);
403  }
404  }
405  }
406  }
407  }
smash::RectangularLattice::identical_to_lattice
bool identical_to_lattice(const L *lat) const
Checks if lattices of possibly different types have identical structure.
Definition: lattice.h:770

◆ compute_four_gradient_lattice()

template<typename T >
void smash::RectangularLattice< T >::compute_four_gradient_lattice ( RectangularLattice< FourVector > &  old_lat,
double  time_step,
RectangularLattice< std::array< FourVector, 4 >> &  grad_lat 
) const
inline

Compute a fourgradient on a lattice of FourVectors jmu via the finite difference method.

Parameters
[in]old_latthe lattice of FourVectors jmu at a previous time step
[in]time_stepthe used time step, needed for the time derivative
[out]grad_lata lattice of 4-arrays of 4-vectors with the following structure: [djmu_dt, djmu_dx, djmu_dy, djmu_dz]

Definition at line 418 of file lattice.h.

420  {
421  if (n_cells_[0] < 2 || n_cells_[1] < 2 || n_cells_[2] < 2) {
422  // Gradient calculation is impossible
423  throw std::runtime_error(
424  "Lattice is too small for gradient calculation"
425  " (should be at least 2x2x2)");
426  }
427  if (!identical_to_lattice(&grad_lat)) {
428  // Lattice for gradient should have identical origin/dims/periodicity
429  throw std::invalid_argument(
430  "Lattice for gradient should have the"
431  " same origin/dims/periodicity as the original one.");
432  }
433  const double inv_2dx = 0.5 / cell_sizes_[0];
434  const double inv_2dy = 0.5 / cell_sizes_[1];
435  const double inv_2dz = 0.5 / cell_sizes_[2];
436  const int dix = 1;
437  const int diy = n_cells_[0];
438  const int diz = n_cells_[0] * n_cells_[1];
439  const int d = diz * n_cells_[2];
440 
441  for (int iz = 0; iz < n_cells_[2]; iz++) {
442  const int z_offset = diz * iz;
443  for (int iy = 0; iy < n_cells_[1]; iy++) {
444  const int y_offset = diy * iy + z_offset;
445  for (int ix = 0; ix < n_cells_[0]; ix++) {
446  const int index = ix + y_offset;
447 
448  // auxiliary vectors used for the calculation of gradients
449  FourVector grad_t_jmu(0.0, 0.0, 0.0, 0.0);
450  FourVector grad_x_jmu(0.0, 0.0, 0.0, 0.0);
451  FourVector grad_y_jmu(0.0, 0.0, 0.0, 0.0);
452  FourVector grad_z_jmu(0.0, 0.0, 0.0, 0.0);
453  // t direction
454  grad_t_jmu = (lattice_[index] - (old_lat)[index]) * (1.0 / time_step);
455  // x direction
456  if (unlikely(ix == 0)) {
457  grad_x_jmu =
458  periodic_
459  ? (lattice_[index + dix] - lattice_[index + diy - dix]) *
460  inv_2dx
461  : (lattice_[index + dix] - lattice_[index]) * 2.0 * inv_2dx;
462  } else if (unlikely(ix == n_cells_[0] - 1)) {
463  grad_x_jmu =
464  periodic_
465  ? (lattice_[index - diy + dix] - lattice_[index - dix]) *
466  inv_2dx
467  : (lattice_[index] - lattice_[index - dix]) * 2.0 * inv_2dx;
468  } else {
469  grad_x_jmu =
470  (lattice_[index + dix] - lattice_[index - dix]) * inv_2dx;
471  }
472  // y direction
473  if (unlikely(iy == 0)) {
474  grad_y_jmu =
475  periodic_
476  ? (lattice_[index + diy] - lattice_[index + diz - diy]) *
477  inv_2dy
478  : (lattice_[index + diy] - lattice_[index]) * 2.0 * inv_2dy;
479  } else if (unlikely(iy == n_cells_[1] - 1)) {
480  grad_y_jmu =
481  periodic_
482  ? (lattice_[index - diz + diy] - lattice_[index - diy]) *
483  inv_2dy
484  : (lattice_[index] - lattice_[index - diy]) * 2.0 * inv_2dy;
485  } else {
486  grad_y_jmu =
487  (lattice_[index + diy] - lattice_[index - diy]) * inv_2dy;
488  }
489  // z direction
490  if (unlikely(iz == 0)) {
491  grad_z_jmu =
492  periodic_
493  ? (lattice_[index + diz] - lattice_[index + d - diz]) *
494  inv_2dz
495  : (lattice_[index + diz] - lattice_[index]) * 2.0 * inv_2dz;
496  } else if (unlikely(iz == n_cells_[2] - 1)) {
497  grad_z_jmu =
498  periodic_
499  ? (lattice_[index - d + diz] - lattice_[index - diz]) *
500  inv_2dz
501  : (lattice_[index] - lattice_[index - diz]) * 2.0 * inv_2dz;
502  } else {
503  grad_z_jmu =
504  (lattice_[index + diz] - lattice_[index - diz]) * inv_2dz;
505  }
506  // fill
507  (grad_lat)[index][0] = grad_t_jmu;
508  (grad_lat)[index][1] = grad_x_jmu;
509  (grad_lat)[index][2] = grad_y_jmu;
510  (grad_lat)[index][3] = grad_z_jmu;
511  }
512  }
513  }
514  }

◆ node()

template<typename T >
T& smash::RectangularLattice< T >::node ( int  ix,
int  iy,
int  iz 
)
inline

Take the value of a cell given its 3-D indices.

Parameters
[in]ixThe index of the cell in x direction.
[in]iyThe index of the cell in y direction.
[in]izThe index of the cell in z direction.
Returns
Physical quantity evaluated at the cell center.

Definition at line 524 of file lattice.h.

524  {
525  return periodic_
526  ? lattice_[positive_modulo(ix, n_cells_[0]) +
527  n_cells_[0] *
528  (positive_modulo(iy, n_cells_[1]) +
529  n_cells_[1] * positive_modulo(iz, n_cells_[2]))]
530  : lattice_[ix + n_cells_[0] * (iy + n_cells_[1] * iz)];
531  }
smash::RectangularLattice::positive_modulo
int positive_modulo(int i, int n) const
Returns division modulo, which is always between 0 and n-1 in is not suitable, because it returns res...
Definition: lattice.h:813

◆ value_at()

template<typename T >
bool smash::RectangularLattice< T >::value_at ( const ThreeVector &  r,
T &  value 
)
inline

Interpolates lattice quantity to coordinate r.

Result is stored in the value variable. Returns true if coordinate r is on the lattice, false if out of the lattice. In the latter case, the value is set to the default value (usually 0).

Parameters
[in]rPosition where the physical quantity would be evaluated.
[out]valuePhysical quantity evaluated at the nearest cell to the given position.
Returns
Boolean indicates whether the position r is located inside the lattice.
Todo:
(oliiny): maybe 1-order interpolation instead of 0-order?

Definition at line 547 of file lattice.h.

547  {
548  const int ix = std::floor((r.x1() - origin_[0]) / cell_sizes_[0]);
549  const int iy = std::floor((r.x2() - origin_[1]) / cell_sizes_[1]);
550  const int iz = std::floor((r.x3() - origin_[2]) / cell_sizes_[2]);
551  if (out_of_bounds(ix, iy, iz)) {
552  value = T();
553  return false;
554  } else {
555  value = node(ix, iy, iz);
556  return true;
557  }
558  }
smash::RectangularLattice::out_of_bounds
bool out_of_bounds(int ix, int iy, int iz) const
Checks if 3D index is out of lattice bounds.
Definition: lattice.h:114
smash::RectangularLattice::node
T & node(int ix, int iy, int iz)
Take the value of a cell given its 3-D indices.
Definition: lattice.h:524

◆ iterate_sublattice()

template<typename T >
template<typename F >
void smash::RectangularLattice< T >::iterate_sublattice ( const std::array< int, 3 > &  lower_bounds,
const std::array< int, 3 > &  upper_bounds,
F &&  func 
)
inline

A sub-lattice iterator, which iterates in a 3D-structured manner and calls a function on every cell.

Template Parameters
FType of the function. Arguments are the current node and the 3 integer indices of the cell.
Parameters
[in]lower_boundsStarting numbers for iterating ix, iy, iz.
[in]upper_boundsEnding numbers for iterating ix, iy, iz.
[in]funcFunction acting on the cells (such as taking value).

Definition at line 571 of file lattice.h.

572  {
573  logg[LLattice].debug(
574  "Iterating sublattice with lower bound index (", lower_bounds[0], ",",
575  lower_bounds[1], ",", lower_bounds[2], "), upper bound index (",
576  upper_bounds[0], ",", upper_bounds[1], ",", upper_bounds[2], ")");
577 
578  if (periodic_) {
579  for (int iz = lower_bounds[2]; iz < upper_bounds[2]; iz++) {
580  const int z_offset = positive_modulo(iz, n_cells_[2]) * n_cells_[1];
581  for (int iy = lower_bounds[1]; iy < upper_bounds[1]; iy++) {
582  const int y_offset =
583  n_cells_[0] * (positive_modulo(iy, n_cells_[1]) + z_offset);
584  for (int ix = lower_bounds[0]; ix < upper_bounds[0]; ix++) {
585  const int index = positive_modulo(ix, n_cells_[0]) + y_offset;
586  func(lattice_[index], ix, iy, iz);
587  }
588  }
589  }
590  } else {
591  for (int iz = lower_bounds[2]; iz < upper_bounds[2]; iz++) {
592  const int z_offset = iz * n_cells_[1];
593  for (int iy = lower_bounds[1]; iy < upper_bounds[1]; iy++) {
594  const int y_offset = n_cells_[0] * (iy + z_offset);
595  for (int ix = lower_bounds[0]; ix < upper_bounds[0]; ix++) {
596  func(lattice_[ix + y_offset], ix, iy, iz);
597  }
598  }
599  }
600  }
601  }

◆ iterate_in_cube()

template<typename T >
template<typename F >
void smash::RectangularLattice< T >::iterate_in_cube ( const ThreeVector &  point,
const double  r_cut,
F &&  func 
)
inline

Iterates only nodes whose cell centers lie not further than r_cut in x, y, z directions from the given point, that is iterates within a cube of side length 2*r_cut, and applies a function to each node.

Useful for adding quantities from one particle to the lattice.

Template Parameters
FType of the function. Arguments are the current node and the 3 integer indices of the cell.
Parameters
[in]pointPosition, usually the position of particle [fm].
[in]r_cutMaximum distance from the cell center to the given position. [fm]
[in]funcFunction acting on the cells (such as taking value).

Definition at line 617 of file lattice.h.

617  {
618  std::array<int, 3> l_bounds, u_bounds;
619 
620  /* Array holds value at the cell center: r_center = r_0 + (i+0.5)cell_size,
621  * where i is index in any direction. Therefore we want cells with condition
622  * (r-r_cut)/csize - 0.5 < i < (r+r_cut)/csize - 0.5, r = r_center - r_0 */
623  for (int i = 0; i < 3; i++) {
624  l_bounds[i] =
625  std::ceil((point[i] - origin_[i] - r_cut) / cell_sizes_[i] - 0.5);
626  u_bounds[i] =
627  std::ceil((point[i] - origin_[i] + r_cut) / cell_sizes_[i] - 0.5);
628  }
629 
630  if (!periodic_) {
631  for (int i = 0; i < 3; i++) {
632  if (l_bounds[i] < 0) {
633  l_bounds[i] = 0;
634  }
635  if (u_bounds[i] > n_cells_[i]) {
636  u_bounds[i] = n_cells_[i];
637  }
638  if (l_bounds[i] > n_cells_[i] || u_bounds[i] < 0) {
639  return;
640  }
641  }
642  }
643  iterate_sublattice(l_bounds, u_bounds, std::forward<F>(func));
644  }
smash::RectangularLattice::iterate_sublattice
void iterate_sublattice(const std::array< int, 3 > &lower_bounds, const std::array< int, 3 > &upper_bounds, F &&func)
A sub-lattice iterator, which iterates in a 3D-structured manner and calls a function on every cell.
Definition: lattice.h:571

◆ integrate_volume()

template<typename T >
template<typename F >
void smash::RectangularLattice< T >::integrate_volume ( F &  integral,
F(*)(ThreeVector, T &, ThreeVector)  integrand,
const double  rcut,
const ThreeVector &  point 
)
inline

Calculate a volume integral with given integrand.

Template Parameters
Freturn type of the integrand
Parameters
[out]integralvariable to store rsult in. Input should be 0.
[in]integrandFunction to be integrated
[in]rcutsize of the integration volume. In total the intgration volume will be a cube with edge length 2*rcut
[in]pointcenter of the integration volume

Definition at line 658 of file lattice.h.

660  {
661  iterate_in_rectangle(
662  point, {rcut, rcut, rcut},
663  [&point, &integral, &integrand, this](T value, int ix, int iy, int iz) {
664  ThreeVector pos = this->cell_center(ix, iy, iz);
665  integral += integrand(pos, value, point) * this->cell_volume_;
666  });
667  }
smash::RectangularLattice::iterate_in_rectangle
void iterate_in_rectangle(const ThreeVector &point, const std::array< double, 3 > &rectangle, F &&func)
Iterates only nodes whose cell centers lie not further than d_x in x-, d_y in y-, and d_z in z-direct...
Definition: lattice.h:683

◆ iterate_in_rectangle()

template<typename T >
template<typename F >
void smash::RectangularLattice< T >::iterate_in_rectangle ( const ThreeVector &  point,
const std::array< double, 3 > &  rectangle,
F &&  func 
)
inline

Iterates only nodes whose cell centers lie not further than d_x in x-, d_y in y-, and d_z in z-direction from the given point, that is iterates within a rectangle of side lengths (2*dx, 2*dy, 2*dz), and applies a function to each node.

Useful for adding quantities from one particle to the lattice.

Template Parameters
FType of the function. Arguments are the current node and the 3 integer indices of the cell.
Parameters
[in]pointPosition, usually the position of particle [fm].
[in]rectangleMaximum distances in the x-, y-, and z-directions from the cell center to the given position. [fm]
[in]funcFunction acting on the cells (such as taking value).

Definition at line 683 of file lattice.h.

684  {
685  std::array<int, 3> l_bounds, u_bounds;
686 
687  /* Array holds value at the cell center: r_center = r_0 + (i+0.5)cell_size,
688  * where i is index in any direction. Therefore we want cells with condition
689  * (r[i]-rectangle[i])/csize - 0.5 < i < (r[i]+rectangle[i])/csize - 0.5,
690  * r[i] = r_center[i] - r_0[i]
691  */
692  for (int i = 0; i < 3; i++) {
693  l_bounds[i] = std::ceil(
694  (point[i] - origin_[i] - rectangle[i]) / cell_sizes_[i] - 0.5);
695  u_bounds[i] = std::ceil(
696  (point[i] - origin_[i] + rectangle[i]) / cell_sizes_[i] - 0.5);
697  }
698 
699  if (!periodic_) {
700  for (int i = 0; i < 3; i++) {
701  if (l_bounds[i] < 0) {
702  l_bounds[i] = 0;
703  }
704  if (u_bounds[i] > n_cells_[i]) {
705  u_bounds[i] = n_cells_[i];
706  }
707  if (l_bounds[i] > n_cells_[i] || u_bounds[i] < 0) {
708  return;
709  }
710  }
711  }
712  iterate_sublattice(l_bounds, u_bounds, std::forward<F>(func));
713  }

◆ iterate_nearest_neighbors()

template<typename T >
template<typename F >
void smash::RectangularLattice< T >::iterate_nearest_neighbors ( const ThreeVector &  point,
F &&  func 
)
inline

Iterates only over nodes corresponding to the center cell (the cell containing the given point) and its nearest neighbors in the -x, +x, -y, +y, -z, +z directions, and applies a function to each node.

Useful for adding quantities from one particle to the lattice.

Template Parameters
FType of the function. Arguments are the current node and the 3 integer indices of the cell.
Parameters
[in]pointPosition, usually the position of particle [fm].
[in]funcFunction acting on the cells (such as taking value).

Definition at line 727 of file lattice.h.

727  {
728  // get the 3D indeces of the cell containing the given point
729  const int ix = std::floor((point.x1() - origin_[0]) / cell_sizes_[0]);
730  const int iy = std::floor((point.x2() - origin_[1]) / cell_sizes_[1]);
731  const int iz = std::floor((point.x3() - origin_[2]) / cell_sizes_[2]);
732 
733  logg[LLattice].debug(
734  "Iterating over nearest neighbors of the cell at ix = ", ix,
735  ", iy = ", iy, ", iz = ", iz);
736 
737  // determine the 1D index of the center cell
738  const int i = index1d(ix, iy, iz);
739  func(lattice_[i], i, i);
740 
741  // determine the indeces of nearby cells, perform function on them
742  int ileft = index_left(ix, iy, iz);
743  func(lattice_[ileft], ileft, i);
744 
745  int iright = index_right(ix, iy, iz);
746  func(lattice_[iright], iright, i);
747 
748  int idown = index_down(ix, iy, iz);
749  func(lattice_[idown], idown, i);
750 
751  int iup = index_up(ix, iy, iz);
752  func(lattice_[iup], iup, i);
753 
754  int ibackward = index_backward(ix, iy, iz);
755  func(lattice_[ibackward], ibackward, i);
756 
757  int iforward = index_forward(ix, iy, iz);
758  func(lattice_[iforward], iforward, i);
759  }
smash::RectangularLattice::index_right
int index_right(int ix, int iy, int iz)
Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +x...
Definition: lattice.h:238
smash::RectangularLattice::index_left
int index_left(int ix, int iy, int iz)
Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -x...
Definition: lattice.h:221
smash::RectangularLattice::index_up
int index_up(int ix, int iy, int iz)
Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +y...
Definition: lattice.h:272
smash::RectangularLattice::index_backward
int index_backward(int ix, int iy, int iz)
Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -z...
Definition: lattice.h:289
smash::RectangularLattice::index_forward
int index_forward(int ix, int iy, int iz)
Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the +z...
Definition: lattice.h:306
smash::RectangularLattice::index_down
int index_down(int ix, int iy, int iz)
Given the indices of a cell in the x, y, and z directions, return index of the nearest cell in the -y...
Definition: lattice.h:255

◆ identical_to_lattice()

template<typename T >
template<typename L >
bool smash::RectangularLattice< T >::identical_to_lattice ( const L *  lat) const
inline

Checks if lattices of possibly different types have identical structure.

Template Parameters
LType of the other lattice.
Parameters
[in]latThe other lattice being compared with the current one
Returns
Whether the two lattices have the same sizes, cell numbers, origins, and boundary conditions.

Definition at line 770 of file lattice.h.

770  {
771  return n_cells_[0] == lat->n_cells()[0] &&
772  n_cells_[1] == lat->n_cells()[1] &&
773  n_cells_[2] == lat->n_cells()[2] &&
774  std::abs(lattice_sizes_[0] - lat->lattice_sizes()[0]) <
775  really_small &&
776  std::abs(lattice_sizes_[1] - lat->lattice_sizes()[1]) <
777  really_small &&
778  std::abs(lattice_sizes_[2] - lat->lattice_sizes()[2]) <
779  really_small &&
780  std::abs(origin_[0] - lat->origin()[0]) < really_small &&
781  std::abs(origin_[1] - lat->origin()[1]) < really_small &&
782  std::abs(origin_[2] - lat->origin()[2]) < really_small &&
783  periodic_ == lat->periodic();
784  }
smash::really_small
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:37

◆ positive_modulo()

template<typename T >
int smash::RectangularLattice< T >::positive_modulo ( int  i,
int  n 
) const
inlineprivate

Returns division modulo, which is always between 0 and n-1 in is not suitable, because it returns results from -(n-1) to n-1.

Parameters
[in]iDividend.
[in]nDivisor.
Returns
Positive remainder.

Definition at line 813 of file lattice.h.

813  {
814  /* (i % n + n) % n would be correct, but slow.
815  * Instead I rely on the fact that i should never go too far
816  * in negative region and replace i%n + n by i + 256 * n = i + (n << 8) */
817  // FIXME: This should use asserts, also checking for under- or overflows.
818  return (i + (n << 8)) % n;
819  }

Member Data Documentation

◆ lattice_

template<typename T >
std::vector<T> smash::RectangularLattice< T >::lattice_
protected

The lattice itself, array containing physical quantities.

Definition at line 788 of file lattice.h.

◆ lattice_sizes_

template<typename T >
const std::array<double, 3> smash::RectangularLattice< T >::lattice_sizes_
protected

Lattice sizes in x, y, z directions.

Definition at line 790 of file lattice.h.

◆ n_cells_

template<typename T >
const std::array<int, 3> smash::RectangularLattice< T >::n_cells_
protected

Number of cells in x,y,z directions.

Definition at line 792 of file lattice.h.

◆ cell_sizes_

template<typename T >
const std::array<double, 3> smash::RectangularLattice< T >::cell_sizes_
protected

Cell sizes in x, y, z directions.

Definition at line 794 of file lattice.h.

◆ cell_volume_

template<typename T >
const double smash::RectangularLattice< T >::cell_volume_
protected

Volume of a cell.

Definition at line 796 of file lattice.h.

◆ origin_

template<typename T >
const std::array<double, 3> smash::RectangularLattice< T >::origin_
protected

Coordinates of the left down nearer corner.

Definition at line 798 of file lattice.h.

◆ periodic_

template<typename T >
const bool smash::RectangularLattice< T >::periodic_
protected

Whether the lattice is periodic.

Definition at line 800 of file lattice.h.

◆ when_update_

template<typename T >
const LatticeUpdate smash::RectangularLattice< T >::when_update_
protected

When the lattice should be recalculated.

Definition at line 802 of file lattice.h.

The documentation for this class was generated from the following file: