← Back to branch summary

~ubuntu-branches/ubuntu/wily/yade/wily

« back to all changes in this revision

Viewing changes to pkg/common/SPHEngine.cpp

  • Committer: Package Import Robot
  • Author(s): Dmitry Shachnev
  • Date: 2014-11-14 12:54:52 UTC
  • mfrom: (20.1.23 sid)
  • Revision ID: package-import@ubuntu.com-20141114125452-t16anreumu4ybg2s
Tags: 1.12.0-2ubuntu1
Add allow-stderr restriction to autopkgtest, to silence a warning
printed by new matplotlib.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/common/SPHEngine.cpp
1
1
#ifdef YADE_SPH
2
2
#include"SPHEngine.hpp"
3
 
#include<yade/core/Scene.hpp>
4
 
#include<yade/pkg/dem/ViscoelasticPM.hpp>
5
 
#include<yade/pkg/common/Sphere.hpp>
 
3
#include<core/Scene.hpp>
 
4
#include<pkg/dem/ViscoelasticPM.hpp>
 
5
#include<pkg/common/Sphere.hpp>
6
6
 
7
 
#include<yade/core/State.hpp>
8
 
#include<yade/core/Omega.hpp>
 
7
#include<core/State.hpp>
 
8
#include<core/Omega.hpp>
9
9
 
10
10
void SPHEngine::action(){
11
11
  {
12
12
    YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
13
13
      if(mask>0 && (b->groupMask & mask)==0) continue;
14
14
      this->calculateSPHRho(b);
15
 
      b->press=k*(b->rho - b->rho0);
16
 
    } YADE_PARALLEL_FOREACH_BODY_END();
17
 
  }
18
 
  /*
19
 
  {
20
 
    YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
21
 
      if(mask>0 && (b->groupMask & mask)==0) continue;
22
 
      this->calculateSPHCs(b);
23
 
    } YADE_PARALLEL_FOREACH_BODY_END();
24
 
  }
25
 
  */
 
15
      b->press=std::max(0.0, k*(b->rho - b->rho0));
 
16
    } YADE_PARALLEL_FOREACH_BODY_END();
 
17
  }
26
18
}
27
19
 
28
20
void SPHEngine::calculateSPHRho(const shared_ptr<Body>& b) {
32
24
  Real rho = 0;
33
25
  
34
26
  // Pointer to kernel function
35
 
  KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, KernFunctionDensity, Norm);
 
27
  KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, Norm);
36
28
  
37
29
  // Calculate rho for every particle
38
30
  for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
49
41
      Real Mass = b2->state->mass;
50
42
      if (Mass == 0) Mass = b->state->mass;
51
43
      
52
 
      const Real SmoothDist = -geom.penetrationDepth + phys.h;
 
44
      const Real SmoothDist = (b2->state->pos - b->state->pos).norm();
53
45
     
54
 
      // [Mueller2003], (3)
55
 
      rho += b2->state->mass*kernelFunctionCurDensity(SmoothDist, phys.h);
 
46
      // [Monaghan1992], (2.7) (3.8) 
 
47
      rho += b2->state->mass*kernelFunctionCurDensity(SmoothDist, h);
56
48
    }
57
 
    // [Mueller2003], (3), we need to consider the density of the current body (?)
58
 
    rho += b->state->mass*kernelFunctionCurDensity(0.0, s->radius);
59
49
  }
 
50
  // Self mass contribution
 
51
  rho += b->state->mass*kernelFunctionCurDensity(0.0, h);
60
52
  b->rho = rho;
61
53
}
62
54
 
63
 
void SPHEngine::calculateSPHCs(const shared_ptr<Body>& b) {
64
 
  if (b->Cs<0) {
65
 
    b->Cs = 0.0;
66
 
  }
67
 
  Real Cs = 0;
68
 
  
69
 
  // Pointer to kernel function
70
 
  KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, KernFunctionDensity, Norm);
71
 
  
72
 
  // Calculate Cs for every particle
73
 
  for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
74
 
    const shared_ptr<Body> b2 = Body::byId((*it).first,scene);
75
 
    Sphere* s=dynamic_cast<Sphere*>(b->shape.get());
76
 
    if(!s) continue;
