MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ChargedBondedInteraction.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 
29 #include "Particles/BaseParticle.h"
30 #include "InteractionHandler.h"
31 #include <iomanip>
32 //#include <cassert>
33 
36 
43  : BaseInteraction(P, I, timeStamp)
44 {
45  //ensuring that, by default, particles are not 'bonded'
46  //i.e. they will not unintentionally 'stick' to any overlapping particles!
47  bonded_ = false;
48 #ifdef DEBUG_CONSTRUCTOR
49  std::cout<<"ChargedBondedInteraction::ChargedBondedInteraction() finished"<<std::endl;
50 #endif
51 }
52 
57  : BaseInteraction(p)
58 {
59  //carrying the history parameter over for copied particles to ensure that any bonded particles
60  //remain bonded!
61  bonded_ = p.bonded_;
62 #ifdef DEBUG_CONSTRUCTOR
63  std::cout<<"ChargedBondedInteraction::ChargedBondedInteraction(const ChargedBondedInteraction &p finished"<<std::endl;
64 #endif
65 }
66 
68 {
69 #ifdef MERCURY_USE_MPI
70  logger(FATAL,"ChargedBondedInteractions are currently not implemented in parallel MercuryDPM");
71 #endif
72 }
73 
78 {
79 #ifdef DEBUG_DESTRUCTOR
80  std::cout<<"ChargedBondedInteraction::ChargedBondedInteractionaction() finished"<<std::endl;
81 #endif
82 }
83 
87 void ChargedBondedInteraction::write(std::ostream& os UNUSED) const
88 {
89  os << " bonded " << bonded_;
90 }
91 
95 void ChargedBondedInteraction::read(std::istream& is UNUSED)
96 {
97  std::string dummy;
98  //logger(INFO,"ChargedBondedSpecies %",dummy);
99  is >> dummy >> bonded_;
100 }
101 
106 {
107 
108  const ChargedBondedSpecies* species = getSpecies();
109  //std::cout << getSpecies()->getCharge() << std::endl;
110 
111  //creating local parameters to store the charges of both particles
112  //involved in the interaction to allow for quick calculation
113  const auto pSpecies = dynamic_cast<const ChargedBondedSpecies*>(getP()->getSpecies());
114  const auto iSpecies = dynamic_cast<const ChargedBondedSpecies*>(getI()->getSpecies());
115  logger.assert_debug(pSpecies,"No ChargedBondedSpecies");
116  logger.assert_debug(iSpecies,"No ChargedBondedSpecies");
117  const int pCharge = pSpecies->getCharge();
118  const int iCharge = iSpecies->getCharge();
119 
120  //similarly, creating local parameters to store the relevant stiffness
121  //and max force values
122  const Mdouble k = species->getAdhesionStiffness();
123  const Mdouble fMax = species->getAdhesionForceMax();
124 
125  const Mdouble kWaals = species->getVanDerWaalsStiffness();
126  const Mdouble fMaxWaals = species->getVanDerWaalsForceMax();
127  const Mdouble rWaals = fMaxWaals / kWaals;
128 
129 
130  //First, adding bonded force if applicable
131  if (bonded_ && getOverlap() >= 0)
132  {
133  addForce(getNormal() * (-species->getBondForceMax()
135  return;
136  }
137 
138 
139  //determining which of the three possible cases for force based on charge -
140  //repulsive (like charges), attractive (unlike charges) or none (1 or more uncharged) -
141  //is relevant for the current combination of particle charges...
142  //(Note that the charge set function contains a safety check that means charge can only be
143  // +/- 1, i.e. the expressions used below should not produce errors!
144 
145  //case 1 - 1 or more particles has no charge, i.e. no EM force between them
146  if ((pCharge == 0) or (iCharge == 0))
147  {
148  //No need to write anything here as nothing needs to be returned!
149  //std::cout << "no charge" << std::endl;
150  }
151  //case 2: unlike charges --> attractive force
152  else if (pCharge == -iCharge)
153  {
154  //std::cout << "dissimilar charge" << std::endl;
155  //in the case of particles with opposing charges, simply implements the
156  //standard attractive force used for adhesive particle interactions
157  if (getOverlap() >= 0)
158  {
159  addForce(getNormal() * -fMax);
160  addForce(getNormal() * -fMaxWaals);
161  }
162  else if (getOverlap() >= -rWaals)
163  {
164  addForce(getNormal() * (-kWaals * getOverlap() - fMaxWaals));
165  addForce(getNormal() * (-k * getOverlap() - fMax));
166  }
167  else
168  {
169  addForce(getNormal() * (-k * getOverlap() - fMax));
170  }
171  }
172  //case 3: like charges --> repulsive force
173  else if (pCharge == iCharge)
174  {
175  //std::cout << "similar charge" << std::endl;
176  //in the case of particles with like charges, reverse the direction of the force applied
177  //such that particles repel one another
178  if (getOverlap() >= 0)
179  {
180  addForce(getNormal() * +fMax);
181  addForce(getNormal() * -fMaxWaals);
183  }
184  else if (getOverlap() >= -rWaals)
185  {
186  addForce(getNormal() * (-kWaals * getOverlap() - fMaxWaals));
187  addForce(getNormal() * (+k * getOverlap() + fMax));
188  //std::cout << "Waals = " << getNormal() * (-kWaals * getOverlap() - fMaxWaals) << std::endl;
189  }
190  else
191  {
192  addForce(getNormal() * (+k * getOverlap() + fMax));
193  }
194  }
195  //if none of the above are satisfied, something must have gone very wrong!
196  else
197  {
198  logger(ERROR, "Particle charge has erroneous value");
199  }
200 }
201 
208 {
209  const ChargedBondedSpecies* species = getSpecies();
210  const auto pSpecies = static_cast<const ChargedBondedSpecies*>(getP()->getSpecies()->getAdhesiveForce());
211  const auto iSpecies = static_cast<const ChargedBondedSpecies*>(getI()->getSpecies()->getAdhesiveForce());
212  logger.assert_debug(pSpecies,"No ChargedBondedSpecies");
213  logger.assert_debug(iSpecies,"No ChargedBondedSpecies");
214  const int pCharge = pSpecies->getCharge();
215  const int iCharge = iSpecies->getCharge();
216 
217  const Mdouble k = species->getAdhesionStiffness();
218  const Mdouble fMax = species->getAdhesionForceMax();
220 
221  const Mdouble kWaals = species->getVanDerWaalsStiffness();
222  const Mdouble fMaxWaals = species->getVanDerWaalsForceMax();
223  const Mdouble rWaals = (fMaxWaals == 0) ? 0 : (fMaxWaals / kWaals);
224 
225 
226  //First, adding bonded force if applicable
227  if (bonded_ && getOverlap() >= 0)
228  {
229  //comment to ignore BondForce
230  Mdouble elasticEnergyAtEquilibrium = getElasticEnergyAtEquilibrium(species->getBondForceMax());
231  return -species->getBondForceMax() * getOverlap() + elasticEnergyAtEquilibrium;
232  }
233 
234  Mdouble elasticEnergy = 0.0;
235  if ((pCharge != 0) && (iCharge != 0))
236  {
237  if (pCharge == -iCharge)
238  {
239  if (getOverlap() >= 0)
240  {
241  elasticEnergy -= (0.5 * rWaals + getOverlap()) * fMaxWaals
242  + (0.5 * r + getOverlap()) * fMax;
243  }
244  else if (getOverlap() >= -rWaals)
245  {
246  elasticEnergy -= (0.5 * kWaals * mathsFunc::square(getOverlap() + rWaals)) +
247  (0.5 * k * mathsFunc::square(getOverlap() + r));
248  }
249  else
250  {
251  elasticEnergy -= (0.5 * k * mathsFunc::square(getOverlap() + r));
252  }
253  }
254  //case 3: like charges --> repulsive force
255  else if (pCharge == iCharge)
256  {
257  if (getOverlap() >= 0)
258  {
259  elasticEnergy -= (0.5 * rWaals + getOverlap()) * fMaxWaals
260  - (0.5 * r + getOverlap()) * fMax;
261  }
262  else if (getOverlap() >= -rWaals)
263  {
264  elasticEnergy -= (0.5 * kWaals * mathsFunc::square(getOverlap() + rWaals)) -
265  (0.5 * k * mathsFunc::square(getOverlap() + r));
266  }
267  else
268  {
269  elasticEnergy += (0.5 * k * mathsFunc::square(getOverlap() + r));
270  }
271  }
272  else
273  {
274  logger(ERROR, "Particle charge has erroneous value");
275  }
276  }
277  return elasticEnergy;
278 }
279 
284 {
285  return static_cast<const ChargedBondedSpecies*> (getBaseSpecies()->getAdhesiveForce()); //downcast
286 }
287 
292 {
293  return "ChargedBonded";
294 }
295 
301 {
302  bonded_ = true;
303 }
304 
311 {
312  bonded_ = false;
313 }
Mdouble getVanDerWaalsForceMax() const
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
BaseAdhesiveForce * getAdhesiveForce() const
Definition: BaseSpecies.h:152
std::string getBaseName() const
Returns the name of the interaction, see Interaction.h.
Mdouble getAdhesionForceMax() const
Allows the spring constant to be accessed.
void write(std::ostream &os) const override
Interaction print function, which accepts an std::ostream as input.
BaseInteractable * getI()
Returns a pointer to the second object involved in the interaction (often a wall or a particle)...
Mdouble getInteractionDistance() const
returns the largest separation distance at which adhesive short-range forces can occur.
Definition: BaseSpecies.h:146
Stores information about interactions between two interactable objects; often particles but could be ...
ChargedBondedSpecies contains the parameters used to describe a linear reversible short-range force...
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
Mdouble getBondDissipation() const
Allows the additional dissipation used to damp oscillations between bondd particles to be accessed...
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
Mdouble getVanDerWaalsStiffness() const
Mdouble getNormalRelativeVelocity() const
Returns a double which is the norm (length) of the relative velocity vector.
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
#define UNUSED
Definition: GeneralDefine.h:39
Mdouble getAdhesionStiffness() const
Allows the spring constant to be accessed.
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
void addForce(Vec3D force)
add an force increment to the total force.
Mdouble getElasticEnergy() const override
Returns the amount of Elastic energy involved in an interaction. Basically used in case you want to w...
const ChargedBondedSpecies * getSpecies() const
Returns a pointer to the adhesive force species ChargedBondedSpecies.
~ChargedBondedInteraction() override
Destructor.
Mdouble getBondForceMax() const
Allows the maximal force for 'bonding' particles together to be accessed.
virtual Mdouble getElasticEnergyAtEquilibrium(Mdouble adhesiveForce) const
void bond()
A pair of functions which can be used to fix or unfix a pair of overlapping particles.
Defines the basic properties that a interactable object can have.
void computeAdhesionForce()
Computes the adhesive forces.
T square(const T val)
squares a number
Definition: ExtendedMath.h:106
bool bonded_
A history parameter to store if the particles were in contact or not. Useful to compute adhesive forc...