particle position

Asked by jamespaul on 2019-05-28

Hi,

I want to export each particle's position in txt.In def addPlotData():,I used x,y,z = b.state.pos and plot.saveDataTxt to save positon,but I found that it only saved just one x/y/z.I think maybe saveDataTxt can only save int type.How can I save all particle's positons in txt?

Thanks

James

#################
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from yade import pack,plot,qt

sphere=O.materials.append(FrictMat(young=6e6,poisson=0.3,density=2600,frictionAngle=float(atan(0.4)),label='sphere'))
wall=O.materials.append(FrictMat(young=6e6,poisson=0.3,density=2600,frictionAngle=float(atan(0.1)),label='wall'))

r=2.2275e-3
sp=pack.SpherePack()
ylength=r*2*45
zlength=r*2*200
sp.makeCloud(minCorner=(0,0,0),maxCorner=(0,ylength,zlength),rMean=r,distributeMass=False)
sp.toSimulation(material='sphere')

for b in O.bodies:
    if isinstance(b.shape,Sphere): b.state.blockedDOFs = 'ZxY'

wallbutton=O.bodies.append(utils.wall(0,axis=2,sense=0,color=(1,0,0),material='wall'))
wallleft=O.bodies.append(utils.wall(0,axis=1,sense=0,color=(1,0,0),material='wall'))
wallright=O.bodies.append(utils.wall(ylength,axis=1,sense=0,color=(1,0,0),material='wall'))

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.2,label='newton'),
   VTKRecorder(fileName='3d-vtk-',recorders=['all'],virtPeriod=0.2),#save data for Paraview
   PyRunner(command='addPlotData()',iterPeriod=1000, label="apd"),
]

O.dt=1*PWaveTimeStep()

O.trackEnergy=True

cut=0.4

x=y=z=0

def addPlotData():
    for b in O.bodies:
         if isinstance(b.shape,Sphere):
   x,y,z = b.state.pos
       print(x,y,z)
    plot.addData(i=O.iter,kineticEnergy=utils.kineticEnergy(),y=y)
    plot.saveDataTxt('f.txt',vars=('i','y'))

plot.plots={'i':(('unbalancedForce','g'),None,'kineticEnergy')}
plot.plot()

qt.Controller()
O.saveTmp()
O.run()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-06-03
Last query:
2019-06-03
Last reply:
2019-05-29

Hi,

Have a look at VTK recorder https://yade-dem.org/doc/yade.wrapper.html?highlight=vtk#yade.wrapper.VTKRecorder.recorders

recorder "sphere" record position and radius of all particles.

William

jamespaul (jamespauljames) said : #2

Thanks William,

I know how to use VTK and paraview.But I want to export in txt for analysizing data in matlab.

James

Robert Caulk (rcaulk) said : #3

plot.AddData is used for plotting data. plot.SaveDataTxt is used for saving the plot data to a text file. These features are not meant to be used for plotting positions of all particles. If you really want to save particle positions to a text file you need to do it using standard python functionality:

f = open("positions.txt","w+")

for b in O.bodies:
    x,y,z = b.state.pos
    f.write("%d %d %d\n" % (x,y,z))

f.close()

Jan Stránský (honzik) said : #4

or standard yade functionality:
###
from yade import export
export.text("positions.txt")
###
Jan

[1] https://yade-dem.org/doc/yade.export.html#yade.export.text

You can read vtk data from matlab, using this function. You need first to compile it using mex:

/* Convert four uint8 to float64 value */

#include "mex.h"
#include "matrix.h"
#include <iostream>
#include <string>
#include <vtkSmartPointer.h>
#include <vtkXMLReader.h>
#include <vtkXMLUnstructuredGridReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkXMLStructuredGridReader.h>
#include <vtkXMLRectilinearGridReader.h>
#include <vtkXMLHyperOctreeReader.h>
#include <vtkXMLCompositeDataReader.h>
#include <vtkXMLStructuredGridReader.h>
#include <vtkXMLImageDataReader.h>
#include <vtkDataSetReader.h>
#include <vtkDataSet.h>
#include <vtkUnstructuredGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkHyperOctree.h>
#include <vtkImageData.h>
#include <vtkPolyData.h>
#include <vtkStructuredGrid.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkFieldData.h>
#include <vtkCellTypes.h>
#include <vtksys/SystemTools.hxx>

