yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10081
[Branch ~yade-pkg/yade/git-trunk] Rev 3698: Added periodic cell VTK export (to VTKRecorder and export.VTKExporter)
------------------------------------------------------------
revno: 3698
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Fri 2013-10-04 17:25:30 +0200
message:
Added periodic cell VTK export (to VTKRecorder and export.VTKExporter)
added:
examples/test/vtkPeriodicCell.py
modified:
pkg/dem/VTKRecorder.cpp
pkg/dem/VTKRecorder.hpp
py/export.py
--
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
=== added file 'examples/test/vtkPeriodicCell.py'
--- examples/test/vtkPeriodicCell.py 1970-01-01 00:00:00 +0000
+++ examples/test/vtkPeriodicCell.py 2013-10-04 15:25:30 +0000
@@ -0,0 +1,23 @@
+######################################################################
+# Simple script to test VTK export of periodic cell
+######################################################################
+# enable periodic cell
+O.periodic=True
+# insert some bodies
+sp = randomPeriPack(radius=1,initSize=(10,20,30),memoizeDb='/tmp/vtkPeriodicCell.sqlite')
+sp.toSimulation()
+# transform the cell a bit
+O.cell.hSize *= Matrix3(1,.1,.1, .1,1,0, .1,0,1) # skew the cell in xy and xz plane
+O.cell.hSize *= Matrix3(1,0,0, 0,.8,.6, 0,-.6,.8) # rotate it along x axis
+
+O.step()
+
+# test of export.VTKExporter
+from yade import export
+vtk1 = export.VTKExporter('/tmp/vtkPeriodicCell-VTKExporter')
+vtk1.exportSpheres()
+vtk1.exportPeriodicCell()
+
+# test of VTKReorder
+vtk2 = VTKRecorder(fileName='/tmp/vtkPeriodicCell-VTKRecorder-',recorders=['spheres','pericell'])
+vtk2() # do the export
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2013-07-23 10:53:20 +0000
+++ pkg/dem/VTKRecorder.cpp 2013-10-04 15:25:30 +0000
@@ -16,6 +16,7 @@
#include<vtkTriangle.h>
#include<vtkLine.h>
#include<vtkQuad.h>
+#include<vtkHexahedron.h>
#ifdef YADE_VTK_MULTIBLOCK
#include<vtkXMLMultiBlockDataWriter.h>
#include<vtkMultiBlockDataSet.h>
@@ -49,6 +50,7 @@
recActive[REC_CLUMPID]=true;
recActive[REC_MATERIALID]=true;
recActive[REC_STRESS]=true;
+ if (scene->isPeriodic) { recActive[REC_PERICELL]=true; }
}
else if(rec=="spheres") recActive[REC_SPHERES]=true;
else if(rec=="velocity") recActive[REC_VELOCITY]=true;
@@ -64,7 +66,8 @@
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, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId). Ignored.");
+ else if(rec=="pericell" && scene->isPeriodic) recActive[REC_PERICELL]=true;
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, pericell). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -189,6 +192,10 @@
intrAbsForceT->SetNumberOfComponents(3);
intrAbsForceT->SetName("absForceT");
+ // pericell
+ vtkSmartPointer<vtkPoints> pericellPoints = vtkSmartPointer<vtkPoints>::New();
+ vtkSmartPointer<vtkCellArray> pericellHexa = vtkSmartPointer<vtkCellArray>::New();
+
// extras for CPM
if(recActive[REC_CPM]){ CpmStateUpdater csu; csu.update(scene); }
vtkSmartPointer<vtkFloatArray> cpmDamage = vtkSmartPointer<vtkFloatArray>::New();
@@ -433,6 +440,32 @@
}
}
+ if (recActive[REC_PERICELL]) {
+ const Matrix3r& hSize = scene->cell->hSize;
+ Vector3r v0 = hSize*Vector3r(0,0,1);
+ Vector3r v1 = hSize*Vector3r(0,1,1);
+ Vector3r v2 = hSize*Vector3r(1,1,1);
+ Vector3r v3 = hSize*Vector3r(1,0,1);
+ Vector3r v4 = hSize*Vector3r(0,0,0);
+ Vector3r v5 = hSize*Vector3r(0,1,0);
+ Vector3r v6 = hSize*Vector3r(1,1,0);
+ Vector3r v7 = hSize*Vector3r(1,0,0);
+ pericellPoints->InsertNextPoint(v0[0],v0[1],v0[2]);
+ pericellPoints->InsertNextPoint(v1[0],v1[1],v1[2]);
+ pericellPoints->InsertNextPoint(v2[0],v2[1],v2[2]);
+ pericellPoints->InsertNextPoint(v3[0],v3[1],v3[2]);
+ pericellPoints->InsertNextPoint(v4[0],v4[1],v4[2]);
+ pericellPoints->InsertNextPoint(v5[0],v5[1],v5[2]);
+ pericellPoints->InsertNextPoint(v6[0],v6[1],v6[2]);
+ pericellPoints->InsertNextPoint(v7[0],v7[1],v7[2]);
+ vtkSmartPointer<vtkHexahedron> h = vtkSmartPointer<vtkHexahedron>::New();
+ vtkIdList* l = h->GetPointIds();
+ for (int i=0; i<8; i++) {
+ l->SetId(i,i);
+ }
+ pericellHexa->InsertNextCell(h);
+ }
+
vtkSmartPointer<vtkDataCompressor> compressor;
if(compress) compressor=vtkSmartPointer<vtkZLibDataCompressor>::New();
@@ -464,7 +497,7 @@
}
if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
-
+
#ifdef YADE_VTK_MULTIBLOCK
if(!multiblock)
#endif
@@ -549,6 +582,23 @@
writer->Write();
}
}
+ vtkSmartPointer<vtkUnstructuredGrid> pericellUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
+ if (recActive[REC_PERICELL]){
+ pericellUg->SetPoints(pericellPoints);
+ pericellUg->SetCells(12,pericellHexa);
+ #ifdef YADE_VTK_MULTIBLOCK
+ if(!multiblock)
+ #endif
+ {
+ vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+ if(compress) writer->SetCompressor(compressor);
+ if(ascii) writer->SetDataModeToAscii();
+ string fn=fileName+"pericell."+lexical_cast<string>(scene->iter)+".vtu";
+ writer->SetFileName(fn.c_str());
+ writer->SetInput(pericellUg);
+ writer->Write();
+ }
+ }
#ifdef YADE_VTK_MULTIBLOCK
if(multiblock){
vtkSmartPointer<vtkMultiBlockDataSet> multiblockDataset = vtkSmartPointer<vtkMultiBlockDataSet>::New();
@@ -556,6 +606,7 @@
if(recActive[REC_SPHERES]) multiblockDataset->SetBlock(i++,spheresUg);
if(recActive[REC_FACETS]) multiblockDataset->SetBlock(i++,facetsUg);
if(recActive[REC_INTR]) multiblockDataset->SetBlock(i++,intrPd);
+ if(recActive[REC_PERICELL]) multiblockDataset->SetBlock(i++,pericellUg);
vtkSmartPointer<vtkXMLMultiBlockDataWriter> writer = vtkSmartPointer<vtkXMLMultiBlockDataWriter>::New();
if(ascii) writer->SetDataModeToAscii();
string fn=fileName+lexical_cast<string>(scene->iter)+".vtm";
@@ -564,6 +615,8 @@
writer->Write();
}
#endif
+
+
};
void VTKRecorder::addWallVTK (vtkSmartPointer<vtkQuad>& boxes, vtkSmartPointer<vtkPoints>& boxesPos, Vector3r& W1, Vector3r& W2, Vector3r& W3, Vector3r& W4) {
=== modified file 'pkg/dem/VTKRecorder.hpp'
--- pkg/dem/VTKRecorder.hpp 2013-08-29 10:30:31 +0000
+++ pkg/dem/VTKRecorder.hpp 2013-10-04 15:25:30 +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_WPM};
+ 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_WPM,REC_PERICELL};
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.",
=== modified file 'py/export.py'
--- py/export.py 2013-10-04 12:04:43 +0000
+++ py/export.py 2013-10-04 15:25:30 +0000
@@ -255,7 +255,10 @@
outFile = open(fName, 'w')
outFile.write("# vtk DataFile Version 3.0.\n%s\nASCII\n\nDATASET POLYDATA\nPOINTS %d double\n"%(comment,n))
for b in bodies:
- pos = b.state.refPos if useRef else b.state.pos
+ if O.periodic:
+ pos = O.cell.wrap(O.cell.b.state.refPos if useRef else b.state.pos)
+ else:
+ pos = O.cell.b.state.refPos if useRef else b.state.pos
outFile.write("%g %g %g\n"%(pos[0],pos[1],pos[2]))
outFile.write("\nPOINT_DATA %d\nSCALARS radius double 1\nLOOKUP_TABLE default\n"%(n))
for b in bodies:
@@ -523,6 +526,36 @@
outFile.close()
self.intrsSnapCount += 1
+ def exportPeriodicCell(self,comment="comment",numLabel=None):
+ """exports spheres (positions and radius) and defined properties.
+
+ :param string comment: comment to add to vtk file
+ :param int numLabel: number of file (e.g. time step), if unspecified, the last used value + 1 will be used
+ """
+ if not O.periodic:
+ print 'WARNING: export.VTKExporter.exportPeriodicCell: scene is not periodic, no export...'
+ return
+ hSize = O.cell.hSize
+ fName = self.baseName+'-periCell-%04d'%(numLabel if numLabel else self.intrsSnapCount)+'.vtk'
+ outFile = open(fName, 'w')
+ outFile.write("# vtk DataFile Version 3.0.\n%s\nASCII\n\nDATASET UNSTRUCTURED_GRID\nPOINTS 8 double\n"%(comment))
+ vertices = [
+ hSize*Vector3(0,0,1),
+ hSize*Vector3(0,1,1),
+ hSize*Vector3(1,1,1),
+ hSize*Vector3(1,0,1),
+ hSize*Vector3(0,0,0),
+ hSize*Vector3(0,1,0),
+ hSize*Vector3(1,1,0),
+ hSize*Vector3(1,0,0),
+ ]
+ for v in vertices:
+ outFile.write('%g %g %g\n'%(v[0],v[1],v[2]))
+ outFile.write('\nCELLS 1 9\n')
+ outFile.write('8 0 1 2 3 4 5 6 7\n')
+ outFile.write('\nCELL_TYPES 1\n12\n')
+ outFile.close()
+
#gmshGeoExport===============================================================