77
 
    
78
 
    if (((*it).second)->geom and ((*it).second)->phys) {
79
 
      const ScGeom geom = *(YADE_PTR_CAST<ScGeom>(((*it).second)->geom));
80
 
      const ViscElPhys phys=*(YADE_PTR_CAST<ViscElPhys>(((*it).second)->phys));
81
 
      
82
 
      if((b2->groupMask & mask)==0)  continue;
83
 
      
84
 
      Real Mass = b2->state->mass;
85
 
      if (Mass == 0) Mass = b->state->mass;
86
 
      
87
 
      const Real SmoothDist = -geom.penetrationDepth + phys.h;
88
 
     
89
 
      // [Mueller2003], (15)
90
 
      Cs += b2->state->mass/b2->rho*kernelFunctionCurDensity(SmoothDist, phys.h);
91
 
    }
92
 
  }
93
 
  b->Cs = Cs;
94
 
}
95
 
 
96
 
Real smoothkernelPoly6(const double & rr, const double & hh) {
97
 
  if (rr<=hh) {
98
 
    const Real h = 1; const Real r = rr/hh;
99
 
    //return 315/(64*M_PI*pow(h,9))*pow((h*h - r*r), 3);
100
 
    return 315/(64*M_PI)*pow((h - r*r), 3);
101
 
  } else {
102
 
    return 0;
103
 
  }
104
 
}
105
 
 
106
 
Real smoothkernelPoly6Grad(const double & rr, const double & hh) {
107
 
  if (rr<=hh) {
108
 
    const Real h = 1; const Real r = rr/hh;
109
 
    //return -315/(64*M_PI*pow(h,9))*(-6*r*pow((h*h-r*r), 2));
110
 
    return -315/(64*M_PI)*(-6*r*pow((h-r*r), 2));
111
 
  } else {
112
 
    return 0;
113
 
  }
114
 
}
115
 
 
116
 
Real smoothkernelPoly6Lapl(const double & rr, const double & hh) {
117
 
  if (rr<=hh) {
118
 
    const Real h = 1; const Real r = rr/hh;
119
 
    //return 315/(64*M_PI*pow(h,9))*(-6*(h*h-r*r)*(h*h - 5*r*r));
120
 
    return 315/(64*M_PI)*(-6*(h*h-r*r)*(h*h - 5*r*r));
121
 
  } else {
122
 
    return 0;
123
 
  }
124
 
}
125
 
 
126
 
Real smoothkernelSpiky(const double & rr, const double & hh) {
127
 
  if (rr<=hh) {
128
 
    const Real h = 1; const Real r = rr/hh;
129
 
    //return 15/(M_PI*pow(h,6))*(pow((h-r), 3));             // [Mueller2003], (21)
130
 
    return 15/(M_PI)*(pow((h-r), 3));             // [Mueller2003], (21)
131
 
  } else {
132
 
    return 0;
133
 
  }
134
 
}
135
 
 
136
 
Real smoothkernelSpikyGrad(const double & rr, const double & hh) {
137
 
  if (rr<=hh) {
138
 
    const Real h = 1; const Real r = rr/hh;
139
 
    //return -15/(M_PI*pow(h,6))*(-3*pow((h-r),2));
140
 
    return -15/(M_PI)*(-3*pow((h-r),2));
141
 
  } else {
142
 
    return 0;
143
 
  }
144
 
}
145
 
 
146
 
Real smoothkernelSpikyLapl(const double & rr, const double & hh) {
147
 
  if (rr<=hh) {
148
 
    const Real h = 1; const Real r = rr/hh;
149
 
    //return 15/(M_PI*pow(h,6))*(6*(h-r));
150
 
    return 15/(M_PI)*(6*(h-r));
151
 
  } else {
152
 
    return 0;
153
 
  }
154
 
}
155
 
 
156
 
Real smoothkernelVisco(const double & rr, const double & hh) {
157
 
  if (rr<=hh and rr!=0 and hh!=0) {
158
 
    const Real h = 1; const Real r = rr/hh;
159
 
    //return 15/(2*M_PI*pow(h,3))*(-r*r*r/(2*h*h*h) + r*r/(h*h) + h/(2*r) -1);   // [Mueller2003], (21)
160
 
    return 15/(2*M_PI)*(-r*r*r/(2) + r*r + h/(2*r) -1);   // [Mueller2003], (21)
161
 
  } else {
162
 
    return 0;
163
 
  }
164
 
}
165
 
 
166
 