#include <map>

mxArray* readData(const char*);

/* The gateway function */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
/* variable declarations here */
    if(nrhs != 1)
    {
        mexErrMsgIdAndTxt("VtkUncompresserToolBox:TooMuchInput","One input required");
    }

    if(!mxIsChar(prhs[0]))
    {
        mexErrMsgIdAndTxt("VtkUncompresserToolBox:BadInput","Input must be a filename");
    }

    int buflen = mxGetNumberOfElements(prhs[0]) + 1;
    char* buf = (char*)mxCalloc(buflen, sizeof(char));

    /* Copy the string data from string_array_ptr and place it into buf. */
    if (mxGetString(prhs[0], buf, buflen) != 0)
    mexErrMsgIdAndTxt( "VtkUncompresserToolBox:invalidStringArray",
    "Could not convert string data.");

    plhs[0] = readData(buf);

/* code here */
}
//
// DumpXMLFile - report on the contents of an XML or legacy vtk file
// Usage: DumpXMLFile XMLFile1 XMLFile2 ...
// where
// XMLFile is a vtk XML file of type .vtu, .vtp, .vts, .vtr,
// .vti, .vto
//

template<class TReader> vtkDataSet *ReadAnXMLFile(const char*fileName)
{
  vtkSmartPointer<TReader> reader =
    vtkSmartPointer<TReader>::New();
  reader->SetFileName(fileName);
  reader->Update();
  reader->GetOutput()->Register(reader);
  return vtkDataSet::SafeDownCast(reader->GetOutput());
}

