Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69409962
fmm_multigrid.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
Mon, Jul 1, 16:28
Size
14 KB
Mime Type
text/x-c++
Expires
Wed, Jul 3, 16:28 (2 d)
Engine
blob
Format
Raw Data
Handle
18705378
Attached To
rLIBMULTISCALE LibMultiScale
fmm_multigrid.hh
View Options
/**
* @file fmm_multigrid.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Oct 11 13:00:16 2013
*
* @brief Part of the fast multipole implementation
*
* @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_FMM_MULTIGRID_HH__
#define __LIBMULTISCALE_FMM_MULTIGRID_HH__
/* -------------------------------------------------------------------------- */
#include "fmm_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
class
MaxNbChildren
{
public
:
// static const UInt max_nb_children;
};
template
<>
class
MaxNbChildren
<
1u
>
{
public
:
enum
{
max_nb_children
=
2
};
};
template
<>
class
MaxNbChildren
<
2u
>
{
public
:
enum
{
max_nb_children
=
4
};
};
template
<>
class
MaxNbChildren
<
3u
>
{
public
:
enum
{
max_nb_children
=
8
};
};
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
class
Block
{
public
:
Block
()
{
nb_children
=
0
;
}
void
addChildren
(
Block
<
Dim
>
&
b
)
{
children
[
nb_children
]
=
&
b
;
++
nb_children
;
};
Real
getSquaredDistanceFromBlock
(
Block
<
Dim
>
&
b
)
{
Real
res
=
0.
;
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
Real
tmp
=
b
.
coordinates
[
i
]
-
this
->
coordinates
[
i
];
res
+=
tmp
*
tmp
;
}
return
res
;
};
UInt
level
;
UInt
cartesianIndex
[
Dim
];
Real
coordinates
[
Dim
];
UInt
index
;
Block
*
father
;
Block
*
children
[
MaxNbChildren
<
Dim
>::
max_nb_children
];
UInt
nb_children
;
};
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
class
FMMMultigrid
:
public
FMMInterface
{
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public
:
typedef
FMM
<
FieldDim
,
Dim
,
order
>
myFMM
;
struct
Level
{
myFMM
*
fmm
;
std
::
vector
<
Block
<
Dim
>>
blocks
;
};
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public
:
inline
FMMMultigrid
(
Cube
&
c
,
UInt
start_level
,
UInt
nlevels
);
inline
virtual
~
FMMMultigrid
();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public
:
//! add a set of points to be taken into account in the grid
template
<
typename
Cont
>
inline
void
addPoints
(
Cont
&
positions
);
//! compute local moments
template
<
typename
Cont
>
void
computeMoments
(
Cont
&
cont
);
//! set the number of pbc replica in a given direction
void
setReplicaNumber
(
UInt
direction
,
UInt
nb_replica
)
{
replica
[
direction
]
=
nb_replica
;
};
//! compute the convolution at a given point
template
<
typename
F
,
typename
Cont
>
inline
void
computeConvolution
(
F
&
myFunc
,
Real
*
position
,
Cont
&
cont
,
Real
*
res
,
Real
rcut
=
-
1
);
//! build the interaction list with a recursive call
void
buildInteractionList
(
Real
rcut
);
private
:
//! shif a position for periodic boundary condition images
void
getReplicatedPosition
(
Real
*
local_pos
,
Real
*
position
,
int
sx
,
int
sy
,
int
sz
);
//! build the tree structure by examining the different levels
void
buildTreeStructure
();
//! build the interaction list: recursive function
void
buildInteractionListRecursive
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
);
//! add a multipole neighbor in the interaction list
void
addMultipoleNeighbor
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
);
//! add a close neighbor in the interaction list
void
addCloseNeighbor
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
);
//! return true if the block is outside the close neighborhood of provided
//! position
bool
isFarAway
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
);
//! return true if this block has no more children
bool
isChildLess
(
Block
<
Dim
>
&
far_block
);
//! simple function which call the right subgrid object for multipole
//! evaluation
template
<
typename
F
>
void
multipoleEvaluation
(
F
&
myFunc
,
Real
*
position
,
Block
<
Dim
>
&
far_block
,
Real
*
res
);
//! simple function which call the right subgrid object for fullEvaluation
template
<
typename
F
,
typename
Cont
>
void
fullEvaluation
(
F
&
myFunc
,
Real
*
position
,
Cont
&
cont
,
Block
<
Dim
>
&
far_block
,
Real
*
res
);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public
:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private
:
std
::
vector
<
Level
>
levels
;
int
replica
[
3
];
Real
domain_size
[
3
];
};
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
bool
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
isFarAway
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
)
{
UInt
level
=
far_block
.
level
;
UInt
indexes
[
Dim
];
levels
[
level
].
fmm
->
grid
.
coordinates2CartesianIndexNoCheck
(
current_block
.
coordinates
,
indexes
);
bool
test
=
false
;
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
int
i1
=
indexes
[
i
];
int
i2
=
far_block
.
cartesianIndex
[
i
];
int
_abs
=
std
::
abs
(
i1
-
i2
);
test
|=
(
_abs
>
1
);
}
return
test
;
};
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
bool
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
isChildLess
(
Block
<
Dim
>
&
far_block
)
{
return
(
!
far_block
.
nb_children
);
};
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
template
<
typename
F
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
multipoleEvaluation
(
F
&
myFunc
,
Real
*
position
,
Block
<
Dim
>
&
far_block
,
Real
*
res
)
{
UInt
level
=
far_block
.
level
;
UInt
index
=
far_block
.
index
;
levels
[
level
].
fmm
->
multipoleEvaluation
(
myFunc
,
position
,
index
,
res
);
};
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
template
<
typename
F
,
typename
Cont
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
fullEvaluation
(
F
&
myFunc
,
Real
*
position
,
Cont
&
cont
,
Block
<
Dim
>
&
far_block
,
Real
*
res
){
// UInt level = far_block.level;
// UInt index = far_block.index;
// levels[level].fmm->fullEvaluation(myFunc,position,cont,index,res);
};
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
FMMMultigrid
(
Cube
&
c
,
UInt
start_level
,
UInt
nlevels
)
{
Real
gridspace
[
Dim
];
UInt
griddivision
[
Dim
];
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
gridspace
[
i
]
=
-
1
;
for
(
UInt
n
=
0
;
n
<
nlevels
;
++
n
)
{
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
griddivision
[
i
]
=
std
::
pow
(
2
,
n
+
start_level
);
Level
level
;
level
.
fmm
=
new
myFMM
(
c
,
gridspace
,
griddivision
);
typename
myFMM
::
MyBlockList
&
blocks
=
level
.
fmm
->
grid
.
getBlocks
();
for
(
UInt
i
=
0
;
i
<
blocks
.
size
();
++
i
)
{
Block
<
Dim
>
b
;
b
.
level
=
n
;
b
.
index
=
i
;
level
.
fmm
->
grid
.
index2CartesianIndex
(
i
,
b
.
cartesianIndex
);
level
.
fmm
->
grid
.
index2Coordinates
(
i
,
b
.
coordinates
);
b
.
father
=
NULL
;
level
.
blocks
.
push_back
(
b
);
}
levels
.
push_back
(
level
);
}
buildTreeStructure
();
Real
*
d_s
=
c
.
getSize
();
for
(
UInt
i
=
0
;
i
<
3
;
++
i
)
domain_size
[
i
]
=
d_s
[
i
];
for
(
UInt
i
=
0
;
i
<
3
;
++
i
)
replica
[
i
]
=
0
;
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
buildTreeStructure
()
{
UInt
nlevels
=
levels
.
size
();
for
(
UInt
n
=
1
;
n
<
nlevels
;
++
n
)
{
Level
&
level
=
levels
[
n
];
Level
&
levelUp
=
levels
[
n
-
1
];
typename
myFMM
::
MyGrid
&
grid
=
level
.
fmm
->
grid
;
typename
myFMM
::
MyGrid
&
gridUp
=
levelUp
.
fmm
->
grid
;
typename
myFMM
::
MyBlockList
&
blocks
=
grid
.
getBlocks
();
for
(
UInt
i
=
0
;
i
<
blocks
.
size
();
++
i
)
{
Real
coords
[
Dim
];
grid
.
index2Coordinates
(
i
,
coords
);
Block
<
Dim
>
&
b
=
level
.
blocks
[
i
];
UInt
index
;
gridUp
.
coordinates2Index
(
coords
,
index
);
levelUp
.
blocks
[
index
].
addChildren
(
b
);
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::~
FMMMultigrid
()
{}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
template
<
typename
Cont
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
addPoints
(
Cont
&
positions
)
{
for
(
UInt
n
=
0
;
n
<
levels
.
size
();
++
n
)
{
levels
[
n
].
fmm
->
addPoints
(
positions
);
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
template
<
typename
Cont
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
computeMoments
(
Cont
&
cont
)
{
for
(
UInt
n
=
0
;
n
<
levels
.
size
();
++
n
)
{
levels
[
n
].
fmm
->
computeMoments
(
cont
);
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
template
<
typename
F
,
typename
Cont
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
computeConvolution
(
F
&
myFunc
,
Real
*
position
,
Cont
&
cont
,
Real
*
res
,
Real
rcut
)
{
for
(
UInt
i
=
0
;
i
<
FieldDim
;
++
i
)
res
[
i
]
=
0.0
;
for
(
UInt
l
=
0
;
l
<
levels
.
size
();
++
l
)
{
levels
[
l
].
fmm
->
computeConvolution
(
myFunc
,
position
,
cont
,
res
);
// for (int sx = -replica[0]; sx < replica[0]+1; ++sx) {
// for (int sy = -replica[1]; sy < replica[1]+1; ++sy) {
// for (int sz = -replica[2]; sz < replica[2]+1; ++sz) {
// Real local_pos[Dim];
// getReplicatedPosition(local_pos,position,sx,sy,sz);
// Real d =
// levels[0].fmm->grid.getSquaredDistanceToBlockCenter(local_pos,i);
// if (rcut > 0 && d > dist_block_limit2) continue;
// computeConvolutionRecursive(myFunc,local_pos,levels[0].blocks[i],cont,res);
// }
// }
// }
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
buildInteractionList
(
Real
rcut
)
{
Real
half_diagonal
=
levels
[
0
].
fmm
->
grid
.
getHalfDiagonal
();
Real
dist_block_limit
=
rcut
+
2
*
half_diagonal
;
Real
dist_block_limit2
=
dist_block_limit
*
dist_block_limit
;
Level
&
last_level
=
levels
.
back
();
Level
&
first_level
=
levels
[
0
];
for
(
UInt
i
=
0
;
i
<
last_level
.
blocks
.
size
();
++
i
)
{
Block
<
Dim
>
&
current_block
=
last_level
.
blocks
[
i
];
for
(
UInt
j
=
0
;
j
<
first_level
.
blocks
.
size
();
++
j
)
{
Block
<
Dim
>
&
far_block
=
first_level
.
blocks
[
j
];
Real
d
=
current_block
.
getSquaredDistanceFromBlock
(
far_block
);
if
(
rcut
>
0
&&
d
>
dist_block_limit2
)
continue
;
buildInteractionListRecursive
(
current_block
,
far_block
);
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
buildInteractionListRecursive
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
)
{
if
(
isFarAway
(
current_block
,
far_block
))
addMultipoleNeighbor
(
current_block
,
far_block
);
else
if
(
isChildLess
(
far_block
))
addCloseNeighbor
(
current_block
,
far_block
);
else
{
for
(
UInt
i
=
0
;
i
<
MaxNbChildren
<
Dim
>::
max_nb_children
;
++
i
)
{
buildInteractionListRecursive
(
current_block
,
*
far_block
.
children
[
i
]);
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
addMultipoleNeighbor
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
)
{
Level
&
level
=
levels
[
far_block
.
level
];
UInt
index
;
level
.
fmm
->
grid
.
coordinates2Index
(
current_block
.
coordinates
,
index
);
level
.
fmm
->
addMultipoleNeighbor
(
index
,
far_block
.
index
);
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
addCloseNeighbor
(
Block
<
Dim
>
&
current_block
,
Block
<
Dim
>
&
far_block
)
{
Level
&
level
=
levels
[
far_block
.
level
];
UInt
index
;
level
.
fmm
->
grid
.
coordinates2Index
(
current_block
.
coordinates
,
index
);
level
.
fmm
->
addCloseNeighbor
(
index
,
far_block
.
index
);
}
/* -------------------------------------------------------------------------- */
template
<
UInt
FieldDim
,
UInt
Dim
,
UInt
order
>
inline
void
FMMMultigrid
<
FieldDim
,
Dim
,
order
>::
getReplicatedPosition
(
Real
*
local_pos
,
Real
*
position
,
int
sx
,
int
sy
,
int
sz
)
{
int
s
[
3
]
=
{
sx
,
sy
,
sz
};
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
local_pos
[
i
]
=
position
[
i
]
+
s
[
i
]
*
domain_size
[
i
];
}
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif
/* __LIBMULTISCALE_FMM_MULTIGRID_HH__ */
Event Timeline
Log In to Comment