Commit 096352d8 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Provided error computation for the AdaptUnifFMM test, provided a dense version...

Provided error computation for the AdaptUnifFMM test, provided a dense version for the UnifKernel, some fixes in AdaptUnifKernel.
parent a46d7dcc
8 8
100 20 0 0 0
0.998298 0.0088693 0.576438 0.396465 14.29748155 -0.7594320072 -0.4684170216 1.050760367
0.107422 -0.2487 -9.62605 0.840485 14.71623438 5.42210735 0.01697689024 22.63840185
-0.784533 0.0389288 6.18863 0.353336 14.01625814 1.342793154 2.541019898 1.023328623
-0.239929 0.958705 1.52707 0.446583 13.82038644 0.3427293794 -0.9177223361 -0.1173689885
0.403898 0.780385 -4.77353 0.318693 18.02659746 -1.8519823 -0.7599168368 -0.4212908058
0.0229265 0.413024 -9.10432 0.886428 12.84480354 -1.476993165 -4.408963236 0.8039101696
0.161035 0.900809 4.03249 0.0155828 17.91858982 0.3674203458 -0.135490968 -0.001703902669
-0.391436 0.31133 -8.6594 0.58409 12.79118476 2.110249652 -0.8154242261 -1.340130298
0.336155 0.892897 -2.99556 0.159369 16.51879858 -0.3413313918 -0.4038742197 -0.1695737321
0.909534 0.318221 2.67363 0.383716 13.63236151 -0.6984426188 -0.08453458298 -0.03184121651
-0.323809 0.118659 -9.38652 0.691004 13.65578397 4.13754236 -3.278345133 2.212893705
0.423993 -0.27039 8.6436 0.0588589 21.37140431 -3.965824309 -5.856611834 -3.175616708
-0.800443 0.592642 -0.898189 0.899854 15.35737508 1.492493525 -1.837674921 -2.119327481
0.487848 -0.782409 3.87093 0.163546 13.25532837 -0.09611802987 0.3182714101 0.082792552
-0.4158 -0.85766 3.02539 0.159072 13.89967484 0.8726699513 0.02768148804 -0.4731774218
0.818692 -0.351912 -4.53763 0.533065 18.53139958 -1.739883317 1.353975794 2.351712185
-0.510301 0.495726 -7.02744 0.604144 12.68431756 0.6622832741 -0.5530973005 1.225103725
0.171849 -0.265669 -9.48624 0.582699 16.70097877 -10.85750933 4.032661984 -16.16935553
-0.410242 0.878534 -2.44705 0.269971 16.36315255 0.1533852258 -0.8053696682 -0.07870660632
-0.00788981 0.954825 2.97065 0.390478 13.80729829 0.6353608641 -0.7302390407 0.6531532648
0.812492 -0.302815 -4.98156 0.293401 17.42390354 -0.9688379933 0.2534905912 1.484205495
-0.846722 -0.267532 -4.59878 0.742377 17.49791845 3.168836804 1.414376324 -0.3656350602
-0.197639 0.654068 7.3016 0.298526 10.97730809 -0.01347553415 -0.5272954632 -0.5457414919
0.531595 -0.522932 6.66295 0.0755381 11.6527424 -0.1081160051 0.08949230156 -0.1083691593
0.418405 -0.309864 -8.53769 0.404983 12.95183706 -0.8191691201 1.590530126 -0.329490809
0.532116 -0.816875 2.22637 0.857378 12.83028575 -1.032696155 1.359460738 0.1149674386
0.0428093 -0.966039 -2.54826 0.941968 16.26675069 -0.3883253602 2.73825496 -0.5316220517
-0.339692 -0.886642 -3.13806 0.662831 16.70765386 1.413716907 1.563798562 0.5779957608
0.596096 -0.158221 7.87169 0.846476 10.26675944 -1.546348649 -0.2492048396 0.7548470147
-0.835556 -0.334092 -4.36152 0.00275508 20.26450745 0.009373692973 0.0136275533 -0.03438362594
-0.730973 0.25523 6.32879 0.462379 14.12117451 0.7259766969 -1.514796122 -0.2433358602
-0.217175 0.974508 -0.56292 0.532596 14.99497744 -0.3727025802 -1.342466841 -0.6883943114
0.782485 0.605231 1.46331 0.787877 17.4939615 4.310470707 -7.543446712 -25.35161807
0.380726 -0.791496 4.78103 0.265612 14.76978812 1.286129406 1.18308357 -2.37147925
-0.0353821 0.866432 -4.9804 0.982752 17.35316672 -1.331900559 -4.291978616 0.5803187573
-0.960368 -0.166482 -2.23557 0.306785 18.02336054 0.7787297365 0.5313155548 1.411507595
0.675061 -0.710818 -1.97559 0.600855 16.33900066 -1.58879407 1.319930849 -0.4905025247
-0.717732 0.562193 -4.10852 0.608716 17.80111434 1.775954064 -1.166262563 -1.242916712
0.407635 -0.831071 -3.78358 0.212439 18.94087426 0.2333054121 1.045851485 0.1422702839
-0.59752 0.611465 -5.1873 0.885895 17.87404217 2.330263656 -2.427078395 5.093755376
0.808438 -0.582506 -0.843452 0.304657 14.7618041 -0.6371379514 0.3717871183 -0.601245541
0.974247 -0.182917 1.31847 0.15186 14.532023 -0.305811479 0.370761115 -0.07421467996
0.426876 0.893875 -1.36987 0.337662 15.95289331 -0.09672938628 -1.156072134 -0.6485779912
-0.906895 0.267389 -3.25646 0.387477 16.80465821 1.171615313 -0.1310559689 -0.2915520221
0.859915 0.35397 3.67765 0.64361 13.56155883 -1.59483203 0.398439636 0.3234847415
-0.106217 0.923875 -3.67658 0.753553 16.91246203 -1.570139229 -2.650253579 0.4046381132
0.456599 0.645136 6.12631 0.603616 11.78750398 -0.776873399 -0.9251797086 -0.9066983702
0.487819 -0.741964 4.59915 0.531628 13.52424296 -1.866601516 0.3345897487 2.451865114
-0.263461 -0.211377 -9.41227 0.45936 14.86271702 2.884619993 3.212359652 1.203027081
-0.804465 0.564279 1.85539 0.652488 12.82868631 1.49687237 -0.1729054568 -0.453300195
-0.350825 0.87264 -3.39738 0.327181 17.98262156 1.189545052 -0.7542418926 -1.556327738
0.719278 0.513161 4.68299 0.94637 13.40021267 -2.510202382 0.2645304858 0.1269231865
0.976244 -0.0912485 -1.96521 0.36804 16.8140201 -1.33151783 0.4012758142 0.3814902231
0.381419 -0.336697 8.60903 0.94389 11.491142 -1.896896355 4.229064696 0.515105202
0.612962 -0.771412 -1.70885 0.00742826 18.28645281 -0.009292540661 0.02696330591 -0.05877513813
-0.971677 0.153174 -1.7995 0.5166 21.38240653 0.4159665734 -8.711526195 10.5521484
0.975193 0.145973 -1.66402 0.272771 17.12632992 -1.120347369 0.1961395631 -0.9029540694
0.0789429 0.896583 4.35783 0.0242992 15.1898529 0.1104537747 -0.06683186503 -0.09279405087
0.0663075 -0.825902 5.59901 0.591955 16.58254059 -9.507211695 2.04394656 17.51168547
0.86939 -0.420831 -2.58964 0.204964 16.87327953 -0.6732864973 0.08038422427 0.22127318
0.902358 -0.0929655 -4.20841 0.877693 17.33314226 -4.043922988 -2.988626443 -1.448524729
-0.310971 -0.220975 -9.24374 0.0593687 16.68404051 0.5356370279 0.4450973984 -0.9499403671
0.58475 0.628044 5.13448 0.260843 15.43499033 -1.676985795 0.6949679164 -1.702324214
-0.728461 -0.387442 -5.65008 0.302829 16.36104328 0.8806147341 0.4319300274 0.9685086786
-0.00839011 -0.818783 5.74042 0.891495 14.5722321 9.706793087 0.8191228752 -18.98136498
0.531829 -0.729876 -4.29463 0.498198 18.29054245 -0.1254323436 2.584365766 0.9601796901
-0.672177 0.315034 6.70024 0.710026 12.10878463 0.6263313778 -0.8626211707 -3.771042809
0.00306144 0.432259 9.01744 0.286414 10.14090971 -0.7109008432 -1.062711173 -0.7829087055
-0.986248 0.0182631 -1.64259 0.864924 22.10581221 2.868524976 7.856483617 15.09145536
-0.391103 -0.911636 1.26331 0.675541 12.60575989 0.6921027527 1.061084013 0.03245341107
-0.703262 -0.319161 -6.35262 0.45849 14.08893904 0.8311435878 0.5870437172 1.330141
0.441326 0.674298 -5.92076 0.959636 16.62361545 -11.12542316 2.465246918 2.661836817
0.194108 0.782199 -5.92019 0.774675 17.52788403 8.229836987 -5.597290683 2.421260513
0.643091 -0.457905 -6.13805 0.376551 15.55764126 -0.5671567363 0.4313177085 5.08606665
-0.611491 0.522452 5.94241 0.228639 13.74001873 -0.02446148821 -1.050253869 0.09030808573
-0.0661106 -0.832189 -5.50537 0.354534 16.23699299 -0.4674788047 1.296576191 1.237277694
0.281414 0.901252 3.29472 0.300318 14.56693558 -0.3879021755 -0.5210172679 0.07814535005
0.811012 0.569713 1.32994 0.669766 18.37889061 -6.27186783 5.938290788 24.85531156
0.472628 0.194104 -8.5962 0.718967 12.53467115 -2.44115358 -1.261721354 -0.7823551363
0.734747 -0.557898 -3.85871 0.565955 18.10639504 -2.127400743 0.9506388968 -1.782062777
-0.402198 0.190946 8.9542 0.824465 9.097815055 2.209256028 -0.6150680949 -1.645466648
-0.772446 0.28541 5.67334 0.390612 12.97935297 1.055492666 -0.06714733981 1.082092007
0.33774 0.850235 4.03773 0.818766 13.70786857 -0.1291249713 -1.71603002 -0.2464575282
-0.472288 -0.72322 -5.03881 0.844008 16.75513341 1.276065346 3.040439577 1.891969281
0.411452 -0.839687 -3.54449 0.180468 18.52831204 -0.1250614432 0.6678103838 -0.9626920093
-0.752738 0.453448 -4.77254 0.943396 18.19673518 4.370174502 -0.5349894537 -1.416239708
0.408807 0.758538 5.07442 0.424887 14.51369479 2.068019917 -2.367422636 -0.6850077683
-0.338412 -0.8306 -4.42246 0.520666 17.5345359 0.4623101451 2.168896619 -0.6770585867
0.740683 0.233553 6.29953 0.0656438 12.48104888 -0.1927865055 0.07022899581 -0.1313581529
-0.00108433 -0.990426 -1.38042 0.913508 15.12869223 -0.1434366029 2.698202115 -1.821630348
0.787219 0.591184 -1.75465 0.882585 15.49192788 -1.821806551 -2.423193199 -0.01620785149
-0.0907875 -0.954644 2.83572 0.761364 12.73059659 0.2559030184 1.519572696 -0.3520593696
-0.0664647 -0.59599 -8.00236 0.398923 12.01965614 0.2073711032 0.8317877456 -0.1493446859
-0.227157 0.973079 0.389491 0.688257 13.71042819 0.3854590296 -1.010951777 -0.4948542166
-0.988393 0.0354439 -1.47727 0.761548 21.57235787 2.585745461 -1.972885298 -27.46026958
-0.296175 0.882092 3.66325 0.405009 13.48236105 1.334825307 -0.5179042414 -0.1181481579
0.986697 -0.145634 0.722434 0.125251 15.55842828 -0.1539269267 0.9463824461 -0.6331477821
-0.858933 -0.485494 -1.62882 0.484634 17.89352969 0.4506722447 3.58065165 -0.4724275055
0.650514 -0.466089 -5.9966 0.222463 16.97435046 -0.7183339394 0.651751984 -3.564899372
0.0480048 -0.51252 8.57332 0.873121 11.03326509 5.236181041 4.194179655 -0.6447565705
8 8
20 10 0 0 0
-0.381596 -0.909478 -1.65028 7.82637e-06 2.541938222 1.557606994e-06 3.489201182e-06 1.9909474e-07
-0.328097 -0.475796 -8.1607 0.131538 1.119259139 0.0008908943508 0.004579383445 0.02555980285
0.773946 -0.488975 4.02381 0.755605 3.011973162 -1.45016502 -1.652190245 -2.411831961
0.841093 0.258549 -4.75094 0.45865 1.98990251 -0.2577510645 0.03655161637 0.1546151569
-0.631449 -0.695168 -3.43532 0.532767 2.15223708 0.07790438363 0.2755104636 -0.006342371017
0.855311 0.420252 3.03037 0.218959 2.376001959 -0.03848747137 -0.1167391388 0.07418806718
-0.0908349 0.632258 7.69414 0.0470446 5.012127274 0.9218738946 -0.04010423188 0.4409636235
-0.837853 0.118 5.3299 0.678865 1.634756411 0.1767291196 -0.03509828097 -0.1226211785
-0.565994 0.690639 -4.50186 0.679296 2.910213541 -1.257188682 -1.850335729 -0.5607101455
0.799983 -0.599417 -0.269652 0.934693 1.907152543 -0.2129228855 0.2885587276 -0.02844632625
-0.42005 0.906916 0.325839 0.383502 2.89458325 0.434894857 0.007264147335 -0.595471763
-0.101258 0.994787 -0.120466 0.519416 2.798499163 -0.3339421311 -0.2683887681 0.3091350594
-0.164731 0.958734 -2.31716 0.830965 2.092163739 -0.003599627918 -0.1920184011 0.01079716782
-0.973424 0.221148 0.594948 0.0345721 2.634378832 0.02483866798 0.01471930828 -0.01397072505
0.532694 0.535145 6.55635 0.0534616 2.001048773 -0.01524642715 -0.004049037322 0.001833595101
0.626587 -0.682149 3.76911 0.5297 3.660159952 1.30712139 1.873943426 2.287116347
0.0663265 0.625817 7.77144 0.671149 1.431005579 -0.9270508856 0.007352727933 -0.6234777785
-0.951257 -0.0795204 -2.9797 0.00769819 2.98995399 0.005367080118 -0.001040743771 -0.003162173107
-0.764777 0.449732 -4.61365 0.383416 3.753458481 1.512786443 1.666239546 0.9383888638
-0.23121 0.970786 -0.641653 0.0668422 3.441324218 0.03394590788 -0.01475826149 0.123436539
......@@ -22,7 +22,8 @@
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
//#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifDenseKernel.hpp"
#include "Kernels/Uniform/FUnifM2LHandler.hpp"
class FTreeCoordinate;
......@@ -54,10 +55,10 @@ class FTreeCoordinate;
*/
template< class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAdaptiveUnifKernel : FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
class FAdaptiveUnifKernel : FUnifDenseKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
typedef FUnifDenseKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
enum {order = ORDER,
nnodes = KernelBaseClass::nnodes};
/// Needed for M2L operator
......@@ -97,7 +98,7 @@ public:
void P2M(CellClass* const pole, const int cellLevel, const ContainerClass* const particles) override {
//pole->setDataUp(pole->getDataUp() + particles->getNbParticles());
//
const FPoint CellCenter(KernelBaseClass::getLeafCellCenter(pole->getCoordinate()));
const FPoint CellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),cellLevel));
const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-cellLevel);
//
......@@ -115,63 +116,94 @@ public:
void M2M(CellClass* const pole, const int poleLevel, const CellClass* const subCell, const int subCellLevel) override {
// pole->setDataUp(pole->getDataUp() + subCell->getDataUp());
const FPoint subCellCenter(KernelBaseClass::getLeafCellCenter(subCell->getCoordinate()));
const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel));
const FReal subCellWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-subCellLevel);
const FPoint poleCellCenter(KernelBaseClass::getLeafCellCenter(pole->getCoordinate()));
const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
const FReal poleCellWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-poleLevel);
// Set sub-child coords
FReal globalChildCoords[3][ORDER];
FUnifTensor<order>::setPolynomialsRoots(subCellCenter, subCellWidth, globalChildCoords);
////////////////////////////////////////////////////////////////////////////
/// p^6 version
// allocate memory
FReal* subChildParentInterpolator = new FReal [nnodes * nnodes];
// set child info
FPoint ChildRoots[nnodes], localChildRoots[nnodes];
FUnifTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots);
// Map global position of sub-child nodes to [-1,1]
FReal localChildCoords[3][ORDER];
// map global position of roots to local position in parent cell
const map_glob_loc map(poleCellCenter, poleCellWidth);
FPoint localChildPoints;
for (unsigned int n=0; n<ORDER; ++n) {
map(FPoint(globalChildCoords[0][n],globalChildCoords[1][n],globalChildCoords[2][n]), localChildPoints);
localChildCoords[0][n] = localChildPoints.getX();
localChildCoords[1][n] = localChildPoints.getY();
localChildCoords[2][n] = localChildPoints.getZ();
}
for (unsigned int n=0; n<nnodes; ++n)
map(ChildRoots[n], localChildRoots[n]);
// assemble interpolator
FReal* subChildParentInterpolator = new FReal [3 * ORDER*ORDER];
KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[0], subChildParentInterpolator);
KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[1], subChildParentInterpolator + 1 * ORDER*ORDER);
KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[2], subChildParentInterpolator + 2 * ORDER*ORDER);
// assemble child - parent - interpolator
KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator);
// get permutation operators
unsigned int perm[3][nnodes];
for (unsigned int i=0;i<3; ++i)
for (unsigned int n=0; n<nnodes; ++n)
perm[i][n] = KernelBaseClass::Interpolator->getPermutationsM2ML2L(i)[n];
// ////////////////////////////////////////////////////////////////////////////
// /// p^4 version
//
// // Set sub-child coords
// FReal globalChildCoords[3][ORDER];
// FUnifTensor<order>::setPolynomialsRoots(subCellCenter, subCellWidth, globalChildCoords);
//
// // Map global position of sub-child nodes to [-1,1]
// FReal localChildCoords[3][ORDER];
// const map_glob_loc map(poleCellCenter, poleCellWidth);
// FPoint localChildPoints;
// for (unsigned int n=0; n<ORDER; ++n) {
// map(FPoint(globalChildCoords[0][n],globalChildCoords[1][n],globalChildCoords[2][n]), localChildPoints);
// localChildCoords[0][n] = localChildPoints.getX();
// localChildCoords[1][n] = localChildPoints.getY();
// localChildCoords[2][n] = localChildPoints.getZ();
// }
//
// // assemble interpolator
// FReal* subChildParentInterpolator = new FReal [3 * ORDER*ORDER];
// KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[0], subChildParentInterpolator);
// KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[1], subChildParentInterpolator + 1 * ORDER*ORDER);
// KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[2], subChildParentInterpolator + 2 * ORDER*ORDER);
//
// // get permutation operators
// unsigned int perm[3][nnodes];
// for (unsigned int i=0;i<3; ++i)
// for (unsigned int n=0; n<nnodes; ++n)
// perm[i][n] = KernelBaseClass::Interpolator->getPermutationsM2ML2L(i)[n];
//
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy (using tensor product M2M with an interpolator computed on the fly)
FBlas::scal(nnodes, FReal(0.), pole->getMultipole(idxRhs));
FReal Exp[nnodes], PermExp[nnodes];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
subChildParentInterpolator, ORDER,
const_cast<FReal*>(subCell->getMultipole(idxRhs)), ORDER, PermExp, ORDER);
// Do we need to reset multipole expansion? Same question for non adap?
FBlas::scal(nnodes, FReal(0.), pole->getMultipole(idxRhs));
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
subChildParentInterpolator + 2 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
/// p^6 version
FBlas::gemtva(nnodes, nnodes, FReal(1.),
subChildParentInterpolator,
const_cast<FReal*>(subCell->getMultipole(idxRhs)), pole->getMultipole(idxRhs));
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
subChildParentInterpolator + 1 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
// /// p^4 version
//
// FReal Exp[nnodes], PermExp[nnodes];
// // ORDER*ORDER*ORDER * (2*ORDER-1)
// FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
// subChildParentInterpolator, ORDER,
// const_cast<FReal*>(subCell->getMultipole(idxRhs)), ORDER, PermExp, ORDER);
//
// for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// // ORDER*ORDER*ORDER * (2*ORDER-1)
// FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
// subChildParentInterpolator + 2 * ORDER*ORDER, ORDER,
// Exp, ORDER, PermExp, ORDER);
//
// for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// // ORDER*ORDER*ORDER * (2*ORDER-1)
// FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
// subChildParentInterpolator + 1 * ORDER*ORDER, ORDER,
// Exp, ORDER, PermExp, ORDER);
//
// for (unsigned int n=0; n<nnodes; ++n) pole->getMultipole(idxRhs)[perm[2][n]] += PermExp[n];
//
for (unsigned int n=0; n<nnodes; ++n) pole->getMultipole(idxRhs)[perm[2][n]] += PermExp[n];
// // 2) Apply Discete Fourier Transform
// M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs),
......@@ -184,28 +216,31 @@ public:
}
void M2L(CellClass* const local, const int localLevel, const CellClass* const pole, const int poleLevel) override {
// local->setDataDown(local->getDataDown() + pole->getDataUp());
//local->setDataDown(local->getDataDown() + pole->getDataUp());
// Source cell: pole
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2, poleLevel)));
const FPoint poleCellCenter(KernelBaseClass::getLeafCellCenter(pole->getCoordinate()));
const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
// Target cell: local
const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2, localLevel)));
const FPoint localCellCenter(KernelBaseClass::getLeafCellCenter(local->getCoordinate()));
const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
// interpolation points of source (Y) and target (X) cell
FPoint X[nnodes], Y[nnodes];
FUnifTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y);
FUnifTensor<order>::setRoots(localCellCenter, localCellWidth, X);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// Dense M2L
const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs);
for (unsigned int m=0; m<nnodes; ++m)
for (unsigned int n=0; n<nnodes; ++n)
local->getLocal(idxRhs)[m]=MatrixKernel->evaluate(X[m], Y[n]) * MultipoleExpansion[n];
for (unsigned int n=0; n<nnodes; ++n){
local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * MultipoleExpansion[n];
}
}
}
......@@ -218,64 +253,90 @@ public:
void L2L(const CellClass* const local, const int localLevel, CellClass* const subCell, const int subCellLevel) override {
// subCell->setDataDown(local->getDataDown() + subCell->getDataDown());
const FPoint subCellCenter(KernelBaseClass::getLeafCellCenter(subCell->getCoordinate()));
const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel));
const FReal subCellWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-subCellLevel);
const FPoint localCenter(KernelBaseClass::getLeafCellCenter(local->getCoordinate()));
const FPoint localCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
const FReal localWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-localLevel);
// Set sub-child coords
FReal globalChildCoords[3][ORDER];
FUnifTensor<order>::setPolynomialsRoots(subCellCenter, subCellWidth, globalChildCoords);
////////////////////////////////////////////////////////////////////////////
/// p^6 version
// allocate memory
FReal* subChildParentInterpolator = new FReal [nnodes * nnodes];
// Map global position of sub-child nodes to [-1,1]
FReal localChildCoords[3][ORDER];
// set child info
FPoint ChildRoots[nnodes], localChildRoots[nnodes];
FUnifTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots);
// map global position of roots to local position in parent cell
const map_glob_loc map(localCenter, localWidth);
FPoint localChildPoints;
for (unsigned int n=0; n<ORDER; ++n) {
map(FPoint(globalChildCoords[0][n],globalChildCoords[1][n],globalChildCoords[2][n]), localChildPoints);
localChildCoords[0][n] = localChildPoints.getX();
localChildCoords[1][n] = localChildPoints.getY();
localChildCoords[2][n] = localChildPoints.getZ();
}
for (unsigned int n=0; n<nnodes; ++n)
map(ChildRoots[n], localChildRoots[n]);
// assemble interpolator
FReal* subChildParentInterpolator = new FReal [3 * ORDER*ORDER];
KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[0], subChildParentInterpolator);
KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[1], subChildParentInterpolator + 1 * ORDER*ORDER);
KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[2], subChildParentInterpolator + 2 * ORDER*ORDER);
// assemble child - parent - interpolator
KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator);
// get permutation operators
unsigned int perm[3][nnodes];
for (unsigned int i=0;i<3; ++i)
for (unsigned int n=0; n<nnodes; ++n)
perm[i][n] = KernelBaseClass::Interpolator->getPermutationsM2ML2L(i)[n];
// ////////////////////////////////////////////////////////////////////////////
// /// p^4 version
// // Set sub-child coords
// FReal globalChildCoords[3][ORDER];
// FUnifTensor<order>::setPolynomialsRoots(subCellCenter, subCellWidth, globalChildCoords);
//
// // Map global position of sub-child nodes to [-1,1]
// FReal localChildCoords[3][ORDER];
// const map_glob_loc map(localCenter, localWidth);
// FPoint localChildPoints;
// for (unsigned int n=0; n<ORDER; ++n) {
// map(FPoint(globalChildCoords[0][n],globalChildCoords[1][n],globalChildCoords[2][n]), localChildPoints);
// localChildCoords[0][n] = localChildPoints.getX();
// localChildCoords[1][n] = localChildPoints.getY();
// localChildCoords[2][n] = localChildPoints.getZ();
// }
//
// // assemble interpolator
// FReal* subChildParentInterpolator = new FReal [3 * ORDER*ORDER];
// KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[0], subChildParentInterpolator);
// KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[1], subChildParentInterpolator + 1 * ORDER*ORDER);
// KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[2], subChildParentInterpolator + 2 * ORDER*ORDER);
//
// // get permutation operators
// unsigned int perm[3][nnodes];
// for (unsigned int i=0;i<3; ++i)
// for (unsigned int n=0; n<nnodes; ++n)
// perm[i][n] = KernelBaseClass::Interpolator->getPermutationsM2ML2L(i)[n];
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // 1) Apply Inverse Discete Fourier Transform
// M2LHandler.unapplyZeroPaddingAndDFT(local->getTransformedLocal(idxRhs),
// const_cast<CellClass*>(local)->getLocal(idxRhs));
// 2) apply Sx
FReal Exp[nnodes], PermExp[nnodes];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
subChildParentInterpolator, ORDER,
const_cast<FReal*>(local->getLocal(idxRhs)), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
subChildParentInterpolator + 2 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
subChildParentInterpolator + 1 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) subCell->getLocal(idxRhs)[perm[2][n]] += PermExp[n];
// total flops count: 3 * ORDER*ORDER*ORDER * (2*ORDER-1)
/// p^6 version
FBlas::gemva(nnodes, nnodes, FReal(1.),
subChildParentInterpolator,
const_cast<FReal*>(local->getLocal(idxRhs)), subCell->getLocal(idxRhs));
// /// p^4 version
//
// // 2) apply Sx
// FReal Exp[nnodes], PermExp[nnodes];
// // ORDER*ORDER*ORDER * (2*ORDER-1)
// FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
// subChildParentInterpolator, ORDER,
// const_cast<FReal*>(local->getLocal(idxRhs)), ORDER, PermExp, ORDER);
//
// for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// // ORDER*ORDER*ORDER * (2*ORDER-1)
// FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
// subChildParentInterpolator + 2 * ORDER*ORDER, ORDER,
// Exp, ORDER, PermExp, ORDER);
//
// for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// // ORDER*ORDER*ORDER * (2*ORDER-1)
// FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
// subChildParentInterpolator + 1 * ORDER*ORDER, ORDER,
// Exp, ORDER, PermExp, ORDER);
//
// for (unsigned int n=0; n<nnodes; ++n) subCell->getLocal(idxRhs)[perm[2][n]] += PermExp[n];
// // total flops count: 3 * ORDER*ORDER*ORDER * (2*ORDER-1)
}
......@@ -286,7 +347,8 @@ public:
// for(int idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += local->getDataDown();
// }
const FPoint CellCenter(KernelBaseClass::getLeafCellCenter(local->getCoordinate()));
const FPoint CellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),cellLevel));
const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-cellLevel);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
......
......@@ -72,6 +72,27 @@ protected:
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
/**
* @brief Return the position of the center of a cell from its tree
* coordinate
* @param FTreeCoordinate
* @param inLevel the current level of Cell
*/
FPoint getCellCenter(const FTreeCoordinate coordinate, int inLevel)
{
//Set the boxes width needed
FReal widthAtCurrentLevel = BoxWidthLeaf*FReal(1 << (TreeHeight-(inLevel+1)));
FReal widthAtCurrentLevelDiv2 = widthAtCurrentLevel/FReal(2.);
//Set the center real coordinates from box corner and widths.
FReal X = BoxCorner.getX() + FReal(coordinate.getX())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
FReal Y = BoxCorner.getY() + FReal(coordinate.getY())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
FReal Z = BoxCorner.getZ() + FReal(coordinate.getZ())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
return FPoint(X,Y,Z);
}
public:
/**
* The constructor initializes all constant attributes and it reads the
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program 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 General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FUNIFDENSEKERNEL_HPP
#define FUNIFDENSEKERNEL_HPP
#include "../../Utils/FGlobal.hpp"
#include "../../Utils/FTrace.hpp"
#include "../../Utils/FSmartPointer.hpp"
#include "./FAbstractUnifKernel.hpp"
#include "./FUnifM2LHandler.hpp"
class FTreeCoordinate;
/**
* @author Pierre Blanchard (pierre.blanchard@inria.fr)
* @class FUnifDenseKernel
* @brief
* Please read the license
*
* This kernels implement the Lagrange interpolation based FMM operators. It
* implements all interfaces (P2P,P2M,M2M,M2L,L2L,L2P) which are required by
* the FFmmAlgorithm and FFmmAlgorithmThread.
* This is the dense version of the kernel. The transfer are done in real space
* and not in Fourier space.
*
* @tparam CellClass Type of cell
* @tparam ContainerClass Type of container to store particles
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Lagrange interpolation order
*/
template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FUnifDenseKernel
: public FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
// private types
typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
// using from
typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
AbstractBaseClass;
/// Needed for P2P and M2L operators
const MatrixKernelClass *const MatrixKernel;
/// Needed for M2L operator
const M2LHandlerClass 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).
*/
FUnifDenseKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const MatrixKernelClass *const inMatrixKernel)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
MatrixKernel(inMatrixKernel),
M2LHandler(MatrixKernel,
inTreeHeight,
inBoxWidth)
{ }
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
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxRhs), SourceParticles);
// // 2) apply Discrete Fourier Transform
// M2LHandler.applyZeroPaddingAndDFT(LeafCell->getMultipole(idxRhs),
// LeafCell->getTransformedMultipole(idxRhs));
}
}
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
// Do we need to reset multipole expansion when adaptive kernel used?
//FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
ParentCell->getMultipole(idxRhs));
}
}
// // 2) Apply Discete Fourier Transform
// M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs),
// ParentCell->getTransformedMultipole(idxRhs));
}
}
// 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)
{
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
// const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
// interpolation points of source (Y) and target (X) cell
FPoint X[AbstractBaseClass::nnodes], Y[AbstractBaseClass::nnodes];
FUnifTensor<ORDER>::setRoots(AbstractBaseClass::getCellCenter(TargetCell->getCoordinate(),TreeLevel), CellWidth, X);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// FComplexe *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxRhs);
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
// M2LHandler.applyFC(idx, TreeLevel, scale,
// SourceCells[idx]->getTransformedMultipole(idxRhs),
// TransformedLocalExpansion);
FUnifTensor<ORDER>::setRoots(AbstractBaseClass::getCellCenter(SourceCells[idx]->getCoordinate(),TreeLevel), CellWidth, Y);
for (unsigned int m=0; m<AbstractBaseClass::nnodes; ++m)
for (unsigned int n=0; n<AbstractBaseClass::nnodes; ++n){
TargetCell->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * SourceCells[idx]->getMultipole(idxRhs)[n];
}
}
}
}
}
// 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);
// }