mxArray* readData(const char* filename)
{

  // Process each file on the command line

    vtkDataSet *dataSet;
    std::string extension =
      vtksys::SystemTools::GetFilenameLastExtension(filename);
    // Dispatch based on the file extension
    if (extension == ".vtu")
    {
        dataSet = ReadAnXMLFile<vtkXMLUnstructuredGridReader> (filename);
    }
    else if (extension == ".vtp")
    {
        dataSet = ReadAnXMLFile<vtkXMLPolyDataReader> (filename);
    }
    else if (extension == ".vts")
    {
        dataSet = ReadAnXMLFile<vtkXMLStructuredGridReader> (filename);
    }
    else if (extension == ".vtr")
    {
        dataSet = ReadAnXMLFile<vtkXMLRectilinearGridReader> (filename);
    }
    else if (extension == ".vti")
    {
        dataSet = ReadAnXMLFile<vtkXMLImageDataReader> (filename);
    }
    else if (extension == ".vto")
    {
        dataSet = ReadAnXMLFile<vtkXMLHyperOctreeReader> (filename);
    }
    else if (extension == ".vtk")
    {
        dataSet = ReadAnXMLFile<vtkDataSetReader> (filename);
    }
    else
    {
          mexErrMsgIdAndTxt("VtkUncompresserToolBox:UnknownExt","Invalid extension");
          return 0;
    }

    const char *fieldnames[5] = {"Points","CellCount","PointData","CellData","FieldData"};

    int numberOfCells = dataSet->GetNumberOfCells();
    int numberOfPoints = dataSet->GetNumberOfPoints();

    mxArray* data = mxCreateStructMatrix(1,1,5,fieldnames);

    mxArray* points = mxCreateDoubleMatrix(numberOfPoints,3,mxREAL);
    double* pointsPtr = mxGetPr(points);
    double pos[3];

    for(int i=0;i<numberOfPoints;i++)
    {
        dataSet->GetPoint(i,pos);
        pointsPtr[i] = pos[0];
        pointsPtr[i+numberOfPoints] = pos[1];
        pointsPtr[i+numberOfPoints*2] = pos[2];
    }

    mxSetFieldByNumber(data,0,0,points);

    //mxSetFieldByNumber(data,0,0,mxCreateDoubleScalar(numberOfPoints));
    mxSetFieldByNumber(data,0,1,mxCreateDoubleScalar(numberOfCells));

    mxArray *PD = mxCreateStructMatrix(1,1,0,0);
    vtkPointData *pd = dataSet->GetPointData();
    if (pd)
    {
        for (int i = 0; i < pd->GetNumberOfArrays(); i++)
        {
            int id = mxAddField(PD,pd->GetArrayName(i));

            vtkAbstractArray* d = pd->GetAbstractArray(i);

            if(d && d->GetDataType() == VTK_DOUBLE)
            {
                mxArray* dat = mxCreateDoubleMatrix(d->GetNumberOfTuples(),d->GetNumberOfComponents(),mxREAL);
                double *datPtr = mxGetPr(dat);
                double *vtkDat = (double*)d->GetVoidPointer(0);

                for(int j=0;j<d->GetNumberOfComponents();j++)
                {
                    for(int k=0;k<d->GetNumberOfTuples();k++)
                    {
                        datPtr[k+j*d->GetNumberOfTuples()] = vtkDat[j+k*d->GetNumberOfComponents()];
                    }
                }

                mxSetFieldByNumber(PD,0,id,dat);
               // mxSetFieldByNumber(PD,0,id,mxCreateDoubleScalar(d->GetDataType()));
            }
        }
    }
    mxSetFieldByNumber(data,0,2,PD);

    // Now check for cell data
    mxArray *CD = mxCreateStructMatrix(1,1,0,0);
    vtkCellData *cd = dataSet->GetCellData();
    if (cd)
    {
        for (int i = 0; i < cd->GetNumberOfArrays(); i++)
        {
            int id = mxAddField(CD,cd->GetArrayName(i));
            vtkAbstractArray* d = cd->GetAbstractArray(i);

            if(d && d->GetDataType() == VTK_DOUBLE)
            {
                mxArray* dat = mxCreateDoubleMatrix(d->GetNumberOfTuples(),d->GetNumberOfComponents(),mxREAL);
                double *datPtr = mxGetPr(dat);
                double *vtkDat = (double*)d->GetVoidPointer(0);

                for(int j=0;j<d->GetNumberOfComponents();j++)
                {
                    for(int k=0;k<d->GetNumberOfTuples();k++)
                    {
                        datPtr[k+j*d->GetNumberOfTuples()] = vtkDat[j+k*d->GetNumberOfComponents()];
                    }
                }

                mxSetFieldByNumber(CD,0,id,dat);
               // mxSetFieldByNumber(PD,0,id,mxCreateDoubleScalar(d->GetDataType()));
            }
        }
    }
    mxSetFieldByNumber(data,0,3,CD);

    // Now check for field data
    mxArray *FD = mxCreateStructMatrix(1,1,0,0);
    if (dataSet->GetFieldData())
    {
        for (int i = 0; i < dataSet->GetFieldData()->GetNumberOfArrays(); i++)
        {
            mxAddField(FD,dataSet->GetFieldData()->GetArray(i)->GetName());
        }
    }
    mxSetFieldByNumber(data,0,4,FD);

    dataSet->Delete();

    return data;
}

Robert Caulk (rcaulk) said : #6

Wow we are all making this way more complicated than it needs to be. Not sure how I forgot about export.text() Jan, thanks.

You could record the 3D view at runtime using screen capture, then use some image analysis tool for segmentation and reconstruction of the displacement field. I know a software running on Atari which could convert the output of this to some ascii that matlab may have a chance to decode.
;)

jamespaul (jamespauljames) said : #8

Thanks Jan and all the people who helped me,

Now I can import the export.text("positions.txt") in matlab.And I programmed a script to subscibe the positions.txt of different iter to plot the heat map to research the particle's displacement.But after I run the script,I think is there a function in yade that calculates the displacement?

Best Jan Stránský (honzik) said : #9

> think is there a function in yade that calculates the displacement?

of course there is :-)

### script saving current position, reference position and displacement
# a simple case to save
from yade import export
s1,s2 = sphs = (
   sphere((1,2,3),4),
   sphere((5,6,7),8),
)
s1.state.vel = (2,3,4)
s2.state.vel = (5,6,7)
O.bodies.append(sphs)
O.step()
# actual save
export.textExt( # [1]
   "dspl.dat", # file name
   format="x_y_z_r_attrs", # if you would like to load it again
   comment="xr yr zr dx dy dz", # 1st line of the file (x y z radius is there by default)
   attrs=["b.state.refPos","b.state.displ()"] # saving reference position and displacement vector
)
###

cheers
Jan

jamespaul (jamespauljames) said : #11

Thanks Jan Stránský, that solved my question.