Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68411141
spatial_grid.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jun 27, 06:54
Size
32 KB
Mime Type
text/x-c++
Expires
Sat, Jun 29, 06:54 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18579201
Attached To
rLIBMULTISCALE LibMultiScale
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__ */
Event Timeline
Log In to Comment