[IC01]

Hello, Foamers

I have processed icoUncoupledKinematicParticle tracking and now I need to extract and process the particle data that landed on a patch
I'm using patchPostProcessing, I need separate data every deltaT

cloudFunctions
{
    patchPostProcessing1
    {
        type            patchPostProcessing;
        maxStoredParcels 1e7;
        patches         ( fixedWalls );
    }

It all works fine but output files have the Barycentric coordinates. The question is how do I convert it? Wither with searate code or in OpenFoam?

I'm okay with OpenFoam but I'm bad in programming in C++

Thank you
Regards
IC01

[Mars409]

From particles.C, line 473 (OF v10), https://cpp.openfoam.org/v10/particl...ce.html#l00473, the particle constructor:

 Foam::particle::particle

 (

     const polyMesh& mesh,

     const barycentric& coordinates,

     const label celli,

     const label tetFacei,

     const label tetPti

  )

The positions file is a list of lines that look like:
(a b c d) cellId faceId vertexInFaceId
where a, b, c, d are the weights (total=1) at the 4 vertices of the tetrahedron, cellId is the cell Id, faceId is the face Id, and vertexInFaceId is the index of a vertex in the face identified by faceId.


The particleIO.C's writePosition() member function (https://cpp.openfoam.org/v10/particl...ce.html#l00087) is what wrote the entries in the positions file:

 void Foam::particle::writePosition(Ostream& os) const

 {

     if (os.format() == IOstream::ASCII)

     {

         os  << coordinates_

             << token::SPACE << celli_

             << token::SPACE << tetFacei_

             << token::SPACE << tetPti_;

     }

     else

     {

         os.write(reinterpret_cast<const char*>(&coordinates_), sizeofPosition_);

     }

 

     // Check state of Ostream

     os.check("particle::writePosition(Ostream& os, bool) const");

 }

The reference tetrahedron has the first vertex at the centroid of the cell ("center").

The 2nd vertex of the tet is the first vertex of the face ("base").

The 3rd vertex is identified by vertexInFaceId.

The 4th vertex is the next vertex in the list of vertices belonging to the face.

With the 4 vertices of the tet all identified and their coordinates found, the cartesian coordinates can be computed.

Using this info, one can write a Python script to go into the owner, faces, and points files found under constant/polymesh to read/compute the coordinates of the 4 vertices of the tet.


A cell id can appear multiple times in the 'owner' and 'neighbour' lists, the sum of which is the number of faces a cell has. The index into the list is the face id. Collect all face ids for the cell from both the 'owner' and 'neighbour' lists.

To find all the vertices of a cell, go into the 'faces' file, which is a list of vertex lists, one for each face. Merge all the vertex lists for all the faceIds collected, removing repeated entries, and there you have all the unique vertex ids of the cell.

Using these vertex ids, get all the vertices' coordinates from the 'points' file.

Average all the x, y, and z values respectively and you have the coordinates of the center (the 1st vertex of the tetrahedron).

Next, we get the coordinates of the base as the coordinates of the 1st vertex of the face identified by the face id given as the 2nd last entry in the current row of the positions file.

Next, we get the coordinates of the 3rd vertex as the coordinates of the vertex whose id is given as the last entry of the current row of the positions file.

Finally, the 4th vertex is the next vertex in the list of vertices of the face.

(Note that the script needs to check the 'owner' file to see if the face is from the 'neighbour' list, the determinant of the linear equation to solve for a, b, c, d needs to be negated.)

[Mars409]

Below code inline function from tetIndicesI.H (https://cpp.openfoam.org/v10/tetIndi...ce.html#l00067) shows how the base point index, and the 3rd and 4th point indices are obtained. (Note, however, that it does NOT assume the base to be the first vertex of the list of vertices that constitutes the face. This, however, contradicts the locate() function (https://cpp.openfoam.org/v10/particl...ource.html#407).)

 inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const

 {

     const Foam::face& f = mesh.faces()[face()];

 

     label faceBasePtI = mesh.tetBasePtIs()[face()];

 

     if (faceBasePtI < 0)

     {

         static labelHashSet badFaces;

         static label badTimeIndex = -1;

 

         if (badTimeIndex != mesh.time().timeIndex())

         {

             badFaces.clear();

             badTimeIndex = mesh.time().timeIndex();

         }

 

         if (!badFaces[face()])

         {

             WarningInFunction

                 << "No base point for face " << face() << ", " << f

                 << ", produces a valid tet decomposition." << endl;

 

             badFaces.set(face());

         }

 

         faceBasePtI = 0;

     }

 

     label facePtI = (tetPt() + faceBasePtI) % f.size();

     label faceOtherPtI = f.fcIndex(facePtI);

 

     if (mesh.faceOwner()[face()] != cell())

     {

         Swap(facePtI, faceOtherPtI);

     }

 

     return triFace(f[faceBasePtI], f[facePtI], f[faceOtherPtI]);

 }

The above faceTriIs() function is called by stationaryTetGeometry() from particleI.H () to determine the cartesian coordinates of the three face triangle vertices by using the point indices to look up the point cartesian coordinates, besides looking up the cell center cartesian coordinates using the cell index in a precomputed cell centers list ccs.

 void Foam::particle::stationaryTetGeometry

 (

     vector& centre,

     vector& base,

     vector& vertex1,

     vector& vertex2

 ) const

 {

     const triFace triIs(currentTetIndices().faceTriIs(mesh_));

     const vectorField& ccs = mesh_.cellCentres();

     const pointField& pts = mesh_.points();

 

     centre = ccs[celli_];

     base = pts[triIs[0]];

     vertex1 = pts[triIs[1]];

     vertex2 = pts[triIs[2]];

 }

In the same particleI.H, as also a private member function, stationaryTetTransform() put together a 3 row by 4 column matrix that when multiplies a 4-by-1 barycentric column vector from the left gives a 3-by-1 cartesian coordinate vector.

 inline Foam::barycentricTensor Foam::particle::stationaryTetTransform() const

 {

     vector centre, base, vertex1, vertex2;

     stationaryTetGeometry(centre, base, vertex1, vertex2);

 

     return barycentricTensor(centre, base, vertex1, vertex2);

 }

If the mesh is moving/stationary, the currentTetTransform() selects movingTetTransform(0)[0]/stationaryTetTransform():

inline Foam::barycentricTensor Foam::particle::currentTetTransform() const
 
{
 
    if (mesh_.moving() && stepFraction_ != 1)
 
    {
 
        return movingTetTransform(0)[0];
 
    }
 
    else
 
    {
 
        return stationaryTetTransform();
 
    }
 
}

And the particle private member function position() indeed multiplies that matrix with the barycentric coordinates to return a cartesian position vector (https://cpp.openfoam.org/v10/particl...ce.html#l00278) (https://www.openfoam.com/documentati...3ca5130a7225):

 inline Foam::vector Foam::particle::position() const

 {

     return currentTetTransform() & coordinates_;

 }

The particle private member function position() is called at the beginning in particle.C's tracking function to convert the particle barycentric coordinates into cartesian position vector (https://cpp.openfoam.org/v10/particl....html#l00658):

Foam::scalar Foam::particle::trackToStationaryTri
 (

     const vector& displacement,

     const scalar fraction,

     label& tetTriI

 )

 {

     const vector x0 = position();

     const vector x1 = displacement;

     const barycentric y0 = coordinates_;

 

     if (debug)

     {

         Info<< "Particle " << origId() << endl << "Tracking from " << x0

             << " along " << x1 << " to " << x0 + x1 << endl;

     }

Is there a way to harness the position() private function indirectly without altering the particle.C file so that we can get the barycentric coordinates and converted into Cartesian position vectors?


There is. The particle public member function deviationFromMeshCentre(): (https://cpp.openfoam.org/v10/particl...ce.html#l01008)

 Foam::vector Foam::particle::deviationFromMeshCentre() const

 {

     if (cmptMin(mesh_.geometricD()) == -1)

     {

         vector pos = position(), posC = pos;

         meshTools::constrainToMeshCentre(mesh_, posC);

         return pos - posC;

     }

     else

     {

         return vector::zero;

     }
 }

The following code snippet looks like will do the trick for each particle p:

    const point& min = mesh().bounds().min();
    const point& max = mesh().bounds().max();         
    vector posC ((min + max)/2);
    
    meshTools::constrainToMeshCentre(mesh(), posC);
    return posC + p.deviationFromMeshCentre();

I imagine a coded function object using this code snippet will do the job, written ala the cloud writeData() public member function in IOPosition.C (https://cpp.openfoam.org/v10/IOPosit...ce.html#l00058), swapping out the iter().writePosition() with the above code snippet wrapped in a function and with p replaced by iter().

 template<class CloudType>

 bool Foam::IOPosition<CloudType>::writeData(Ostream& os) const

 {

     os  << cloud_.size() << nl << token::BEGIN_LIST << nl;

 

     forAllConstIter(typename CloudType, cloud_, iter)

     {

         iter().writePosition(os);

         os  << nl;

     }

 

     os  << token::END_LIST << endl;

 

     return os.good();

 }

[Mars409]

An even simpler & more direct way to get the cartesian coordinates of the particles in a cloud seems to be provided by file Cloud.C's Cloud class' writePositions() public member function: (https://cpp.openfoam.org/v10/Cloud_8...ce.html#l00441)

template<class ParticleType>
void Foam::Cloud<ParticleType>::writePositions() const
{
    OFstream pObj
    (
        this->db().time().path()/this->name() + "_positions.obj"
    );

    forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
    {
        const ParticleType& p = pIter();
        pObj<< "v " << p.position().x() << " " << p.position().y() << " "
            << p.position().z() << nl;
    }

    pObj.flush();
}

[Mars409]

Modify foamFormatConvert by inserting call to Cloud::writePositions() will write out cartesian coordinates in <cloud>_positions.obj files. All write function calls in the original code may be commented out.

#include <cstdio>        // for rename(), strerror()
#include <cstring>        // for c_string()
#include <cerrno>        // for errno



        parcels.writePositions();  // Write <cloud>_positions.obj file

                    /*
                     * Cloud::writePositions() wrote the <cloud>_positions.obj 
                     *   file in the case directory. Now, move it to the <time>/
                     *   lagrangian/<cloud name>/ directory where lagrangian 
                     *   field files of the named cloud belongs.
                     */
        fileName    fullToPath{lagrangianDir / parcels.name() };

        fileName    fullFromPath(parcels.time().path());

        word        objFileName(parcels.name() + "_positions.obj");

        if(std::rename( \
            (fullFromPath / objFileName).c_str(), \
            (fullToPath / objFileName).c_str()) < 0) {
            
            Info << "Failed to move " + \
            cloudName + "_positions.obj file from " + \
            fullFromPath + " to " + \
            fullToPath
            << strerror(errno) << endl;
        }

[sarina]

Hi Mars409, if I understand well, one of the Face IDs I find using the cellId must correspond to the faceId in the position file, right?

[Mars409]

Originally Posted by sarina

Hi Mars409, if I understand well, one of the Face IDs I find using the cellId must correspond to the faceId in the position file, right?

Yes, that's correct.

[sarina]

Thank you! Another question: if vertexInFaceID is equal to 1, then the 2th and the 3th vertex of the set will be the same, right?

[Mars409]

Originally Posted by sarina

Thank you! Another question: if vertexInFaceID is equal to 1, then the 2th and the 3th vertex of the set will be the same, right?

No, b/c the 2nd vertex of the tet is the 1st vertex of the face and this vertex has index 0 (index values begin at 0). If vertexInFaceID=1, the 3rd vertex is the 2nd vertex of the face, and the 4th vertex is the 3rd vertex of the face.