libhoomd/data_structures/ParticleGroup.cc

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2009-2014 The Regents of
   4 the University of Michigan All rights reserved.
   5
   6 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   7 copyright is held, by various Contributors who have granted The Regents of the
   8 University of Michigan the right to modify and/or distribute such Contributions.
   9
  10 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  11 and binary forms, provided you abide by the following conditions:
  12
  13 * Redistributions of source code must retain the above copyright notice, this
  14 list of conditions, and the following disclaimer both in the code and
  15 prominently in any materials provided with the distribution.
  16
  17 * Redistributions in binary form must reproduce the above copyright notice, this
  18 list of conditions, and the following disclaimer in the documentation and/or
  19 other materials provided with the distribution.
  20
  21 * All publications and presentations based on HOOMD-blue, including any reports
  22 or published results obtained, in whole or in part, with HOOMD-blue, will
  23 acknowledge its use according to the terms posted at the time of submission on:
  24 http://codeblue.umich.edu/hoomd-blue/citations.html
  25
  26 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  27 http://codeblue.umich.edu/hoomd-blue/
  28
  29 * Apart from the above required attributions, neither the name of the copyright
  30 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  31 promote products derived from this software without specific prior written
  32 permission.
  33
  34 Disclaimer
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  37 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  38 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  39 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  40
  41 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  42 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  43 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  44 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  46 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  47 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48 */
  49
  50 // Maintainer: joaander
  51
  52 #ifdef WIN32
  53 #pragma warning( push )
  54 #pragma warning( disable : 4267 4244 )
  55 #endif
  56
  57 #include "ParticleGroup.h"
  58
  59 #ifdef ENABLE_CUDA
  60 #include "ParticleGroup.cuh"
  61 #include "CachedAllocator.h"
  62 #endif
  63
  64 #include <boost/python.hpp>
  65 #include <boost/bind.hpp>
  66 using namespace boost::python;
  67 using namespace boost;
  68
  69 #include <algorithm>
  70 #include <iostream>
  71 using namespace std;
  72
  73 /*! \file ParticleGroup.cc
  74     \brief Defines the ParticleGroup and related classes
  75 */
  76
  77 //////////////////////////////////////////////////////////////////////////////
  78 // ParticleSelector
  79
  80 /*! \param sysdef System the particles are to be selected from
  81 */
  82 ParticleSelector::ParticleSelector(boost::shared_ptr<SystemDefinition> sysdef)
  83     : m_sysdef(sysdef), m_pdata(sysdef->getParticleData())
  84     {
  85     assert(m_sysdef);
  86     assert(m_pdata);
  87     m_exec_conf = m_pdata->getExecConf();
  88     }
  89
  90 /*! \param tag Tag of the particle to check
  91     \returns true if the particle is selected
  92     \returns false if it is not
  93 */
  94 bool ParticleSelector::isSelected(unsigned int tag) const
  95     {
  96     // base class doesn't do anything useful
  97     return false;
  98     }
  99
 100 //////////////////////////////////////////////////////////////////////////////
 101 // ParticleSelectorTag
 102
 103 /*! \param sysdef System the particles are to be selected from
 104     \param tag_min Minimum tag to select (inclusive)
 105     \param tag_max Maximum tag to select (inclusive)
 106 */
 107 ParticleSelectorTag::ParticleSelectorTag(boost::shared_ptr<SystemDefinition> sysdef,
 108                                          unsigned int tag_min,
 109                                          unsigned int tag_max)
 110     : ParticleSelector(sysdef), m_tag_min(tag_min), m_tag_max(tag_max)
 111     {
 112     // make a quick check on the sanity of the input data
 113     if (m_tag_max < m_tag_min)
 114         m_exec_conf->msg->warning() << "group: max < min specified when selecting particle tags" << endl;
 115
 116     if (m_tag_max >= m_pdata->getNGlobal())
 117         {
 118         m_exec_conf->msg->error() << "Cannot select particles with tags larger than the number of particles "
 119              << endl;
 120         throw runtime_error("Error selecting particles");
 121         }
 122     }
 123
 124 /*! \param tag Tag of the particle to check
 125     \returns true if \a m_tag_min <= \a tag <= \a m_tag_max
 126 */
 127 bool ParticleSelectorTag::isSelected(unsigned int tag) const
 128     {
 129     assert(tag < m_pdata->getNGlobal());
 130     return (m_tag_min <= tag && tag <= m_tag_max);
 131     }
 132
 133 //////////////////////////////////////////////////////////////////////////////
 134 // ParticleSelectorType
 135
 136 /*! \param sysdef System the particles are to be selected from
 137     \param typ_min Minimum type id to select (inclusive)
 138     \param typ_max Maximum type id to select (inclusive)
 139 */
 140 ParticleSelectorType::ParticleSelectorType(boost::shared_ptr<SystemDefinition> sysdef,
 141                                            unsigned int typ_min,
 142                                            unsigned int typ_max)
 143     : ParticleSelector(sysdef), m_typ_min(typ_min), m_typ_max(typ_max)
 144     {
 145     // make a quick check on the sanity of the input data
 146     if (m_typ_max < m_typ_min)
 147         m_exec_conf->msg->warning() << "group: max < min specified when selecting particle types" << endl;
 148
 149     if (m_typ_max >= m_pdata->getNTypes())
 150         m_exec_conf->msg->warning() << "group: Requesting the selection of a non-existant particle type" << endl;
 151     }
 152
 153 /*! \param tag Tag of the particle to check
 154     \returns true if the type of particle \a tag is in the inclusive range [ \a m_typ_min, \a m_typ_max ]
 155 */
 156 bool ParticleSelectorType::isSelected(unsigned int tag) const
 157     {
 158     assert(tag < m_pdata->getNGlobal());
 159     unsigned int typ = m_pdata->getType(tag);
 160
 161     // see if it matches the criteria
 162     bool result = (m_typ_min <= typ && typ <= m_typ_max);
 163
 164     return result;
 165     }
 166
 167 //////////////////////////////////////////////////////////////////////////////
 168 // ParticleSelectorRigid
 169
 170 /*! \param sysdef System the particles are to be selected from
 171     \param rigid true selects particles that are in rigid bodies, false selects particles that are not part of a body
 172 */
 173 ParticleSelectorRigid::ParticleSelectorRigid(boost::shared_ptr<SystemDefinition> sysdef,
 174                                              bool rigid)
 175     : ParticleSelector(sysdef), m_rigid(rigid)
 176     {
 177     }
 178
 179 /*! \param tag Tag of the particle to check
 180     \returns true if the type of particle \a tag meets the rigid critera selected
 181 */
 182 bool ParticleSelectorRigid::isSelected(unsigned int tag) const
 183     {
 184     assert(tag < m_pdata->getNGlobal());
 185
 186     // get body id of current particle tag
 187     unsigned int body = m_pdata->getBody(tag);
 188
 189     // see if it matches the criteria
 190     bool result = false;
 191     if (m_rigid && body != NO_BODY)
 192         result = true;
 193     if (!m_rigid && body == NO_BODY)
 194         result = true;
 195
 196     return result;
 197     }
 198
 199 //////////////////////////////////////////////////////////////////////////////
 200 // ParticleSelectorCuboid
 201
 202 ParticleSelectorCuboid::ParticleSelectorCuboid(boost::shared_ptr<SystemDefinition> sysdef, Scalar3 min, Scalar3 max)
 203     :ParticleSelector(sysdef), m_min(min), m_max(max)
 204     {
 205     // make a quick check on the sanity of the input data
 206     if (m_min.x >= m_max.x || m_min.y >= m_max.y || m_min.z >= m_max.z)
 207         m_exec_conf->msg->warning() << "group: max < min specified when selecting particle in a cuboid" << endl;
 208     }
 209
 210 /*! \param tag Tag of the particle to check
 211     \returns true if the type of particle \a tag is in the cuboid
 212
 213     Evaluation is performed by \a m_min.x <= x < \a m_max.x so that multiple cuboids stacked next to each other
 214     do not have overlapping sets of particles.
 215 */
 216 bool ParticleSelectorCuboid::isSelected(unsigned int tag) const
 217     {
 218     assert(tag < m_pdata->getNGlobal());
 219
 220     // identify the index of the current particle tag
 221     Scalar3 pos = m_pdata->getPosition(tag);
 222
 223     // see if it matches the criteria
 224     bool result = (m_min.x <= pos.x && pos.x < m_max.x &&
 225                    m_min.y <= pos.y && pos.y < m_max.y &&
 226                    m_min.z <= pos.z && pos.z < m_max.z);
 227
 228     return result;
 229     }
 230
 231 //////////////////////////////////////////////////////////////////////////////
 232 // ParticleGroup
 233
 234 /*! \param sysdef System definition to build the group from
 235     \param selector ParticleSelector used to choose the group members
 236
 237     Particles where criteria falls within the range [min,max] (inclusive) are added to the group.
 238 */
 239 ParticleGroup::ParticleGroup(boost::shared_ptr<SystemDefinition> sysdef, boost::shared_ptr<ParticleSelector> selector)
 240     : m_sysdef(sysdef),
 241       m_pdata(sysdef->getParticleData()),
 242       m_exec_conf(m_pdata->getExecConf()),
 243       m_num_local_members(0),
 244       m_particles_sorted(true)
 245     {
 246     #ifdef ENABLE_CUDA
 247     if (m_pdata->getExecConf()->isCUDAEnabled())
 248         {
 249         // create a ModernGPU context
 250         m_mgpu_context = mgpu::CreateCudaDeviceAttachStream(0);
 251         }
 252     #endif
 253
 254     // update member tag arrays
 255     updateMemberTags(selector);
 256
 257     // connect to the particle sort signal
 258     m_sort_connection = m_pdata->connectParticleSort(bind(&ParticleGroup::slotParticleSort, this));
 259
 260     // connect reallocate() method to maximum particle number change signal
 261     m_max_particle_num_change_connection = m_pdata->connectMaxParticleNumberChange(bind(&ParticleGroup::reallocate, this));
 262     }
 263
 264 /*! \param sysdef System definition to build the group from
 265     \param member_tags List of particle tags that belong to the group
 266
 267     All particles specified in \a member_tags will be added to the group.
 268 */
 269 ParticleGroup::ParticleGroup(boost::shared_ptr<SystemDefinition> sysdef, const std::vector<unsigned int>& member_tags)
 270     : m_sysdef(sysdef),
 271       m_pdata(sysdef->getParticleData()),
 272       m_exec_conf(m_pdata->getExecConf()),
 273       m_num_local_members(0),
 274       m_particles_sorted(true)
 275     {
 276     // let's make absolutely sure that the tag order given from outside is sorted
 277     std::vector<unsigned int> sorted_member_tags =  member_tags;
 278     sort(sorted_member_tags.begin(), sorted_member_tags.end());
 279
 280     // store member tags
 281     GPUArray<unsigned int> member_tags_array(member_tags.size(), m_exec_conf);
 282     m_member_tags.swap(member_tags_array);
 283
 284         {
 285         ArrayHandle<unsigned int> h_member_tags(m_member_tags, access_location::host, access_mode::overwrite);
 286         std::copy(sorted_member_tags.begin(), sorted_member_tags.end(), h_member_tags.data);
 287         }
 288
 289     // one byte per particle to indicate membership in the group, initialize with current number of local particles
 290     GPUArray<unsigned char> is_member(m_pdata->getMaxN(), m_pdata->getExecConf());
 291     m_is_member.swap(is_member);
 292
 293     GPUArray<unsigned char> is_member_tag(m_pdata->getNGlobal(), m_pdata->getExecConf());
 294     m_is_member_tag.swap(is_member_tag);
 295
 296     // build the reverse lookup table for tags
 297     buildTagHash();
 298
 299     GPUArray<unsigned int> member_idx(member_tags.size(), m_pdata->getExecConf());
 300     m_member_idx.swap(member_idx);
 301
 302     #ifdef ENABLE_CUDA
 303     if (m_pdata->getExecConf()->isCUDAEnabled())
 304         {
 305         // create a ModernGPU context
 306         m_mgpu_context = mgpu::CreateCudaDeviceAttachStream(0);
 307         }
 308     #endif
 309
 310     // now that the tag list is completely set up and all memory is allocated, rebuild the index list
 311     rebuildIndexList();
 312
 313     // connect to the particle sort signal
 314     m_sort_connection = m_pdata->connectParticleSort(bind(&ParticleGroup::slotParticleSort, this));
 315
 316     // connect reallocate() method to maximum particle number change signal
 317     m_max_particle_num_change_connection = m_pdata->connectMaxParticleNumberChange(bind(&ParticleGroup::reallocate, this));
 318     }
 319
 320 ParticleGroup::~ParticleGroup()
 321     {
 322     // disconnect the sort connection, but only if there was a particle data to connect it to in the first place
 323     if (m_pdata)
 324         {
 325         m_sort_connection.disconnect();
 326         m_max_particle_num_change_connection.disconnect();
 327         }
 328     }
 329
 330 void ParticleGroup::updateMemberTags(boost::shared_ptr<ParticleSelector> selector)
 331     {
 332     // assign all of the particles that belong to the group
 333     // for each particle in the (global) data
 334     vector<unsigned int> member_tags;
 335     for (unsigned int tag = 0; tag < m_pdata->getNGlobal(); tag++)
 336         {
 337         // add the tag to the list if it matches the selection
 338         if (selector->isSelected(tag))
 339             member_tags.push_back(tag);
 340         }
 341
 342     // store member tags
 343     GPUArray<unsigned int> member_tags_array(member_tags.size(), m_pdata->getExecConf());
 344     m_member_tags.swap(member_tags_array);
 345
 346         {
 347         ArrayHandle<unsigned int> h_member_tags(m_member_tags, access_location::host, access_mode::overwrite);
 348         std::copy(member_tags.begin(), member_tags.end(), h_member_tags.data);
 349         }
 350
 351
 352     // one byte per particle to indicate membership in the group, initialize with current number of local particles
 353     GPUArray<unsigned char> is_member(m_pdata->getMaxN(), m_pdata->getExecConf());
 354     m_is_member.swap(is_member);
 355
 356     GPUArray<unsigned char> is_member_tag(m_pdata->getNGlobal(), m_pdata->getExecConf());
 357     m_is_member_tag.swap(is_member_tag);
 358
 359     // build the reverse lookup table for tags
 360     buildTagHash();
 361
 362     GPUArray<unsigned int> member_idx(member_tags.size(), m_pdata->getExecConf());
 363     m_member_idx.swap(member_idx);
 364
 365     // now that the tag list is completely set up and all memory is allocated, rebuild the index list
 366     rebuildIndexList();
 367     }
 368
 369 void ParticleGroup::reallocate()
 370     {
 371     m_is_member.resize(m_pdata->getMaxN());
 372
 373     if (m_is_member_tag.getNumElements() != m_pdata->getNGlobal())
 374         {
 375         // reallocate if necessary
 376         GPUArray<unsigned char> is_member_tag(m_pdata->getNGlobal(), m_exec_conf);
 377         m_is_member_tag.swap(is_member_tag);
 378
 379         buildTagHash();
 380         }
 381     }
 382
 383 /*! \returns Total mass of all particles in the group
 384     \note This method aquires the ParticleData internally
 385 */
 386 Scalar ParticleGroup::getTotalMass() const
 387     {
 388     // grab the particle data
 389     ArrayHandle< Scalar4 > h_vel(m_pdata->getVelocities(), access_location::host, access_mode::read);
 390
 391     // loop  through all indices in the group and total the mass
 392     Scalar total_mass = 0.0;
 393     for (unsigned int i = 0; i < getNumMembers(); i++)
 394         {
 395         unsigned int idx = getMemberIndex(i);
 396         total_mass += h_vel.data[idx].w;
 397         }
 398     return total_mass;
 399     }
 400
 401 /*! \returns The center of mass of the group, in unwrapped coordinates
 402     \note This method aquires the ParticleData internally
 403 */
 404 Scalar3 ParticleGroup::getCenterOfMass() const
 405     {
 406     // grab the particle data
 407     ArrayHandle< Scalar4 > h_vel(m_pdata->getVelocities(), access_location::host, access_mode::read);
 408     ArrayHandle< Scalar4 > h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 409     ArrayHandle< int3 > h_image(m_pdata->getImages(), access_location::host, access_mode::read);
 410
 411     // grab the box dimensions
 412     BoxDim box = m_pdata->getBox();
 413
 414     // loop  through all indices in the group and compute the weighted average of the positions
 415     Scalar total_mass = 0.0;
 416     Scalar3 center_of_mass = make_scalar3(Scalar(0.0), Scalar(0.0), Scalar(0.0));
 417     for (unsigned int i = 0; i < getNumMembers(); i++)
 418         {
 419         unsigned int idx = getMemberIndex(i);
 420         Scalar mass = h_vel.data[idx].w;
 421         total_mass += mass;
 422         Scalar3 pos = make_scalar3(h_pos.data[idx].x, h_pos.data[idx].y, h_pos.data[idx].z);
 423         Scalar3 unwrapped = box.shift(pos, h_image.data[idx]);
 424         center_of_mass.x += mass * unwrapped.x;
 425         center_of_mass.y += mass * unwrapped.y;
 426         center_of_mass.z += mass * unwrapped.z;
 427         }
 428     center_of_mass.x /= total_mass;
 429     center_of_mass.y /= total_mass;
 430     center_of_mass.z /= total_mass;
 431
 432     return center_of_mass;
 433     }
 434
 435 /*! \param a First particle group
 436     \param b Second particle group
 437
 438     \returns A shared pointer to a newly created particle group that contains all the elements present in \a a and
 439     \a b
 440 */
 441 boost::shared_ptr<ParticleGroup> ParticleGroup::groupUnion(boost::shared_ptr<ParticleGroup> a,
 442                                                            boost::shared_ptr<ParticleGroup> b)
 443     {
 444     // vector to store the new list of tags
 445     vector<unsigned int> member_tags;
 446
 447     if (a != b)
 448         {
 449         // make the union
 450         ArrayHandle<unsigned int> h_members_a(a->m_member_tags, access_location::host, access_mode::read);
 451         ArrayHandle<unsigned int> h_members_b(b->m_member_tags, access_location::host, access_mode::read);
 452
 453         insert_iterator< vector<unsigned int> > ii(member_tags, member_tags.begin());
 454         set_union(h_members_a.data,
 455                   h_members_a.data + a->getNumMembersGlobal(),
 456                   h_members_b.data,
 457                   h_members_b.data + b->getNumMembersGlobal(),
 458                   ii);
 459         }
 460     else
 461         {
 462         // If the two arguments are the same, just return a copy of the whole group (we cannot
 463         // acquire the member_tags array twice)
 464         ArrayHandle<unsigned int> h_members_a(a->m_member_tags, access_location::host, access_mode::read);
 465
 466         insert_iterator< vector<unsigned int> > ii(member_tags, member_tags.begin());
 467         std::copy(h_members_a.data,
 468                   h_members_a.data + a->getNumMembersGlobal(),
 469                   ii);
 470         }
 471
 472
 473
 474     // create the new particle group
 475     boost::shared_ptr<ParticleGroup> new_group(new ParticleGroup(a->m_sysdef, member_tags));
 476
 477     // return the newly created group
 478     return new_group;
 479     }
 480
 481 /*! \param a First particle group
 482     \param b Second particle group
 483
 484     \returns A shared pointer to a newly created particle group that contains only the elements present in both \a a and
 485     \a b
 486 */
 487 boost::shared_ptr<ParticleGroup> ParticleGroup::groupIntersection(boost::shared_ptr<ParticleGroup> a,
 488                                                                   boost::shared_ptr<ParticleGroup> b)
 489     {
 490     // vector to store the new list of tags
 491     vector<unsigned int> member_tags;
 492
 493     if (a != b)
 494         {
 495         // make the intersection
 496         ArrayHandle<unsigned int> h_members_a(a->m_member_tags, access_location::host, access_mode::read);
 497         ArrayHandle<unsigned int> h_members_b(b->m_member_tags, access_location::host, access_mode::read);
 498
 499         insert_iterator< vector<unsigned int> > ii(member_tags, member_tags.begin());
 500         set_intersection(h_members_a.data,
 501                          h_members_a.data + a->getNumMembersGlobal(),
 502                          h_members_b.data,
 503                          h_members_b.data + b->getNumMembersGlobal(),
 504                          ii);
 505         }
 506     else
 507         {
 508         // If the two arguments are the same, just return a copy of the whole group (we cannot
 509         // acquire the member_tags array twice)
 510         ArrayHandle<unsigned int> h_members_a(a->m_member_tags, access_location::host, access_mode::read);
 511
 512         insert_iterator< vector<unsigned int> > ii(member_tags, member_tags.begin());
 513         std::copy(h_members_a.data,
 514                   h_members_a.data + a->getNumMembersGlobal(),
 515                   ii);
 516         }
 517
 518     // create the new particle group
 519     boost::shared_ptr<ParticleGroup> new_group(new ParticleGroup(a->m_sysdef, member_tags));
 520
 521     // return the newly created group
 522     return new_group;
 523     }
 524
 525 /*! \param a First particle group
 526     \param b Second particle group
 527
 528     \returns A shared pointer to a newly created particle group that contains only the elements present in \a a, and
 529     not any present in \a b
 530 */
 531 boost::shared_ptr<ParticleGroup> ParticleGroup::groupDifference(boost::shared_ptr<ParticleGroup> a,
 532                                                                 boost::shared_ptr<ParticleGroup> b)
 533     {
 534     // vector to store the new list of tags
 535     vector<unsigned int> member_tags;
 536
 537     if (a != b)
 538         {
 539         // make the difference
 540         ArrayHandle<unsigned int> h_members_a(a->m_member_tags, access_location::host, access_mode::read);
 541         ArrayHandle<unsigned int> h_members_b(b->m_member_tags, access_location::host, access_mode::read);
 542
 543         insert_iterator< vector<unsigned int> > ii(member_tags, member_tags.begin());
 544         set_difference(h_members_a.data,
 545                   h_members_a.data + a->getNumMembersGlobal(),
 546                   h_members_b.data,
 547                   h_members_b.data + b->getNumMembersGlobal(),
 548                   ii);
 549         }
 550     else
 551         {
 552         // If the two arguments are the same, just return an empty group (we cannot
 553         // acquire the member_tags array twice)
 554         }
 555
 556
 557     // create the new particle group
 558     boost::shared_ptr<ParticleGroup> new_group(new ParticleGroup(a->m_sysdef, member_tags));
 559
 560     // return the newly created group
 561     return new_group;
 562     }
 563
 564 /*! Builds the by-tag-lookup table for group membership
 565  */
 566 void ParticleGroup::buildTagHash()
 567     {
 568     ArrayHandle<unsigned char> h_is_member_tag(m_is_member_tag, access_location::host, access_mode::overwrite);
 569     ArrayHandle<unsigned int> h_member_tags(m_member_tags, access_location::host, access_mode::read);
 570
 571     // reset member ship flags
 572     memset(h_is_member_tag.data, 0, sizeof(unsigned char)*m_pdata->getNGlobal());
 573
 574     unsigned int num_members = m_member_tags.getNumElements();
 575     for (unsigned int member = 0; member < num_members; member++)
 576         h_is_member_tag.data[h_member_tags.data[member]] = 1;
 577     }
 578
 579 /*! \pre m_member_tags has been filled out, listing all particle tags in the group
 580     \pre memory has been allocated for m_is_member and m_member_idx
 581     \post m_is_member is updated so that it reflects the current indices of the particles in the group
 582     \post m_member_idx is updated listing all particle indices belonging to the group, in index order
 583 */
 584 void ParticleGroup::rebuildIndexList() const
 585     {
 586     // notice message
 587     m_pdata->getExecConf()->msg->notice(10) << "ParticleGroup: rebuilding index" << std::endl;
 588
 589     #ifdef ENABLE_CUDA
 590     if (m_pdata->getExecConf()->isCUDAEnabled() )
 591         {
 592         rebuildIndexListGPU();
 593         }
 594     else
 595     #endif
 596         {
 597
 598         // rebuild the membership flags for the  indices in the group and construct member list
 599         ArrayHandle<unsigned char> h_is_member(m_is_member, access_location::host, access_mode::readwrite);
 600         ArrayHandle<unsigned char> h_is_member_tag(m_is_member_tag, access_location::host, access_mode::read);
 601         ArrayHandle<unsigned int> h_tag(m_pdata->getTags(), access_location::host, access_mode::read);
 602         ArrayHandle<unsigned int> h_member_idx(m_member_idx, access_location::host, access_mode::readwrite);
 603         unsigned int nparticles = m_pdata->getN();
 604         unsigned int cur_member = 0;
 605         for (unsigned int idx = 0; idx < nparticles; idx ++)
 606             {
 607             assert(h_tag.data[idx] < m_pdata->getNGlobal());
 608             unsigned char is_member = h_is_member_tag.data[h_tag.data[idx]];
 609             h_is_member.data[idx] =  is_member;
 610             if (is_member)
 611                 {
 612                 h_member_idx.data[cur_member] = idx;
 613                 cur_member++;
 614                 }
 615             }
 616
 617         m_num_local_members = cur_member;
 618         assert(m_num_local_members <= m_member_tags.getNumElements());
 619         }
 620
 621     // index has been rebuilt
 622     m_particles_sorted = false;
 623     }
 624
 625 #ifdef ENABLE_CUDA
 626 //! rebuild index list on the GPU
 627 void ParticleGroup::rebuildIndexListGPU() const
 628     {
 629     ArrayHandle<unsigned char> d_is_member(m_is_member, access_location::device, access_mode::overwrite);
 630     ArrayHandle<unsigned char> d_is_member_tag(m_is_member_tag, access_location::device, access_mode::read);
 631     ArrayHandle<unsigned int> d_member_idx(m_member_idx, access_location::device, access_mode::overwrite);
 632     ArrayHandle<unsigned int> d_tag(m_pdata->getTags(), access_location::device, access_mode::read);
 633
 634     // get temporary buffer
 635     ScopedAllocation<unsigned int> d_tmp(m_pdata->getExecConf()->getCachedAllocator(), m_pdata->getN());
 636
 637     // reset membership properties
 638     if (m_member_tags.getNumElements() > 0)
 639         {
 640         gpu_rebuild_index_list(m_pdata->getN(),
 641                            d_is_member_tag.data,
 642                            d_is_member.data,
 643                            d_member_idx.data,
 644                            d_tag.data,
 645                            m_num_local_members,
 646                            d_tmp.data,
 647                            m_mgpu_context);
 648         if (m_exec_conf->isCUDAErrorCheckingEnabled())
 649             CHECK_CUDA_ERROR();
 650         }
 651     else
 652         m_num_local_members = 0;
 653     }
 654 #endif
 655
 656 void export_ParticleGroup()
 657     {
 658     class_<ParticleGroup, boost::shared_ptr<ParticleGroup>, boost::noncopyable>
 659             ("ParticleGroup", init< boost::shared_ptr<SystemDefinition>, boost::shared_ptr<ParticleSelector> >())
 660             .def(init<boost::shared_ptr<SystemDefinition>, const std::vector<unsigned int>& >())
 661             .def(init<>())
 662             .def("getNumMembersGlobal", &ParticleGroup::getNumMembersGlobal)
 663             .def("getMemberTag", &ParticleGroup::getMemberTag)
 664             .def("getTotalMass", &ParticleGroup::getTotalMass)
 665             .def("getCenterOfMass", &ParticleGroup::getCenterOfMass)
 666             .def("groupUnion", &ParticleGroup::groupUnion)
 667             .def("groupIntersection", &ParticleGroup::groupIntersection)
 668             .def("groupDifference", &ParticleGroup::groupDifference)
 669             .def("updateMemberTags", &ParticleGroup::updateMemberTags)
 670             ;
 671
 672     class_<ParticleSelector, boost::shared_ptr<ParticleSelector>, boost::noncopyable>
 673             ("ParticleSelector", init< boost::shared_ptr<SystemDefinition> >())
 674             .def("isSelected", &ParticleSelector::isSelected)
 675             ;
 676
 677     class_<ParticleSelectorTag, boost::shared_ptr<ParticleSelectorTag>, bases<ParticleSelector>, boost::noncopyable>
 678         ("ParticleSelectorTag", init< boost::shared_ptr<SystemDefinition>, unsigned int, unsigned int >())
 679         ;
 680
 681     class_<ParticleSelectorType, boost::shared_ptr<ParticleSelectorType>, bases<ParticleSelector>, boost::noncopyable>
 682         ("ParticleSelectorType", init< boost::shared_ptr<SystemDefinition>, unsigned int, unsigned int >())
 683         ;
 684
 685     class_<ParticleSelectorRigid, boost::shared_ptr<ParticleSelectorRigid>, bases<ParticleSelector>, boost::noncopyable>
 686         ("ParticleSelectorRigid", init< boost::shared_ptr<SystemDefinition>, bool >())
 687         ;
 688
 689     class_<ParticleSelectorCuboid, boost::shared_ptr<ParticleSelectorCuboid>, bases<ParticleSelector>, boost::noncopyable>
 690         ("ParticleSelectorCuboid", init< boost::shared_ptr<SystemDefinition>, Scalar3, Scalar3 >())
 691         ;
 692     }
 693
 694 #ifdef WIN32
 695 #pragma warning( pop )
 696 #endif