Real smoothkernelViscoGrad(const double & rr, const double & hh) {
167
 
  if (rr<=hh and rr!=0 and hh!=0) {
168
 
    const Real h = 1; const Real r = rr/hh;
169
 
    //return -15/(2*M_PI*pow(h,3))*(-3*r*r/(2*h*h*h) + 2*r/(h*h) - h/(2*r*r));
170
 
    return -15/(2*M_PI)*(-3*r*r/(2) + 2*r - h/(2*r*r));
171
 
  } else {
172
 
    return 0;
173
 
  }
174
 
}
175
 
 
176
 
Real smoothkernelViscoLapl(const double & rr, const double & hh) {
177
 
  if (rr<=hh and rr!=0 and hh!=0) {
178
 
    const Real h = 1; const Real r = rr/hh;
179
 
    //return 45/(M_PI*pow(h,6))*(h - rrj);                     // [Mueller2003], (22+)
180
 
    //return 15/(2*M_PI*pow(h,3))*(-3*r*r/(2*h*h*h) + 2*r/(h*h) - h/(2*r*r));
181
 
    return 15/(2*M_PI)*(-3*r*r/(2) + 2*r - h/(2*r*r));
182
 
  } else {
183
 
    return 0;
184
 
  }
185
 
}
186
 
 
187
 
Real smoothkernelLucy(const double & rr, const double & hh) {
188
 
  if (rr<=hh and rr!=0 and hh!=0) {
189
 
    //const Real h = 1; 
190
 
    //return 5/(9*M_PI*pow(h,2))*(1+3*r/h)*pow((1-r/h),3);
191
 
    const Real r = rr/hh;
192
 
    return 5/(9*M_PI)*(1+3*r)*pow((1-r),3);
193
 
  } else {
194
 
    return 0;
195
 
  }
196
 
}
197
 
 
198
 
Real smoothkernelLucyGrad(const double & rr, const double & hh) {
199
 
  if (rr<=hh and rr!=0 and hh!=0) {
200
 
    //const Real h = 1; 
201
 
    //return -5/(9*M_PI*pow(h,2))*(-12*r/(h*h))*pow((1-r/h),2);
202
 
    const Real r = rr/hh;
203
 
    return -5/(9*M_PI)*(-12*r)*pow((1-r),2);
204
 
  } else {
205
 
    return 0;
206
 
  }
207
 
}
208
 
 
209
 
Real smoothkernelLucyLapl(const double & rr, const double & hh) {
210
 
  if (rr<=hh and rr!=0 and hh!=0) {
211
 
    //const Real h = 1; 
212
 
    //return  5/(9*M_PI*pow(h,2))*(-12/(h*h))*(1-r/h)*(1-3*r/h);
213
 
    const Real r = rr/hh;
214
 
    return  5/(9*M_PI)*(-12)*(1-r)*(1-3*r);
215
 
  } else {
216
 
    return 0;
217
 
  }
218
 
}
219
 
 
220
 
Real smoothkernelMonaghan(const double & rr, const double & hh) {
221
 
  Real ret = 0.0;
222
 
  if (hh!=0) {
223
 
    //const Real h = 1; 
224
 
    const Real r = rr/hh;
225
 
    if (rr/hh<0.5) {
226
 
      //ret = 40/(7*M_PI*h*h)*(1 - 6*pow((r/h),2) + 6*pow((r/h),3));
227
 
      ret = 40/(7*M_PI)*(1 - 6*pow((r),2) + 6*pow((r),3));
228
 
    } else {
229
 
      //ret = 80/(7*M_PI*h*h)*pow((1 - (r/h)), 3);
230
 
      ret = 80/(7*M_PI)*pow((1 - (r)), 3);
231
 
    }
232
 
  }
233
 
  return ret;
234
 
}
235
 
 
236
 
Real smoothkernelMonaghanGrad(const double & rr, const double & hh) {
237
 
  Real ret = 0.0;
238
 
  if (hh!=0) {
239
 
    const Real h = 1; const Real r = rr/hh;
240
 
    if (rr/hh<0.5) {
241
 
      //ret = -40/(7*M_PI*h*h)*( - 6*r/(h*h))*(2 - 3 * r/(h*h*h));
242
 
      ret = -40/(7*M_PI)*( - 6*r)*(2 - 3 * r);
243
 
    } else {
244
 
      //ret = -80/(7*M_PI*h*h)*( -3/h)*pow((1 -r/h), 2);
245
 
      ret = -80/(7*M_PI)*( -3/h)*pow((1 -r), 2);
246
 
    }
247
 
  }
248
 
  return ret;
249
 
}
250
 
 
251
 
