MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StandardFields.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 #include "StandardFields.h"
27 #include <Particles/BaseParticle.h>
28 
29 namespace CGFields
30 {
31 
33 {
34  setZero();
35 #ifdef DEBUG_CONSTRUCTOR
36  std::cerr << "StandardFields::StandardFields() finished" << std::endl;
37 #endif
38 }
39 
44 void StandardFields::writeNames(std::ostream& os, const unsigned countVariables)
45 {
46  os << countVariables + 1 << ":VolumeFraction "; //volume
47  os << countVariables + 2 << ":Density "; //mass
48  os << countVariables + 3 << '-' << countVariables + 5 << ":Momentum ";
49  os << countVariables + 6 << '-' << countVariables + 11 << ":MomentumFlux ";
50  os << countVariables + 12 << '-' << countVariables + 20 << ":ContactStress ";
51  os << countVariables + 21 << '-' << countVariables + 23 << ":InteractionForce ";
52  os << countVariables + 24 << '-' << countVariables + 29 << ":ParticleSize ";
53 }
54 
58 void StandardFields::write(std::ostream& os) const
59 {
60  os << volumeFraction_;
61  os << ' ' << density_;
62  os << ' ' << momentum_;
63  os << ' ' << momentumFlux_;
64  os << ' ' << contactStress_;
65  os << ' ' << interactionForceDensity_;
66  for (const auto p : particleSizeDensity_) os << ' ' << p;
67 }
68 
72 void StandardFields::output(std::ostream& os) const
73 {
74  os << "VolumeFraction " << volumeFraction_;
75  os << " Density " << density_;
76  os << " Momentum " << momentum_;
77  os << " MomentumFlux " << momentumFlux_;
78  os << " ContactStress " << contactStress_;
79  os << " InteractionForce " << interactionForceDensity_;
80  os << " ParticleSize";
81  for (const auto p : particleSizeDensity_) os << ' ' << p;
82 }
83 
85 {
86  volumeFraction_ = 0.0;
87  density_ = 0.0;
92  for (auto& p : particleSizeDensity_) p = 0;
93 }
94 
99 {
100  StandardFields P;
108  for (auto& p : P.particleSizeDensity_) p *= p;
109  return P;
110 }
111 
117 = default;
118 
124 {
126  density_ += P.density_;
127  momentum_ += P.momentum_;
131  for (size_t i = 0; i < particleSizeDensity_.size(); ++i) particleSizeDensity_[i] += P.particleSizeDensity_[i];
132  return *this;
133 }
134 
140 {
142  density_ -= P.density_;
143  momentum_ -= P.momentum_;
147  for (size_t i = 0; i < particleSizeDensity_.size(); ++i) particleSizeDensity_[i] -= P.particleSizeDensity_[i];
148  return *this;
149 }
150 
156 {
157  StandardFields p;
159  p.density_ = density_ * a;
160  p.momentum_ = momentum_ * a;
164  for (size_t i = 0; i < particleSizeDensity_.size(); ++i)
166  return p;
167 }
168 
174 {
175  volumeFraction_ /= a;
176  density_ /= a;
177  momentum_ /= a;
178  momentumFlux_ /= a;
179  contactStress_ /= a;
181  for (auto& p : particleSizeDensity_) p /= a;
182  return *this;
183 }
184 
190 {
191  volumeFraction_ += currentInteraction.getVolumeFraction() * phi;
192  density_ += currentInteraction.getDensity() * phi;
193  momentum_ += currentInteraction.getMomentum() * phi;
194  momentumFlux_ += currentInteraction.getMomentumFlux() * phi;
195  for (size_t i = 0; i < particleSizeDensity_.size(); ++i)
196  particleSizeDensity_[i] += currentInteraction.getParticleSizeDensity(i) * phi;
197 }
198 
200 {
201 }
202 
208 {
209  contactStress_ += currentInteraction.getContactStress() * psi;
210 }
211 
218 {
219  interactionForceDensity_ += currentInteraction.getInteractionForceDensity() * phi;
220 }
221 
229 {
230  return true;
231 }
232 
234 {
235  //P is always real, but I not
236  if (type == IntegralType::I_TO_P)
237  {
239  }
240  else if (type == IntegralType::I_TO_CONTACT)
241  {
243  }
244  else
245  {
247  }
249 }
250 
252 {
254  density_ = p.getMass();
255  momentum_ = p.getVelocity() * p.getMass();
257  Mdouble particleSize = 1.0;
258  for (auto& ps : particleSizeDensity_)
259  {
260  ps = particleSize;
261  particleSize *= p.getRadius();
262  }
263 }
264 
266 {
267  setFields(c, type);
270 }
271 
273 {
274  setFields(p);
277 }
278 
280 {
281  //https://en.wikipedia.org/wiki/Moment_%28mathematics%29#Higher_moments
282  auto momenta = getStandardisedParticleSizeMomenta();
283  os << "Particle number: " << momenta[0] << '\n';
284  os << "Mean radius: " << momenta[1] << '\n';
285  os << "Standard deviation: " << std::sqrt(momenta[2]) << '\n';
286  os << "Skewness: " << momenta[3] << '\n'; //left- or right-sided
287  os << "Kurtosis: " << momenta[4] << '\n'; //heavy or light tailed (3 for normal distribution)
288  os << "5-th moment: " << momenta[5] << '\n'; //hyper skewness (0 for normal dist)
289 }
290 
291 std::array<Mdouble, 6> StandardFields::getParticleSizeMomenta() const
292 {
293  //https://en.wikipedia.org/wiki/Moment_(mathematics)
294  auto momenta = particleSizeDensity_;
295  for (size_t i = 1; i < particleSizeDensity_.size(); ++i)
296  {
297  momenta[i] /= particleSizeDensity_[0];
298  }
299  return momenta;
300 }
301 
302 std::array<Mdouble, 6> StandardFields::getCentralParticleSizeMomenta() const
303 {
304  //https://en.wikipedia.org/wiki/Central_moment
305  auto momenta = getParticleSizeMomenta();
306  Mdouble mean = momenta[1];
307  momenta[5] += -5 * mean * momenta[4] + 10 * mean * mean * momenta[3]
308  - 10 * mean * mean * mean * momenta[2] + 4 * mean * mean * mean * mean * mean;
309  momenta[4] += -4 * mean * momenta[3] + 6 * mean * mean * momenta[2] - 3 * mean * mean * mean * mean;
310  momenta[3] += -3 * mean * momenta[2] + 2 * mean * mean * mean;
311  momenta[2] += -mean * mean;
312  return momenta;
313 }
314 
316 {
317  auto momenta = getCentralParticleSizeMomenta();
318  Mdouble std = std::sqrt(momenta[2]);
319  momenta[3] /= std * std * std;
320  momenta[4] /= std * std * std * std;
321  momenta[5] /= std * std * std * std * std;
322  return momenta;
323 }
324 
325 }
static MatrixSymmetric3D square(const MatrixSymmetric3D &A)
Calculates the pointwise square.
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
void write(std::ostream &os) const
Writes class content into an output stream, typically a stat file.
virtual Mdouble getVolume() const
Get Particle volume function, which required a reference to the Species vector. It returns the volume...
Mdouble getDensity() const
double Mdouble
Definition: GeneralDefine.h:34
std::array< Mdouble, 6 > getCentralParticleSizeMomenta() const
Vec3D getInteractionForceDensity() const
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
StandardFields & operator/=(Mdouble a)
Divides the field values on the LHS by the RHS of the equation.
Vec3D getCylindricalTensorField(const Vec3D &position) const
Returns this vector field at point p to cylindrical coordinates.
Definition: Vector.cc:261
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
Matrix3D getContactStress() const
static Vec3D square(const Vec3D &a)
Calculates the pointwise square of a Vec3D.
Definition: Vector.cc:114
Mdouble getParticleSizeDensity(size_t i) const
Stores information about interactions between two interactable objects; often particles but could be ...
static Matrix3D dyadic(const Vec3D &a, const Vec3D &b)
Calculates the dyadic product of a two Vec3D: .
Definition: Matrix.cc:323
MatrixSymmetric3D momentumFlux_
StandardFields & operator=(const StandardFields &P)
Copies all field values.
IntegralType
Specifies the two points between which a line integral is computed.
Definition: IntegralType.h:33
void addParticleStatistics(Mdouble phi, const StandardFields &currentInteraction)
This function should be called from within a loop over all particles to compute all the fields that a...
Mdouble getMass() const
Returns the particle's mass.
Definition: BaseParticle.h:322
Mdouble getVolumeFraction() const
std::array< Mdouble, 6 > getParticleSizeMomenta() const
StandardFields operator*(Mdouble a) const
Multiplies the field values on the left of the '*' by the scalar value on the right of the '*' and re...
Vec3D getMomentum() const
void setZero()
Sets all elements to zero.
void addInteractionStatistics(Mdouble psi, const StandardFields &currentInteraction)
This function should be called from within a loop over all Interactions to compute all the fields tha...
const Vec3D & getForce() const
Gets the current force (vector) between the two interacting objects.
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
Vec3D getIP() const
void outputStandardisedParticleSizeMomenta(std::ostream &os) const
Matrix3D getCylindricalTensorField(const Vec3D &p) const
Returns the matrix in cylindrical coordinates.
Definition: Matrix.cc:394
Vec3D getIC() const
static MatrixSymmetric3D selfDyadic(const Vec3D &a)
Calculates the dyadic product of a Vec3D with itself: .
Vec3D getCP() const
void setZero()
Sets all fields to zero.
StandardFields & operator+=(const StandardFields &P)
Adds the field values on the RHS to the LHS of the equation.
std::array< Mdouble, 6 > particleSizeDensity_
void setZero()
Sets all elements to zero.
Definition: Matrix.cc:75
void addParticleDifferentialStatistics(Vec3D &dphi, const StandardFields &currentInteraction)
std::array< Mdouble, 6 > getStandardisedParticleSizeMomenta() const
StandardFields & operator-=(const StandardFields &P)
Subtracts the field values on the RHS from the LHS of the equation.
void output(std::ostream &os) const
Writes human-readable class content into an output stream, typically a stat file. ...
MatrixSymmetric3D getCylindricalTensorField(const Vec3D &p) const
Returns the matrix in cylindrical coordinates.
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Definition: Vector.h:49
void setCylindricalFields(const BaseInteraction &c, IntegralType type)
MatrixSymmetric3D getMomentumFlux() const
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
Contains the computed field values, like density, momentum and stress.
void addContactPointStatistics(Mdouble phi, const StandardFields &currentInteraction)
This function should be called from within a loop over all Interactions to compute all the fields tha...
static bool doInteractionStatistics()
Returns true if the class contains fields that are defined as a sum over all Interactions (e...
static void writeNames(std::ostream &os, unsigned countVariables)
static Matrix3D square(const Matrix3D &A)
Calculates the pointwise square.
Definition: Matrix.cc:298
StandardFields getSquared() const
Returns the square of all field values (to calculate standard deviation).
StandardFields()
Default constructor, sets all field values to zero.
void setFields(const BaseInteraction &c, IntegralType type)