Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
solverstack
ScalFMM
Commits
7d30cb59
Commit
7d30cb59
authored
Jan 09, 2015
by
BLANCHARD Pierre
Browse files
Massive cleanup in interpolation based kernels (mostly reindentation and comment removal).
parent
7d114481
Changes
12
Expand all
Hide whitespace changes
Inline
Side-by-side
Src/Kernels/Chebyshev/FChebKernel.hpp
View file @
7d30cb59
...
...
@@ -40,51 +40,51 @@ class FTreeCoordinate;
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template
<
class
CellClass
,
class
ContainerClass
,
class
MatrixKernelClass
,
int
ORDER
,
int
NVALS
=
1
>
template
<
class
CellClass
,
class
ContainerClass
,
class
MatrixKernelClass
,
int
ORDER
,
int
NVALS
=
1
>
class
FChebKernel
:
public
FAbstractChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVALS
>
{
// private types
typedef
FChebM2LHandler
<
ORDER
,
MatrixKernelClass
>
M2LHandlerClass
;
// private types
typedef
FChebM2LHandler
<
ORDER
,
MatrixKernelClass
>
M2LHandlerClass
;
// using from
// using from
typedef
FAbstractChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVALS
>
AbstractBaseClass
;
AbstractBaseClass
;
/// Needed for P2P and M2L operators
const
MatrixKernelClass
*
const
MatrixKernel
;
/// Needed for P2P and M2L operators
const
MatrixKernelClass
*
const
MatrixKernel
;
/// Needed for M2L operator
FSmartPointer
<
M2LHandlerClass
,
FSmartPointerMemory
>
M2LHandler
;
/// Needed for M2L operator
FSmartPointer
<
M2LHandlerClass
,
FSmartPointerMemory
>
M2LHandler
;
public:
/**
* The constructor initializes all constant attributes and it reads the
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FChebKernel
(
const
int
inTreeHeight
,
const
FReal
inBoxWidth
,
const
FPoint
&
inBoxCenter
,
const
MatrixKernelClass
*
const
inMatrixKernel
,
const
FReal
Epsilon
)
/**
* The constructor initializes all constant attributes and it reads the
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FChebKernel
(
const
int
inTreeHeight
,
const
FReal
inBoxWidth
,
const
FPoint
&
inBoxCenter
,
const
MatrixKernelClass
*
const
inMatrixKernel
,
const
FReal
Epsilon
)
:
FAbstractChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVALS
>
(
inTreeHeight
,
inBoxWidth
,
inBoxCenter
),
MatrixKernel
(
inMatrixKernel
),
M2LHandler
(
new
M2LHandlerClass
(
MatrixKernel
,
Epsilon
))
{
// read precomputed compressed m2l operators from binary file
//M2LHandler->ReadFromBinaryFileAndSet();
M2LHandler
->
ComputeAndCompressAndSet
();
}
{
// read precomputed compressed m2l operators from binary file
//M2LHandler->ReadFromBinaryFileAndSet();
M2LHandler
->
ComputeAndCompressAndSet
();
}
FChebKernel
(
const
int
inTreeHeight
,
const
FReal
inBoxWidth
,
const
FPoint
&
inBoxCenter
,
const
MatrixKernelClass
*
const
inMatrixKernel
)
:
FChebKernel
(
inTreeHeight
,
inBoxWidth
,
inBoxCenter
,
inMatrixKernel
,
FMath
::
pow
(
10.0
,
static_cast
<
FReal
>
(
-
ORDER
)))
{
FChebKernel
(
const
int
inTreeHeight
,
const
FReal
inBoxWidth
,
const
FPoint
&
inBoxCenter
,
const
MatrixKernelClass
*
const
inMatrixKernel
)
:
FChebKernel
(
inTreeHeight
,
inBoxWidth
,
inBoxCenter
,
inMatrixKernel
,
FMath
::
pow
(
10.0
,
static_cast
<
FReal
>
(
-
ORDER
)))
{
}
}
void
P2M
(
CellClass
*
const
LeafCell
,
const
ContainerClass
*
const
SourceParticles
)
{
void
P2M
(
CellClass
*
const
LeafCell
,
const
ContainerClass
*
const
SourceParticles
)
{
const
FPoint
LeafCellCenter
(
AbstractBaseClass
::
getLeafCellCenter
(
LeafCell
->
getCoordinate
()));
for
(
int
idxRhs
=
0
;
idxRhs
<
NVALS
;
++
idxRhs
){
// 1) apply Sy
...
...
@@ -93,13 +93,13 @@ public:
// 2) apply B
M2LHandler
->
applyB
(
LeafCell
->
getMultipole
(
idxRhs
),
LeafCell
->
getMultipole
(
idxRhs
)
+
AbstractBaseClass
::
nnodes
);
}
}
}
void
M2M
(
CellClass
*
const
FRestrict
ParentCell
,
const
CellClass
*
const
FRestrict
*
const
FRestrict
ChildCells
,
const
int
/*TreeLevel*/
)
{
void
M2M
(
CellClass
*
const
FRestrict
ParentCell
,
const
CellClass
*
const
FRestrict
*
const
FRestrict
ChildCells
,
const
int
/*TreeLevel*/
)
{
for
(
int
idxRhs
=
0
;
idxRhs
<
NVALS
;
++
idxRhs
){
// 1) apply Sy
FBlas
::
scal
(
AbstractBaseClass
::
nnodes
*
2
,
FReal
(
0.
),
ParentCell
->
getMultipole
(
idxRhs
));
...
...
@@ -112,32 +112,32 @@ public:
// 2) apply B
M2LHandler
->
applyB
(
ParentCell
->
getMultipole
(
idxRhs
),
ParentCell
->
getMultipole
(
idxRhs
)
+
AbstractBaseClass
::
nnodes
);
}
}
//
void M2L(CellClass* const FRestrict TargetCell,
//
const CellClass* SourceCells[343],
//
const int NumSourceCells,
//
const int TreeLevel) const
//
{
//
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
//
const FTreeCoordinate& cx = TargetCell->getCoordinate();
//
for (int idx=0; idx<NumSourceCells; ++idx) {
//
const FTreeCoordinate& cy = SourceCells[idx]->getCoordinate();
//
const int transfer[3] = {cy.getX()-cx.getX(),
//
cy.getY()-cx.getY(),
//
cy.getZ()-cx.getZ()};
//
M2LHandler->applyC(transfer, CellWidth,
//
SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
//
TargetCell->getLocal() + AbstractBaseClass::nnodes);
//
}
//
}
void
M2L
(
CellClass
*
const
FRestrict
TargetCell
,
const
CellClass
*
SourceCells
[
343
],
const
int
/*NumSourceCells*/
,
const
int
TreeLevel
)
{
}
//
void M2L(CellClass* const FRestrict TargetCell,
//
const CellClass* SourceCells[343],
//
const int NumSourceCells,
//
const int TreeLevel) const
//
{
//
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
//
const FTreeCoordinate& cx = TargetCell->getCoordinate();
//
for (int idx=0; idx<NumSourceCells; ++idx) {
//
const FTreeCoordinate& cy = SourceCells[idx]->getCoordinate();
//
const int transfer[3] = {cy.getX()-cx.getX(),
//
cy.getY()-cx.getY(),
//
cy.getZ()-cx.getZ()};
//
M2LHandler->applyC(transfer, CellWidth,
//
SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
//
TargetCell->getLocal() + AbstractBaseClass::nnodes);
//
}
//
}
void
M2L
(
CellClass
*
const
FRestrict
TargetCell
,
const
CellClass
*
SourceCells
[
343
],
const
int
/*NumSourceCells*/
,
const
int
TreeLevel
)
{
for
(
int
idxRhs
=
0
;
idxRhs
<
NVALS
;
++
idxRhs
){
FReal
*
const
CompressedLocalExpansion
=
TargetCell
->
getLocal
(
idxRhs
)
+
AbstractBaseClass
::
nnodes
;
const
FReal
CellWidth
(
AbstractBaseClass
::
BoxWidth
/
FReal
(
FMath
::
pow
(
2
,
TreeLevel
)));
...
...
@@ -148,33 +148,33 @@ public:
}
}
}
}
//
void M2L(CellClass* const FRestrict TargetCell,
//
const CellClass* SourceCells[343],
//
const int NumSourceCells,
//
const int TreeLevel) const
//
{
//
const unsigned int rank = M2LHandler.getRank();
//
FBlas::scal(343*rank, FReal(0.), MultipoleExpansion);
//
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
//
for (int idx=0; idx<343; ++idx)
//
if (SourceCells[idx])
//
FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idx]->getMultipole())+AbstractBaseClass::nnodes,
//
MultipoleExpansion+idx*rank);
//
//
M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + AbstractBaseClass::nnodes);
//
}
void
L2L
(
const
CellClass
*
const
FRestrict
ParentCell
,
CellClass
*
FRestrict
*
const
FRestrict
ChildCells
,
const
int
/*TreeLevel*/
)
{
}
//
void M2L(CellClass* const FRestrict TargetCell,
//
const CellClass* SourceCells[343],
//
const int NumSourceCells,
//
const int TreeLevel) const
//
{
//
const unsigned int rank = M2LHandler.getRank();
//
FBlas::scal(343*rank, FReal(0.), MultipoleExpansion);
//
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
//
for (int idx=0; idx<343; ++idx)
//
if (SourceCells[idx])
//
FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idx]->getMultipole())+AbstractBaseClass::nnodes,
//
MultipoleExpansion+idx*rank);
//
//
M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + AbstractBaseClass::nnodes);
//
}
void
L2L
(
const
CellClass
*
const
FRestrict
ParentCell
,
CellClass
*
FRestrict
*
const
FRestrict
ChildCells
,
const
int
/*TreeLevel*/
)
{
for
(
int
idxRhs
=
0
;
idxRhs
<
NVALS
;
++
idxRhs
){
// 1) apply U
M2LHandler
->
applyU
(
ParentCell
->
getLocal
(
idxRhs
)
+
AbstractBaseClass
::
nnodes
,
const_cast
<
CellClass
*>
(
ParentCell
)
->
getLocal
(
idxRhs
));
const_cast
<
CellClass
*>
(
ParentCell
)
->
getLocal
(
idxRhs
));
// 2) apply Sx
for
(
unsigned
int
ChildIndex
=
0
;
ChildIndex
<
8
;
++
ChildIndex
){
if
(
ChildCells
[
ChildIndex
]){
...
...
@@ -182,11 +182,11 @@ public:
}
}
}
}
}
void
L2P
(
const
CellClass
*
const
LeafCell
,
ContainerClass
*
const
TargetParticles
)
{
void
L2P
(
const
CellClass
*
const
LeafCell
,
ContainerClass
*
const
TargetParticles
)
{
const
FPoint
LeafCellCenter
(
AbstractBaseClass
::
getLeafCellCenter
(
LeafCell
->
getCoordinate
()));
for
(
int
idxRhs
=
0
;
idxRhs
<
NVALS
;
++
idxRhs
){
...
...
@@ -195,26 +195,26 @@ public:
//// 2.a) apply Sx
//AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
//
AbstractBaseClass::BoxWidthLeaf,
//
LeafCell->getLocal(),
//
TargetParticles);
//
AbstractBaseClass::BoxWidthLeaf,
//
LeafCell->getLocal(),
//
TargetParticles);
//// 2.b) apply Px (grad Sx)
//AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
//
AbstractBaseClass::BoxWidthLeaf,
//
LeafCell->getLocal(),
//
TargetParticles);
//
AbstractBaseClass::BoxWidthLeaf,
//
LeafCell->getLocal(),
//
TargetParticles);
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass
::
Interpolator
->
applyL2PTotal
(
LeafCellCenter
,
AbstractBaseClass
::
BoxWidthLeaf
,
LeafCell
->
getLocal
(
idxRhs
),
TargetParticles
);
}
}
}
void
P2P
(
const
FTreeCoordinate
&
/* LeafCellCoordinate */
,
// needed for periodic boundary conditions
ContainerClass
*
const
FRestrict
TargetParticles
,
const
ContainerClass
*
const
FRestrict
/*SourceParticles*/
,
ContainerClass
*
const
NeighborSourceParticles
[
27
],
const
int
/* size */
)
ContainerClass
*
const
FRestrict
TargetParticles
,
const
ContainerClass
*
const
FRestrict
/*SourceParticles*/
,
ContainerClass
*
const
NeighborSourceParticles
[
27
],
const
int
/* size */
)
{
DirectInteractionComputer
<
MatrixKernelClass
::
NCMP
,
NVALS
>::
P2P
(
TargetParticles
,
NeighborSourceParticles
,
MatrixKernel
);
}
...
...
@@ -222,7 +222,8 @@ public:
void
P2PRemote
(
const
FTreeCoordinate
&
/*inPosition*/
,
ContainerClass
*
const
FRestrict
inTargets
,
const
ContainerClass
*
const
FRestrict
/*inSources*/
,
ContainerClass
*
const
inNeighbors
[
27
],
const
int
/*inSize*/
){
ContainerClass
*
const
inNeighbors
[
27
],
const
int
/*inSize*/
)
{
DirectInteractionComputer
<
MatrixKernelClass
::
NCMP
,
NVALS
>::
P2PRemote
(
inTargets
,
inNeighbors
,
27
,
MatrixKernel
);
}
...
...
Src/Kernels/Chebyshev/FChebSymKernel.hpp
View file @
7d30cb59
This diff is collapsed.
Click to expand it.
Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp
View file @
7d30cb59
This diff is collapsed.
Click to expand it.
Src/Kernels/Chebyshev/FChebTensorialKernel.hpp
View file @
7d30cb59
...
...
@@ -41,204 +41,193 @@ class FTreeCoordinate;
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template
<
class
CellClass
,
class
ContainerClass
,
class
MatrixKernelClass
,
int
ORDER
,
int
NVALS
=
1
>
template
<
class
CellClass
,
class
ContainerClass
,
class
MatrixKernelClass
,
int
ORDER
,
int
NVALS
=
1
>
class
FChebTensorialKernel
:
public
FAbstractChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVALS
>
{
enum
{
nRhs
=
MatrixKernelClass
::
NRHS
,
nLhs
=
MatrixKernelClass
::
NLHS
,
nPot
=
MatrixKernelClass
::
NPOT
,
nPv
=
MatrixKernelClass
::
NPV
};
enum
{
nRhs
=
MatrixKernelClass
::
NRHS
,
nLhs
=
MatrixKernelClass
::
NLHS
,
nPot
=
MatrixKernelClass
::
NPOT
,
nPv
=
MatrixKernelClass
::
NPV
};
protected:
//PB: for OptiDis
// private types
typedef
FChebTensorialM2LHandler
<
ORDER
,
MatrixKernelClass
,
MatrixKernelClass
::
Type
>
M2LHandlerClass
;
// private types
typedef
FChebTensorialM2LHandler
<
ORDER
,
MatrixKernelClass
,
MatrixKernelClass
::
Type
>
M2LHandlerClass
;
// using from
// using from
typedef
FAbstractChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVALS
>
AbstractBaseClass
;
AbstractBaseClass
;
/// Needed for P2P and M2L operators
const
MatrixKernelClass
*
const
MatrixKernel
;
/// Needed for P2P and M2L operators
const
MatrixKernelClass
*
const
MatrixKernel
;
/// Needed for M2L operator
FSmartPointer
<
M2LHandlerClass
,
FSmartPointerMemory
>
M2LHandler
;
/// Needed for M2L operator
FSmartPointer
<
M2LHandlerClass
,
FSmartPointerMemory
>
M2LHandler
;
public:
/**
* The constructor initializes all constant attributes and it reads the
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FChebTensorialKernel
(
const
int
inTreeHeight
,
const
FReal
inBoxWidth
,
const
FPoint
&
inBoxCenter
,
const
MatrixKernelClass
*
const
inMatrixKernel
,
const
FReal
inBoxWidthExtension
,
const
FReal
Epsilon
)
/**
* The constructor initializes all constant attributes and it reads the
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FChebTensorialKernel
(
const
int
inTreeHeight
,
const
FReal
inBoxWidth
,
const
FPoint
&
inBoxCenter
,
const
MatrixKernelClass
*
const
inMatrixKernel
,
const
FReal
inBoxWidthExtension
,
const
FReal
Epsilon
)
:
FAbstractChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVALS
>
(
inTreeHeight
,
inBoxWidth
,
inBoxCenter
,
inBoxWidthExtension
),
MatrixKernel
(
inMatrixKernel
),
M2LHandler
(
new
M2LHandlerClass
(
MatrixKernel
,
inTreeHeight
,
inBoxWidth
,
inBoxWidthExtension
,
Epsilon
))
{
}
M2LHandler
(
new
M2LHandlerClass
(
MatrixKernel
,
inTreeHeight
,
inBoxWidth
,
inBoxWidthExtension
,
Epsilon
))
{
}
void
P2M
(
CellClass
*
const
LeafCell
,
const
ContainerClass
*
const
SourceParticles
)
{
const
FPoint
LeafCellCenter
(
AbstractBaseClass
::
getLeafCellCenter
(
LeafCell
->
getCoordinate
()));
const
FReal
ExtendedLeafCellWidth
(
AbstractBaseClass
::
BoxWidthLeaf
+
AbstractBaseClass
::
BoxWidthExtension
);
void
P2M
(
CellClass
*
const
LeafCell
,
const
ContainerClass
*
const
SourceParticles
)
{
const
FPoint
LeafCellCenter
(
AbstractBaseClass
::
getLeafCellCenter
(
LeafCell
->
getCoordinate
()));
const
FReal
ExtendedLeafCellWidth
(
AbstractBaseClass
::
BoxWidthLeaf
+
AbstractBaseClass
::
BoxWidthExtension
);
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
// 1) apply Sy
AbstractBaseClass
::
Interpolator
->
applyP2M
(
LeafCellCenter
,
ExtendedLeafCellWidth
,
LeafCell
->
getMultipole
(
idxV
*
nRhs
),
SourceParticles
);
// 1) apply Sy
AbstractBaseClass
::
Interpolator
->
applyP2M
(
LeafCellCenter
,
ExtendedLeafCellWidth
,
LeafCell
->
getMultipole
(
idxV
*
nRhs
),
SourceParticles
);
for
(
int
idxRhs
=
0
;
idxRhs
<
nRhs
;
++
idxRhs
){
// update multipole index
int
idxMul
=
idxV
*
nRhs
+
idxRhs
;
for
(
int
idxRhs
=
0
;
idxRhs
<
nRhs
;
++
idxRhs
){
// update multipole index
int
idxMul
=
idxV
*
nRhs
+
idxRhs
;
// 2) apply B (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyB
(
LeafCell
->
getMultipole
(
idxMul
),
LeafCell
->
getMultipole
(
idxMul
)
+
AbstractBaseClass
::
nnodes
);
// 2) apply B (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyB
(
LeafCell
->
getMultipole
(
idxMul
),
LeafCell
->
getMultipole
(
idxMul
)
+
AbstractBaseClass
::
nnodes
);
}
}
}
}
}
void
M2M
(
CellClass
*
const
FRestrict
ParentCell
,
const
CellClass
*
const
FRestrict
*
const
FRestrict
ChildCells
,
const
int
TreeLevel
)
{
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxRhs
=
0
;
idxRhs
<
nRhs
;
++
idxRhs
){
// update multipole index
int
idxMul
=
idxV
*
nRhs
+
idxRhs
;
// 1) apply Sy
FBlas
::
scal
(
AbstractBaseClass
::
nnodes
*
2
,
FReal
(
0.
),
ParentCell
->
getMultipole
(
idxMul
));
for
(
unsigned
int
ChildIndex
=
0
;
ChildIndex
<
8
;
++
ChildIndex
){
if
(
ChildCells
[
ChildIndex
]){
AbstractBaseClass
::
Interpolator
->
applyM2M
(
ChildIndex
,
ChildCells
[
ChildIndex
]
->
getMultipole
(
idxMul
),
ParentCell
->
getMultipole
(
idxMul
),
TreeLevel
/*Cell width extension specific*/
);
}
void
M2M
(
CellClass
*
const
FRestrict
ParentCell
,
const
CellClass
*
const
FRestrict
*
const
FRestrict
ChildCells
,
const
int
TreeLevel
)
{
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxRhs
=
0
;
idxRhs
<
nRhs
;
++
idxRhs
){
// update multipole index
int
idxMul
=
idxV
*
nRhs
+
idxRhs
;
// 1) apply Sy
FBlas
::
scal
(
AbstractBaseClass
::
nnodes
*
2
,
FReal
(
0.
),
ParentCell
->
getMultipole
(
idxMul
));
for
(
unsigned
int
ChildIndex
=
0
;
ChildIndex
<
8
;
++
ChildIndex
){
if
(
ChildCells
[
ChildIndex
]){
AbstractBaseClass
::
Interpolator
->
applyM2M
(
ChildIndex
,
ChildCells
[
ChildIndex
]
->
getMultipole
(
idxMul
),
ParentCell
->
getMultipole
(
idxMul
),
TreeLevel
/*Cell width extension specific*/
);
}
}
// 2) apply B (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyB
(
ParentCell
->
getMultipole
(
idxMul
),
ParentCell
->
getMultipole
(
idxMul
)
+
AbstractBaseClass
::
nnodes
);
}
}
}
// 2) apply B (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyB
(
ParentCell
->
getMultipole
(
idxMul
),
ParentCell
->
getMultipole
(
idxMul
)
+
AbstractBaseClass
::
nnodes
);
void
M2L
(
CellClass
*
const
FRestrict
TargetCell
,
const
CellClass
*
SourceCells
[
343
],
const
int
/*NumSourceCells*/
,
const
int
TreeLevel
)
{
// scale factor for homogeneous case
const
FReal
CellWidth
(
AbstractBaseClass
::
BoxWidth
/
FReal
(
FMath
::
pow
(
2
,
TreeLevel
)));
const
FReal
ExtendedCellWidth
(
CellWidth
+
AbstractBaseClass
::
BoxWidthExtension
);
const
FReal
scale
(
MatrixKernel
->
getScaleFactor
(
ExtendedCellWidth
));
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxLhs
=
0
;
idxLhs
<
nLhs
;
++
idxLhs
){
// update local index
const
int
idxLoc
=
idxV
*
nLhs
+
idxLhs
;
FReal
*
const
CompressedLocalExpansion
=
TargetCell
->
getLocal
(
idxLoc
)
+
AbstractBaseClass
::
nnodes
;
// update idxRhs
const
int
idxRhs
=
idxLhs
%
nPv
;
// update multipole index
const
int
idxMul
=
idxV
*
nRhs
+
idxRhs
;
// get index in matrix kernel
const
unsigned
int
d
=
MatrixKernel
->
getPosition
(
idxLhs
);
for
(
int
idx
=
0
;
idx
<
343
;
++
idx
){
if
(
SourceCells
[
idx
]){
M2LHandler
->
applyC
(
idx
,
TreeLevel
,
scale
,
d
,
SourceCells
[
idx
]
->
getMultipole
(
idxMul
)
+
AbstractBaseClass
::
nnodes
,
CompressedLocalExpansion
);
}
}
}
// NLHS=NPOT*NPV
}
// NVALS
}
}
}
void
M2L
(
CellClass
*
const
FRestrict
TargetCell
,
const
CellClass
*
SourceCells
[
343
],
const
int
/*NumSourceCells*/
,
const
int
TreeLevel
)
{
// scale factor for homogeneous case
const
FReal
CellWidth
(
AbstractBaseClass
::
BoxWidth
/
FReal
(
FMath
::
pow
(
2
,
TreeLevel
)));
const
FReal
ExtendedCellWidth
(
CellWidth
+
AbstractBaseClass
::
BoxWidthExtension
);
const
FReal
scale
(
MatrixKernel
->
getScaleFactor
(
ExtendedCellWidth
));
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxLhs
=
0
;
idxLhs
<
nLhs
;
++
idxLhs
){
// update local index
const
int
idxLoc
=
idxV
*
nLhs
+
idxLhs
;
FReal
*
const
CompressedLocalExpansion
=
TargetCell
->
getLocal
(
idxLoc
)
+
AbstractBaseClass
::
nnodes
;
// update idxRhs
const
int
idxRhs
=
idxLhs
%
nPv
;
// update multipole index
const
int
idxMul
=
idxV
*
nRhs
+
idxRhs
;
// get index in matrix kernel
const
unsigned
int
d
=
MatrixKernel
->
getPosition
(
idxLhs
);
for
(
int
idx
=
0
;
idx
<
343
;
++
idx
){
if
(
SourceCells
[
idx
]){
M2LHandler
->
applyC
(
idx
,
TreeLevel
,
scale
,
d
,
SourceCells
[
idx
]
->
getMultipole
(
idxMul
)
+
AbstractBaseClass
::
nnodes
,
CompressedLocalExpansion
);
}
}
}
// NLHS=NPOT*NPV
}
// NVALS
}
void
L2L
(
const
CellClass
*
const
FRestrict
ParentCell
,
CellClass
*
FRestrict
*
const
FRestrict
ChildCells
,
void
L2L
(
const
CellClass
*
const
FRestrict
ParentCell
,
CellClass
*
FRestrict
*
const
FRestrict
ChildCells
,
const
int
TreeLevel
)
{
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxLhs
=
0
;
idxLhs
<
nLhs
;
++
idxLhs
){
int
idxLoc
=
idxV
*
nLhs
+
idxLhs
;
// 1) apply U (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyU
(
ParentCell
->
getLocal
(
idxLoc
)
+
AbstractBaseClass
::
nnodes
,
const_cast
<
CellClass
*>
(
ParentCell
)
->
getLocal
(
idxLoc
));
// 2) apply Sx
for
(
unsigned
int
ChildIndex
=
0
;
ChildIndex
<
8
;
++
ChildIndex
){
if
(
ChildCells
[
ChildIndex
]){
AbstractBaseClass
::
Interpolator
->
applyL2L
(
ChildIndex
,
ParentCell
->
getLocal
(
idxLoc
),
ChildCells
[
ChildIndex
]
->
getLocal
(
idxLoc
),
TreeLevel
/*Cell width extension specific*/
);
}
{
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxLhs
=
0
;
idxLhs
<
nLhs
;
++
idxLhs
){
int
idxLoc
=
idxV
*
nLhs
+
idxLhs
;
// 1) apply U (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyU
(
ParentCell
->
getLocal
(
idxLoc
)
+
AbstractBaseClass
::
nnodes
,
const_cast
<
CellClass
*>
(
ParentCell
)
->
getLocal
(
idxLoc
));
// 2) apply Sx
for
(
unsigned
int
ChildIndex
=
0
;
ChildIndex
<
8
;
++
ChildIndex
){
if
(
ChildCells
[
ChildIndex
]){
AbstractBaseClass
::
Interpolator
->
applyL2L
(
ChildIndex
,
ParentCell
->
getLocal
(
idxLoc
),
ChildCells
[
ChildIndex
]
->
getLocal
(
idxLoc
),
TreeLevel
/*Cell width extension specific*/
);
}
}
}
//NLHS
}
// NVALS
}
void
L2P
(
const
CellClass
*
const
LeafCell
,
ContainerClass
*
const
TargetParticles
)
{
const
FPoint
LeafCellCenter
(
AbstractBaseClass
::
getLeafCellCenter
(
LeafCell
->
getCoordinate
()));
const
FReal
ExtendedLeafCellWidth
(
AbstractBaseClass
::
BoxWidthLeaf
+
AbstractBaseClass
::
BoxWidthExtension
);
for
(
int
idxV
=
0
;
idxV
<
NVALS
;
++
idxV
){
for
(
int
idxLhs
=
0
;
idxLhs
<
nLhs
;
++
idxLhs
){
int
idxLoc
=
idxV
*
nLhs
+
idxLhs
;
// 1) apply U (PB: Tensorial version is just a basic copy)
M2LHandler
->
applyU
(
LeafCell
->
getLocal
(
idxLoc
)
+
AbstractBaseClass
::
nnodes
,
const_cast
<
CellClass
*>
(
LeafCell
)
->
getLocal
(
idxLoc
));
}
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass
::
Interpolator
->
applyL2PTotal
(
LeafCellCenter
,
ExtendedLeafCellWidth
,
LeafCell
->
getLocal
(
idxV
*
nLhs
),
TargetParticles
);