Real smoothkernelMonaghanLapl(const double & rr, const double & hh) {
252
 
  Real ret = 0.0;
253
 
  if (hh!=0) {
254
 
    //const Real h = 1; 
255
 
    const Real r = rr/hh;
256
 
    if (rr/hh<0.5) {
257
 
      //ret = 40/(7*M_PI*h*h)*( - 12/(h*h))*(1 - 3 * r/(h*h*h*h*h));
258
 
      ret = 40/(7*M_PI)*( - 12)*(1 - 3 * r);
259
 
    } else {
260
 
      //ret = 80/(7*M_PI*h*h)*( 6/(h*h))*(1 -r/h);
261
 
      ret = 80/(7*M_PI)*( 6)*(1 -r);
262
 
    }
263
 
  }
264
 
  return ret;
 
55
Real smoothkernelLucy(const double & r, const double & h) {
 
56
  if (r<=h) {
 
57
    // Lucy Kernel function, 
 
58
    return 105./(16.*M_PI*h*h*h) * (1 + 3 * r / h) * pow((1.0 - r / h), 3);
 
59
  } else {
 
60
    return 0;
 
61
  }
 
62
}
 
63
 
 
64
Real smoothkernelLucyGrad(const double & r, const double & h) {
 
65
  if (r<=h) {
 
66
    // 1st derivative of Lucy Kernel function, 
 
67
    return 105./(16.*M_PI*h*h*h) * (-12 * r) * (1/( h * h ) - 2 * r / (h * h * h ) + r * r / (h * h * h * h));
 
68
  } else {
 
69
    return 0;
 
70
  }
 
71
}
 
72
 
 
73
Real smoothkernelLucyLapl(const double & r, const double & h) {
 
74
  if (r<=h) {
 
75
    // 2nd derivative of Lucy Kernel function, 
 
76
    return 105./(16.*M_PI*h*h*h) * (-12) / (h * h * h * h) * (h * h - 2 * r * h + 3 * r * r);
 
77
  } else {
 
78
    return 0;
 
79
  }
 
80
}
 
81
 
 
82
KernelFunction returnKernelFunction(const int a, const typeKernFunctions typeF) {
 
83
  return returnKernelFunction(a, a, typeF);
265
84
}
266
85
 
267
86
KernelFunction returnKernelFunction(const int a, const int b, const typeKernFunctions typeF) {
268
87
  if (a != b) {
269
88
    throw runtime_error("Kernel types should be equal!");
270
89
  }
271
 
  if (a==Poly6) {
272
 
    if (typeF==Norm) {
273
 
      return smoothkernelPoly6;
274
 
    } else if (typeF==Grad) {
275
 
      return smoothkernelPoly6Grad;
276
 
    } else if (typeF==Lapl) {
277
 
      return smoothkernelPoly6Lapl;
278
 
    } else {
279
 
      KERNELFUNCDESCR
280
 
    }
281
 
  } else if (a==Spiky) {
282
 
    if (typeF==Norm) {
283
 
      return smoothkernelSpiky;
284
 
    } else if (typeF==Grad) {
285
 
      return smoothkernelSpikyGrad;
286
 
    } else if (typeF==Lapl) {
287
 
      return smoothkernelSpikyLapl;
288
 
    } else {
289
 
      KERNELFUNCDESCR
290
 
    }
291
 
  } else if (a==Visco) {
292
 
    if (typeF==Norm) {
293
 
      return smoothkernelVisco;
294
 
    } else if (typeF==Grad) {
295
 
      return smoothkernelViscoGrad;
296
 
    } else if (typeF==Lapl) {
297
 
      return smoothkernelViscoLapl;
298
 
    } else {
299
 
    }
300
 
  } else if (a==Lucy) {
 
90
  if (a==Lucy) {
301
91
    if (typeF==Norm) {
302
92
      return smoothkernelLucy;
303
93
    } else if (typeF==Grad) {
307
97
    } else {
308
98
      KERNELFUNCDESCR
309
99
    }
310
 
  } else if (a==Monaghan) {
311
 
    if (typeF==Norm) {
312
 
      return smoothkernelMonaghan;
313
 
    } else if (typeF==Grad) {
314
 
      return smoothkernelMonaghanGrad;
315
 
    } else if (typeF==Lapl) {
316
 
      return smoothkernelMonaghanLapl;
317
 
    } else {
318
 
      KERNELFUNCDESCR
319
 
    }
320
100
  } else {
321
101
    KERNELFUNCDESCR
322
102
  }
323
103
}
324
104
 
325
105
bool computeForceSPH(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force) {
326
 
  
327
106
  const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
328
107
  Scene* scene=Omega::instance().getScene().get();
329
108
  ViscElPhys& phys=*static_cast<ViscElPhys*>(_phys.get());
335
114
  
336
115
  //////////////////////////////////////////////////////////////////
337
116
  // Copy-paste
 
117
  // Handle periodicity.
 
118
  const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); 
 
119
  const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); 
338
120
  
339
121
  const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
340
122
  const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
341
 
  
342
 
    // Handle periodicity.
343
 
  const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); 
