yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11548
[Branch ~yade-pkg/yade/git-trunk] Rev 3488: SPH-code refactoring.
------------------------------------------------------------
revno: 3488
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Fri 2014-10-17 19:14:29 +0200
message:
SPH-code refactoring.
added:
examples/sph/dam_break.py
modified:
doc/references.bib
examples/sph/sph_box.py
examples/sph/toystar.py
examples/sph/watercolumn.py
pkg/common/SPHEngine.cpp
pkg/common/SPHEngine.hpp
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'doc/references.bib'
--- doc/references.bib 2014-10-16 12:19:03 +0000
+++ doc/references.bib 2014-10-17 17:14:29 +0000
@@ -775,6 +775,7 @@
pages = {7341--7354},
number = {C4}
}
+
@book{Dallavalle1948,
title = {Micrometrics : The technology of fine particles},
publisher = {Pitman Pub. Corp},
@@ -815,3 +816,13 @@
doi="10.1007/s10035-008-0116-0",
url="http://www.springerlink.com/content/w0x307g110421035";
}
+
+@article{Monaghan1992,
+ author = {{Monaghan}, J.~J.},
+ title = "{Smoothed particle hydrodynamics}",
+ journal = {\araa},
+ year = 1992,
+ volume = 30,
+ pages = {543-574},
+ doi = {10.1146/annurev.aa.30.090192.002551}
+}
=== added file 'examples/sph/dam_break.py'
--- examples/sph/dam_break.py 1970-01-01 00:00:00 +0000
+++ examples/sph/dam_break.py 2014-10-17 17:14:29 +0000
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+from yade import utils, plot, qt
+o = Omega()
+
+# Physical parameters
+fr = 0.5;
+rho= 1000.0
+
+k = 1000.0
+mu = 10.0
+tc = 0.01; en = 0.7; et = 0.7;
+vel = 0.0
+
+#Rad = 0.006
+#h = 0.011
+
+Rad = 0.015
+h = 0.03
+
+#Rad = 0.02
+#h = 0.04
+
+o.dt = 0.0005
+
+X = 4.0
+Z = 3.0
+yCoeff = 4.0
+
+SpheresX = 1.0
+SpheresZ = 2.0
+
+# Add material
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True, h=h, mu=mu,tc=tc, en=en, et=et, KernFunctionPressure = 1, KernFunctionVisco = 1))
+idBox = pack.regularHexa(pack.inAlignedBox((0.0,-Rad*yCoeff,0.0),(X,Rad*yCoeff,Z)),radius=Rad,gap=0.0,color=(0,1,1), material=mat1, mask = 1, fixed=True)
+
+idBoxAdd = []
+idSpheresAdd = []
+
+for i in range(len(idBox)):
+ if (((idBox[i].state.pos[0])<yCoeff*Rad) or
+ ((idBox[i].state.pos[0])>(X-yCoeff*Rad - 2.0*Rad)) or
+ ((idBox[i].state.pos[2])<(yCoeff*Rad)) or
+ ((idBox[i].state.pos[2])>(Z - 2.0*Rad))):
+ idBoxAdd.append(idBox[i])
+ elif (((idBox[i].state.pos[0])<SpheresX) and
+ ((idBox[i].state.pos[2])<SpheresZ) and
+ ((idBox[i].state.pos[1])>-Rad) and
+ ((idBox[i].state.pos[1])<2.0*Rad)):
+ idBox[i].shape.color = Vector3(1,0,0)
+ idBox[i].state.fixed = False
+ idBox[i].state.blockedDOFs = 'y'
+ idSpheresAdd.append(idBox[i])
+
+
+
+O.bodies.append(idBoxAdd)
+idSpheres = O.bodies.append(idSpheresAdd)
+
+
+# Add engines
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(label='is2aabb')],ompThreads=1),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(label='ss2sc')],
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0.01,gravity=[0,0,-9.81]),
+ SPHEngine(mask=1, k=k, rho0 = rho, h=h, KernFunctionDensity= 1),
+ PyRunner(command='addPlotData()',iterPeriod=10,initRun=True,dead=False),
+]
+
+enlargeF = h/Rad*1.1
+print "enlargeF = %g"%enlargeF
+is2aabb.aabbEnlargeFactor = enlargeF
+ss2sc.interactionDetectionFactor = enlargeF
+
+
+# Function to add data to plot
+def addPlotData():
+ plot.addData(t = O.time, eKin = utils.kineticEnergy())
+
+plot.plots={'t':('eKin')};
+plot.plot()
+
+qt.View()
=== modified file 'examples/sph/sph_box.py'
--- examples/sph/sph_box.py 2014-04-09 09:11:49 +0000
+++ examples/sph/sph_box.py 2014-10-17 17:14:29 +0000
@@ -11,28 +11,29 @@
fr = 0.5;
rho=1000.0
-k = 1.0
-mu = 0.001
-tc = 0.0001; en = 0.7; et = 0.7;
-Rad = 10.0e-3
+k = 2000.0
+tc = 0.001; en = 0.5; et = 0.5;
+Rad = 7.0e-3
+h = 2*7.0e-3
o.dt = 0.00025
# Add material
-mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True,mu=mu,tc=tc, en=en, et=et))
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True, h=h, tc=tc, en=en, et=et, KernFunctionPressure = 1, KernFunctionVisco = 1))
+mat2 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=False,tc=tc, en=en, et=et))
# Add spheres
-d = 0.18
+d = 0.17
idSpheres = O.bodies.append(
pack.regularHexa(
pack.inAlignedBox(
- (-d,-d,-d), # lower angle
- (d,d,d)), # upper angle
- radius=Rad,gap=0.2*Rad, material=mat1,color=(0,1,1)))
+ (-d*0.5,-d*0.5,-d), # lower angle
+ (d*0.5,d*0.5,d)), # upper angle
+ radius=Rad,gap=0.01*Rad, mask=3, material=mat1,color=(0,1,1)))
id1 = O.bodies.append(geom.facetBox((Vector3(0.0,0,0.0)),
(Vector3(0.2, 0.2, 0.2)),
- material=mat1, mask=1, color=(1,0,0), wire=True))
+ material=mat2, mask=5, color=(1,0,0), wire=True))
# Add engines
o.engines = [
@@ -44,9 +45,13 @@
[Law2_ScGeom_ViscElPhys_Basic()],
),
NewtonIntegrator(damping=0.0,gravity=[0,0,-9.81]),
- SPHEngine(mask=1, k=k, rho0 = rho),
+ SPHEngine(mask=3, k=k, rho0 = rho, h=h, KernFunctionDensity= 1),
]
+enlargeF = h/Rad*1.1
+print "enlargeF = %g"%enlargeF
+is2aabb.aabbEnlargeFactor = enlargeF
+ss2sc.interactionDetectionFactor = enlargeF
O.step()
qt.View()
=== modified file 'examples/sph/toystar.py'
--- examples/sph/toystar.py 2014-04-10 06:21:57 +0000
+++ examples/sph/toystar.py 2014-10-17 17:14:29 +0000
@@ -12,13 +12,13 @@
rho=1000.0
k = 1.0
-mu = 0.01
tc = 0.0001; en = 0.7; et = 0.7;
Rad = 10.0e-3
+h = 2*Rad
o.dt = 0.00001
# Add material
-mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True,mu=mu,tc=tc, en=en, et=et))
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True,h=h,tc=tc, en=en, et=et, KernFunctionPressure = 1, KernFunctionVisco = 1))
# Add spheres
d = 0.8
@@ -42,10 +42,15 @@
),
CentralGravityEngine(accel=50.0, label='gr', centralBody=idCentralBody),
NewtonIntegrator(damping=0.1),
- SPHEngine(mask=1, k=k, rho0 = rho),
+ SPHEngine(mask=1, k=k, rho0 = rho, h=h, KernFunctionDensity= 1),
VTKRecorder(iterPeriod=1000, mask=1, fileName='./cpt/spheres-', recorders=['spheres','velocity','colors','intr','ids','mask','materialId','stress']),
]
+enlargeF = h/Rad*1.1
+print "enlargeF = %g"%enlargeF
+is2aabb.aabbEnlargeFactor = enlargeF
+ss2sc.interactionDetectionFactor = enlargeF
+
O.step()
qt.View()
=== modified file 'examples/sph/watercolumn.py'
--- examples/sph/watercolumn.py 2014-04-09 14:03:16 +0000
+++ examples/sph/watercolumn.py 2014-10-17 17:14:29 +0000
@@ -8,19 +8,21 @@
fr = 0.0;
rho=1000.0
-k = 5.0
-mu = 0.001
-tc = 0.0001; en = 0.7; et = 0.7;
+k = 5000.0
+tc = 0.001; en = 0.7; et = 0.7;
vel = 0.05
-Rad = 15.0e-3
+Rad = 12.0e-3
+h = 2*Rad
o.dt = 0.0002
scaleF = 0.001
# Add material
-mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True,mu=mu,tc=tc, en=en, et=et))
-id1 = O.bodies.append(ymport.gmsh("box.mesh", scale=scaleF, material=mat1, color=(1,0,0), mask = 1, wire=True))
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=True,h=h,tc=tc, en=en, et=et, KernFunctionPressure = 1, KernFunctionVisco = 1))
+mat2 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho, SPHmode=False,h=h,tc=tc, en=en, et=et))
+
+id1 = O.bodies.append(ymport.gmsh("box.mesh", scale=scaleF, material=mat2, color=(1,0,0), mask = 5, wire=True))
d = 15.0*scaleF
@@ -30,8 +32,7 @@
pack.inAlignedBox(
(0, -2000.0*scaleF, 0.0),
(1000*scaleF, 200*scaleF, 50.0*scaleF)),
- radius=d,gap=0.01*d, material=mat1, mask=1, color=(0,1,1)))
-
+ radius=d,gap=0.001*d, material=mat1, mask=3, color=(0,1,1)))
print len(idSpheres)
@@ -45,7 +46,7 @@
[Law2_ScGeom_ViscElPhys_Basic()],
),
NewtonIntegrator(damping=0.05,gravity=[0,-9.81,0]),
- SPHEngine(mask=1, k=k, rho0 = rho),
+ SPHEngine(mask=3, k=k, rho0 = rho, h = h, KernFunctionDensity= 1),
VTKRecorder(iterPeriod=100,fileName='./cpt/spheres-', recorders=['spheres','velocity','colors','intr','ids','mask','materialId','stress']),
VTKRecorder(iterPeriod=100,fileName='./cpt/facet-', recorders=['facets'],label='VTK_box2'),
PyRunner(command='addPlotData()',iterPeriod=50,dead=False),
@@ -56,6 +57,11 @@
plot.plots={'t':('Ekin')}; plot.plot()
+enlargeF = h/Rad*1.1
+print "enlargeF = %g"%enlargeF
+is2aabb.aabbEnlargeFactor = enlargeF
+ss2sc.interactionDetectionFactor = enlargeF
+
O.run(1, True)
qt.View()
=== modified file 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp 2014-10-15 06:44:01 +0000
+++ pkg/common/SPHEngine.cpp 2014-10-17 17:14:29 +0000
@@ -12,17 +12,9 @@
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
if(mask>0 && (b->groupMask & mask)==0) continue;
this->calculateSPHRho(b);
- b->press=k*(b->rho - b->rho0);
- } YADE_PARALLEL_FOREACH_BODY_END();
- }
- /*
- {
- YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
- if(mask>0 && (b->groupMask & mask)==0) continue;
- this->calculateSPHCs(b);
- } YADE_PARALLEL_FOREACH_BODY_END();
- }
- */
+ b->press=std::max(0.0, k*(b->rho - b->rho0));
+ } YADE_PARALLEL_FOREACH_BODY_END();
+ }
}
void SPHEngine::calculateSPHRho(const shared_ptr<Body>& b) {
@@ -32,7 +24,7 @@
Real rho = 0;
// Pointer to kernel function
- KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, KernFunctionDensity, Norm);
+ KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, Norm);
// Calculate rho for every particle
for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
@@ -49,255 +41,53 @@
Real Mass = b2->state->mass;
if (Mass == 0) Mass = b->state->mass;
- const Real SmoothDist = -geom.penetrationDepth + phys.h;
+ const Real SmoothDist = (b2->state->pos - b->state->pos).norm();
- // [Mueller2003], (3)
- rho += b2->state->mass*kernelFunctionCurDensity(SmoothDist, phys.h);
+ // [Monaghan1992], (2.7) (3.8)
+ rho += b2->state->mass*kernelFunctionCurDensity(SmoothDist, h);
}
- // [Mueller2003], (3), we need to consider the density of the current body (?)
- rho += b->state->mass*kernelFunctionCurDensity(0.0, s->radius);
}
+ // Self mass contribution
+ rho += b->state->mass*kernelFunctionCurDensity(0.0, h);
b->rho = rho;
}
-void SPHEngine::calculateSPHCs(const shared_ptr<Body>& b) {
- if (b->Cs<0) {
- b->Cs = 0.0;
- }
- Real Cs = 0;
-
- // Pointer to kernel function
- KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, KernFunctionDensity, Norm);
-
- // Calculate Cs for every particle
- for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
- const shared_ptr<Body> b2 = Body::byId((*it).first,scene);
- Sphere* s=dynamic_cast<Sphere*>(b->shape.get());
- if(!s) continue;
-
- if (((*it).second)->geom and ((*it).second)->phys) {
- const ScGeom geom = *(YADE_PTR_CAST<ScGeom>(((*it).second)->geom));
- const ViscElPhys phys=*(YADE_PTR_CAST<ViscElPhys>(((*it).second)->phys));
-
- if((b2->groupMask & mask)==0) continue;
-
- Real Mass = b2->state->mass;
- if (Mass == 0) Mass = b->state->mass;
-
- const Real SmoothDist = -geom.penetrationDepth + phys.h;
-
- // [Mueller2003], (15)
- Cs += b2->state->mass/b2->rho*kernelFunctionCurDensity(SmoothDist, phys.h);
- }
- }
- b->Cs = Cs;
-}
-
-Real smoothkernelPoly6(const double & rr, const double & hh) {
- if (rr<=hh) {
- const Real h = 1; const Real r = rr/hh;
- //return 315/(64*M_PI*pow(h,9))*pow((h*h - r*r), 3);
- return 315/(64*M_PI)*pow((h - r*r), 3);
- } else {
- return 0;
- }
-}
-
-Real smoothkernelPoly6Grad(const double & rr, const double & hh) {
- if (rr<=hh) {
- const Real h = 1; const Real r = rr/hh;
- //return -315/(64*M_PI*pow(h,9))*(-6*r*pow((h*h-r*r), 2));
- return -315/(64*M_PI)*(-6*r*pow((h-r*r), 2));
- } else {
- return 0;
- }
-}
-
-Real smoothkernelPoly6Lapl(const double & rr, const double & hh) {
- if (rr<=hh) {
- const Real h = 1; const Real r = rr/hh;
- //return 315/(64*M_PI*pow(h,9))*(-6*(h*h-r*r)*(h*h - 5*r*r));
- return 315/(64*M_PI)*(-6*(h*h-r*r)*(h*h - 5*r*r));
- } else {
- return 0;
- }
-}
-
-Real smoothkernelSpiky(const double & rr, const double & hh) {
- if (rr<=hh) {
- const Real h = 1; const Real r = rr/hh;
- //return 15/(M_PI*pow(h,6))*(pow((h-r), 3)); // [Mueller2003], (21)
- return 15/(M_PI)*(pow((h-r), 3)); // [Mueller2003], (21)
- } else {
- return 0;
- }
-}
-
-Real smoothkernelSpikyGrad(const double & rr, const double & hh) {
- if (rr<=hh) {
- const Real h = 1; const Real r = rr/hh;
- //return -15/(M_PI*pow(h,6))*(-3*pow((h-r),2));
- return -15/(M_PI)*(-3*pow((h-r),2));
- } else {
- return 0;
- }
-}
-
-Real smoothkernelSpikyLapl(const double & rr, const double & hh) {
- if (rr<=hh) {
- const Real h = 1; const Real r = rr/hh;
- //return 15/(M_PI*pow(h,6))*(6*(h-r));
- return 15/(M_PI)*(6*(h-r));
- } else {
- return 0;
- }
-}
-
-Real smoothkernelVisco(const double & rr, const double & hh) {
- if (rr<=hh and rr!=0 and hh!=0) {
- const Real h = 1; const Real r = rr/hh;
- //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)
- return 15/(2*M_PI)*(-r*r*r/(2) + r*r + h/(2*r) -1); // [Mueller2003], (21)
- } else {
- return 0;
- }
-}
-
-Real smoothkernelViscoGrad(const double & rr, const double & hh) {
- if (rr<=hh and rr!=0 and hh!=0) {
- const Real h = 1; const Real r = rr/hh;
- //return -15/(2*M_PI*pow(h,3))*(-3*r*r/(2*h*h*h) + 2*r/(h*h) - h/(2*r*r));
- return -15/(2*M_PI)*(-3*r*r/(2) + 2*r - h/(2*r*r));
- } else {
- return 0;
- }
-}
-
-Real smoothkernelViscoLapl(const double & rr, const double & hh) {
- if (rr<=hh and rr!=0 and hh!=0) {
- const Real h = 1; const Real r = rr/hh;
- //return 45/(M_PI*pow(h,6))*(h - rrj); // [Mueller2003], (22+)
- //return 15/(2*M_PI*pow(h,3))*(-3*r*r/(2*h*h*h) + 2*r/(h*h) - h/(2*r*r));
- return 15/(2*M_PI)*(-3*r*r/(2) + 2*r - h/(2*r*r));
- } else {
- return 0;
- }
-}
-
-Real smoothkernelLucy(const double & rr, const double & hh) {
- if (rr<=hh and rr!=0 and hh!=0) {
- //const Real h = 1;
- //return 5/(9*M_PI*pow(h,2))*(1+3*r/h)*pow((1-r/h),3);
- const Real r = rr/hh;
- return 5/(9*M_PI)*(1+3*r)*pow((1-r),3);
- } else {
- return 0;
- }
-}
-
-Real smoothkernelLucyGrad(const double & rr, const double & hh) {
- if (rr<=hh and rr!=0 and hh!=0) {
- //const Real h = 1;
- //return -5/(9*M_PI*pow(h,2))*(-12*r/(h*h))*pow((1-r/h),2);
- const Real r = rr/hh;
- return -5/(9*M_PI)*(-12*r)*pow((1-r),2);
- } else {
- return 0;
- }
-}
-
-Real smoothkernelLucyLapl(const double & rr, const double & hh) {
- if (rr<=hh and rr!=0 and hh!=0) {
- //const Real h = 1;
- //return 5/(9*M_PI*pow(h,2))*(-12/(h*h))*(1-r/h)*(1-3*r/h);
- const Real r = rr/hh;
- return 5/(9*M_PI)*(-12)*(1-r)*(1-3*r);
- } else {
- return 0;
- }
-}
-
-Real smoothkernelMonaghan(const double & rr, const double & hh) {
- Real ret = 0.0;
- if (hh!=0) {
- //const Real h = 1;
- const Real r = rr/hh;
- if (rr/hh<0.5) {
- //ret = 40/(7*M_PI*h*h)*(1 - 6*pow((r/h),2) + 6*pow((r/h),3));
- ret = 40/(7*M_PI)*(1 - 6*pow((r),2) + 6*pow((r),3));
- } else {
- //ret = 80/(7*M_PI*h*h)*pow((1 - (r/h)), 3);
- ret = 80/(7*M_PI)*pow((1 - (r)), 3);
- }
- }
- return ret;
-}
-
-Real smoothkernelMonaghanGrad(const double & rr, const double & hh) {
- Real ret = 0.0;
- if (hh!=0) {
- const Real h = 1; const Real r = rr/hh;
- if (rr/hh<0.5) {
- //ret = -40/(7*M_PI*h*h)*( - 6*r/(h*h))*(2 - 3 * r/(h*h*h));
- ret = -40/(7*M_PI)*( - 6*r)*(2 - 3 * r);
- } else {
- //ret = -80/(7*M_PI*h*h)*( -3/h)*pow((1 -r/h), 2);
- ret = -80/(7*M_PI)*( -3/h)*pow((1 -r), 2);
- }
- }
- return ret;
-}
-
-Real smoothkernelMonaghanLapl(const double & rr, const double & hh) {
- Real ret = 0.0;
- if (hh!=0) {
- //const Real h = 1;
- const Real r = rr/hh;
- if (rr/hh<0.5) {
- //ret = 40/(7*M_PI*h*h)*( - 12/(h*h))*(1 - 3 * r/(h*h*h*h*h));
- ret = 40/(7*M_PI)*( - 12)*(1 - 3 * r);
- } else {
- //ret = 80/(7*M_PI*h*h)*( 6/(h*h))*(1 -r/h);
- ret = 80/(7*M_PI)*( 6)*(1 -r);
- }
- }
- return ret;
+Real smoothkernelLucy(const double & r, const double & h) {
+ if (r<=h) {
+ // Lucy Kernel function,
+ return 105./(16.*M_PI*h*h*h) * (1 + 3 * r / h) * pow((1.0 - r / h), 3);
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelLucyGrad(const double & r, const double & h) {
+ if (r<=h) {
+ // 1st derivative of Lucy Kernel function,
+ return 105./(16.*M_PI*h*h*h) * (-12 * r) * (1/( h * h ) - 2 * r / (h * h * h ) + r * r / (h * h * h * h));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelLucyLapl(const double & r, const double & h) {
+ if (r<=h) {
+ // 2nd derivative of Lucy Kernel function,
+ return 105./(16.*M_PI*h*h*h) * (-12) / (h * h * h * h) * (h * h - 2 * r * h + 3 * r * r);
+ } else {
+ return 0;
+ }
+}
+
+KernelFunction returnKernelFunction(const int a, const typeKernFunctions typeF) {
+ return returnKernelFunction(a, a, typeF);
}
KernelFunction returnKernelFunction(const int a, const int b, const typeKernFunctions typeF) {
if (a != b) {
throw runtime_error("Kernel types should be equal!");
}
- if (a==Poly6) {
- if (typeF==Norm) {
- return smoothkernelPoly6;
- } else if (typeF==Grad) {
- return smoothkernelPoly6Grad;
- } else if (typeF==Lapl) {
- return smoothkernelPoly6Lapl;
- } else {
- KERNELFUNCDESCR
- }
- } else if (a==Spiky) {
- if (typeF==Norm) {
- return smoothkernelSpiky;
- } else if (typeF==Grad) {
- return smoothkernelSpikyGrad;
- } else if (typeF==Lapl) {
- return smoothkernelSpikyLapl;
- } else {
- KERNELFUNCDESCR
- }
- } else if (a==Visco) {
- if (typeF==Norm) {
- return smoothkernelVisco;
- } else if (typeF==Grad) {
- return smoothkernelViscoGrad;
- } else if (typeF==Lapl) {
- return smoothkernelViscoLapl;
- } else {
- }
- } else if (a==Lucy) {
+ if (a==Lucy) {
if (typeF==Norm) {
return smoothkernelLucy;
} else if (typeF==Grad) {
@@ -307,23 +97,12 @@
} else {
KERNELFUNCDESCR
}
- } else if (a==Monaghan) {
- if (typeF==Norm) {
- return smoothkernelMonaghan;
- } else if (typeF==Grad) {
- return smoothkernelMonaghanGrad;
- } else if (typeF==Lapl) {
- return smoothkernelMonaghanLapl;
- } else {
- KERNELFUNCDESCR
- }
} else {
KERNELFUNCDESCR
}
}
bool computeForceSPH(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force) {
-
const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
Scene* scene=Omega::instance().getScene().get();
ViscElPhys& phys=*static_cast<ViscElPhys*>(_phys.get());
@@ -333,66 +112,29 @@
const BodyContainer& bodies = *scene->bodies;
- //////////////////////////////////////////////////////////////////
- // Copy-paste
-
const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
- // Handle periodicity.
- const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero();
- const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero();
-
- const Vector3r c1x = (geom.contactPoint - de1.pos);
- const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
-
- const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
- const Real normalVelocity = geom.normal.dot(relativeVelocity);
- //const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom.normal;
-
- // Copy-paste
- //////////////////////////////////////////////////////////////////
-
- if (phys.h<0) {
- Sphere* s1=dynamic_cast<Sphere*>(bodies[id1]->shape.get());
- Sphere* s2=dynamic_cast<Sphere*>(bodies[id2]->shape.get());
- if (s1 and s2) {
- phys.h = s1->radius + s2->radius;
- } else if (s1 and not(s2)) {
- phys.h = s1->radius;
- } else {
- phys.h = s2->radius;
- }
- }
-
- Real Mass = bodies[id2]->state->mass;
- if (Mass==0.0 and bodies[id1]->state->mass!= 0.0) {
- Mass = bodies[id1]->state->mass;
- }
-
- Real Rho = bodies[id2]->rho;
- if (Rho==0.0 and bodies[id1]->rho!=0.0) {
- Rho = bodies[id1]->rho;
- }
-
- const Vector3r xixj = (-geom.penetrationDepth + phys.h)*geom.normal;
+ const Real Mass1 = bodies[id1]->state->mass;
+ const Real Mass2 = bodies[id2]->state->mass;
+
+ const Real Rho1 = bodies[id1]->rho;
+ const Real Rho2 = bodies[id2]->rho;
+
+ const Vector3r xixj = de2.pos - de1.pos;
if (xixj.norm() < phys.h) {
Real fpressure = 0.0;
- if (Rho!=0.0) {
- // [Mueller2003], (10)
- fpressure = -Mass *
- (bodies[id1]->press + bodies[id2]->press)/(2.0*Rho) *
- phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
+ if (bodies[id1]->rho!=0.0 and bodies[id2]->rho!=0.0) {
+ // from [Monaghan1992], (3.3), multiply by Mass2, because we need a force, not du/dt
+ fpressure = - Mass1 * Mass2 * (
+ bodies[id1]->press/(Rho1*Rho1) +
+ bodies[id2]->press/(Rho2*Rho2)
+ )
+ * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
}
- Vector3r fvisc = Vector3r::Zero();
- if (Rho!=0.0) {
- fvisc = phys.mu * Mass *
- normalVelocity*geom.normal/Rho *
- phys.kernelFunctionCurrentVisco(xixj.norm(), phys.h);
- }
- force = fpressure*geom.normal + fvisc;
+ force = fpressure*geom.normal;
return true;
} else {
return false;
=== modified file 'pkg/common/SPHEngine.hpp'
--- pkg/common/SPHEngine.hpp 2014-10-15 06:44:01 +0000
+++ pkg/common/SPHEngine.hpp 2014-10-17 17:14:29 +0000
@@ -6,40 +6,29 @@
typedef Real (* KernelFunction)(const double & r, const double & h);
-enum KernFunctions {Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5};
+enum KernFunctions {Lucy=1};
#define KERNELFUNCDESCR throw runtime_error("Type of kernel function undefined! The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5.");
enum typeKernFunctions {Norm, Grad, Lapl};
class SPHEngine: public PartialEngine{
public:
void calculateSPHRho(const shared_ptr<Body>& b);
- void calculateSPHCs (const shared_ptr<Body>& b);
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS(SPHEngine,PartialEngine,"Apply given torque (momentum) value at every subscribed particle, at every step. ",
((int, mask,-1,, "Bitmask for SPH-particles."))
((Real,k,-1,, "Gas constant for SPH-interactions (only for SPH-model). See Mueller [Mueller2003]_ .")) // [Mueller2003], (11)
((Real,rho0,-1,, "Rest density. See Mueller [Mueller2003]_ .")) // [Mueller2003], (1)
- ((int,KernFunctionDensity, Poly6,, "Kernel function for density calculation (by default - Poly6). The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."))
+ ((Real,h,-1,, "Core radius. See Mueller [Mueller2003]_ .")) // [Mueller2003], (1)
+ ((int,KernFunctionDensity, Lucy,, "Kernel function for density calculation (by default - Lucy). The following kernel functions are available: Lucy=1."))
);
};
REGISTER_SERIALIZABLE(SPHEngine);
-Real smoothkernelPoly6(const double & r, const double & h); // [Mueller2003] (20)
-Real smoothkernelPoly6Grad(const double & r, const double & h);
-Real smoothkernelPoly6Lapl(const double & r, const double & h);
-Real smoothkernelSpiky(const double & r, const double & h); // [Mueller2003] (21)
-Real smoothkernelSpikyGrad(const double & r, const double & h);
-Real smoothkernelSpikyLapl(const double & r, const double & h);
-Real smoothkernelVisco(const double & r, const double & h); // [Mueller2003] (22)
-Real smoothkernelViscoGrad(const double & r, const double & h); // [Mueller2003] (22)
-Real smoothkernelViscoLapl(const double & r, const double & h);
Real smoothkernelLucy(const double & r, const double & h);
Real smoothkernelLucyGrad(const double & r, const double & h);
Real smoothkernelLucyLapl(const double & r, const double & h);
-Real smoothkernelMonaghan(const double & r, const double & h);
-Real smoothkernelMonaghanGrad(const double & r, const double & h);
-Real smoothkernelMonaghanLapl(const double & r, const double & h);
KernelFunction returnKernelFunction(const int a, const int b, const typeKernFunctions typeF);
+KernelFunction returnKernelFunction(const int a, const typeKernFunctions typeF);
bool computeForceSPH(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force);
#endif
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2014-10-17 17:14:29 +0000
@@ -256,6 +256,7 @@
if (mat1->SPHmode and mat2->SPHmode) {
phys->SPHmode=true;
phys->mu=(mat1->mu+mat2->mu)/2.0;
+ phys->h=(mat1->h+mat2->h)/2.0;
}
phys->kernelFunctionCurrentPressure = returnKernelFunction (mat1->KernFunctionPressure, mat2->KernFunctionPressure, Grad);
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-10-17 17:14:29 +0000
@@ -35,8 +35,9 @@
#ifdef YADE_SPH
((bool,SPHmode,false,,"True, if SPH-mode is enabled."))
((Real,mu,-1,, "Viscosity. See Mueller [Mueller2003]_ .")) // [Mueller2003], (14)
- ((int,KernFunctionPressure,Spiky,, "Kernel function for pressure calculation (by default - Spiky). The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."))
- ((int,KernFunctionVisco, Visco,, "Kernel function for viscosity calculation (by default - Visco). The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."))
+ ((Real,h,-1,, "Core radius. See Mueller [Mueller2003]_ .")) // [Mueller2003], (1)
+ ((int,KernFunctionPressure,Lucy,, "Kernel function for pressure calculation (by default - Lucy). The following kernel functions are available: Lucy=1."))
+ ((int,KernFunctionVisco, Lucy,, "Kernel function for viscosity calculation (by default - Lucy). The following kernel functions are available: Lucy=1."))
#endif
((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_.")),
createIndex();
