MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StressStrainControlBoundary.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 provid->d 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 provid->d 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 
27 #include "ParticleHandler.h"
28 #include "Particles/BaseParticle.h"
29 #include "MpiDataClass.h"
30 #include "MpiContainer.h"
31 #include "DPMBase.h"
32 #include "PeriodicBoundary.h"
33 #include "LeesEdwardsBoundary.h"
34 
39  : BaseBoundary()
40 {
48  integratedShift_ = 0.0;
49  //
50  logger(DEBUG, "StressStrainControlBoundary::StressStrainControlBoundary() finished");
51 }
52 
57 {
58  return new StressStrainControlBoundary(*this);
59 }
60 
65 void StressStrainControlBoundary::write(std::ostream& os) const
66 {
68  os << " stressGoal " << stressGoal_;
69  os << " strainRate " << strainRate_;
70  os << " gainFactor " << gainFactor_;
71  os << " isStrainRateControlled " << isStrainRateControlled_;
72  os << " integratedShift " << integratedShift_;
73 // os << leesEdwardsBoundaries_.size() << ' ';
74 // for (const LeesEdwardsBoundary& b : leesEdwardsBoundaries_)
75 // {
76 // os << b << ' ';
77 // }
78 //
79 // os << periodicBoundaries_.size() << ' ';
80 // for (const PeriodicBoundary& b : periodicBoundaries_)
81 // {
82 // os << b << ' ';
83 // }
84 }
85 
90 void StressStrainControlBoundary::read(std::istream& is)
91 {
93  std::string dummy;
94  is >> dummy >> stressGoal_;
95  is >> dummy >> strainRate_;
96  is >> dummy >> gainFactor_;
97  is >> dummy >> isStrainRateControlled_;
98  set(stressGoal_, strainRate_, gainFactor_, isStrainRateControlled_);
99  is >> dummy >> integratedShift_;
100 }
101 
107 {
108  return "StressStrainControlBoundary";
109 }
110 
117 {
119 
120  //Real Stress and StrainRate Control every time step
121  logger.assert(getHandler() != nullptr,
122  "you need to set the handler of this boundary before the parameters can be set");
123 
125 
126  //this checks if stressGoal matrix is non-zero and then activate only the stress control
127  if (stressGoal_.XX != 0 || stressGoal_.YY != 0 || stressGoal_.ZZ != 0 ||
128  stressGoal_.XY != 0)
129  {
130  //the strainrate in the corresponding direction will be calculated
132  }
133  activateStrainRateControl(particleHandler);
134 }
135 
141 {
142  // Call Boundaries
144  {
145  b.checkBoundaryAfterParticlesMove(particleHandler);
146  }
148  {
149  b.checkBoundaryAfterParticlesMove(particleHandler);
150  }
151 }
152 
158 {
159  const DPMBase* const dpm = getHandler()->getDPMBase();
160  lengthBox_.X = (dpm->getXMax() - dpm->getXMin());
161  lengthBox_.Y = (dpm->getYMax() - dpm->getYMin());
162  lengthBox_.Z = (dpm->getZMax() - dpm->getZMin());
163  // Box length Lx, Lyand Lz, and center point C of the box
164  centerBox_.X = (dpm->getXMax() + dpm->getXMin()) / 2.0;
165  centerBox_.Y = (dpm->getYMax() + dpm->getYMin()) / 2.0;
166  centerBox_.Z = (dpm->getZMax() + dpm->getZMin()) / 2.0;
167 
168 }
169 
175 {
176  //Get the timestep dt
177  const Mdouble timeStep = getHandler()->getDPMBase()->getTimeStep();
178  // calculate the stress total and average over the volume
179  Matrix3D stressTotal = getHandler()->getDPMBase()->getTotalStress();
180 
181  // integral controller
182  // amount by which the strainrate tensor has to be changed
183  if (stressGoal_.XX != 0)
184  {
185  strainRate_.XX += gainFactor_.XX * timeStep * (stressTotal.XX - stressGoal_.XX);
186  logger(VERBOSE, "StressXX = %",stressTotal.XX);
187  }
188  if (stressGoal_.YY != 0)
189  {
190  strainRate_.YY += gainFactor_.YY * timeStep * (stressTotal.YY - stressGoal_.YY);
191  logger(VERBOSE, "StressYY = %",stressTotal.YY);
192  }
193  if (stressGoal_.ZZ != 0)
194  {
195  strainRate_.ZZ += gainFactor_.ZZ * timeStep * (stressTotal.ZZ - stressGoal_.ZZ);
196  logger(VERBOSE, "StressZZ = %",stressTotal.ZZ);
197  }
198  if (stressGoal_.XY != 0)
199  {
200  strainRate_.XY += gainFactor_.XY * timeStep * (stressTotal.XY - stressGoal_.XY);
201  logger(VERBOSE, "StressXY = %",stressTotal.XY);
202  logger(VERBOSE, "StressYX = %",stressTotal.YX);
203  }
204 }
205 
212 {
213  const DPMBase* const dpm = getHandler()->getDPMBase();
214  // update the domain size based on the new strainrate tensor
216 
217  // Move the boundaries to next time step
218  if (strainRate_.XY != 0 || stressGoal_.XY != 0)
219  {
221  }
222  else
223  {
224  // Move the boundary in x direction to next time step
225  periodicBoundaries_[0].set(Vec3D(1.0, 0.0, 0.0), dpm->getXMin(), dpm->getXMax());
226  // Move the boundary in y direction to next time step
227  periodicBoundaries_[1].set(Vec3D(0.0, 1.0, 0.0), dpm->getYMin(), dpm->getYMax());
228  // Move the boundary in z direction to next time step
229  periodicBoundaries_[2].set(Vec3D(0.0, 0.0, 1.0), dpm->getZMin(), dpm->getZMax());
230  }
231 
232  // Give the strain-rate for all particles and move them to next timestep before integration
234  {
235  for (auto& p : particleHandler)
236  {
237  relativeToCenter_.X = p->getPosition().X - centerBox_.X;
238  relativeToCenter_.Y = p->getPosition().Y - centerBox_.Y;
239  relativeToCenter_.Z = p->getPosition().Z - centerBox_.Z;
240  p->move(Vec3D(strainRate_.XX * dpm->getTimeStep() * relativeToCenter_.X +
244  }
245  }
246 }
247 
253 {
254  DPMBase* const dpm = getHandler()->getDPMBase();
255  // Change the system size according to next time step
256  dpm->setXMax(dpm->getXMax() + 0.5 * lengthBox_.X * strainRate_.XX * dpm->getTimeStep());
257  dpm->setXMin(dpm->getXMin() - 0.5 * lengthBox_.X * strainRate_.XX * dpm->getTimeStep());
258  dpm->setYMax(dpm->getYMax() + 0.5 * lengthBox_.Y * strainRate_.YY * dpm->getTimeStep());
259  dpm->setYMin(dpm->getYMin() - 0.5 * lengthBox_.Y * strainRate_.YY * dpm->getTimeStep());
260  dpm->setZMax(dpm->getZMax() + 0.5 * lengthBox_.Z * strainRate_.ZZ * dpm->getTimeStep());
261  dpm->setZMin(dpm->getZMin() - 0.5 * lengthBox_.Z * strainRate_.ZZ * dpm->getTimeStep());
262 
263  // Box length Lx, Lyand Lz for next time step
264  lengthBox_.X = (dpm->getXMax() - dpm->getXMin());
265  lengthBox_.Y = (dpm->getYMax() - dpm->getYMin());
266  lengthBox_.Z = (dpm->getZMax() - dpm->getZMin());
267 }
268 
275 {
276  const DPMBase* const dpm = getHandler()->getDPMBase();
277  //Determine the current shear velocity and shift of the Lees-Edwards boundary
278  const Mdouble velocity_xy = strainRate_.XY * lengthBox_.Y;
279  integratedShift_ += velocity_xy * dpm->getTimeStep();
280 
281  // Determine how to move boundaries based on strainrate control or boundary control
283  {
284  // Move the Lees-Edwards boundary in z direction to next time step
285  leesEdwardsBoundaries_[0].set(
286  [this](Mdouble time)
287  { return integratedShift_; },
288  [](Mdouble time UNUSED)
289  { return 0; },
290  dpm->getXMin(), dpm->getXMax(), dpm->getYMin(), dpm->getYMax());
291  }
292  else
293  {
294  // Update the velocity of Lees-Edwards boundary in x-y direction to next time step
295  leesEdwardsBoundaries_[0].set(
296  [this](Mdouble time)
297  { return integratedShift_; },
298  [velocity_xy](Mdouble time UNUSED)
299  { return velocity_xy; },
300  dpm->getXMin(), dpm->getXMax(), dpm->getYMin(), dpm->getYMax());
301  }
302  // Move the boundary in z direction to next time step
303  periodicBoundaries_[0].set(Vec3D(0.0, 0.0, 1.0), dpm->getZMin(), dpm->getZMax());
304 }
305 
316 void StressStrainControlBoundary::set(const Matrix3D& stressGoal, const Matrix3D& strainRate, const Matrix3D& gainFactor, bool isStrainRateControlled)
317 {
318  periodicBoundaries_.clear();
319  leesEdwardsBoundaries_.clear();
320 
321  isStrainRateControlled_ = isStrainRateControlled;
322  stressGoal_ = stressGoal;
323  strainRate_ = strainRate;
324  gainFactor_ = gainFactor;
325 
326  logger.assert_always(stressGoal.XZ == 0, "Shear stress in XZ cannot be controlled; use shear stress in XY instead");
327  logger.assert_always(stressGoal.ZX == 0, "Shear stress in ZX cannot be controlled; use shear stress in XY instead");
328  logger.assert_always(stressGoal.YZ == 0, "Shear stress in YZ cannot be controlled; use shear stress in XY instead");
329  logger.assert_always(stressGoal.ZY == 0, "Shear stress in ZY cannot be controlled; use shear stress in XY instead");
330  logger.assert_always(strainRate.XZ == 0, "Strain rate in XZ cannot be controlled; use strain rate in XY instead");
331  logger.assert_always(strainRate.ZX == 0, "Strain rate in ZX cannot be controlled; use strain rate in XY instead");
332  logger.assert_always(strainRate.YZ == 0, "Strain rate in YZ cannot be controlled; use strain rate in XY instead");
333  logger.assert_always(strainRate.ZY == 0, "Strain rate in ZY cannot be controlled; use strain rate in XY instead");
334  logger.assert_always(stressGoal.XZ != 0 ? gainFactor.XZ != 0 : true,
335  "You need to set a gain factor in XZ in order to control stress");
336  logger.assert_always(stressGoal.XX != 0 ? gainFactor.XX != 0 : true,
337  "You need to set a gain factor in XX in order to control stress");
338  logger.assert_always(stressGoal.YY != 0 ? gainFactor.YY != 0 : true,
339  "You need to set a gain factor in YY in order to control stress");
340  logger.assert_always(stressGoal.ZZ != 0 ? gainFactor.ZZ != 0 : true,
341  "You need to set a gain factor in ZZ in order to control stress");
342 
343  logger.assert_always(getHandler() != nullptr, "You need to set the handler of this boundary");
344  const DPMBase *dpm = getHandler()->getDPMBase();
345 
346  // if there is no shear (compression only)
347  if (stressGoal_.XY == 0 && strainRate_.XY == 0)
348  {
349  logger(INFO, "Shear rate is zero, setting up three periodic boundaries");
350  // Set up new box periodic boundaries, no simple shear activated
351  PeriodicBoundary boundary;
352  boundary.setHandler(getHandler());
353  boundary.set(Vec3D(1.0, 0.0, 0.0), dpm->getXMin(), dpm->getXMax());
354  periodicBoundaries_.push_back(boundary);
355  boundary.set(Vec3D(0.0, 1.0, 0.0), dpm->getYMin(), dpm->getYMax());
356  periodicBoundaries_.push_back(boundary);
357  boundary.set(Vec3D(0.0, 0.0, 1.0), dpm->getZMin(), dpm->getZMax());
358  periodicBoundaries_.push_back(boundary);
359  }
360  else
361  {
362  logger(INFO, "Shear rate is not zero, setting up Lees-Edwards in xy and a periodic boundary in z");
363 
364  PeriodicBoundary boundary;
365  boundary.setHandler(getHandler());
366  boundary.set(Vec3D(0.0, 0.0, 1.0), dpm->getZMin(), dpm->getZMax());
367  periodicBoundaries_.push_back(boundary);
368 
369  // Lees Edwards bc in y direction & periodic boundary in x direction, simple shear boundary activated
370  LeesEdwardsBoundary leesEdwardsBoundary;
371  leesEdwardsBoundary.setHandler(getHandler());
372  double integratedShift = integratedShift_;
373  leesEdwardsBoundary.set(
374  [&integratedShift](Mdouble time UNUSED) { return integratedShift; },
375  [](Mdouble time UNUSED) { return 0; },
376  dpm->getXMin(), dpm->getXMax(), dpm->getYMin(), dpm->getYMax());
377  leesEdwardsBoundaries_.push_back(leesEdwardsBoundary);
378  }
379 }
380 
386 {
387  // Call Boundaries
389  {
390  b.createPeriodicParticles(particleHandler);
391  }
393  {
394  b.createPeriodicParticles(particleHandler);
395  }
396 }
Mdouble XY
Definition: Matrix.h:43
void set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
Defines a PeriodicBoundary by its normal and positions.
void setXMax(Mdouble newXMax)
Sets the value of XMax, the upper bound of the problem domain in the x-direction. ...
Definition: DPMBase.cc:1126
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble integratedShift_
Shift integrated for all the time when using Lees-Edwards Boundary.
Class which creates a boundary with Lees-Edwards type periodic boundary conditions.
void checkBoundaryAfterParticlesMove(ParticleHandler &particleHandler) override
Virtual function that does things to particles, each timestep after particles have moved...
Mdouble ZY
Definition: Matrix.h:43
Vec3D lengthBox_
Box length in x-y-z.
void activateStrainRateControl(const ParticleHandler &particleHandler)
Activate the strainrate control for particle movement based on user's boolean input.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
std::vector< PeriodicBoundary > periodicBoundaries_
double Mdouble
Definition: GeneralDefine.h:34
void determineStressControlledShearBoundaries()
Determines stress-controlled shear Lees-Edwards boundary in x-y direction and normal periodic in z di...
void setYMin(Mdouble newYMin)
Sets the value of YMin, the lower bound of the problem domain in the y-direction. ...
Definition: DPMBase.cc:995
Mdouble XZ
Definition: Matrix.h:43
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.h:617
StressStrainControlBoundary()
default constructor.
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.h:586
std::vector< LeesEdwardsBoundary > leesEdwardsBoundaries_
Store boundaries into a vector for the pushback.
Mdouble ZX
Definition: Matrix.h:43
void updateDomainSize()
Update the domain to new sizes.
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.h:599
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
void setZMax(Mdouble newZMax)
Sets the value of ZMax, the upper bound of the problem domain in the z-direction. ...
Definition: DPMBase.cc:1178
void read(std::istream &is) override=0
Reads the object's id_ from given istream NB: purely virtual function, overriding the version of Base...
Definition: BaseBoundary.cc:61
Mdouble ZZ
Definition: Matrix.h:43
Vec3D relativeToCenter_
Particle position relative to the center of domain.
Mdouble YZ
Definition: Matrix.h:43
Defines a pair of periodic walls. Inherits from BaseBoundary.
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.h:593
Matrix3D getTotalStress() const
Calculate the total stress tensor in the system averaged over the whole volume.
Definition: DPMBase.cc:5166
void setYMax(Mdouble newYMax)
Sets the value of YMax, the upper bound of the problem domain in the y-direction. ...
Definition: DPMBase.cc:1152
void read(std::istream &is) override
Reads the object's id_ from given istream.
Matrix3D stressGoal_
Stores the stress value the boundary should attain.
void set(const Matrix3D &stressGoal, const Matrix3D &strainRate, const Matrix3D &gainFactor, bool isStrainRateControlled)
Sets all boundary inputs at once and determines which deformation mode it is, then combine the right ...
void computeStrainRate()
Compute the change of strainrate tensor based on the stress difference and then update the tensor...
void setHandler(BoundaryHandler *handler)
Sets the boundary's BoundaryHandler.
void setXMin(Mdouble newXMin)
Sets the value of XMin, the lower bound of the problem domain in the x-direction. ...
Definition: DPMBase.cc:971
Vec3D centerBox_
Center position of the domain.
Mdouble YX
Definition: Matrix.h:43
#define UNUSED
Definition: GeneralDefine.h:39
A cuboid box consists of periodic boundaries that can be strain/stress controlled and achieve differe...
void write(std::ostream &os) const override=0
Adds object's id_ to given ostream NB: purely virtual function, overriding the version of BaseObject...
Definition: BaseBoundary.cc:70
BoundaryHandler * getHandler() const
Returns the boundary's BoundaryHandler.
void setZMin(Mdouble newZMin)
Sets the value of ZMin, the lower bound of the problem domain in the z-direction. ...
Definition: DPMBase.cc:1019
StressStrainControlBoundary * copy() const override
Used to create a copy of the object.
Container to store all BaseParticle.
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.h:605
Mdouble Y
Definition: Vector.h:65
bool isStrainRateControlled_
The boolean input, true means switch on the strain rate control for particles affine movements...
std::string getName() const override
Sets the name of the boundary.
Mdouble XX
all nine matrix elements
Definition: Matrix.h:43
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:611
void setZero()
Sets all elements to zero.
Definition: Matrix.cc:75
void write(std::ostream &os) const override
Adds object's id_ to given ostream.
void set(std::function< Mdouble(Mdouble)> shift, std::function< Mdouble(Mdouble)> velocity, Mdouble left, Mdouble right, Mdouble down, Mdouble up)
Sets all boundary properties.
Mdouble YY
Definition: Matrix.h:43
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1211
void createPeriodicParticles(ParticleHandler &particleHandler) override
Create the periodic particles after read in from a restart file to attain right information.
Mdouble Z
Definition: Vector.h:65
void checkPeriodicLeesEdwardsBoundariesAfterParticlesMove(ParticleHandler &particleHandler)
Call the boundary and update them based on the new domain size after the stress-control movement...
void determineLengthAndCentre()
Determine the length in x,y,z and center location of domain.