CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Post-Processing

Lagrangian PatchPost Processing in Global coordinates

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 18, 2021, 00:44
Default Lagrangian PatchPost Processing in Global coordinates
  #1
New Member
 
Illia
Join Date: Aug 2020
Posts: 10
Rep Power: 4
IC01 is on a distinguished road
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
Code:
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
IC01 is offline   Reply With Quote

Old   September 18, 2022, 19:28
Default
  #2
New Member
 
Join Date: May 2020
Posts: 27
Blog Entries: 1
Rep Power: 4
Mars409 is on a distinguished road
From particles.C, line 473 (OF v10), https://cpp.openfoam.org/v10/particl...ce.html#l00473, the particle constructor:

Code:
 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:

Code:
 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.)

Last edited by Mars409; September 20, 2022 at 13:23. Reason: Corrections and additional details.
Mars409 is offline   Reply With Quote

Old   September 20, 2022, 07:40
Default
  #3
New Member
 
Join Date: May 2020
Posts: 27
Blog Entries: 1
Rep Power: 4
Mars409 is on a distinguished road
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).)
Code:
 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.
Code:
 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.
Code:
 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():
Code:
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):
Code:
 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):
Code:
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)

Code:
 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:
Code:
    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().

Code:
 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();

 }

Last edited by Mars409; September 20, 2022 at 13:57. Reason: Corrections
Mars409 is offline   Reply With Quote

Old   September 21, 2022, 12:52
Default
  #4
New Member
 
Join Date: May 2020
Posts: 27
Blog Entries: 1
Rep Power: 4
Mars409 is on a distinguished road
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)
Code:
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 is offline   Reply With Quote

Reply

Tags
barycentric, lagrangian particle, patchpostprocessing

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lagrangian data - particle velocity post processing Ake OpenFOAM Post-Processing 12 August 16, 2020 08:54
Lagrangian post processing in OpenFOAM 5 + jkeep OpenFOAM Post-Processing 0 June 17, 2019 02:07
time step continuity problem in VAWT simulation lpz_michele OpenFOAM Running, Solving & CFD 5 February 22, 2018 20:50
Help for the small implementation in turbulence model shipman OpenFOAM Programming & Development 25 March 19, 2014 11:08
Conversion of computed velocities in Lagrangian coordinates to Eulerian coordinates SMM STAR-CD 4 June 5, 2011 10:20


All times are GMT -4. The time now is 18:15.