particle position

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

Thanks William,

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

James

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()

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;
}

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.
;)

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?

forgotten reference:
[1] https://yade-dem.org/doc/yade.export.html#yade.export.textExt

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