Homec4science
Diffusion LibMultiScale (features/components) Browse

libmultiscale/src/geometry/spatial_grid.hh
b05d5f2aba78features/components

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
/* -------------------------------------------------------------------------- */
__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 {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
//! constructor provided the number of celles in each direction
SpatialGrid(const Real * xmin, const Real * xmax,const UInt * s);
//! constructor provided the sizes of the cells in each direction
SpatialGrid(const Real * xmin, const Real * xmax,const Real * sizes);
//! constructor provided the sizes of the cells in each direction
SpatialGrid(const Real * xmin, const Real * xmax,const Real * sizes, const UInt * divisions);
//! constructor provided the size for all directions
SpatialGrid(const Real * xmin, const Real * xmax,Real p);
//! destructor
virtual ~SpatialGrid();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! add a finite element given
void addFiniteElement(T & el,std::vector<Real> & nodal_coord);
//! add a single element at the block given its index
void addElement(T & el,UInt 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,Real * 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,Real * x);
//! add a single element at the block given by its cartesian indexes
void addElement(T & el,UInt * indexes);
//! return a vector with the block index of the current neighborhood
void getNeighborhood(Real * x, std::vector<UInt> & neighbors, UInt * pbc = NULL);
//! return a vector with the block index of the current neighborhood
void getNeighborhood(UInt * indexes, std::vector<UInt> & neighbors, UInt * pbc = NULL);
//! return a vector with the block index of the current neighborhood
void getNeighborhood(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
void index2Coordinates(const UInt index,Real & x,Real & y,Real & z);
//! convert block index to cartesian index
void index2CartesianIndex(const UInt index,UInt * indexes);
//! convert block index to coordinates
void index2Coordinates(const UInt index,Real * x);
//! convert spatial coordinates to block index
void coordinates2Index(const Real * Xcoord, UInt & index);
//! convert spatial coordinates to block index
void cartesianIndex2Index(const UInt * indexes, UInt & index);
//! convert spatial coordinates to cartesian block index
void coordinates2CartesianIndex(const Real * Xcoord,UInt * indexes);
//! convert spatial coordinates to cartesian block index (no check if inside the grid)
void coordinates2CartesianIndexNoCheck(const Real * x, UInt * index);
private:
//! add an element to all grid cells intersecting a bounding box
void addToBoundingBox(T & el,Real * min, Real * max);
//! initialize the structure provided number of cells in each direction
void init(UInt dim,const Real * xmin, const Real * xmax, const UInt * s);
//! initialize the structure provided size of cells in each direction
void init(UInt dim,const Real * xmin, const Real * xmax, const 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(const Real * 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 Real * xmin, const Real * xmax,const UInt * s){
for (UInt i = 0; i < Dim; ++i) {
init(i,xmin,xmax,s);
}
allocate();
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
SpatialGrid<T,Dim>::SpatialGrid(const Real * xmin, const Real * xmax,const Real * p, const UInt * s){
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 Real * xmin, const Real * xmax,const Real * p){
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 Real * xmin, const Real * xmax,Real p){
Real _p[3];
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 Real * xmi, const Real * xma,const 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 Real * xmi, const Real * xma, const 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>::addElement(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,Real * min, Real * max){
UInt index_min[3];
UInt index_max[3];
UInt indexes[3];
coordinates2CartesianIndex(min,index_min);
coordinates2CartesianIndex(max,index_max);
if (Dim==1){
for (UInt i = index_min[0]; i <= index_max[0] ; ++i){
addElement(el,&i);
}
}
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;
addElement(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;
addElement(el,indexes);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::addFiniteElement(T & el,std::vector<Real> & nodal_coord){
Real max[3] = {0.,0.,0.};
Real min[3] = {0.,0.,0.};
for (UInt i = 0; i < nodal_coord.size(); ++i) {
UInt coord = i%Dim;
if (i/Dim == 0) {
max[coord] = - lm_real_max;
min[coord] = lm_real_max;
}
max[coord] = std::max(nodal_coord[i],max[coord]);
min[coord] = std::min(nodal_coord[i],min[coord]);
}
for (UInt i = 0; i < Dim; ++i) {
min[i] = std::max(xmin[i],min[i]);
max[i] = std::min(xmax[i],max[i]);
DUMP("actual box(x) : " << min[i] << "-" << max[i],DBG_DETAIL);
}
//I add the element at the intersection with grid cell of this bounding box
addToBoundingBox(el,min,max);
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::getNeighborhood(Real * x, std::vector<UInt> & neighbors,
UInt * pbc){
UInt indexes[Dim];
coordinates2CartesianIndex(x,indexes);
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;
getNeighborhood(indexes,neighbors,pbc);
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::getNeighborhood(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>::getNeighborhood(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];
UInt ii[3] = {0,0,0};
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,index);
#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);
neighbors.push_back(index);
}
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::addAtom(T & value,Real * x){
getNeighborhood(x,neighborhood);
for (UInt i = 0; i < neighborhood.size(); ++i) {
addElement(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>
void SpatialGrid<T,Dim>::addElement(T & el,Real * x){
UInt index;
coordinates2Index(x,index);
DUMP("adding element " << &el << " at position "
<< x[0] << " " << x[1] << " " << x[2] << " : in box "
<< index,DBG_DETAIL);
addElement(el,index);
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::addElement(T & el,UInt * indexes){
UInt index;
cartesianIndex2Index(indexes,index);
addElement(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(const Real * x){
UInt index;
coordinates2Index(x,index);
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 >= 0 && 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>
void SpatialGrid<T,Dim>::cartesianIndex2Index(const UInt * indexes, UInt & index){
for (UInt i = 0; i < Dim; ++i) {
if(indexes[i] == UINT_MAX) {
index = UINT_MAX;
return;
}
}
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);
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::coordinates2Index(const Real * Xcoord, UInt & index){
UInt norm[Dim];
coordinates2CartesianIndex(Xcoord,norm);
cartesianIndex2Index(norm,index);
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::coordinates2CartesianIndex(const Real * x, UInt * index){
for (UInt i = 0; i < Dim; ++i) {
if(x[i] < xmin[i] || x[i] > xmax[i]){
for (UInt j = 0; j < Dim; ++j)
index[j] = UINT_MAX;
return;
}
}
coordinates2CartesianIndexNoCheck(x,index);
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::coordinates2CartesianIndexNoCheck(const Real * x, 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);
}
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::index2Coordinates(const UInt index,Real * x){
UInt indexes[Dim];
index2CartesianIndex(index,indexes);
for (UInt i = 0; i < Dim; ++i) {
x[i] = xmin[i] + cell_size[i]/2. + cell_size[i]*indexes[i];
}
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::index2Coordinates(const UInt index,Real & x, Real & y, Real & z){
Real Xcoord[Dim];
index2Coordinates(index,Xcoord);
x = Xcoord[0];
if (Dim > 1) y = Xcoord[1];
if (Dim == 3) z = Xcoord[2];
}
/* -------------------------------------------------------------------------- */
template <typename T,UInt Dim>
void SpatialGrid<T,Dim>::index2CartesianIndex(const UInt index,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];
}
}
/* -------------------------------------------------------------------------- */
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