yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10700
[Branch ~yade-pkg/yade/git-trunk] Rev 3899: Add Cs parameters to particles, visualize them in VTK.
------------------------------------------------------------
revno: 3899
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-09 16:03:16 +0200
message:
Add Cs parameters to particles, visualize them in VTK.
modified:
core/Body.hpp
pkg/common/SPHEngine.cpp
pkg/common/SPHEngine.hpp
pkg/dem/VTKRecorder.cpp
--
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 'core/Body.hpp'
--- core/Body.hpp 2014-04-09 09:11:49 +0000
+++ core/Body.hpp 2014-04-09 14:03:16 +0000
@@ -88,9 +88,10 @@
((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)
+ ((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)
+ ((Real,Cs,0.0,, "Color field (only for SPH-model)")) // [Mueller2003], (15)
#endif
,
/* ctor */,
=== modified file 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp 2014-04-09 14:03:16 +0000
+++ pkg/common/SPHEngine.cpp 2014-04-09 14:03:16 +0000
@@ -8,11 +8,19 @@
#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();
+ {
+ 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();
+ }
}
void SPHEngine::calculateSPHRho(const shared_ptr<Body>& b) {
@@ -23,7 +31,8 @@
// Pointer to kernel function
KernelFunction kernelFunctionCurDensity = returnKernelFunction (KernFunctionDensity, KernFunctionDensity, Norm);
-
+
+ // Calculate rho 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());
@@ -49,6 +58,39 @@
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;
=== modified file 'pkg/common/SPHEngine.hpp'
--- pkg/common/SPHEngine.hpp 2014-04-09 14:03:16 +0000
+++ pkg/common/SPHEngine.hpp 2014-04-09 14:03:16 +0000
@@ -13,6 +13,7 @@
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."))
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2013-11-04 14:55:45 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-04-09 14:03:16 +0000
@@ -97,7 +97,13 @@
vtkSmartPointer<vtkFloatArray> spheresId = vtkSmartPointer<vtkFloatArray>::New();
spheresId->SetNumberOfComponents(1);
spheresId->SetName("id");
-
+
+#ifdef YADE_SPH
+ vtkSmartPointer<vtkFloatArray> spheresCs = vtkSmartPointer<vtkFloatArray>::New();
+ spheresCs->SetNumberOfComponents(1);
+ spheresCs->SetName("cs");
+#endif
+
vtkSmartPointer<vtkFloatArray> spheresMask = vtkSmartPointer<vtkFloatArray>::New();
spheresMask->SetNumberOfComponents(1);
spheresMask->SetName("mask");
@@ -414,7 +420,9 @@
damage->InsertNextValue(YADE_PTR_CAST<JCFpmState>(b->state)->tensBreak + YADE_PTR_CAST<JCFpmState>(b->state)->shearBreak);
damageRel->InsertNextValue(YADE_PTR_CAST<JCFpmState>(b->state)->tensBreakRel + YADE_PTR_CAST<JCFpmState>(b->state)->shearBreakRel);
}
-
+#ifdef YADE_SPH
+ spheresCs->InsertNextValue(b->Cs);
+#endif
if (recActive[REC_MATERIALID]) spheresMaterialId->InsertNextValue(b->material->id);
continue;
}
@@ -550,6 +558,9 @@
spheresUg->GetPointData()->AddArray(spheresLinVelLen);
spheresUg->GetPointData()->AddArray(spheresAngVelLen);
}
+#ifdef YADE_SPH
+ spheresUg->GetPointData()->AddArray(spheresCs);
+#endif
if (recActive[REC_STRESS]){
spheresUg->GetPointData()->AddArray(spheresNormalStressVec);
spheresUg->GetPointData()->AddArray(spheresShearStressVec);
