Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/geometry/spatial_grid.hh
4fcaa8b3c306lammps_2020

Last Change
Raw File

spatial_grid.hh
View Options

/**
* @file spatial_grid.hh
*
* @author Till Junge <till.junge@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 18 22:30:00 2013
*
* @brief Spatial grid used to sort objects
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_SPATIAL_GRID_HH__
#define __LIBMULTISCALE_SPATIAL_GRID_HH__
/* -------------------------------------------------------------------------- */
#include <cmath>
#include <vector>
#ifndef LM_OPTIMIZED
#include <sstream>
#endif // LM_OPTIMIZED
#include "cube.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
// neighbors coords
/* -------------------------------------------------------------------------- */
template <UInt Dim> struct NeighborIndexes {};
template <> struct NeighborIndexes<3> {
static const int topologyCoords[27][3];
static const UInt n_neighs = 27;
};
template <> struct NeighborIndexes<2> {
const static int topologyCoords[9][2];
const static UInt n_neighs = 9;
};
template <> struct NeighborIndexes<1> {
const static int topologyCoords[3][1];
const static UInt n_neighs = 3;
};
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> class SpatialGrid : public virtual LMObject {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
//! constructor provided the number of celles in each direction
SpatialGrid(const Vector<Dim, Real> &xmin, const Vector<Dim, Real> &xmax,
const Vector<Dim, UInt> &s);
//! constructor provided the sizes of the cells in each direction
SpatialGrid(const Vector<Dim, Real> &xmin, const Vector<Dim, Real> &xmax,
const Vector<Dim, Real> &sizes);
//! constructor provided the sizes of the cells in each direction
SpatialGrid(const Vector<Dim, Real> &xmin, const Vector<Dim, Real> &xmax,
const Vector<Dim, Real> &sizes,
const Vector<Dim, UInt> &divisions);
//! constructor provided the size for all directions
SpatialGrid(const Vector<Dim, Real> &xmin, const Vector<Dim, Real> &xmax,
Real p);
//! destructor
virtual ~SpatialGrid();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! add a finite element given
void addFiniteElement(T &el, std::vector<Vector<Dim>> &nodal_coord);
//! add a single element at the block given its index
void addElementFromIndex(T &el, UInt index);
template <typename Vec> void addElement(T &el, Vec &&index);
//! flag to say wether the elements outside the grid boundaries should be
//! considered
void doNotSearchOutSide();
//! add element including to neighbors
void addAtom(T &value, VectorView<Dim> x);
// //! add a single element at the block given by its position
// void addElement(T &el, Real x, Real y, Real z);
// //! add element including to neighbors
// void addAtom(T &value, Real x, Real y, Real z);
//! add a single element at the block given by its position
void addElement(T &el, Vector<Dim> &x);
//! add a single element at the block given by its cartesian indexes
void addElementFromIndexes(T &el, const Vector<Dim, UInt> &indexes);
//! return a vector with the block index of the current neighborhood
void getNeighborhood(const Vector<Dim> &x, std::vector<UInt> &neighbors,
UInt *pbc = NULL);
//! return a vector with the block index of the current neighborhood
void getNeighborhoodFromIndexes(const Vector<Dim, UInt> &indexes,
std::vector<UInt> &neighbors,
UInt *pbc = NULL);
//! return a vector with the block index of the current neighborhood
void getNeighborhoodFromIndex(UInt index, std::vector<UInt> &neighbors,
UInt *pbc = NULL);
//! return the half diagonal distance of blocks
Real getHalfDiagonal();
//! get the squared (to avoid computing the square root) distance between two
//! blocks
// Real getSquaredDistanceToBlockCenter(Real * pos, UInt j);
//! get the min and max coords of the grid
void getMinMax(Real *min, Real *max);
private:
/* ------------------------------------------------------------------------ */
/* conversions functions */
/* ------------------------------------------------------------------------ */
public:
//! convert spatial coordinates to block index
Vector<Dim> index2Coordinates(const UInt index);
//! convert block index to cartesian index
Vector<Dim, UInt> index2CartesianIndex(const UInt index);
//! convert spatial coordinates to block index
template <typename Vec> UInt coordinates2Index(Vec &&Xcoord);
//! convert spatial coordinates to block index
UInt cartesianIndex2Index(const Vector<Dim, UInt> &indexes);
//! convert spatial coordinates to cartesian block index
template <typename Vec>
Vector<Dim, UInt> coordinates2CartesianIndex(Vec &&Xcoord);
//! convert spatial coordinates to cartesian block index (no check if inside
//! the grid)
template <typename Vec>
Vector<Dim, UInt> coordinates2CartesianIndexNoCheck(Vec &&x);
private:
//! add an element to all grid cells intersecting a bounding box
void addToBoundingBox(T &el, Cube &bbox);
//! initialize the structure provided number of cells in each direction
void init(UInt dim, const Vector<Dim> &xmin, const Vector<Dim> &xmax,
const Vector<Dim, UInt> &s);
//! initialize the structure provided size of cells in each direction
void init(UInt dim, const Vector<Dim> &xmin, const Vector<Dim> &xmax,
const Vector<Dim, Real> &p);
//! allocate the structure
void allocate();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! return the total number of cells
UInt getSize();
//! return the number of cells in a provided direction
UInt getCellNumber(UInt i);
//! return the size of cells in a provided direction
Real getCellSize(UInt i);
//! return the average number of elements per block
UInt getAverageEltByBlock();
//! return the number of elements in block index
UInt getNumberElt(UInt index);
//! return the block vectors for direct use
std::vector<std::vector<T> *> &getBlocks();
//! return the outside block
std::vector<T> &getBlockOutside();
//! return a vector of elements provided the block coordinates
std::vector<T> &findSet(Real x, Real y, Real z);
//! return a vector of elements provided the block coordinates
std::vector<T> &findSet(Vector<Dim> &x);
//! function which applies an operator to reduce the values per block
template <typename Op> T extractBoxValuePerOperator(UInt index, Op &op);
//! function which returns the connectivity
std::vector<UInt> &getConnectivity();
//! build the nodes of the grid
std::vector<Real> &getNodes();
//! build a vector of data using a reduce operator
// template <typename Op>
// std::vector<T> & getCellDataFromOperator(Op & op);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! max bounds of the grid
Real xmax[Dim];
//! min bounds of the grid
Real xmin[Dim];
//! inverse of cell size
Real inv_cell_size[Dim];
//! number of grid elements per direction
UInt cell_number[Dim];
//! size of grid elements per direction
Real cell_size[Dim];
//! total number of elements
UInt size;
//! grid of blocks
std::vector<std::vector<T> *> blocks;
//! the block containing outside elements
std::vector<T> block_outside;
//! the empty block to be returned
std::vector<T> empty_block;
//! flag to know if we need to search/store elements outside of the grid
bool search_out_flag;
//! vector to store the connectivity
std::vector<UInt> connectivity;
//! vector to store the node coordinates
std::vector<Real> nodes;
private:
//! storage of neighbor indexes (internal use only)
std::vector<UInt> neighborhood;
};
/* -------------------------------------------------------------------------- */
// marco to convert three coordinates to a small array
#define CONVERT_xyz_2_X \
Real Xcoord[3]; \
Xcoord[0] = x; \
if (Dim > 1) \
Xcoord[1] = y; \
if (Dim == 3) \
Xcoord[2] = z;
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
SpatialGrid<T, Dim>::SpatialGrid(const Vector<Dim> &xmin,
const Vector<Dim> &xmax,
const Vector<Dim, UInt> &s)
: LMObject("spatial_grid") {
for (UInt i = 0; i < Dim; ++i) {
init(i, xmin, xmax, s);
}
allocate();
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
SpatialGrid<T, Dim>::SpatialGrid(const Vector<Dim> &xmin,
const Vector<Dim> &xmax, const Vector<Dim> &p,
const Vector<Dim, UInt> &s)
: LMObject("spatial_grid") {
for (UInt i = 0; i < Dim; ++i) {
if (p[i] == -1 && s[i] == UINT_MAX)
LM_FATAL("cannot initialize spatialgrid object");
if (p[i] != -1)
init(i, xmin, xmax, p);
else
init(i, xmin, xmax, s);
}
allocate();
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
SpatialGrid<T, Dim>::SpatialGrid(const Vector<Dim> &xmin,
const Vector<Dim> &xmax, const Vector<Dim> &p)
: LMObject("spatial_grid") {
for (UInt i = 0; i < Dim; ++i) {
if (i > 1 && p[i] <= -1)
cell_size[i] = cell_size[0];
init(i, xmin, xmax, p);
}
allocate();
}
/* ------------------------------------------------------------------------ */
template <typename T, UInt Dim>
SpatialGrid<T, Dim>::SpatialGrid(const Vector<Dim> &xmin,
const Vector<Dim> &xmax, Real p)
: LMObject("spatial_grid") {
Vector<Dim, Real> _p;
for (UInt i = 0; i < Dim; ++i) {
_p[i] = p;
}
for (UInt i = 0; i < Dim; ++i) {
init(i, xmin, xmax, _p);
}
allocate();
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::init(UInt dim, const Vector<Dim> &xmi,
const Vector<Dim> &xma,
const Vector<Dim, UInt> &cell_n) {
xmax[dim] = xma[dim];
xmin[dim] = xmi[dim];
cell_number[dim] = cell_n[dim];
cell_size[dim] = (xmax[dim] - xmin[dim]) / cell_number[dim];
search_out_flag = true;
if (cell_size[dim] == 0)
LM_FATAL("grid size is zero along axis " << dim << " : Check geometries - "
<< xmax[dim] << " " << xmin[dim]);
inv_cell_size[dim] = 1. / cell_size[dim];
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::init(UInt dim, const Vector<Dim> &xmi,
const Vector<Dim> &xma,
const Vector<Dim, Real> &cell_s) {
xmax[dim] = xma[dim];
xmin[dim] = xmi[dim];
cell_size[dim] = cell_s[dim];
cell_number[dim] = (int)((xmax[dim] - xmin[dim]) / cell_size[dim]);
if (cell_number[dim] <= 0)
LM_FATAL("grid size is <= 0 (" << cell_number[dim] << ") along axis " << dim
<< " cell_size is " << cell_size[dim]
<< " range " << xmin[dim] << " " << xmax[dim]
<< " : Check geometries");
cell_size[dim] = (xmax[dim] - xmin[dim]) / cell_number[dim];
inv_cell_size[dim] = 1. / cell_size[dim];
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> void SpatialGrid<T, Dim>::allocate() {
// compute the total number of cells
size = 1;
for (UInt k = 0; k < Dim; ++k)
size *= cell_number[k];
for (UInt i = 0; i < Dim; ++i) {
DUMP("Creating grid of size[" << i << "] = " << cell_number[i] << " - "
<< size << " : " << xmin[i] << " " << xmax[i]
<< " "
<< " cell_size " << cell_size[i],
DBG_INFO_STARTUP);
}
blocks.resize(size);
for (UInt i = 0; i < size; ++i) {
blocks[i] = NULL;
}
if (Dim == 1)
neighborhood.reserve(3);
if (Dim == 2)
neighborhood.reserve(9);
if (Dim == 3)
neighborhood.reserve(27);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> SpatialGrid<T, Dim>::~SpatialGrid() {
for (UInt i = 0; i < size; ++i) {
if (blocks[i]) {
blocks[i]->clear();
delete blocks[i];
}
}
block_outside.clear();
empty_block.clear();
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> void SpatialGrid<T, Dim>::doNotSearchOutSide() {
search_out_flag = false;
}
/* -------------------------------------------------------------------------- */
// insertion methods
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::addElementFromIndex(T &el, UInt index) {
DUMP("adding element " << &el << " : in box " << index, DBG_DETAIL);
if (index == UINT_MAX) {
block_outside.push_back(el);
DUMP("accessing outside block ", DBG_DETAIL);
return;
}
if (!blocks[index])
blocks[index] = new std::vector<T>();
LM_ASSERT(index < size, "overflow access " << index << " >= " << size);
DUMP("accessing block " << index << " p= " << blocks[index], DBG_DETAIL);
blocks[index]->push_back(el);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::addToBoundingBox(T &el, Cube &bbox) {
// make the intersection with the grid bounding box
auto min = bbox.getXmin();
auto max = bbox.getXmax();
for (UInt i = 0; i < Dim; ++i) {
min[i] = std::max(this->xmin[i], min[i]);
max[i] = std::min(this->xmax[i], max[i]);
DUMP("actual box(x) : " << min[i] << "-" << max[i], DBG_DETAIL);
}
Vector<Dim, UInt> indexes;
auto index_min = coordinates2CartesianIndex(min);
auto index_max = coordinates2CartesianIndex(max);
if (Dim == 1) {
for (UInt i = index_min[0]; i <= index_max[0]; ++i) {
indexes[0] = i;
addElementFromIndexes(el, indexes);
}
}
if (Dim == 2) {
for (UInt i = index_min[0]; i <= index_max[0]; ++i) {
indexes[0] = i;
for (UInt j = index_min[1]; j <= index_max[1]; ++j) {
indexes[1] = j;
addElementFromIndexes(el, indexes);
}
}
}
if (Dim == 3) {
for (UInt i = index_min[0]; i <= index_max[0]; ++i) {
indexes[0] = i;
for (UInt j = index_min[1]; j <= index_max[1]; ++j) {
indexes[1] = j;
for (UInt k = index_min[2]; k <= index_max[2]; ++k) {
indexes[2] = k;
addElementFromIndexes(el, indexes);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::addFiniteElement(
T &el, std::vector<Vector<Dim>> &nodal_coords) {
if (nodal_coords.size() == 0)
return;
Vector<Dim> min = nodal_coords[0];
Vector<Dim> max = nodal_coords[0];
for (auto &&coord : nodal_coords) {
for (UInt d = 0; d < Dim; ++d) {
max[d] = std::max(coord[d], max[d]);
min[d] = std::min(coord[d], min[d]);
}
}
// I add the element at the intersection with grid cell of this bounding box
Cube c;
c.setXmin(min);
c.setXmax(max);
addToBoundingBox(el, c);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::getNeighborhood(const Vector<Dim> &x,
std::vector<UInt> &neighbors,
UInt *pbc) {
const auto indexes = coordinates2CartesianIndex(x);
bool flag = false;
for (unsigned int i = 0; i < Dim; ++i) {
if (indexes[i] == UINT_MAX) {
// DUMP("atom at position " << x[i]
// << " out of the grid (dim " << i << ") not in ["
// << xmin[i] << "," << xmax[i] << "]",DBG_WARNING);
flag = true;
}
}
if (flag)
return;
getNeighborhoodFromIndexes(indexes, neighbors, pbc);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::getNeighborhoodFromIndex(UInt index,
std::vector<UInt> &neighbors,
UInt *pbc) {
UInt indexes[Dim];
index2CartesianIndex(index, indexes);
getNeighborhood(indexes, neighbors, pbc);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
inline void SpatialGrid<T, Dim>::getNeighborhoodFromIndexes(
const Vector<Dim, UInt> &indexes, std::vector<UInt> &neighbors, UInt *pbc) {
UInt pbc_flag[Dim];
if (pbc == NULL)
for (UInt i = 0; i < Dim; ++i)
pbc_flag[i] = 0;
else
for (UInt i = 0; i < Dim; ++i)
pbc_flag[i] = pbc[i];
Vector<Dim, UInt> ii;
ii.setZero();
neighbors.clear();
for (UInt n = 0; n < NeighborIndexes<Dim>::n_neighs; ++n) {
const int *neigh_coord = NeighborIndexes<Dim>::topologyCoords[n];
UInt k = 0;
for (k = 0; k < Dim; ++k) {
if (indexes[k] == 0 && neigh_coord[k] == -1) {
if (pbc_flag[k])
ii[k] = cell_number[k] - 1;
else
break;
} else if (indexes[k] + neigh_coord[k] == cell_number[k]) {
if (pbc_flag[k])
ii[k] = 0;
else
break;
} else
ii[k] = indexes[k] + neigh_coord[k];
}
// one of the dimensions is outside of the domain: go to next neighbor
if (k < Dim)
continue;
UInt index = cartesianIndex2Index(ii);
#ifndef LM_OPTIMIZED
std::stringstream sstr_msg;
sstr_msg << index << " ";
for (UInt k = 0; k < Dim; ++k)
sstr_msg << indexes[k] << " ";
sstr_msg << std::endl;
for (UInt k = 0; k < Dim; ++k)
sstr_msg << " in [" << xmin[k] << " - " << xmax[k] << "]" << std::endl;
#endif // LM_OPTIMIZED
LM_ASSERT(index != UINT_MAX,
"invalid neighbor block number " << sstr_msg.str());
neighbors.push_back(index);
}
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::addAtom(T &value, VectorView<Dim> x [[gnu::unused]]) {
LM_TOIMPLEMENT;
// getNeighborhood(x.data(), neighborhood);
for (UInt i = 0; i < neighborhood.size(); ++i) {
addElementFromIndex(value, neighborhood[i]);
}
}
/* -------------------------------------------------------------------------- */
// template <typename T, UInt Dim>
// void SpatialGrid<T, Dim>::addAtom(T &value, Real x, Real y, Real z) {
// CONVERT_xyz_2_X;
// addAtom(value, Xcoord);
// }
// /* --------------------------------------------------------------------------
// */
// template <typename T, UInt Dim>
// void SpatialGrid<T, Dim>::addElement(T &el, Real x, Real y, Real z) {
// CONVERT_xyz_2_X;
// addElement(el, Xcoord);
// }
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
template <typename Vec>
void SpatialGrid<T, Dim>::addElement(T &el, Vec &&x) {
UInt index = coordinates2Index(x);
DUMP("adding element " << &el << " at position " << x[0] << " " << x[1] << " "
<< x[2] << " : in box " << index,
DBG_DETAIL);
addElementFromIndex(el, index);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::addElementFromIndexes(
T &el, const Vector<Dim, UInt> &indexes) {
UInt index = cartesianIndex2Index(indexes);
addElementFromIndex(el, index);
}
/* -------------------------------------------------------------------------- */
// accessor functions
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
template <typename Op>
T SpatialGrid<T, Dim>::extractBoxValuePerOperator(UInt index, Op &op) {
std::vector<T> &block = *blocks[index];
op.init();
for (UInt i = 0; i < block.size(); ++i) {
op.account(block[i]);
}
return op.result();
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
std::vector<T> &SpatialGrid<T, Dim>::findSet(Real x, Real y, Real z) {
CONVERT_xyz_2_X;
return findSet(Xcoord);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
std::vector<T> &SpatialGrid<T, Dim>::findSet(Vector<Dim> &x) {
UInt index = coordinates2Index(x);
if (index == UINT_MAX) {
for (UInt i = 0; i < Dim; ++i) {
DUMP("findIndex(" << i << ") " << x[i] << " " << xmin[i] << " "
<< xmax[i],
DBG_MESSAGE);
}
DUMP("Warning : searching outside of the spatial grid", DBG_MESSAGE);
if (search_out_flag)
return block_outside;
else
return empty_block;
}
if (!blocks[index]) {
return empty_block;
}
return *blocks[index];
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> UInt SpatialGrid<T, Dim>::getSize() {
return size;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
UInt SpatialGrid<T, Dim>::getCellNumber(UInt i) {
LM_ASSERT(i < Dim, "out of bounds access");
return cell_number[i];
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> Real SpatialGrid<T, Dim>::getCellSize(UInt i) {
switch (i) {
case 0:
return cell_size[0];
break;
case 1:
return cell_size[1];
break;
case 2:
return cell_size[2];
break;
default:
LM_FATAL("this should not happend :"
<< " getSizeArray take 0,1 or 2 as argument->exit");
}
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
std::vector<std::vector<T> *> &SpatialGrid<T, Dim>::getBlocks() {
return blocks;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
std::vector<T> &SpatialGrid<T, Dim>::getBlockOutside() {
return block_outside;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
UInt SpatialGrid<T, Dim>::getAverageEltByBlock() {
UInt result = 0, nb = 0;
for (UInt i = 0; i < blocks.size(); ++i)
if (blocks[i]) {
result += blocks[i]->size();
++nb;
}
result += block_outside.size();
++nb;
return result / nb;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
UInt SpatialGrid<T, Dim>::getNumberElt(UInt index) {
if (blocks[index])
return blocks[index]->size();
return 0;
}
/* -------------------------------------------------------------------------- */
// conversion functions
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
UInt SpatialGrid<T, Dim>::cartesianIndex2Index(
const Vector<Dim, UInt> &indexes) {
for (UInt i = 0; i < Dim; ++i) {
if (indexes[i] == UINT_MAX) {
return UINT_MAX;
}
}
UInt index;
index = indexes[0];
if (Dim >= 2)
index += cell_number[0] * indexes[1];
if (Dim == 3)
index += cell_number[0] * cell_number[1] * indexes[2];
DUMP("found index = " << index, DBG_DETAIL);
return index;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
template <typename Vec>
UInt SpatialGrid<T, Dim>::coordinates2Index(Vec &&Xcoord) {
auto norm = coordinates2CartesianIndex(Xcoord);
return cartesianIndex2Index(norm);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
template <typename Vec>
Vector<Dim, UInt> SpatialGrid<T, Dim>::coordinates2CartesianIndex(Vec &&x) {
for (UInt i = 0; i < Dim; ++i) {
if (x[i] < xmin[i] || x[i] > xmax[i]) {
Vector<Dim, UInt> index;
for (UInt j = 0; j < Dim; ++j)
index[j] = UINT_MAX;
return index;
}
}
return coordinates2CartesianIndexNoCheck(x);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
template <typename Vec>
Vector<Dim, UInt>
SpatialGrid<T, Dim>::coordinates2CartesianIndexNoCheck(Vec &&x) {
Vector<Dim, UInt> index;
for (UInt i = 0; i < Dim; ++i) {
index[i] = (int)(inv_cell_size[i] * (x[i] - xmin[i]));
if (index[i] == cell_number[i])
--index[i];
DUMP("found index[" << i << "] = " << index[i], DBG_DETAIL);
}
return index;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
Vector<Dim> SpatialGrid<T, Dim>::index2Coordinates(const UInt index) {
auto indexes = index2CartesianIndex(index);
Vector<Dim> x;
for (UInt i = 0; i < Dim; ++i) {
x[i] = xmin[i] + cell_size[i] / 2. + cell_size[i] * indexes[i];
}
return x;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
Vector<Dim, UInt> SpatialGrid<T, Dim>::index2CartesianIndex(const UInt index) {
Vector<Dim, UInt> indexes;
if (Dim == 1)
indexes[0] = index;
if (Dim == 2) {
indexes[1] = index / cell_number[0];
indexes[0] = index - indexes[1] * cell_number[0];
}
if (Dim == 3) {
indexes[2] = index / cell_number[0] / cell_number[1];
indexes[1] =
(index - indexes[2] * cell_number[0] * cell_number[1]) / cell_number[0];
indexes[0] = index - indexes[2] * cell_number[0] * cell_number[1] -
indexes[1] * cell_number[0];
}
return indexes;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
std::vector<UInt> &SpatialGrid<T, Dim>::getConnectivity() {
if (connectivity.size())
return connectivity;
if (Dim == 1)
for (UInt i = 0; i < cell_number[0]; ++i) {
connectivity.push_back(i);
connectivity.push_back(i + 1);
}
if (Dim == 2)
for (UInt i = 0; i < cell_number[0]; ++i)
for (UInt j = 0; j < cell_number[1]; ++j) {
UInt p1, p2, p3, p4;
p1 = i + j * (cell_number[0] + 1);
p2 = i + 1 + j * (cell_number[0] + 1);
p3 = i + (j + 1) * (cell_number[0] + 1);
p4 = i + 1 + (j + 1) * (cell_number[0] + 1);
connectivity.push_back(p1);
connectivity.push_back(p2);
connectivity.push_back(p4);
connectivity.push_back(p3);
}
if (Dim == 3)
for (UInt i = 0; i < cell_number[0]; ++i)
for (UInt j = 0; j < cell_number[1]; ++j)
for (UInt k = 0; k < cell_number[2]; ++k) {
UInt p1, p2, p3, p4, p5, p6, p7, p8;
p1 = i + j * (cell_number[0] + 1) +
k * (cell_number[0] + 1) * (cell_number[1] + 1);
p2 = i + 1 + j * (cell_number[0] + 1) +
k * (cell_number[0] + 1) * (cell_number[1] + 1);
p3 = i + (j + 1) * (cell_number[0] + 1) +
k * (cell_number[0] + 1) * (cell_number[1] + 1);
p4 = i + 1 + (j + 1) * (cell_number[0] + 1) +
k * (cell_number[0] + 1) * (cell_number[1] + 1);
p5 = i + j * (cell_number[0] + 1) +
(k + 1) * (cell_number[0] + 1) * (cell_number[1] + 1);
p6 = i + 1 + j * (cell_number[0] + 1) +
(k + 1) * (cell_number[0] + 1) * (cell_number[1] + 1);
p7 = i + (j + 1) * (cell_number[0] + 1) +
(k + 1) * (cell_number[0] + 1) * (cell_number[1] + 1);
p8 = i + 1 + (j + 1) * (cell_number[0] + 1) +
(k + 1) * (cell_number[0] + 1) * (cell_number[1] + 1);
connectivity.push_back(p1);
connectivity.push_back(p2);
connectivity.push_back(p4);
connectivity.push_back(p3);
connectivity.push_back(p5);
connectivity.push_back(p6);
connectivity.push_back(p8);
connectivity.push_back(p7);
}
return connectivity;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
std::vector<Real> &SpatialGrid<T, Dim>::getNodes() {
if (nodes.size())
return nodes;
Real x[3];
for (UInt i = 0; i < Dim; ++i) {
x[i] = xmin[0];
}
if (Dim == 1)
for (UInt i = 0; i < cell_number[0] + 1; ++i, x[0] += cell_size[0]) {
nodes.push_back(x[0]);
}
if (Dim == 2)
for (UInt j = 0; j < cell_number[1] + 1; ++j, x[1] += cell_size[1]) {
x[0] = xmin[0];
for (UInt i = 0; i < cell_number[0] + 1; ++i, x[0] += cell_size[0]) {
nodes.push_back(x[0]);
nodes.push_back(x[1]);
}
}
if (Dim == 3)
for (UInt k = 0; k < cell_number[2] + 1; ++k, x[2] += cell_size[2]) {
x[1] = xmin[1];
for (UInt j = 0; j < cell_number[1] + 1; ++j, x[1] += cell_size[1]) {
x[0] = xmin[0];
for (UInt i = 0; i < cell_number[0] + 1; ++i, x[0] += cell_size[0]) {
nodes.push_back(x[0]);
nodes.push_back(x[1]);
nodes.push_back(x[2]);
}
}
}
return nodes;
}
/* -------------------------------------------------------------------------- */
// template <typename T,UInt Dim>
// template <typename Op>
// std::vector<T> & SpatialGrid<T,Dim>::getCellDataFromOperator(Op & op){
// for (UInt i = 0; i < size ; ++i){
// T v = extractBoxValuePerOperator(i, op);
// }
// }
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim> Real SpatialGrid<T, Dim>::getHalfDiagonal() {
Real res = 0;
for (UInt i = 0; i < Dim; ++i) {
res += 0.25 * cell_size[i] * cell_size[i];
}
return std::sqrt(res);
}
/* -------------------------------------------------------------------------- */
// template <typename T,UInt Dim>
// Real SpatialGrid<T,Dim>::getSquaredDistanceToBlockCenter(Real * position,
// UInt index){
// Real pos[Dim];
// index2Coordinates(index,pos);
// Real res = 0.;
// for (UInt i = 0; i < Dim; ++i) {
// Real tmp = position[i]-pos[i];
// res += tmp*tmp;
// }
// return res;
// }
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim>
void SpatialGrid<T, Dim>::getMinMax(Real *min, Real *max) {
for (UInt i = 0; i < Dim; ++i) {
min[i] = xmin[i];
max[i] = xmax[i];
}
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_SPATIAL_GRID_HH__ */
c4science · Help