yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11060
[Branch ~yade-pkg/yade/git-trunk] Rev 4083: Export principal stresses and directions in VTKRecorder
------------------------------------------------------------
revno: 4083
committer: Jerome Duriez <jerome.duriez@xxxxxxxxxxxxxxx>
timestamp: Wed 2014-07-16 10:34:55 +0200
message:
Export principal stresses and directions in VTKRecorder
From bodyStressTensors(), see https://answers.launchpad.net/yade/+question/251712
modified:
pkg/dem/VTKRecorder.cpp
pkg/dem/VTKRecorder.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 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2014-07-11 17:39:21 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-07-16 08:34:55 +0000
@@ -79,7 +79,8 @@
else if(rec=="cracks") recActive[REC_CRACKS]=true;
else if(rec=="pericell" && scene->isPeriodic) recActive[REC_PERICELL]=true;
else if(rec=="liquidcontrol") recActive[REC_LIQ]=true;
- else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, jcfpm, cracks, pericell). Ignored.");
+ else if(rec=="bstresses") recActive[REC_BSTRESS]=true;
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, jcfpm, cracks, pericell, liquidcontrol, bstresses). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -102,6 +103,30 @@
radii->SetNumberOfComponents(1);
radii->SetName("radii");
+ vtkSmartPointer<vtkDoubleArray> spheresSigI = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresSigI->SetNumberOfComponents(1);
+ spheresSigI->SetName("sigI");
+
+ vtkSmartPointer<vtkDoubleArray> spheresSigII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresSigII->SetNumberOfComponents(1);
+ spheresSigII->SetName("sigII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresSigIII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresSigIII->SetNumberOfComponents(1);
+ spheresSigIII->SetName("sigIII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresDirI = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresDirI->SetNumberOfComponents(3);
+ spheresDirI->SetName("dirI");
+
+ vtkSmartPointer<vtkDoubleArray> spheresDirII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresDirII->SetNumberOfComponents(3);
+ spheresDirII->SetName("dirII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresDirIII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresDirIII->SetNumberOfComponents(3);
+ spheresDirIII->SetName("dirIII");
+
vtkSmartPointer<vtkDoubleArray> spheresMass = vtkSmartPointer<vtkDoubleArray>::New();
spheresMass->SetNumberOfComponents(1);
spheresMass->SetName("mass");
@@ -418,6 +443,13 @@
vector<Shop::bodyState> bodyStates;
if(recActive[REC_STRESS]) Shop::getStressForEachBody(bodyStates);
+
+ vector<Matrix3r> bStresses;
+ if (recActive[REC_BSTRESS])
+ {
+ Shop::getStressLWForEachBody(bStresses);
+ }
+
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
if (!b) continue;
if(mask!=0 && !b->maskCompatible(mask)) continue;
@@ -430,6 +462,73 @@
pid[0] = spheresPos->InsertNextPoint(pos[0], pos[1], pos[2]);
spheresCells->InsertNextCell(1,pid);
radii->InsertNextValue(sphere->radius);
+
+ if (recActive[REC_BSTRESS])
+ {
+ const Matrix3r& bStress = - bStresses[b->getId()]; // compressive states are negativ for getStressLWForEachBody; I want them as positiv
+ Eigen::SelfAdjointEigenSolver<Matrix3r> solver(bStress); // bStress is probably not symmetric (= self-adjoint for real matrices), but the solver hopefully works (considering only one half of bStress). And, moreover, existence of (real) eigenvalues is not sure for not symmetric bStress..
+ Matrix3r dirAll = solver.eigenvectors();
+ Vector3r valPropres = solver.eigenvalues(); // cf http://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html#a30caf3c3884a7f4a46b8ec94efd23c5e to be sure that valPropres[i] * dirAll.col(i) = bStress * dirAll.col(i)
+ int whereSigI(-1), whereSigII(-1), whereSigIII(-1); // all whereSig_i are in [0;2] : whereSigI = 2 => sigI=valPropres[2]
+ if ( valPropres[0] > std::max(valPropres[1],valPropres[2]) )
+ {
+ whereSigI = 0;
+ if (valPropres[1]>valPropres[2])
+ {
+ whereSigII=1;
+ whereSigIII=2;
+ }
+ else //valPropres[0] > valPropres[2] >= valPropres[1]
+ {
+ whereSigII = 2;
+ whereSigIII=1;
+ }
+ }
+ else // max(lambda1,lambda2) >= lambda0
+ {
+ if (valPropres[1]>=valPropres[2]) //max(lambda1,lambda2) = lambda1 : lambda 1 >= lambda2
+ {
+ whereSigI = 1;
+ if (valPropres[2]>=valPropres[0]) //lambda1 >= lambda2 >= lambda0
+ {
+ whereSigII=2;
+ whereSigIII=0;
+ }
+ else //lambda1 >= lambda0 > lambda2
+ {
+ whereSigII=0;
+ whereSigIII=2;
+ }
+ }
+ else //max(lambda1,lambda2) = lambda2 : lambda2 > lambda1
+ {
+ whereSigI = 2;
+ if (valPropres[1] > valPropres[0])
+ {
+ whereSigII = 1;
+ whereSigIII = 0;
+ }
+ else
+ {
+ whereSigIII = 1;
+ whereSigII = 0;
+ }
+ }
+ }
+
+ spheresSigI->InsertNextValue(valPropres[whereSigI]);
+ spheresSigII->InsertNextValue(valPropres[whereSigII]);
+ spheresSigIII->InsertNextValue(valPropres[whereSigIII]);
+ Real dirI[3] { (Real) dirAll(0,whereSigI), (Real) dirAll(1,whereSigI), (Real) dirAll(2,whereSigI) };
+ spheresDirI->InsertNextTupleValue(dirI);
+
+ Real dirII[3] { (Real) dirAll(0,whereSigII), (Real) dirAll(1,whereSigII), (Real) dirAll(2,whereSigII) };
+ spheresDirII->InsertNextTupleValue(dirII);
+
+ Real dirIII[3] { (Real) dirAll(0,whereSigIII), (Real) dirAll(1,whereSigIII), (Real) dirAll(2,whereSigIII) };
+ spheresDirIII->InsertNextTupleValue(dirIII);
+ }
+
if (recActive[REC_ID]) spheresId->InsertNextValue(b->getId());
if (recActive[REC_MASK]) spheresMask->InsertNextValue(GET_MASK(b));
if (recActive[REC_MASS]) spheresMass->InsertNextValue(b->state->mass);
@@ -647,7 +746,15 @@
if (recActive[REC_JCFPM]) {
spheresUg->GetPointData()->AddArray(damage);
}
-
+ if (recActive[REC_BSTRESS])
+ {
+ spheresUg->GetPointData()->AddArray(spheresSigI);
+ spheresUg->GetPointData()->AddArray(spheresSigII);
+ spheresUg->GetPointData()->AddArray(spheresSigIII);
+ spheresUg->GetPointData()->AddArray(spheresDirI);
+ spheresUg->GetPointData()->AddArray(spheresDirII);
+ spheresUg->GetPointData()->AddArray(spheresDirIII);
+ }
if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
#ifdef YADE_VTK_MULTIBLOCK
=== modified file 'pkg/dem/VTKRecorder.hpp'
--- pkg/dem/VTKRecorder.hpp 2014-06-12 15:11:52 +0000
+++ pkg/dem/VTKRecorder.hpp 2014-07-16 08:34:55 +0000
@@ -11,7 +11,7 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL,REC_LIQ};
+ enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL,REC_LIQ,REC_BSTRESS};
virtual void action();
void addWallVTK (vtkSmartPointer<vtkQuad>& boxes, vtkSmartPointer<vtkPoints>& boxesPos, Vector3r& W1, Vector3r& W2, Vector3r& W3, Vector3r& W4);
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.",
@@ -23,7 +23,7 @@
((bool,multiblock,false,,"Use multi-block (``.vtm``) files to store data, rather than separate ``.vtu`` files."))
#endif
((string,fileName,"",,"Base file name; it will be appended with {spheres,intrs,facets}-243100.vtu (unless *multiblock* is ``True``) depending on active recorders and step number (243100 in this case). It can contain slashes, but the directory must exist already."))
- ((vector<string>,recorders,vector<string>(1,string("all")),,"List of active recorders (as strings). ``all`` (the default value) enables all base and generic recorders.\n\n.. admonition:: Base recorders\n\n\tBase recorders save the geometry (unstructured grids) on which other data is defined. They are implicitly activated by many of the other recorders. Each of them creates a new file (or a block, if :yref:`multiblock <VTKRecorder.multiblock>` is set).\n\n\t``spheres``\n\t\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n\t``facets``\n\t\tSave :yref:`facets<Facet>` positions (vertices).\n\t``boxes``\n\t\tSave :yref:`boxes<Box>` positions (edges).\n\t``intr``\n\t\tStore interactions as lines between nodes at respective particles positions. Additionally stores magnitude of normal (``forceN``) and shear (``absForceT``) forces on interactions (the :yref:`geom<Interaction.geom> must be of type :yref:`NormShearPhys`). \n\n.. admonition:: Generic recorders\n\n\tGeneric recorders do not depend on specific model being used and save commonly useful data.\n\n\t``id``\n\t\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n\t``mass``\n\t\tSaves masses (field ``mass``) of spheres; active only if ``spheres`` is active.\n\t``clumpId``\n\t\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n\t``colors``\n\t\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n\t``mask``\n\t\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n\t``materialId``\n\t\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n\t``velocity``\n\t\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n\t``stress``\n\t\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n\t``pericell``\n\t\tSaves the shape of the cell (simulation has to be periodic).\n\n.. admonition:: Specific recorders\n\n\tThe following should only be activated in appropriate cases, otherwise crashes can occur due to violation of type presuppositions.\n\n\t``cpm``\n\t\tSaves data pertaining to the :yref:`concrete model<Law2_ScGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmStress`` (stress on particle); ``intr`` is activated automatically by ``cpm``\n\t``wpm``\n\t\tSaves data pertaining to the :yref:`wire particle model<Law2_ScGeom_WirePhys_WirePM>`: ``wpmForceNFactor`` shows the loading factor for the wire, e.g. normal force divided by threshold normal force.\n\t``jcfpm``\n\t\tSaves data pertaining to the :yref:`rock (smooth)-jointed model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``damage`` is defined by :yref:`JCFpmState.tensBreak` + :yref:`JCFpmState.shearBreak`; ``intr`` is activated automatically by ``jcfpm``, and :yref:`on joint<JCFpmPhys.isOnJoint>` or :yref:`cohesive<JCFpmPhys.isCohesive>` interactions can be vizualized.\n\t``cracks``\n\t\tSaves other data pertaining to the :yref:`rock model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``cracks`` shows locations where cohesive bonds failed during the simulation, with their types (0/1 for tensile/shear breakages), their sizes (0.5*(R1+R2)), and their normal directions. The :yref:`corresponding attribute<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` has to be activated, and Key attributes have to be consistent.\n\n"))
+ ((vector<string>,recorders,vector<string>(1,string("all")),,"List of active recorders (as strings). ``all`` (the default value) enables all base and generic recorders.\n\n.. admonition:: Base recorders\n\n\tBase recorders save the geometry (unstructured grids) on which other data is defined. They are implicitly activated by many of the other recorders. Each of them creates a new file (or a block, if :yref:`multiblock <VTKRecorder.multiblock>` is set).\n\n\t``spheres``\n\t\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n\t``facets``\n\t\tSave :yref:`facets<Facet>` positions (vertices).\n\t``boxes``\n\t\tSave :yref:`boxes<Box>` positions (edges).\n\t``intr``\n\t\tStore interactions as lines between nodes at respective particles positions. Additionally stores magnitude of normal (``forceN``) and shear (``absForceT``) forces on interactions (the :yref:`geom<Interaction.geom> must be of type :yref:`NormShearPhys`). \n\n.. admonition:: Generic recorders\n\n\tGeneric recorders do not depend on specific model being used and save commonly useful data.\n\n\t``id``\n\t\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n\t``mass``\n\t\tSaves masses (field ``mass``) of spheres; active only if ``spheres`` is active.\n\t``clumpId``\n\t\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n\t``colors``\n\t\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n\t``mask``\n\t\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n\t``materialId``\n\t\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n\t``velocity``\n\t\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n\t``stress``\n\t\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n\t``pericell``\n\t\tSaves the shape of the cell (simulation has to be periodic).\n\t``bstresses``\n\t\tSaves per-particle principal stresses (sigI >= sigII >= sigIII) and associated principal directions (dirI/II/III). Per-particle stress tensor is computed from :yref:`bodyStressTensors<yade.utils.bodyStressTensors>` considering positiv values for compressive states.\n\n.. admonition:: Specific recorders\n\n\tThe following should only be activated in appropriate cases, otherwise crashes can occur due to violation of type presuppositions.\n\n\t``cpm``\n\t\tSaves data pertaining to the :yref:`concrete model<Law2_ScGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmStress`` (stress on particle); ``intr`` is activated automatically by ``cpm``\n\t``wpm``\n\t\tSaves data pertaining to the :yref:`wire particle model<Law2_ScGeom_WirePhys_WirePM>`: ``wpmForceNFactor`` shows the loading factor for the wire, e.g. normal force divided by threshold normal force.\n\t``jcfpm``\n\t\tSaves data pertaining to the :yref:`rock (smooth)-jointed model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``damage`` is defined by :yref:`JCFpmState.tensBreak` + :yref:`JCFpmState.shearBreak`; ``intr`` is activated automatically by ``jcfpm``, and :yref:`on joint<JCFpmPhys.isOnJoint>` or :yref:`cohesive<JCFpmPhys.isCohesive>` interactions can be vizualized.\n\t``cracks``\n\t\tSaves other data pertaining to the :yref:`rock model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``cracks`` shows locations where cohesive bonds failed during the simulation, with their types (0/1 for tensile/shear breakages), their sizes (0.5*(R1+R2)), and their normal directions. The :yref:`corresponding attribute<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` has to be activated, and Key attributes have to be consistent.\n\n"))
((string,Key,"",,"Necessary if :yref:`recorders<VTKRecorder.recorders>` contains 'cracks'. A string specifying the name of file 'cracks___.txt' that is considered in this case (see :yref:`corresponding attribute<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.Key>`)."))
((int,mask,0,,"If mask defined, only bodies with corresponding groupMask will be exported. If 0, all bodies will be exported.")),
/*ctor*/