344
 
  const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); 
345
123
 
346
124
  const Vector3r c1x = (geom.contactPoint - de1.pos);
347
125
  const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
348
126
  
349
127
  const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
350
128
  const Real normalVelocity     = geom.normal.dot(relativeVelocity);
351
 
  //const Vector3r shearVelocity        = relativeVelocity-normalVelocity*geom.normal;
352
129
  
353
 
  // Copy-paste
 
130
  // Copy-paste  
354
131
  //////////////////////////////////////////////////////////////////
355
132
  
356
 
  if (phys.h<0) {
357
 
    Sphere* s1=dynamic_cast<Sphere*>(bodies[id1]->shape.get());
358
 
    Sphere* s2=dynamic_cast<Sphere*>(bodies[id2]->shape.get());
359
 
    if (s1 and s2) {
360
 
      phys.h = s1->radius + s2->radius;
361
 
    } else if (s1 and not(s2)) {
362
 
      phys.h = s1->radius;
363
 
    } else {
364
 
      phys.h = s2->radius;
365
 
    }
366
 
  }
367
 
  
368
 
  Real Mass = bodies[id2]->state->mass;
369
 
  if (Mass==0.0 and bodies[id1]->state->mass!= 0.0) {
370
 
    Mass = bodies[id1]->state->mass;
371
 
  }
372
 
  
373
 
  Real Rho = bodies[id2]->rho;
374
 
  if (Rho==0.0 and bodies[id1]->rho!=0.0) {
375
 
    Rho = bodies[id1]->rho;
376
 
  }
377
 
  
378
 
  const Vector3r xixj = (-geom.penetrationDepth + phys.h)*geom.normal;
 
133
  
 
134
  const Real Mass1 = bodies[id1]->state->mass;
 
135
  const Real Mass2 = bodies[id2]->state->mass;
 
136
  
 
137
  const Real Rho1 = bodies[id1]->rho;
 
138
  const Real Rho2 = bodies[id2]->rho;
 
139
  
 
140
  const Vector3r xixj = de2.pos - de1.pos;
379
141
  
380
142
  if (xixj.norm() < phys.h) {
381
143
    Real fpressure = 0.0;
382
 
    if (Rho!=0.0) {
383
 
      // [Mueller2003], (10)
384
 
      fpressure = -Mass * 
385
 
                  (bodies[id1]->press + bodies[id2]->press)/(2.0*Rho) *
386
 
                  phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
 
144
    if (Rho1!=0.0 and Rho2!=0.0) {
 
145
      // from [Monaghan1992], (3.3), multiply by Mass2, because we need a force, not du/dt
 
146
      fpressure = - Mass1 * Mass2 * (
 
147
                  bodies[id1]->press/(Rho1*Rho1) + 
 
148
                  bodies[id2]->press/(Rho2*Rho2) 
 
149
                  ) 
 
150
                  * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
387
151
    }
388
152
    
389
153
    Vector3r fvisc = Vector3r::Zero();
390
 
    if (Rho!=0.0) {
391
 
      fvisc = phys.mu * Mass * 
392
 
                  normalVelocity*geom.normal/Rho *
393
 
                  phys.kernelFunctionCurrentVisco(xixj.norm(), phys.h);
 
154
    if (Rho1!=0.0 and Rho2!=0.0) {
 
155
      // from [Morris1997], (22), multiply by Mass2, because we need a force, not du/dt
 
156
      fvisc = phys.mu * Mass1 * Mass2 *
 
157
      (-normalVelocity*geom.normal)/(Rho1*Rho1) *
 
158
      1 / (xixj.norm()) *
 
159
      phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
 
160
      //phys.kernelFunctionCurrentVisco(xixj.norm(), phys.h);
394
161
    }
395
162
    force = fpressure*geom.normal + fvisc;
396
163
    return true;