ESyS-Particle  2.3.2
BodyForceGroup.hpp
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 #include "ppa/src/pp_array.h"
15 #include <math.h>
16 
17 namespace esys
18 {
19  namespace lsm
20  {
21  template <class TmplParticle>
23  const BodyForceIGP &prms,
24  ParticleArray &particleArray
25  )
26  : m_acceleration(prms.getAcceleration()),
27  m_pParticleArray(&particleArray)
28  {
29  }
30 
31  template <class TmplParticle>
33  {
34  }
35 
36  template <class TmplParticle>
38  {
39  }
40 
41  template <class TmplParticle>
43  {
44  return m_acceleration*mass;
45  }
46 
47  template <class TmplParticle>
48  void BodyForceGroup<TmplParticle>::applyForce(TmplParticle &particle) const
49  {
50  particle.applyForce(getForce(particle.getMass()), particle.getPos());
51  }
52 
53  template <class TmplParticle>
55  {
56  typename ParticleArray::ParticleListHandle plh = m_pParticleArray->getAllParticles();
57  for (
58  ParticleIterator it = plh->begin();
59  it != plh->end();
60  it++
61  )
62  {
63  applyForce(*(*it));
64  }
65  }
66 
67  template <class TmplParticle>
69  const BuoyancyIGP &prms,
70  ParticleArray &particleArray
71  )
72  : m_acceleration(prms.getAcceleration()),
73  m_pParticleArray(&particleArray) ,
74  m_fluidDensity(prms.getFluidDensity()),
75  m_fluidHeight(prms.getFluidHeight())
76  {
77  }
78 
79  template <class TmplParticle>
81  {
82  }
83 
84  template <class TmplParticle>
86  {
87  Vec3 force;
88 
89  force = -1.0*m_acceleration*m_fluidDensity*volume;
90 
91  return force;
92  }
93 
94  template <class TmplParticle>
95  void BuoyancyForceGroup<TmplParticle>::applyForce(TmplParticle &particle) const
96  {
97  double volume, radius, height;
98  double theta, d;
99 
100  radius = particle.getRad();
101  height = -1.0*particle.getPos() * m_acceleration.unit();
102  if(particle.getDo2dCalculations()){
103  if (m_fluidHeight >= height + radius) {
104  // completely submerged
105  volume = M_PI*radius*radius;
106  }
107  else if (m_fluidHeight >= height - radius) {
108  // partially submerged
109  if (m_fluidHeight < height) {
110  d = height - m_fluidHeight;
111  theta = 2.0 * acos (d / radius);
112  volume = 0.5*radius*radius*(theta - sin (theta));
113  }
114  else {
115  d = m_fluidHeight - height;
116  theta = 2.0 * acos (d / radius);
117  volume = M_PI*radius*radius*(M_PI - 0.5*(theta - sin (theta)));
118  }
119  }
120  else {
121  // not submerged
122  volume = 0.0;
123  }
124  }
125  else {
126  if (m_fluidHeight >= height + radius) {
127  // completely submerged
128  volume = 4.0*M_PI*radius*radius*radius/3.0;
129  }
130  else if (m_fluidHeight >= height - radius) {
131  // partially submerged
132  if (m_fluidHeight < height) {
133  d = radius - (height - m_fluidHeight);
134  volume = M_PI*d*d*(3.0*radius - d)/3.0;
135  }
136  else {
137  d = radius - (m_fluidHeight - height);
138  volume = M_PI*(4.0*radius*radius*radius - d*d*(3.0*radius - d))/3.0;
139  }
140  }
141  else {
142  // not submerged
143  volume = 0.0;
144  }
145  }
146  particle.applyForce(getForce(volume), particle.getPos());
147  }
148 
149  template <class TmplParticle>
151  {
152  typename ParticleArray::ParticleListHandle plh = m_pParticleArray->getAllParticles();
153  for (
154  ParticleIterator it = plh->begin();
155  it != plh->end();
156  it++
157  )
158  {
159  applyForce(*(*it));
160  }
161  }
162 
163  template <class TmplParticle>
165  {
166  }
167  }
168 }