1 /* This file is part of StepCore library.
2    Copyright (C) 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
3 
4    StepCore library is free software; you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation; either version 2 of the License, or
7    (at your option) any later version.
8 
9    StepCore library is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13 
14    You should have received a copy of the GNU General Public License
15    along with StepCore; if not, write to the Free Software
16    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
17 */
18 
19 /** \file gas.h
20  *  \brief Gas-related classes
21  */
22 
23 #ifndef STEPCORE_GAS_H
24 #define STEPCORE_GAS_H
25 
26 #include "particle.h"
27 #include "world.h"
28 #include <cmath>
29 
30 namespace StepCore {
31 
32 class GasParticle;
33 class GasLJForce;
34 class Gas;
35 
36 /** \ingroup bodies
37  *  \brief Gas particle
38  */
39 class GasParticle: public Particle
40 {
STEPCORE_OBJECT(GasParticle)41     STEPCORE_OBJECT(GasParticle)
42 
43 public:
44     /** Constructs a GasParticle */
45     explicit GasParticle(const Vector2d &position = Vector2d::Zero(), const Vector2d &velocity = Vector2d::Zero(), double mass = 1)
46         : Particle(position, velocity, mass) {}
47 };
48 
49 /** \ingroup errors
50  *  \brief Errors object for GasLJForce
51  */
52 class GasLJForceErrors: public ObjectErrors
53 {
STEPCORE_OBJECT(GasLJForceErrors)54     STEPCORE_OBJECT(GasLJForceErrors)
55 
56 public:
57     /** Constructs GasLJForceErrors */
58     explicit GasLJForceErrors(Item* owner = 0)
59         : ObjectErrors(owner), _depthVariance(0), _rminVariance(0) {}
60 
61     /** Get owner as GasLJForce */
62     GasLJForce* gasLJForce() const;
63 
64     /** Get depth variance */
depthVariance()65     double depthVariance() const { return _depthVariance; }
66     /** Set depth variance */
setDepthVariance(double depthVariance)67     void setDepthVariance(double depthVariance) { _depthVariance = depthVariance; }
68 
69     /** Get rmin variance */
rminVariance()70     double rminVariance() const { return _rminVariance; }
71     /** Set rmin variance */
setRminVariance(double rminVariance)72     void setRminVariance(double rminVariance) { _rminVariance = rminVariance; }
73 
74 protected:
75     double _depthVariance;
76     double _rminVariance;
77     friend class GasLJForce;
78 };
79 
80 /** \ingroup forces
81  *  \brief Lennard-Jones force with cut-off which acts between particles in the Gas
82  *
83  *  The force acts between pairs of GasParticle and equals:
84  *  \f{eqnarray*}
85  *      \overrightarrow{f} = & 12 \epsilon \left(
86  *             \frac{ r_{min}^{12} }{ r^{13} } -
87  *             \frac{ r_{min}^{6} }{ r^{7} }
88  *         \right) \frac{\overrightarrow{r}}{r} & \mbox{  if  } r<\mbox{cutoff} \\
89  *     \overrightarrow{f} = & 0 & \mbox{  if  } r \ge \mbox{cutoff}
90  *  \f}
91  *  where:\n
92  *  \f$\epsilon\f$ is the depth of the potential\n
93  *  \f$r_{min}\f$ is the distance at which the interparticle force is zero\n
94  *  \f$\overrightarrow{r}\f$ is difference of GasParticle::position
95                              of the first and second particle\n
96  *  \f$\mbox{cutoff}\f$ is a cut-off distance (can be set to infinity)
97  *
98  */
99 class GasLJForce : public Force
100 {
101     STEPCORE_OBJECT(GasLJForce)
102 
103 public:
104     /** Constructs GasLJForce */
105     explicit GasLJForce(double depth = 1, double rmin = 1, double cutoff = HUGE_VAL);
106 
107     void calcForce(bool calcVariances) override;
108 
109     /** Get depth of the potential */
depth()110     double depth() const { return _depth; }
111     /** Set depth of the potential */
setDepth(double depth)112     void setDepth(double depth) { _depth = depth; calcABC(); }
113 
114     /** Get distance at which the interparticle force is zero */
rmin()115     double rmin() const { return _rmin; }
116     /** Set distance at which the interparticle force is zero */
setRmin(double rmin)117     void setRmin(double rmin) { _rmin = rmin; calcABC(); }
118 
119     /** Get cut-off distance */
cutoff()120     double cutoff() const { return _cutoff; }
121     /** Set cut-off distance */
setCutoff(double cutoff)122     void setCutoff(double cutoff) { _cutoff = cutoff; calcABC(); }
123 
124     /** Get (and possibly create) GasLJForceErrors object */
gasLJForceErrors()125     GasLJForceErrors* gasLJForceErrors() {
126         return static_cast<GasLJForceErrors*>(objectErrors()); }
127 
128 protected:
createObjectErrors()129     ObjectErrors* createObjectErrors() override { return new GasLJForceErrors(this); }
130     void calcABC();
131 
132     double _depth;
133     double _rmin;
134     double _cutoff;
135     double _a, _b, _c;
136     double _rmin6, _rmin12;
137 };
138 
139 typedef std::vector<GasParticle*> GasParticleList;
140 
141 /** \ingroup errors
142  *  \brief Errors object for Gas
143  */
144 class GasErrors: public ObjectErrors
145 {
STEPCORE_OBJECT(GasErrors)146     STEPCORE_OBJECT(GasErrors)
147 
148 public:
149     /** Constructs GasErrors */
150     explicit GasErrors(Item* owner = 0)
151         : ObjectErrors(owner) {}
152 
153     /** Get owner as Gas */
154     Gas* gas() const;
155 
156     double rectTemperatureVariance() const;
157     double rectPressureVariance() const;
158     Vector2d rectMeanVelocityVariance() const;
159     double rectMeanKineticEnergyVariance() const;
160     double rectMeanParticleMassVariance() const;
161     double rectMassVariance() const;
162 
163 protected:
164     friend class Gas;
165 };
166 
167 /** \ingroup bodies
168  *  \brief Gas - a group of several GasParticle and a force
169  */
170 class Gas: public ItemGroup
171 {
STEPCORE_OBJECT(Gas)172     STEPCORE_OBJECT(Gas)
173 
174 public:
175     Gas() : _measureRectCenter(0,0), _measureRectSize(1,1) {
176         setColor(0xffff0000); objectErrors();
177     }
178 
179     /** Creates particles with given temperature
180      *  \todo XXX Normalize temperature after particle creation */
181     GasParticleList rectCreateParticles(int count,
182                                 double mass, double temperature,
183                                 const Vector2d& meanVelocity);
184 
185     void addParticles(const GasParticleList& particles);
186 
187     double rectVolume() const;
188     double rectParticleCount() const;
189     double rectConcentration() const;
190     double rectTemperature() const;
191     double rectPressure() const;
192     Vector2d rectMeanVelocity() const;
193     double rectMeanKineticEnergy() const;
194     double rectMeanParticleMass() const;
195     double rectMass() const;
196 
measureRectCenter()197     const Vector2d& measureRectCenter() const { return _measureRectCenter; }
setMeasureRectCenter(const Vector2d & measureRectCenter)198     void setMeasureRectCenter(const Vector2d& measureRectCenter) { _measureRectCenter = measureRectCenter; }
199 
measureRectSize()200     const Vector2d& measureRectSize() const { return _measureRectSize; }
setMeasureRectSize(const Vector2d & measureRectSize)201     void setMeasureRectSize(const Vector2d& measureRectSize) { _measureRectSize = measureRectSize.array().abs().matrix(); }
202 
203     /** Get (and possibly create) GasErrors object */
gasErrors()204     GasErrors* gasErrors() {
205         return static_cast<GasErrors*>(objectErrors()); }
206 
207 protected:
createObjectErrors()208     ObjectErrors* createObjectErrors() override { return new GasErrors(this); }
209 
210     double randomUniform(double min=0, double max=1);
211     double randomGauss(double mean=0, double deviation=1);
212 
213     Vector2d _measureRectCenter;
214     Vector2d _measureRectSize;
215 
216     friend class GasErrors;
217 };
218 
219 } // namespace StepCore
220 
221 #endif
222 
223