00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "Particle/ParticleInteractLayout.h"
00028 #include "Particle/ParticleBConds.h"
00029 #include "Particle/ParticleBase.h"
00030 #include "Region/RegionLayout.h"
00031 #include "FieldLayout/FieldLayout.h"
00032 #include "Utility/IpplInfo.h"
00033 #include "Message/Communicate.h"
00034 #include "Message/Message.h"
00035 #include "Profile/Profiler.h"
00036
00037 #ifdef IPPL_STDSTL
00038 #include <algorithm>
00039 using namespace std;
00040 #else
00041 #include <algo.h>
00042 #endif // IPPL_STDSTL
00043
00045
00046 template < class T, unsigned Dim, class Mesh >
00047 ParticleInteractLayout<T,Dim,Mesh>::ParticleInteractLayout(FieldLayout<Dim>&
00048 fl)
00049 : ParticleSpatialLayout<T,Dim,Mesh>(fl) {
00050 TAU_TYPE_STRING(taustr, CT(*this) + " void (" + CT(fl) + " )" );
00051 TAU_PROFILE("ParticleInteractLayout::ParticleInteractLayout()", taustr,
00052 TAU_PARTICLE);
00053 setup();
00054 }
00055
00056
00058
00059 template < class T, unsigned Dim, class Mesh >
00060 ParticleInteractLayout<T,Dim,Mesh>::ParticleInteractLayout(FieldLayout<Dim>&
00061 fl, Mesh& mesh)
00062 : ParticleSpatialLayout<T,Dim,Mesh>(fl,mesh) {
00063 TAU_TYPE_STRING(taustr, "void (" + CT(fl) + ", " + CT(mesh) + " )" );
00064 TAU_PROFILE("ParticleInteractLayout::ParticleInteractLayout()", taustr,
00065 TAU_PARTICLE);
00066 setup();
00067 }
00068
00069
00071
00072 template < class T, unsigned Dim, class Mesh >
00073 ParticleInteractLayout<T,Dim,Mesh>::ParticleInteractLayout(const
00074 RegionLayout<T,Dim,Mesh>& rl) : ParticleSpatialLayout<T,Dim,Mesh>(rl) {
00075 TAU_TYPE_STRING(taustr, "void (" + CT(rl) + " )" );
00076 TAU_PROFILE("ParticleInteractLayout::ParticleInteractLayout()", taustr,
00077 TAU_PARTICLE);
00078 setup();
00079 }
00080
00081
00083
00084
00085 template < class T, unsigned Dim, class Mesh >
00086 ParticleInteractLayout<T,Dim,Mesh>::ParticleInteractLayout()
00087 : ParticleSpatialLayout<T,Dim,Mesh>() {
00088 TAU_TYPE_STRING(taustr, CT(*this) + " void ()");
00089 TAU_PROFILE("ParticleInteractLayout::ParticleInteractLayout()", taustr,
00090 TAU_PARTICLE);
00091 setup();
00092 }
00093
00094
00096
00097 template < class T, unsigned Dim, class Mesh >
00098 void ParticleInteractLayout<T,Dim,Mesh>::setup() {
00099 TAU_TYPE_STRING(taustr, CT(*this) + " void ()");
00100 TAU_PROFILE("ParticleInteractLayout::setup()", taustr, TAU_PARTICLE);
00101
00102
00103 unsigned N = Ippl::getNodes();
00104 InterNodeList = new bool[N];
00105 SentToNodeList = new bool[N];
00106 InteractionNodes = 0;
00107
00108
00109 InterRadius = 0;
00110 InterRadiusArray = 0;
00111 MaxGlobalInterRadius = 0;
00112 }
00113
00114
00116
00117 template < class T, unsigned Dim, class Mesh >
00118 ParticleInteractLayout<T,Dim,Mesh>::~ParticleInteractLayout() {
00119 TAU_TYPE_STRING(taustr, CT(*this) + " void ()");
00120 TAU_PROFILE("ParticleInteractLayout::~ParticleInteractLayout()", taustr,
00121 TAU_PARTICLE);
00122
00123 delete [] InterNodeList;
00124 delete [] SentToNodeList;
00125
00126 for (int i=(PairList.size() - 1); i >= 0; --i)
00127 delete (PairList[i]);
00128 }
00129
00130
00132
00133 template < class T, unsigned Dim, class Mesh >
00134 T ParticleInteractLayout<T,Dim,Mesh>::getMaxLocalInteractionRadius() {
00135 TAU_TYPE_STRING(taustr, CT(InterRadius) + " ()");
00136 TAU_PROFILE("ParticleInteractLayout::getMaxLocalInteractionRadius()",
00137 taustr, TAU_PARTICLE);
00138
00139 if (InterRadiusArray != 0) {
00140 if (InterRadiusArray->size() > 0)
00141 return *(max_element(InterRadiusArray->begin(),
00142 InterRadiusArray->end()));
00143 else
00144 return 0.0;
00145 } else {
00146 return InterRadius;
00147 }
00148 }
00149
00150
00152
00153
00154
00155
00156
00157 template < class T, unsigned Dim, class Mesh >
00158 void ParticleInteractLayout<T,Dim,Mesh>::getPairlist(unsigned n,
00159 pair_iterator& bpi, pair_iterator& epi,
00160 ParticleBase< ParticleInteractLayout<T,Dim,Mesh> >& PData) {
00161
00162 TAU_TYPE_STRING(taustr, "void (unsigned, pair_iterator, pair_iterator, "
00163 + CT(PData) + " )");
00164 TAU_PROFILE("ParticleInteractLayout::getPairlist()", taustr, TAU_PARTICLE);
00165
00166
00167 if (getUpdateFlag(ParticleLayout<T,Dim>::BCONDS)) {
00168
00169 ParticleBConds<T,Dim>& pBConds = this->getBConds();
00170 bool periodicBC[2*Dim];
00171 unsigned numPeriodicBC = 0;
00172 typename ParticleBConds<T,Dim>::ParticleBCond periodicBCond = ParticlePeriodicBCond;
00173 for (unsigned d=0; d<2*Dim; ++d) {
00174 periodicBC[d] = (pBConds[d] == periodicBCond);
00175 if (periodicBC[d]) ++numPeriodicBC;
00176 }
00177 if (numPeriodicBC>0) {
00178
00179
00180 swap_ghost_particles(PData.getLocalNum(), PData, periodicBC);
00181 }
00182 else {
00183 swap_ghost_particles(PData.getLocalNum(), PData);
00184 }
00185 }
00186 else {
00187
00188
00189 swap_ghost_particles(PData.getLocalNum(), PData);
00190 }
00191
00192
00193 bpi = PairList[n]->begin();
00194 epi = PairList[n]->end();
00195
00196 return;
00197 }
00198
00199
00201
00202
00203
00204 template < class T, unsigned Dim, class Mesh >
00205 void ParticleInteractLayout<T,Dim,Mesh>::rebuild_interaction_data() {
00206 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
00207 TAU_PROFILE("ParticleInteractLayout::rebuild_interaction_data()",
00208 taustr, TAU_PARTICLE);
00209
00210 int j;
00211 unsigned d;
00212
00213 DEBUGMSG("ParticleInteractLayout: rebuilding interaction node data."<<endl);
00214
00215
00216 InteractionNodes = 0;
00217 T interRad = 2.0 * getMaxInteractionRadius();
00218
00219
00220 unsigned N = Ippl::getNodes();
00221 for (j=0; j < N; ++j)
00222 InterNodeList[j] = false;
00223
00224
00225 if (interRad <= 0.0)
00226 return;
00227
00228
00229 typename RegionLayout<T,Dim,Mesh>::iterator_iv localVN, endLocalVN = this->RLayout.end_iv();
00230
00231 for (localVN = this->RLayout.begin_iv(); localVN != endLocalVN; ++localVN) {
00232
00233
00234 NDRegion<T,Dim> chkDom((*localVN).second->getDomain());
00235 for (d=0; d < Dim; ++d) {
00236 chkDom[d] = PRegion<T>(chkDom[d].first() - interRad,
00237 chkDom[d].last() + interRad);
00238 }
00239
00240
00241
00242 DEBUGMSG(" Checking domain " << chkDom << endl);
00243 typename RegionLayout<T,Dim,Mesh>::touch_range_dv touchingVN =
00244 this->RLayout.touch_range_rdv(chkDom);
00245 typename RegionLayout<T,Dim,Mesh>::touch_iterator_dv tVN = touchingVN.first;
00246 DEBUGMSG(" Touching vnodes:" << endl);
00247 for ( ; tVN != touchingVN.second; ++tVN) {
00248
00249
00250 DEBUGMSG(" vnode: node = " << ((*tVN).second)->getNode());
00251 DEBUGMSG(", domain = " << ((*tVN).second)->getDomain() << endl);
00252 unsigned vn = ((*tVN).second)->getNode();
00253 if ( ! InterNodeList[vn] ) {
00254 InterNodeList[vn] = true;
00255 InteractionNodes++;
00256 }
00257 }
00258 }
00259
00260 DEBUGMSG(" There are " << InteractionNodes << " inter. nodes." << endl);
00261
00262
00263
00264
00265 NeedGhostSwap = true;
00266 return;
00267 }
00268
00269
00271
00272
00273
00274
00275 template < class T, unsigned Dim, class Mesh >
00276 void ParticleInteractLayout<T,Dim,Mesh>::rebuild_interaction_data(
00277 const bool periodicBC[2*Dim])
00278 {
00279 TAU_TYPE_STRING(taustr, CT(*this) + " void (const bool [])" );
00280 TAU_PROFILE("ParticleInteractLayout::rebuild_interaction_data()",
00281 taustr, TAU_PARTICLE);
00282
00283 int j;
00284 unsigned d;
00285 unsigned pe = Ippl::myNode();
00286
00287 DEBUGMSG("ParticleInteractLayout: rebuilding interaction node data."<<endl);
00288
00289
00290 InteractionNodes = 0;
00291 T interRad = 2.0 * getMaxInteractionRadius();
00292
00293
00294 unsigned N = Ippl::getNodes();
00295 for (j=0; j < N; ++j)
00296 InterNodeList[j] = false;
00297
00298
00299 if (interRad <= 0.0)
00300 return;
00301
00302
00303 const NDRegion<T,Dim>& globalDom = this->RLayout.getDomain();
00304
00305
00306 T offset[Dim];
00307 unsigned numRef;
00308 bool flipBit, activeBit[Dim], refBit[Dim];
00309 NDRegion<T,Dim> chkDom, refDom;
00310
00311
00312 typename RegionLayout<T,Dim,Mesh>::iterator_iv localVN, endLocalVN = this->RLayout.end_iv();
00313 typename RegionLayout<T,Dim,Mesh>::iterator_iv localVN2;
00314 typename RegionLayout<T,Dim,Mesh>::touch_range_dv touchingVN;
00315 typename RegionLayout<T,Dim,Mesh>::touch_iterator_dv tVN;
00316
00317
00318 for (localVN = this->RLayout.begin_iv(); localVN != endLocalVN; ++localVN) {
00319
00320
00321 chkDom = (*localVN).second->getDomain();
00322 for (d=0; d<Dim; ++d) {
00323 chkDom[d] = PRegion<T>(chkDom[d].first() - interRad,
00324 chkDom[d].last() + interRad);
00325 }
00326
00327
00328
00329 DEBUGMSG(" Checking domain " << chkDom << endl);
00330 touchingVN = this->RLayout.touch_range_rdv(chkDom);
00331 DEBUGMSG(" Touching vnodes:" << endl);
00332 for (tVN = touchingVN.first; tVN != touchingVN.second; ++tVN) {
00333
00334
00335 DEBUGMSG(" vnode: node = " << ((*tVN).second)->getNode());
00336 DEBUGMSG(", domain = " << ((*tVN).second)->getDomain() << endl);
00337 unsigned vn = ((*tVN).second)->getNode();
00338 if ( ! InterNodeList[vn] ) {
00339 InterNodeList[vn] = true;
00340 InteractionNodes++;
00341 }
00342 }
00343
00344
00345 numRef = 0;
00346 for (d=0; d<Dim; ++d) {
00347 if (periodicBC[2*d] && chkDom[d].first()<globalDom[d].first()) {
00348
00349 offset[d] = globalDom[d].length();
00350 activeBit[d] = true;
00351 numRef = 2 * numRef + 1;
00352 }
00353 else if (periodicBC[2*d+1] && chkDom[d].last()>globalDom[d].last()) {
00354
00355 offset[d] = -globalDom[d].length();
00356 activeBit[d] = true;
00357 numRef = 2 * numRef + 1;
00358 }
00359 else {
00360 offset[d] = 0.0;
00361 activeBit[d] = false;
00362 }
00363 refBit[d] = false;
00364 }
00365
00366
00367 for (j=0; j<numRef; ++j) {
00368
00369 refDom = chkDom;
00370
00371 d = 0;
00372 flipBit = false;
00373 while (d<Dim && !flipBit) {
00374
00375 if (activeBit[d]) {
00376
00377 if (refBit[d]) {
00378
00379 refBit[d] = false;
00380 }
00381 else {
00382
00383 refBit[d] = true;
00384 flipBit = true;
00385 }
00386 }
00387 ++d;
00388 }
00389 PAssert(flipBit);
00390
00391
00392 for (d=0; d<Dim; ++d) {
00393 if (refBit[d]) refDom[d] = refDom[d] + offset[d];
00394 }
00395
00396
00397
00398 DEBUGMSG(" Checking domain " << refDom << endl);
00399 touchingVN = this->RLayout.touch_range_rdv(refDom);
00400 DEBUGMSG(" Touching vnodes:" << endl);
00401 for (tVN = touchingVN.first; tVN != touchingVN.second; ++tVN) {
00402
00403
00404 DEBUGMSG(" vnode: node = " << ((*tVN).second)->getNode());
00405 DEBUGMSG(", domain = " << ((*tVN).second)->getDomain() << endl);
00406 unsigned vn = ((*tVN).second)->getNode();
00407 if ( ! InterNodeList[vn] ) {
00408 InterNodeList[vn] = true;
00409 InteractionNodes++;
00410 }
00411 }
00412
00413 if (!InterNodeList[pe]) {
00414
00415 bool interact = false;
00416 localVN2 = this->RLayout.begin_iv();
00417 while (localVN2 != endLocalVN && !interact) {
00418 interact = refDom.touches((*localVN2).second->getDomain());
00419 ++localVN2;
00420 }
00421 if (interact) {
00422 InterNodeList[pe] = true;
00423 InteractionNodes++;
00424 }
00425 }
00426 }
00427
00428 }
00429
00430 DEBUGMSG(" There are " << InteractionNodes << " inter. nodes." << endl);
00431
00432
00433
00434
00435 NeedGhostSwap = true;
00436 return;
00437 }
00438
00439
00441
00442
00443
00444
00445 template < class T, unsigned Dim, class Mesh >
00446 void ParticleInteractLayout<T,Dim,Mesh>::update(
00447 ParticleBase< ParticleInteractLayout<T,Dim,Mesh> >& PData,
00448 const ParticleAttrib<char>* canSwap) {
00449
00450 TAU_TYPE_STRING(taustr, "void (" + CT(PData) + ", ParticleAttrib<char>* )");
00451 TAU_PROFILE("ParticleInteractLayout::update()", taustr, TAU_PARTICLE);
00452
00453 unsigned N = Ippl::getNodes();
00454 unsigned myN = Ippl::myNode();
00455 unsigned LocalNum = PData.getLocalNum();
00456 unsigned DestroyNum = PData.getDestroyNum();
00457 unsigned TotalNum;
00458 T maxrad = getMaxLocalInteractionRadius();
00459 int node;
00460
00461
00462 PData.performDestroy();
00463 LocalNum -= DestroyNum;
00464
00465
00466
00467 if ( ! this->RLayout.initialized())
00468 rebuild_layout(LocalNum,PData);
00469
00470
00471 if (getUpdateFlag(ParticleLayout<T,Dim>::BCONDS))
00472 apply_bconds(LocalNum, PData.R, this->getBConds(), this->RLayout.getDomain());
00473
00474
00475
00476
00477 if (N > 1 && getUpdateFlag(this->SWAP)) {
00478 if (canSwap==0)
00479 LocalNum = swap_particles(LocalNum, PData);
00480 else
00481 LocalNum = swap_particles(LocalNum, PData, *canSwap);
00482 }
00483
00484
00485 NeedGhostSwap = true;
00486
00487
00488 TotalNum = this->NodeCount[myN] = LocalNum;
00489
00490
00491
00492 if (N > 1) {
00493
00494
00495 int tag1 = Ippl::Comm->next_tag(P_SPATIAL_LAYOUT_TAG, P_LAYOUT_CYCLE);
00496 int tag2 = Ippl::Comm->next_tag(P_SPATIAL_RETURN_TAG, P_LAYOUT_CYCLE);
00497 if (myN != 0) {
00498 Message *msg = new Message;
00499
00500
00501 msg->put(LocalNum);
00502
00503
00504 msg->put(maxrad);
00505
00506
00507 Ippl::Comm->send(msg, 0, tag1);
00508
00509
00510
00511 node = 0;
00512 msg = Ippl::Comm->receive_block(node, tag2);
00513 msg->get(this->NodeCount);
00514 msg->get(maxrad);
00515 msg->get(TotalNum);
00516 } else {
00517
00518 int notrecvd = N - 1;
00519 TotalNum = LocalNum;
00520 while (notrecvd > 0) {
00521
00522 node = Communicate::COMM_ANY_NODE;
00523 Message *msg = Ippl::Comm->receive_block(node, tag1);
00524 int remNodeCount = 0;
00525 T remMaxRad = 0.0;
00526 msg->get(remNodeCount);
00527 msg->get(remMaxRad);
00528 delete msg;
00529 notrecvd--;
00530
00531
00532 TotalNum += remNodeCount;
00533 this->NodeCount[node] = remNodeCount;
00534 if (remMaxRad > maxrad)
00535 maxrad = remMaxRad;
00536 }
00537
00538
00539 Message *msg = new Message;
00540 msg->put(this->NodeCount, this->NodeCount + N);
00541 msg->put(maxrad);
00542 msg->put(TotalNum);
00543 Ippl::Comm->broadcast_others(msg, tag2);
00544 }
00545 }
00546
00547
00548 PData.setTotalNum(TotalNum);
00549 PData.setLocalNum(LocalNum);
00550
00551
00552 if (maxrad != getMaxInteractionRadius()) {
00553 setMaxInteractionRadius(maxrad);
00554
00555 if (getUpdateFlag(ParticleLayout<T,Dim>::BCONDS)) {
00556
00557 ParticleBConds<T,Dim>& pBConds = this->getBConds();
00558 bool periodicBC[2*Dim];
00559 unsigned numPeriodicBC = 0;
00560 typename ParticleBConds<T,Dim>::ParticleBCond periodicBCond=ParticlePeriodicBCond;
00561 for (unsigned d=0; d<2*Dim; ++d) {
00562 periodicBC[d] = (pBConds[d] == periodicBCond);
00563 if (periodicBC[d]) ++numPeriodicBC;
00564 }
00565 if (numPeriodicBC>0) {
00566
00567
00568 rebuild_interaction_data(periodicBC);
00569 }
00570 else {
00571 rebuild_interaction_data();
00572 }
00573 }
00574 else {
00575 rebuild_interaction_data();
00576 }
00577 }
00578 return;
00579 }
00580
00581
00583
00584
00585
00586
00587
00588 template < class T, unsigned Dim, class Mesh >
00589 void ParticleInteractLayout<T,Dim,Mesh>::swap_ghost_particles(unsigned
00590 LocalNum,
00591 ParticleBase< ParticleInteractLayout<T,Dim,Mesh> >& PData) {
00592
00593 TAU_TYPE_STRING(taustr, "void (unsigned, " + CT(PData) + " )");
00594 TAU_PROFILE("ParticleInteractLayout::swap_ghost_particles()", taustr,
00595 TAU_PARTICLE);
00596
00597 int i, j;
00598
00599
00600 if ( ! NeedGhostSwap ) return;
00601
00602
00603 NeedGhostSwap = false;
00604
00605
00606
00607
00608 PData.ghostDestroy(PData.getGhostNum(), 0);
00609
00610
00611 unsigned N = Ippl::getNodes();
00612 unsigned sendnum = 0;
00613 for (i=0; i < N; i++)
00614 if (InterNodeList[i] && this->NodeCount[i] > 0)
00615 sendnum++;
00616
00617
00618
00619 if (sendnum == 0 || LocalNum == 0) {
00620 find_pairs(LocalNum, 0, LocalNum, true, PData);
00621 return;
00622 }
00623
00624
00625
00626 T interRad = 2.0 * getMaxInteractionRadius();
00627
00628
00629 NDRegion<T,Dim> pLoc;
00630
00631
00632
00633
00634
00635
00636 for (i=0; i < N; i++)
00637 if (InterNodeList[i] && this->NodeCount[i] > 0)
00638 this->SwapMsgList[i] = new Message;
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698 find_pairs(LocalNum, 0, LocalNum, true, PData);
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717 }
00718
00719
00721
00722
00723
00724
00725
00726
00727 template < class T, unsigned Dim, class Mesh >
00728 void ParticleInteractLayout<T,Dim,Mesh>::swap_ghost_particles(
00729 unsigned LocalNum,
00730 ParticleBase< ParticleInteractLayout<T,Dim,Mesh> >& PData,
00731 const bool periodicBC[2*Dim])
00732 {
00733 TAU_TYPE_STRING(taustr, "void (unsigned, " + CT(PData) + ", const bool [])");
00734 TAU_PROFILE("ParticleInteractLayout::swap_ghost_particles()", taustr,
00735 TAU_PARTICLE);
00736
00737 int i, j;
00738 unsigned d;
00739
00740
00741 if ( ! NeedGhostSwap ) return;
00742
00743
00744 NeedGhostSwap = false;
00745
00746
00747
00748
00749 PData.ghostDestroy(PData.getGhostNum(), 0);
00750
00751
00752 unsigned N = Ippl::getNodes();
00753 unsigned pe = Ippl::myNode();
00754 unsigned sendnum = 0;
00755 for (i=0; i < N; i++)
00756 if (InterNodeList[i] && this->NodeCount[i] > 0)
00757 sendnum++;
00758
00759
00760
00761 if (sendnum == 0 || LocalNum == 0) {
00762 find_pairs(LocalNum, 0, LocalNum, true, PData);
00763 return;
00764 }
00765
00766
00767
00768 T interRad = 2.0 * getMaxInteractionRadius();
00769
00770
00771 const NDRegion<T,Dim>& globalDom = this->RLayout.getDomain();
00772
00773
00774 T offset[Dim];
00775 unsigned numRef;
00776 bool flipBit, activeBit[Dim], refBit[Dim];
00777 NDRegion<T,Dim> pLoc, refLoc;
00778
00779
00780 typename RegionLayout<T,Dim,Mesh>::iterator_iv localVN, endLocalVN = this->RLayout.end_iv();
00781 typename RegionLayout<T,Dim,Mesh>::touch_range_dv touchingVN;
00782 typename RegionLayout<T,Dim,Mesh>::touch_iterator_dv tVNit;
00783 SingleParticlePos_t savePos;
00784
00785
00786
00787
00788
00789
00790 for (i=0; i < N; i++)
00791 if (InterNodeList[i] && this->NodeCount[i] > 0)
00792 this->SwapMsgList[i] = new Message;
00793
00794
00795
00796
00797
00798 for (i=0; i < LocalNum; ++i) {
00799
00800
00801 memset((void *)SentToNodeList, 0, N * sizeof(bool));
00802
00803
00804
00805 for (j=0; j < Dim; ++j)
00806 pLoc[j] = PRegion<T>(PData.R[i][j] - interRad,
00807 PData.R[i][j] + interRad);
00808
00809
00810 touchingVN = this->RLayout.touch_range_rdv(pLoc);
00811 for (tVNit = touchingVN.first; tVNit != touchingVN.second; ++tVNit) {
00812 Rnode<T,Dim> *tVN = (*tVNit).second;
00813 unsigned node = tVN->getNode();
00814
00815
00816 if (this->NodeCount[node] > 0 && ! SentToNodeList[node]) {
00817 if (! InterNodeList[node]) {
00818 ERRORMSG("ParticleInteractLayout: Cannot send ghost " << i);
00819 ERRORMSG(" to node " << node << " ... skipping." << endl);
00820 }
00821 else {
00822 PData.ghostPutMessage(*(this->SwapMsgList[node]), 1, i);
00823 SentToNodeList[node] = true;
00824 }
00825 }
00826 }
00827
00828
00829 numRef = 0;
00830 for (d=0; d<Dim; ++d) {
00831 if (periodicBC[2*d] && pLoc[d].first()<globalDom[d].first()) {
00832
00833 offset[d] = globalDom[d].length();
00834 activeBit[d] = true;
00835 numRef = 2 * numRef + 1;
00836 }
00837 else if (periodicBC[2*d+1] && pLoc[d].last()>globalDom[d].last()) {
00838
00839 offset[d] = -globalDom[d].length();
00840 activeBit[d] = true;
00841 numRef = 2 * numRef + 1;
00842 }
00843 else {
00844 offset[d] = 0.0;
00845 activeBit[d] = false;
00846 }
00847 refBit[d] = false;
00848 }
00849
00850 if (numRef>0) savePos = PData.R[i];
00851
00852
00853 for (j=0; j<numRef; ++j) {
00854
00855 refLoc = pLoc;
00856 PData.R[i] = savePos;
00857
00858 d = 0;
00859 flipBit = false;
00860 while (d<Dim && !flipBit) {
00861
00862 if (activeBit[d]) {
00863
00864 if (refBit[d]) {
00865
00866 refBit[d] = false;
00867 }
00868 else {
00869
00870 refBit[d] = true;
00871 flipBit = true;
00872 }
00873 }
00874 ++d;
00875 }
00876 PAssert(flipBit);
00877
00878
00879 for (d=0; d<Dim; ++d) {
00880 if (refBit[d]) {
00881 refLoc[d] = refLoc[d] + offset[d];
00882 PData.R[i][d] = PData.R[i][d] + offset[d];
00883 }
00884 }
00885
00886
00887 memset((void *)SentToNodeList, 0, N * sizeof(bool));
00888
00889
00890 touchingVN = this->RLayout.touch_range_rdv(refLoc);
00891 for (tVNit = touchingVN.first; tVNit != touchingVN.second; ++tVNit) {
00892 Rnode<T,Dim> *tVN = (*tVNit).second;
00893 unsigned node = tVN->getNode();
00894
00895
00896 if (this->NodeCount[node] > 0 && ! SentToNodeList[node]) {
00897 if (! InterNodeList[node]) {
00898 ERRORMSG("ParticleInteractLayout: Cannot send ghost " << i);
00899 ERRORMSG(" to node " << node << " ... skipping." << endl);
00900 }
00901 else {
00902 PData.ghostPutMessage(*(this->SwapMsgList[node]), 1, i);
00903 SentToNodeList[node] = true;
00904 }
00905 }
00906 }
00907
00908 if (InterNodeList[pe]) {
00909
00910 bool interact = false;
00911 localVN = this->RLayout.begin_iv();
00912 while (localVN != endLocalVN && !interact) {
00913 interact = refLoc.touches((*localVN).second->getDomain());
00914 ++localVN;
00915 }
00916 if (interact) {
00917 PData.ghostPutMessage(*(this->SwapMsgList[pe]), 1, i);
00918 }
00919 }
00920 }
00921 if (numRef>0) PData.R[i] = savePos;
00922
00923 }
00924
00925
00926 int tag = Ippl::Comm->next_tag(P_SPATIAL_GHOST_TAG, P_LAYOUT_CYCLE);
00927 for (i=0; i < N; ++i) {
00928 if (InterNodeList[i] && this->NodeCount[i] > 0) {
00929
00930 PData.ghostPutMessage(*(this->SwapMsgList[i]), (unsigned)0, (unsigned)0);
00931
00932
00933 Ippl::Comm->send(this->SwapMsgList[i], i, tag);
00934 }
00935 }
00936
00937
00938 find_pairs(LocalNum, 0, LocalNum, true, PData);
00939
00940
00941 while (sendnum-- > 0) {
00942 int node = Communicate::COMM_ANY_NODE;
00943 unsigned oldGN = PData.getGhostNum();
00944 Message *recmsg = Ippl::Comm->receive_block(node, tag);
00945
00946 while (PData.ghostGetMessage(*recmsg, node) > 0);
00947 delete recmsg;
00948
00949
00950 find_pairs(LocalNum, LocalNum + oldGN, LocalNum + PData.getGhostNum(),
00951 false, PData);
00952 }
00953 }
00954
00955
00957
00958
00959 template < class T, unsigned Dim, class Mesh >
00960 void ParticleInteractLayout<T,Dim,Mesh>::find_pairs(const unsigned LocalNum,
00961 const unsigned a1, const unsigned a2, const bool initLists,
00962 ParticleBase< ParticleInteractLayout<T,Dim,Mesh> >& PData) {
00963
00964 TAU_TYPE_STRING(taustr, "void (unsigned, unsigned, unsigned, bool, "
00965 + CT(PData) + " )");
00966 TAU_PROFILE("ParticleInteractLayout::find_pairs()", taustr,
00967 TAU_PARTICLE);
00968
00969 unsigned i, j;
00970
00971
00972 if (initLists) {
00973 unsigned vlen = PairList.size();
00974 if (vlen > LocalNum)
00975 vlen = LocalNum;
00976 for (i=0; i < vlen; ++i)
00977 PairList[i]->erase(PairList[i]->begin(), PairList[i]->end());
00978
00979
00980 if (PairList.size() < LocalNum) {
00981 int newamt = LocalNum - PairList.size();
00982 PairList.reserve(newamt);
00983 for (i=0; i < newamt; ++i)
00984 PairList.push_back(new vector<pair_t>);
00985 }
00986 }
00987
00988
00989 if (a2 <= a1) return;
00990
00991
00992 for (i=0; i < LocalNum; ++i) {
00993
00994 T intrad1 = getInteractionRadius(i);
00995
00996
00997 j = (a1 > i ? a1 : i + 1);
00998
00999
01000
01001 for (; j < LocalNum; ++j) {
01002
01003 T intrad2 = intrad1 + getInteractionRadius(j);
01004 intrad2 *= intrad2;
01005
01006
01007 Vektor<T,Dim> rsep = PData.R[j];
01008 rsep -= PData.R[i];
01009 T sep2 = dot(rsep, rsep);
01010
01011
01012
01013 if (sep2 < intrad2) {
01014 PairList[i]->push_back(pair_t(j, sep2));
01015 PairList[j]->push_back(pair_t(i, sep2));
01016 }
01017 }
01018
01019
01020
01021 for (; j < a2; ++j) {
01022
01023 T intrad2 = intrad1 + getInteractionRadius(j);
01024 intrad2 *= intrad2;
01025
01026
01027 Vektor<T,Dim> rsep = PData.R[j];
01028 rsep -= PData.R[i];
01029 T sep2 = dot(rsep, rsep);
01030
01031
01032
01033
01034 if (sep2 < intrad2) {
01035 PairList[i]->push_back(pair_t(j, sep2));
01036 }
01037 }
01038 }
01039 }
01040
01041
01043
01044 template < class T, unsigned Dim, class Mesh >
01045 ostream& operator<<(ostream& out, const ParticleInteractLayout<T,Dim,Mesh>& L) {
01046 TAU_TYPE_STRING(taustr, "ostream (ostream, " + CT(L) + " )");
01047 TAU_PROFILE("ParticleInteractLayout::operator<<()", taustr,
01048 TAU_PARTICLE | TAU_IO);
01049
01050 out << "ParticleInteractLayout, with particle distribution:\n ";
01051 for (unsigned int i=0; i < Ippl::getNodes(); ++i)
01052 out << L.getNodeCount(i) << " ";
01053 out << "\nInteractLayout decomposition = " << L.getLayout();
01054 return out;
01055 }
01056
01057
01059
01060
01061
01062 template < class T, unsigned Dim, class Mesh >
01063 void ParticleInteractLayout<T,Dim,Mesh>::Repartition(UserList* userlist) {
01064 TAU_TYPE_STRING(taustr, "void (UserList *)");
01065 TAU_PROFILE("ParticleInteractLayout::Repartition()", taustr, TAU_PARTICLE);
01066
01067
01068
01069 if (userlist->getUserListID() == this->RLayout.get_Id()) {
01070
01071 this->rebuild_neighbor_data();
01072
01073
01074
01075
01076 InteractionNodes = 0;
01077 setMaxInteractionRadius(0);
01078 NeedGhostSwap = true;
01079 }
01080 }
01081
01082
01083
01084
01085
01086
01087