Permeability example customization


I am new to Palabos and am already a big fan. A big thank you and kudos to the development team!

I am trying to feed in my own voxelized geometry of random packed spheres into the examples/tutorials/permeability tutorial. I have a periodic sphere pack in a cube-shaped domain which I have voxelized into an (NX,NY,NZ) = (100,100,100) voxel pack. For better or worse, I produced my voxel inputs by treating the X and Y directions as periodic, and Z as non-periodic (the pack itself is periodic in all directions); setting sphere surfaces to be value = 1, interior sphere voxels = 2, and otherwise 0. This all checks out OK. Then I write the voxel pack to file by writing NZ planes, each with NY rows of NX values. This written file is the one I use for input into the permeability tutorial example code.

I suspect that there is a mismatch between the orientation of my geometry on disk vs. how it gets read in, but if anyone can provide any insight on the following questions, it would be deeply appreciated.

  1. New (seemingly extraneous) geometry I/O code? - Starting in version 1.4, a new section of code in the I/O section attempts to read a bunch of slices into the geometry:


    MultiScalarField3D<int> geometry(nx,ny,nz);
    MultiScalarField3D<int> slice(1,ny,nz);
    for (iX=0; iX<nx-1; ++iX) {
        string fname = createFileName("slice_", iX, 4)+"_truc.dat";
        pcout << "Reading slice " << fname;
        plb_ifstream geometryFile(fNameIn.c_str());
        if (!geometryFile.is_open()) {
     	    pcout << "Error: could not open geometry file " << fNameIn << endl;
                return -1;
        geometryFile >> slice;
        copy(slice, slice.getBoundingBox(), geometry, Box3D(iX,iX,0,ny-1,0,nz-1));


I have been unable to find any documentation on this section of code, nor do I find any utility or procedure for generating the slices that this section of code is trying to read.   Looking back at earlier versions of Palabos, I see that this code block is new since version 1.4 and I assume it can be safely removed from the example/tutorial code.  Is that true?  Actually, it won't even compile with this block of code in there since there is a duplicate declaration of the geometry object.

2. I have tried (to no avail) to understand how the voxel data is parsed from the *.dat file. My data is written in rows with X,Y,Z in order of most to least quickly increment. [in other words, i write NZ X/Y planes, each with NY rows of NX values] Can anyone let me know how the data should be ordered in the *.dat file?  

I have tried to parse the createDAT.m file, but it is a little confusing to me.  It actually seems to write things in Y/Z planes or planes normal to the flow direction.  I don't have access to any software that can actually use the *.m file to do the preprocessing.

Thanks (in advance) for your help!