yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05173
[Branch ~yade-dev/yade/trunk] Rev 2331: 1. Engine for PSD-analyze is finished, need testing.
------------------------------------------------------------
revno: 2331
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
branch nick: trunk
timestamp: Thu 2010-07-08 16:58:21 +0200
message:
1. Engine for PSD-analyze is finished, need testing.
modified:
pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.cpp
pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.hpp
pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp
pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp
pkg/dem/meta/RockPM.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.cpp'
--- pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.cpp 2010-07-07 15:42:59 +0000
+++ pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.cpp 2010-07-08 14:58:21 +0000
@@ -12,6 +12,7 @@
FOREACH(const shared_ptr<Body>& b, *scene->bodies){ // Annulling specimenNumber
if (!b) continue;
YADE_PTR_CAST<RpmState>(b->state)->specimenNumber = 0;
+ YADE_PTR_CAST<RpmState>(b->state)->specimenMass = 0;
}
//Check all interactions
@@ -21,56 +22,156 @@
body_id_t id1 = i->getId1(); //Get bodies ids from interaction
body_id_t id2 = i->getId2();
- int& specimenNumberId1 = YADE_PTR_CAST<RpmState>(Body::byId(id1)->state)->specimenNumber;
- int& specimenNumberId2 = YADE_PTR_CAST<RpmState>(Body::byId(id2)->state)->specimenNumber;
-
- bool cohesiveState = contPhys->isCohesive;
- if (cohesiveState==true){
- if ((specimenNumberId1 == 0) and (specimenNumberId2 == 0)){ //Both bodies are "untouched"
- specimenNumberId1 = curBin;
- specimenNumberId2 = curBin;
- curBin++;
- } else if ((specimenNumberId1 != 0) and (specimenNumberId2 == 0)){ //On of bodies is 0, another one has number specimen
- specimenNumberId2 = specimenNumberId1;
- } else if ((specimenNumberId1 == 0) and (specimenNumberId2 != 0)){ //On of bodies is 0, another one has number specimen
- specimenNumberId1 = specimenNumberId2;
- } else if ((specimenNumberId1 != 0) and (specimenNumberId2 != 0) and (specimenNumberId1 != specimenNumberId2) ){ //Bodies have different specimen number, but they are cohesive! Update it
- int minIdR = std::min(specimenNumberId1, specimenNumberId2);
- int maxIdR = std::max(specimenNumberId1, specimenNumberId2);
- specimenNumberId1 = minIdR;
- specimenNumberId2 = minIdR;
- identicalIds tempVar(minIdR, maxIdR);
- arrayIdentIds.push_back(tempVar);
- }
- }
-
-
+
+ const Sphere* sphere1 = dynamic_cast<Sphere*>(Body::byId(id1)->shape.get());
+ const Sphere* sphere2 = dynamic_cast<Sphere*>(Body::byId(id2)->shape.get());
+
+ if ((sphere1)&&(sphere2)) { //This class is ONLY for spheres
+ int& specimenNumberId1 = YADE_PTR_CAST<RpmState>(Body::byId(id1)->state)->specimenNumber;
+ int& specimenNumberId2 = YADE_PTR_CAST<RpmState>(Body::byId(id2)->state)->specimenNumber;
+
+ bool cohesiveState = contPhys->isCohesive;
+ if (cohesiveState==true){
+ if ((specimenNumberId1 == 0) and (specimenNumberId2 == 0)){ //Both bodies are "untouched"
+ specimenNumberId1 = curBin;
+ specimenNumberId2 = curBin;
+ curBin++;
+ } else if ((specimenNumberId1 != 0) and (specimenNumberId2 == 0)){ //On of bodies is 0, another one has number specimen
+ specimenNumberId2 = specimenNumberId1;
+ } else if ((specimenNumberId1 == 0) and (specimenNumberId2 != 0)){ //On of bodies is 0, another one has number specimen
+ specimenNumberId1 = specimenNumberId2;
+ } else if ((specimenNumberId1 != 0) and (specimenNumberId2 != 0) and (specimenNumberId1 != specimenNumberId2) ){ //Bodies have different specimen number, but they are cohesive! Update it
+ int minIdR = std::min(specimenNumberId1, specimenNumberId2);
+ int maxIdR = std::max(specimenNumberId1, specimenNumberId2);
+ specimenNumberId1 = minIdR;
+ specimenNumberId2 = minIdR;
+ identicalIds tempVar(minIdR, maxIdR, 0);
+ arrayIdentIds.push_back(tempVar);
+ }
+ }
+ if (contPhys->isCohesive==true) { //Check whether they are cohesive
+ numberCohesiveContacts++; //If yes - calculate them
+ }
+ }
+ }
+
+ //Clean dublicates
+ for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
+ for (unsigned int d=0; d<arrayIdentIds.size(); d++) {
+ if ((i!=d)&&(arrayIdentIds[i].id1 == arrayIdentIds[d].id1)&&(arrayIdentIds[i].id2 == arrayIdentIds[d].id2)) {
+ arrayIdentIds.erase(arrayIdentIds.begin()+d);
+ }
+ }
+ }
+
+ //Check for identical Ids and drop them off.
+ bool flagChange=true;
+ while (flagChange){
+ flagChange=false;
for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
for (unsigned int d=0; d<arrayIdentIds.size(); d++) {
if (i!=d){
identicalIds& tempAr1 = arrayIdentIds[i];
identicalIds& tempAr2 = arrayIdentIds[d];
- if (tempAr1.id2 == tempAr2.id1) {
+ if ((tempAr1.id2 == tempAr2.id1) and (tempAr2.id1>tempAr1.id1)) {
tempAr2.id1 = tempAr1.id1;
+ flagChange=true;
+ }
+ if (tempAr1.id2 == tempAr2.id2) {
+ if (tempAr2.id1>tempAr1.id1) {
+ tempAr2.id2 = tempAr2.id1;
+ tempAr2.id1 = tempAr1.id1;
+ flagChange=true;
+ } else if (tempAr2.id1<tempAr1.id1) {
+ tempAr1.id2 = tempAr1.id1;
+ tempAr1.id1 = tempAr2.id1;
+ flagChange=true;
+ } else {
+ arrayIdentIds.erase(arrayIdentIds.begin()+d);
+ }
}
}
}
}
-
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if (!b) continue;
- if ((YADE_PTR_CAST<RpmState>(b->state)->specimenNumber)==arrayIdentIds[i].id2) {
+ }
+
+ //Update RpmState, specimenIds
+ for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if (!b) continue;
+ const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
+ if (sphere) {
+ int tempSpecimenNum = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
+ if ((tempSpecimenNum==arrayIdentIds[i].id1) or (tempSpecimenNum==arrayIdentIds[i].id2)){
YADE_PTR_CAST<RpmState>(b->state)->specimenNumber=arrayIdentIds[i].id1;
}
}
}
-
- if (contPhys->isCohesive==true) { //Check whether they are cohesive
- numberCohesiveContacts++; //If yes - calculate them
- }
- }
+ }
+
+ int maximalSpecimenId = curBin;
+ arrayIdentIds.clear();
+ Real totalMass = 0;
+ //Calculate specimen masses, create vector for storing it
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if (!b) continue;
+ const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
+ if (sphere) {
+ Real massTemp = b->state->mass;
+ totalMass += massTemp;
+ int specimenNumberId = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
+
+ if (specimenNumberId != 0) {
+ bool foundItemInArray = false;
+ for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
+ if (arrayIdentIds[i].id1 == specimenNumberId) {
+ arrayIdentIds[i].mass+=massTemp;
+ foundItemInArray = true;
+ }
+ if (foundItemInArray) break;
+ }
+ if (!foundItemInArray) {
+ identicalIds tempVar(specimenNumberId, specimenNumberId+1, massTemp);
+ arrayIdentIds.push_back(tempVar);
+ }
+ } else {
+ YADE_PTR_CAST<RpmState>(b->state)->specimenNumber = maximalSpecimenId;
+ identicalIds tempVar(maximalSpecimenId, maximalSpecimenId+1, massTemp);
+ arrayIdentIds.push_back(tempVar);
+ maximalSpecimenId++;
+ }
+ }
+ }
+
+ //Update specimen masses
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if (!b) continue;
+ const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
+ if (sphere) {
+ int specimenNumberId = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
+ for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
+ if (arrayIdentIds[i].id1 == specimenNumberId) {
+ YADE_PTR_CAST<RpmState>(b->state)->specimenMass = arrayIdentIds[i].mass;
+ break;
+ }
+ }
+ }
+ }
+
+ //Check Masses in bodies loop and in dedicated vector. ONLY for debugging.
+ Real totalMass2 = 0;
+ for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
+ totalMass2+=arrayIdentIds[i].mass;
+ }
+ std::sort (arrayIdentIds.begin(), arrayIdentIds.end(), identicalIds::sortArrayIdentIds);
+
//Save data to a file
- out<<Omega::instance().getCurrentIteration()<<" "<<curBin<<"\n";
+ out<<"**********\n";
+ out<<"iter totalMass1 totalMass2 numSpecimen\n";
+ out<<Omega::instance().getCurrentIteration()<<" "<<totalMass<<" "<<totalMass2<<" "<<arrayIdentIds.size()<<"\n";
+ out<<"id mass\n";
+ for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
+ out<<arrayIdentIds[i].id1<<" "<<arrayIdentIds[i].mass<<"\n";
+ }
out.close();
}
=== modified file 'pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.hpp'
--- pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.hpp 2010-07-07 15:42:59 +0000
+++ pkg/dem/Engine/GlobalEngine/ParticleSizeDistrbutionRPMRecorder.hpp 2010-07-08 14:58:21 +0000
@@ -20,18 +20,13 @@
struct identicalIds{
int id1, id2;
- identicalIds (int id1r, int id2r){
+ Real mass;
+ identicalIds (int id1r, int id2r, Real massr){
assert(id1r<id2r);
id1 = id1r;
id2 = id2r;
- }
-
- static bool checkIdentical(identicalIds param1, identicalIds param2) {
- if ((param1.id1 == param2.id1) and (param1.id2 == param2.id2)) {
- return true;
- } else {
- return false;
- }
- }
+ mass = massr;
+ }
+ static bool sortArrayIdentIds (identicalIds i, identicalIds d) {return i.mass>d.mass;}
};
REGISTER_SERIALIZABLE(ParticleSizeDistrbutionRPMRecorder);
=== modified file 'pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp 2010-06-08 21:22:07 +0000
+++ pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp 2010-07-08 14:58:21 +0000
@@ -16,6 +16,7 @@
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-common/Facet.hpp>
#include<yade/pkg-dem/ConcretePM.hpp>
+#include<yade/pkg-dem/RockPM.hpp>
#include<yade/pkg-dem/Shop.hpp>
@@ -43,13 +44,14 @@
else if(rec=="facets") recActive[REC_FACETS]=true;
else if((rec=="colors") || (rec=="color"))recActive[REC_COLORS]=true;
else if(rec=="cpm") recActive[REC_CPM]=true;
+ else if(rec=="rpm") recActive[REC_RPM]=true;
else if(rec=="intr") recActive[REC_INTR]=true;
else if((rec=="ids") || (rec=="id")) recActive[REC_ID]=true;
else if(rec=="mask") recActive[REC_MASK]=true;
else if((rec=="clumpids") || (rec=="clumpId")) recActive[REC_CLUMPID]=true;
else if(rec=="materialId") recActive[REC_MATERIALID]=true;
else if(rec=="stress") recActive[REC_STRESS]=true;
- else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, color, stress, cpm, intr, id, clumpId, materialId). Ignored.");
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, color, stress, cpm, rpm, intr, id, clumpId, materialId). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -153,6 +155,15 @@
vtkSmartPointer<vtkFloatArray> cpmTau = vtkSmartPointer<vtkFloatArray>::New();
cpmTau->SetNumberOfComponents(3);
cpmTau->SetName("cpmTau");
+
+ // extras for RPM
+ vtkSmartPointer<vtkFloatArray> rpmSpecNum = vtkSmartPointer<vtkFloatArray>::New();
+ rpmSpecNum->SetNumberOfComponents(1);
+ rpmSpecNum->SetName("rpmSpecNum");
+ vtkSmartPointer<vtkFloatArray> rpmSpecMass = vtkSmartPointer<vtkFloatArray>::New();
+ rpmSpecMass->SetNumberOfComponents(1);
+ rpmSpecMass->SetName("rpmSpecMass");
+
if(recActive[REC_INTR]){
// save body positions, referenced by ids by vtkLine
@@ -238,6 +249,11 @@
cpmSigmaM->InsertNextValue((ss[0]+ss[1]+ss[2])/3.);
cpmTau->InsertNextTupleValue(t);
}
+ if (recActive[REC_RPM]){
+ rpmSpecNum->InsertNextValue(YADE_PTR_CAST<RpmState>(b->state)->specimenNumber);
+ rpmSpecMass->InsertNextValue(YADE_PTR_CAST<RpmState>(b->state)->specimenMass);
+ }
+
if (recActive[REC_MATERIALID]) spheresMaterialId->InsertNextValue(b->material->id);
continue;
}
@@ -303,6 +319,11 @@
spheresUg->GetPointData()->AddArray(cpmSigmaM);
spheresUg->GetPointData()->AddArray(cpmTau);
}
+ if (recActive[REC_RPM]){
+ spheresUg->GetPointData()->AddArray(rpmSpecNum);
+ spheresUg->GetPointData()->AddArray(rpmSpecMass);
+ }
+
if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
=== modified file 'pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp 2010-05-25 08:59:56 +0000
+++ pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp 2010-07-08 14:58:21 +0000
@@ -4,14 +4,14 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK};
+ enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM};
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(VTKRecorder,PeriodicEngine,"Engine recording snapshots of simulation into series of *.vtu files, readable by VTK-based postprocessing programs such as Paraview. Both bodies (spheres and facets) and interactions can be recorded, with various vector/scalar quantities that are defined on them.\n\n:yref:`PeriodicEngine.initRun` is initialized to ``True`` automatically.",
((bool,compress,false,"Compress output XML files [experimental]."))
((bool,skipFacetIntr,true,"Skip interactions with facets, when saving interactions"))
((bool,skipNondynamic,false,"Skip non-dynamic spheres (but not facets)."))
((string,fileName,"","Base file name; it will be appended with {spheres,intrs,facets}-243100.vtu depending on active recorders and step number (243100 in this case). It can contain slashes, but the directory must exist already."))
- ((vector<string>,recorders,,"List of active recorders (as strings). Accepted recorders are:\n\n``all``\n\tSaves all possible parameters, except of specific. Default value.\n``spheres``\n\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n``id``\n\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n``clumpId``\n\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n``colors``\n\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n``mask``\n\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n``materialId``\n\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n``velocity``\n\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n``facets``\n\tSave :yref:`facets<Facet>` positions (vertices).\n``stress``\n\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n``cpm``\n\tSaves data pertaining to the :yref:`concrete model<Law2_Dem3DofGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmSigma`` (stress on particle, normal components), ``cpmTau`` (shear components of stress on particle), ``cpmSigmaM`` (mean stress around particle); ``intr`` is activated automatically by ``cpm``.\n``intr``\n\tWhen ``cpm`` is used, it saves magnitude of normal (``forceN``) and shear (``absForceT``) forces.\n\n\tWithout ``cpm``, saves [TODO]\n\n"))
+ ((vector<string>,recorders,,"List of active recorders (as strings). Accepted recorders are:\n\n``all``\n\tSaves all possible parameters, except of specific. Default value.\n``spheres``\n\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n``id``\n\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n``clumpId``\n\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n``colors``\n\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n``mask``\n\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n``materialId``\n\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n``velocity``\n\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n``facets``\n\tSave :yref:`facets<Facet>` positions (vertices).\n``stress``\n\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n``cpm``\n\tSaves data pertaining to the :yref:`concrete model<Law2_Dem3DofGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmSigma`` (stress on particle, normal components), ``cpmTau`` (shear components of stress on particle), ``cpmSigmaM`` (mean stress around particle); ``intr`` is activated automatically by ``cpm``.\n``intr``\n\tWhen ``cpm`` is used, it saves magnitude of normal (``forceN``) and shear (``absForceT``) forces.\n\n\tWithout ``cpm``, saves [TODO]\n``rpm``\n\tSaves data pertaining to the :yref:`rock particle model<Law2_Dem3DofGeom_RockPMPhys_Rpm>`: ``rpmSpecNum`` shows different pieces of separated stones, only ids.\n ``rpmSpecMass`` shows masses of separated stones.\n\n"))
((int,mask,0,"If mask defined, only bodies with corresponding groupMask will be exported. If 0 - all bodies will be exported.")),
/*ctor*/
initRun=true;
=== modified file 'pkg/dem/meta/RockPM.hpp'
--- pkg/dem/meta/RockPM.hpp 2010-07-07 15:42:59 +0000
+++ pkg/dem/meta/RockPM.hpp 2010-07-08 14:58:21 +0000
@@ -19,6 +19,7 @@
class RpmState: public State {
YADE_CLASS_BASE_DOC_ATTRS(RpmState,State,"State information about Rpm body.",
((int,specimenNumber,0,"The variable is used for particle size distribution analyze. Indicates, to which part of specimen belongs para of particles."))
+ ((Real,specimenMass,0,"Indicates the mass of the whole stone, which owns the particle."))
);
};
REGISTER_SERIALIZABLE(RpmState);
Follow ups
