yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10693
[Branch ~yade-pkg/yade/git-trunk] Rev 3893: Implement experimental SPH-model
------------------------------------------------------------
revno: 3893
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-09 11:11:49 +0200
message:
Implement experimental SPH-model
Basic implementation of smoothed particles hydrodynamics
Squashed commit of the following:
commit 2e49d31e1f10c6d18636dc0de9295e7107592a89
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Apr 4 10:58:11 2014 +0200
Add SPH-example (experimental).
commit 628ff6f9f2bdf130b0e69e5f791d76a8e881bf13
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Apr 4 10:14:41 2014 +0200
Add experimental comment on SPH.
commit 35cc070c4bc21917096f46aa92880e6c070924a3
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Apr 4 10:10:56 2014 +0200
Move force SPH calculations into the separate function.
commit 330ccb3f96fa0e9da89c1fb0861b917a356e823c
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Apr 3 20:05:54 2014 +0200
Fix forces calculation.
commit fad16b5beddf3390633aaf8cefa92f4b40062127
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Apr 3 15:32:02 2014 +0200
Re-enable SPH-functions.
commit 8829b10bcdda2ddeec7e2c869694639c4f33af42
Merge: 0cc1f70 c00c469
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Apr 3 13:55:11 2014 +0200
Merge branch 'master' into SPH
commit 0cc1f700c93ed389fef49e2388a03500ea1e495e
Merge: d06a92f fccbe5d
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Apr 3 10:51:00 2014 +0200
Merge with master.
commit d06a92fca289cbafb36076662c8d753d6e2a0c4d
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Mar 7 13:55:48 2014 +0100
Minor fix.
commit 7281ddf414486b72dd1320c7adc424467d10daed
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Mar 6 09:43:28 2014 +0100
Fix calculation.
commit 8a101259776a07aae32763feab3274f239ccf38d
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Mar 5 10:35:51 2014 +0100
Fix calculation of rho.
commit b215b6a2ed0172edba4266a5d303d3014339c4d2
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Tue Mar 4 17:16:21 2014 +0100
Fix smoothkernels .
commit f34b8bbe57570861dbb9c49f8d07bd04b1c3781f
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Tue Mar 4 15:57:58 2014 +0100
Fix compilation.
commit e71d0cb6e89b8d5a3cb96e3aeffcf27d6427f102
Merge: c470a00 5f72cd2
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Tue Mar 4 15:37:01 2014 +0100
Merge with newer trunk version.
commit c470a0033f4568b7cde9f2d96812199b0bbf4726
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Tue Mar 4 15:33:51 2014 +0100
Add some more smoothing kernels.
commit fd15132ae2091fcfa75422359c32aab6e6aec316
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Jan 17 15:35:41 2014 +0100
Update SPH-modell.
commit 96facbc435743bada75a3e35b28f82f7d734e3e2
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Jan 17 09:02:56 2014 +0100
Change velocity calculation.
commit 9efc074ae303e962da1ea757defba223a210f4e9
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Jan 16 17:33:47 2014 +0100
Change the press function in SPH.
commit 2ea260102cf98b05d5c43cabc4273122117e2a51
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 15 16:34:02 2014 +0100
Add vescous forces.
commit 23e7d58a61d5b2e872b52bc2788f9ff310ee83e8
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 15 15:50:38 2014 +0100
Prevent crash, if rho==0.0
commit 3806a200512169fa4f76a1675cc67ecdd2fec664
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 15 15:46:25 2014 +0100
Cosmetic change in fpress calculation.
commit d87f5c13f5ee2560da1a6e0cdcdb22894bf7d2d6
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 15 15:34:53 2014 +0100
Fix Fpress calculation in SPH-modell.
commit 91a082bb37a60d6531c1d1e20db9fe60c4ca1175
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 15 11:34:28 2014 +0100
Make smoothkernel function available for all classes.
commit cb54b34811e855a886d30e20b468f352c53fa404
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Fri Jan 10 17:09:58 2014 +0100
Add pressure force calculation into SPH-model.
commit a5220663df167a0d78a4ceefe8a8e296f1014fb2
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Jan 9 17:04:34 2014 +0100
Calculate press in SPH-model.
commit 5fd35c52942670704ccd51c98976b790c9006e85
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Jan 9 16:37:22 2014 +0100
Add basic files for SPH-model.
commit b8e808700c85f35afe12608411227a84b189e734
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Thu Jan 9 14:45:27 2014 +0100
Add a reference to SPH-article.
commit 48a39516e00f6e56d83b45bb1aedd7a834ccd732
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 8 13:27:06 2014 +0100
Add some additional variables for SPH into body.hpp
commit 68b8da73a16ab4aa0bda0bb1cceec51fd295e954
Author: Anton Gladky <gladky.anton@xxxxxxxxx>
Date: Wed Jan 8 11:24:28 2014 +0100
Add SPH-option to cmake.
added:
examples/sph/
examples/sph/sph_box.py
pkg/common/SPHEngine.cpp
pkg/common/SPHEngine.hpp
modified:
CMakeLists.txt
core/Body.hpp
doc/references.bib
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 'CMakeLists.txt'
--- CMakeLists.txt 2014-03-13 13:09:40 +0000
+++ CMakeLists.txt 2014-04-09 09:11:49 +0000
@@ -14,6 +14,7 @@
# ENABLE_GL2PS: enable GL2PS-option (ON by default)
# ENABLE_LINSOLV: enable LINSOLV-option (ON by default)
# ENABLE_PFVFLOW: enable PFVFLOW-option, FlowEngine (ON by default)
+# ENABLE_SPH: enable SPH-option, Smoothed Particle Hydrodynamics (OFF by default, experimental)
# runtimePREFIX: used for packaging, when install directory is not the same is runtime directory (/usr/local by default)
# CHUNKSIZE: set >1, if you want several sources to be compiled at once. Increases compilation speed and RAM-consumption during it (1 by default).
@@ -116,6 +117,7 @@
OPTION(ENABLE_GL2PS "Enable GL2PS" ${DEFAULT})
OPTION(ENABLE_LINSOLV "Enable direct solver for the flow engines (experimental)" ${DEFAULT})
OPTION(ENABLE_PFVFLOW "Enable flow engine (experimental)" ${DEFAULT})
+OPTION(ENABLE_SPH "Enable SPH" OFF)
#===========================================================
# Use Eigen3 by default
@@ -286,6 +288,13 @@
SET(DISABLED_FEATS "${DISABLED_FEATS} LinSolv")
ENDIF(ENABLE_LINSOLV)
#===============================================
+IF(ENABLE_SPH)
+ ADD_DEFINITIONS("-DYADE_SPH")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} SPH")
+ELSE(ENABLE_SPH)
+ SET(DISABLED_FEATS "${DISABLED_FEATS} SPH")
+ENDIF(ENABLE_SPH)
+#===============================================
IF(ENABLE_GL2PS)
FIND_PACKAGE(GL2PS)
IF(GL2PS_FOUND)
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2012-09-21 19:03:18 +0000
+++ core/Body.hpp 2014-04-09 09:11:49 +0000
@@ -87,6 +87,11 @@
((long,chain,-1,,"Id of chain to which the body belongs."))
((long,iterBorn,-1,,"Step number at which the body was added to simulation."))
((Real,timeBorn,-1,,"Time at which the body was added to simulation."))
+#ifdef YADE_SPH
+ ((Real,rho, -1.0,, "Current density (only for SPH-model)")) // [Mueller2003], (12)
+ ((Real,rho0,-1.0,, "Rest density (only for SPH-model)")) // [Mueller2003], (12)
+ ((Real,press,0.0,,"Pressure (only for SPH-model)")) // [Mueller2003], (12)
+#endif
,
/* ctor */,
/* py */
@@ -103,6 +108,11 @@
.add_property("timeBorn",&Body::timeBorn,"Returns time at which the body was added to simulation.")
.def_readwrite("chain",&Body::chain,"Returns Id of chain to which the body belongs.")
.def("intrs",&Body::py_intrs,"Return all interactions in which this body participates.")
+#ifdef YADE_SPH
+ .add_property("rho", &Body::rho, "Returns the current density (only for SPH-model).")
+ .add_property("rho0", &Body::rho0,"Returns the rest density (only for SPH-model).")
+ .add_property("press",&Body::press,"Returns the pressure (only for SPH-model).")
+#endif
);
};
REGISTER_SERIALIZABLE(Body);
=== modified file 'doc/references.bib'
--- doc/references.bib 2014-04-02 15:33:41 +0000
+++ doc/references.bib 2014-04-09 09:11:49 +0000
@@ -668,6 +668,22 @@
publisher={ACS Publications}
}
+@inproceedings{Mueller2003,
+ author = {M\"{u}ller, Matthias and Charypar, David and Gross, Markus},
+ title = {Particle-based Fluid Simulation for Interactive Applications},
+ booktitle = {Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation},
+ series = {SCA '03},
+ year = {2003},
+ isbn = {1-58113-659-5},
+ location = {San Diego, California},
+ pages = {154--159},
+ numpages = {6},
+ url = {http://dl.acm.org/citation.cfm?id=846276.846298},
+ acmid = {846298},
+ publisher = {Eurographics Association},
+ address = {Aire-la-Ville, Switzerland, Switzerland},
+}
+
@article {Soulie2006,
author = {Soulié, F. and Cherblanc, F. and El Youssoufi, M.S. and Saix, C.},
title = {Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials},
=== added directory 'examples/sph'
=== added file 'examples/sph/sph_box.py'
--- examples/sph/sph_box.py 1970-01-01 00:00:00 +0000
+++ examples/sph/sph_box.py 2014-04-09 09:11:49 +0000
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+# This script demonstrates SPH-engine in Yade
+# !!! Very experimental at the moment!!!
+
+from yade import utils, plot, qt
+o = Omega()
+
+# Physical parameters
+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
+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))
+
+# Add spheres
+d = 0.18
+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)))
+
+
+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))
+
+# Add engines
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(label='is2aabb'),Bo1_Facet_Aabb(label='is3aabb')]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(label='ss2sc'),Ig2_Facet_Sphere_ScGeom()],
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0.0,gravity=[0,0,-9.81]),
+ SPHEngine(mask=1, k=k, rho0 = rho),
+]
+
+
+O.step()
+qt.View()
=== added file 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp 1970-01-01 00:00:00 +0000
+++ pkg/common/SPHEngine.cpp 2014-04-09 09:11:49 +0000
@@ -0,0 +1,155 @@
+#ifdef YADE_SPH
+#include"SPHEngine.hpp"
+#include<yade/core/Scene.hpp>
+#include<yade/pkg/dem/ViscoelasticPM.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+
+#include<yade/core/State.hpp>
+#include<yade/core/Omega.hpp>
+
+void SPHEngine::action(){
+ 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();
+}
+
+void SPHEngine::calculateSPHRho(const shared_ptr<Body>& b) {
+ if (b->rho0<0) {
+ b->rho0 = rho0;
+ }
+ Real rho = 0;
+ 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;
+
+ rho += b2->state->mass*smoothkernelPoly6(SmoothDist, phys.h);
+ }
+ rho += b->state->mass*smoothkernelPoly6(0.0, s->radius);
+ }
+ b->rho = rho;
+}
+
+Real smoothkernelPoly6(const double & rrj, const double & h) {
+ if (rrj<=h) {
+ return 315/(64*M_PI*pow(h,9))*pow((h*h - rrj*rrj), 3); // [Mueller2003], (20)
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelSpiky(const double & rrj, const double & h) {
+ if (rrj<=h) {
+ return 15/(M_PI*pow(h,6))*pow((h - rrj), 3); // [Mueller2003], (21)
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelVisco(const double & rrj, const double & h) {
+ if (rrj<=h) {
+ return 15/(2*M_PI*pow(h,3))*(-rrj*rrj*rrj/(2*h*h*h) + rrj*rrj/(h*h) + h/(2*rrj) - 1); // [Mueller2003], (21)
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelViscoLapl(const double & rrj, const double & h) {
+ if (rrj<=h) {
+ return 45/(M_PI*pow(h,6))*(h - rrj); // [Mueller2003], (22+)
+ } else {
+ return 0;
+ }
+}
+
+void 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());
+
+ const int id1 = I->getId1();
+ const int id2 = I->getId2();
+
+ 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!=0.0;
+ }
+
+ const Vector3r xixj = (-geom.penetrationDepth + phys.h)*geom.normal;
+
+ if (xixj.norm() < phys.h) {
+ Real fpressure = 0.0;
+ if (Rho!=0.0) {
+ fpressure = Mass *
+ (bodies[id1]->press + bodies[id2]->press)/(2.0*Rho) *
+ smoothkernelSpiky(xixj.norm(), phys.h);
+ }
+
+ Vector3r fvisc = Vector3r::Zero();
+ if (Rho!=0.0) {
+ fvisc = phys.mu * Mass *
+ normalVelocity*geom.normal/Rho *
+ smoothkernelViscoLapl(xixj.norm(), phys.h);
+ }
+ force = fpressure*geom.normal - fvisc;
+ return;
+ } else {
+ return;
+ }
+}
+YADE_PLUGIN((SPHEngine));
+#endif
+
=== added file 'pkg/common/SPHEngine.hpp'
--- pkg/common/SPHEngine.hpp 1970-01-01 00:00:00 +0000
+++ pkg/common/SPHEngine.hpp 2014-04-09 09:11:49 +0000
@@ -0,0 +1,26 @@
+#ifdef YADE_SPH
+#pragma once
+
+#include<yade/core/PartialEngine.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+
+class SPHEngine: public PartialEngine{
+ public:
+ void calculateSPHRho(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)
+ );
+};
+REGISTER_SERIALIZABLE(SPHEngine);
+Real smoothkernelPoly6(const double & rrj, const double & h); // [Mueller2003] (20)
+Real smoothkernelSpiky(const double & rrj, const double & h); // [Mueller2003] (21)
+Real smoothkernelVisco(const double & rrj, const double & h); // [Mueller2003] (22)
+Real smoothkernelViscoLapl(const double & rrj, const double & h); // [Mueller2003] (22+)
+
+void 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-04-03 09:06:23 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2014-04-09 09:11:49 +0000
@@ -6,6 +6,10 @@
#include<yade/core/Scene.hpp>
#include<yade/pkg/common/Sphere.hpp>
+#ifdef YADE_SPH
+#include<yade/pkg/common/SPHEngine.hpp>
+#endif
+
YADE_PLUGIN((ViscElMat)(ViscElPhys)(Ip2_ViscElMat_ViscElMat_ViscElPhys)(Law2_ScGeom_ViscElPhys_Basic));
/* ViscElMat */
@@ -41,13 +45,22 @@
}
void computeForceTorqueViscEl(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2) {
+ ViscElPhys& phys=*static_cast<ViscElPhys*>(_phys.get());
const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
Scene* scene=Omega::instance().getScene().get();
- ViscElPhys& phys=*static_cast<ViscElPhys*>(_phys.get());
+
+#ifdef YADE_SPH
+//=======================================================================================================
+ if (phys.SPHmode) {
+ computeForceSPH(_geom, _phys, I, force);
+ return;
+ }
+//=======================================================================================================
+#endif
const int id1 = I->getId1();
const int id2 = I->getId2();
-
+
const BodyContainer& bodies = *scene->bodies;
const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
@@ -222,6 +235,12 @@
} else {
phys->mRtype = mRtype1;
}
+#ifdef YADE_SPH
+ if (mat1->SPHmode and mat2->SPHmode) {
+ phys->SPHmode=true;
+ phys->mu=(mat1->mu+mat2->mu)/2.0;
+ }
+#endif
}
/* Contact parameter calculation function */
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-04-03 09:27:51 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-04-09 09:11:49 +0000
@@ -11,6 +11,11 @@
#include<yade/pkg/dem/DemXDofGeom.hpp>
#include<yade/pkg/common/MatchMaker.hpp>
+#ifdef YADE_SPH
+#include<yade/pkg/common/SPHEngine.hpp>
+#endif
+
+
/* Simple viscoelastic model */
/// Material
@@ -27,6 +32,10 @@
((Real,ks,NaN,,"Shear elastic stiffness. Attention, this parameter cannot be set if tc, en or es is defined!"))
((Real,cs,NaN,,"Shear viscous constant. Attention, this parameter cannot be set if tc, en or es is defined!"))
((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
+#ifdef YADE_SPH
+ ((bool,SPHmode,false,,"True, if SPH-mode is enabled."))
+ ((Real,mu,-1,, "Viscosity. See Mueller [Mueller2003]_ .")) // [Mueller2003], (14)
+#endif
((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_.")),
createIndex();
);
@@ -43,6 +52,11 @@
((Real,cn,NaN,,"Normal viscous constant"))
((Real,cs,NaN,,"Shear viscous constant"))
((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
+#ifdef YADE_SPH
+ ((bool,SPHmode,false,,"True, if SPH-mode is enabled."))
+ ((Real,h,-1,, "Core radius. See Mueller [Mueller2003]_ .")) // [Mueller2003], (1)
+ ((Real,mu,-1,, "Viscosity. See Mueller [Mueller2003]_ .")) // [Mueller2003], (14)
+#endif
((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_")),
createIndex();
)
