libhoomd/computes/ConstraintSphere.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 : 4103 4244 )
  55 #endif
  56
  57 #include <boost/python.hpp>
  58 using namespace boost::python;
  59
  60 #include "ConstraintSphere.h"
  61 #include "EvaluatorConstraint.h"
  62 #include "EvaluatorConstraintSphere.h"
  63
  64 using namespace std;
  65
  66 /*! \file ConstraintSphere.cc
  67     \brief Contains code for the ConstraintSphere class
  68 */
  69
  70 /*! \param sysdef SystemDefinition containing the ParticleData to compute forces on
  71     \param group Group of particles on which to apply this constraint
  72     \param P position of the sphere
  73     \param r radius of the sphere
  74 */
  75 ConstraintSphere::ConstraintSphere(boost::shared_ptr<SystemDefinition> sysdef,
  76                                    boost::shared_ptr<ParticleGroup> group,
  77                                    Scalar3 P,
  78                                    Scalar r)
  79         : ForceConstraint(sysdef), m_group(group), m_P(P), m_r(r)
  80     {
  81     m_exec_conf->msg->notice(5) << "Constructing ConstraintSphere" << endl;
  82
  83     validate();
  84     }
  85
  86 ConstraintSphere::~ConstraintSphere()
  87     {
  88     m_exec_conf->msg->notice(5) << "Destroying ConstraintSphere" << endl;
  89     }
  90
  91 /*!
  92     \param P position of the sphere
  93     \param r radius of the sphere
  94 */
  95 void ConstraintSphere::setSphere(Scalar3 P, Scalar r)
  96     {
  97     m_P = P;
  98     m_r = r;
  99     validate();
 100     }
 101
 102 /*! ConstraintSphere removes 1 degree of freedom per particle in the group
 103 */
 104 unsigned int ConstraintSphere::getNDOFRemoved()
 105     {
 106     return m_group->getNumMembers();
 107     }
 108
 109 /*! Computes the specified constraint forces
 110     \param timestep Current timestep
 111 */
 112 void ConstraintSphere::computeForces(unsigned int timestep)
 113     {
 114     unsigned int group_size = m_group->getNumMembers();
 115     if (group_size == 0)
 116         return;
 117
 118     if (m_prof) m_prof->push("ConstraintSphere");
 119
 120     assert(m_pdata);
 121     // access the particle data arrays
 122     ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 123     ArrayHandle<Scalar4> h_vel(m_pdata->getVelocities(), access_location::host, access_mode::read);
 124
 125     const GPUArray< Scalar4 >& net_force = m_pdata->getNetForce();
 126     ArrayHandle<Scalar4> h_net_force(net_force, access_location::host, access_mode::read);
 127
 128
 129     ArrayHandle<Scalar4> h_force(m_force,access_location::host, access_mode::overwrite);
 130     ArrayHandle<Scalar> h_virial(m_virial,access_location::host, access_mode::overwrite);
 131     unsigned int virial_pitch = m_virial.getPitch();
 132
 133     // Zero data for force calculation.
 134     memset((void*)h_force.data,0,sizeof(Scalar4)*m_force.getNumElements());
 135     memset((void*)h_virial.data,0,sizeof(Scalar)*m_virial.getNumElements());
 136
 137    // there are enough other checks on the input data: but it doesn't hurt to be safe
 138     assert(h_force.data);
 139     assert(h_virial.data);
 140
 141     // for each of the particles in the group
 142     for (unsigned int group_idx = 0; group_idx < group_size; group_idx++)
 143         {
 144         // get the current particle properties
 145         unsigned int j = m_group->getMemberIndex(group_idx);
 146         Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z);
 147         Scalar3 V = make_scalar3(h_vel.data[j].x, h_vel.data[j].y, h_vel.data[j].z);
 148         Scalar3 F = make_scalar3(h_net_force.data[j].x, h_net_force.data[j].y, h_net_force.data[j].z);
 149         Scalar m = h_vel.data[j].w;
 150
 151         // evaluate the constraint position
 152         EvaluatorConstraint constraint(X, V, F, m, m_deltaT);
 153         EvaluatorConstraintSphere sphere(m_P, m_r);
 154         Scalar3 C = sphere.evalClosest(constraint.evalU());
 155
 156         // evaluate the constraint force
 157         Scalar3 FC;
 158         Scalar virial[6];
 159         constraint.evalConstraintForce(FC, virial, C);
 160
 161         // apply the constraint force
 162         h_force.data[j].x = FC.x;
 163         h_force.data[j].y = FC.y;
 164         h_force.data[j].z = FC.z;
 165         for (int k = 0; k < 6; k++)
 166             h_virial.data[k*virial_pitch+j]  = virial[k];
 167         }
 168
 169     if (m_prof)
 170         m_prof->pop();
 171     }
 172
 173 /*! Print warning messages if the sphere is outside the box.
 174     Generate an error if any particle in the group is not near the sphere.
 175 */
 176 void ConstraintSphere::validate()
 177     {
 178     BoxDim box = m_pdata->getBox();
 179     Scalar3 lo = box.getLo();
 180     Scalar3 hi = box.getHi();
 181
 182     if (m_P.x + m_r > hi.x || m_P.x - m_r < lo.x ||
 183         m_P.y + m_r > hi.y || m_P.y - m_r < lo.y ||
 184         m_P.z + m_r > hi.z || m_P.z - m_r < lo.z)
 185         {
 186         m_exec_conf->msg->warning() << "constrain.sphere: Sphere constraint is outside of the box. Constrained particle positions may be incorrect"
 187              << endl;
 188         }
 189
 190     unsigned int group_size = m_group->getNumMembers();
 191     if (group_size == 0)
 192         return;
 193
 194     ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 195     ArrayHandle<unsigned int> h_tag(m_pdata->getTags(), access_location::host, access_mode::read);
 196     ArrayHandle<unsigned int> h_body(m_pdata->getBodies(), access_location::host, access_mode::read);
 197
 198     // for each of the particles in the group
 199     bool errors = false;
 200     for (unsigned int group_idx = 0; group_idx < group_size; group_idx++)
 201         {
 202         // get the current particle properties
 203         unsigned int j = m_group->getMemberIndex(group_idx);
 204         Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z);
 205
 206         // evaluate the constraint position
 207         EvaluatorConstraintSphere sphere(m_P, m_r);
 208         Scalar3 C = sphere.evalClosest(X);
 209         Scalar3 V;
 210         V.x = C.x - X.x;
 211         V.y = C.y - X.y;
 212         V.z = C.z - X.z;
 213         Scalar dist = sqrt(V.x*V.x + V.y*V.y + V.z*V.z);
 214
 215         if (dist > Scalar(1.0))
 216             {
 217             m_exec_conf->msg->error() << "constrain.sphere: Particle " << h_tag.data[j] << " is more than 1 unit of"
 218                                       << " distance away from the closest point on the sphere constraint" << endl;
 219             errors = true;
 220             }
 221
 222         if (h_body.data[j] != NO_BODY)
 223             {
 224             m_exec_conf->msg->error() << "constrain.sphere: Particle " << h_tag.data[j] << " belongs to a rigid body"
 225                                       << " - cannot constrain" << endl;
 226             errors = true;
 227             }
 228         }
 229
 230     if (errors)
 231         {
 232         throw std::runtime_error("Invalid constraint specified");
 233         }
 234     }
 235
 236
 237 void export_ConstraintSphere()
 238     {
 239     class_< ConstraintSphere, boost::shared_ptr<ConstraintSphere>, bases<ForceConstraint>, boost::noncopyable >
 240     ("ConstraintSphere", init< boost::shared_ptr<SystemDefinition>,
 241                                                  boost::shared_ptr<ParticleGroup>,
 242                                                  Scalar3,
 243                                                  Scalar >())
 244     .def("setSphere", &ConstraintSphere::setSphere)
 245     .def("getNDOFRemoved", &ConstraintSphere::getNDOFRemoved)
 246     ;
 247     }
 248
 249 #ifdef WIN32
 250 #pragma warning( pop )
 251 #endif