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

write cell volumes

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

Like Tree6Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   December 8, 2010, 01:17
Default write cell volumes
  #1
Member
 
Daniel
Join Date: Jul 2010
Location: California
Posts: 39
Rep Power: 7
hyperion is on a distinguished road
Hi there - I am trying to modify the writeCellCentres application to get something that will write cell volumes. The heart of the code is as follows:

// Main program:

int main(int argc, char *argv[])

{
timeSelector::addOptions();

# include "setRootCase.H"

# include "createTime.H"

instantList timeDirs = timeSelector::select0(runTime, args);

# include "createMesh.H"

forAll(timeDirs, timeI)
{

runTime.setTime(timeDirs[timeI], timeI);

Info<< "Time = " << runTime.timeName() << endl;

// Check for new mesh

mesh.readUpdate();

scalarField cv
(
IOobject
(
"cellVolumes",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE

),
mesh.V()
);

Info<< "Writing cellCentre positions to " << cv.name() << " in "

<< runTime.timeName() << endl;

cv.write();



This code mainly differs from the writeCellCentres in that I replaced mesh.C() with mesh.V(). Also, where I define scalarField cv, the original code used volVectorField cc.

mesh.V() is used in a variety of places in openFOAM, but apparently it is not a scalarField or a volScalarField because I get a variety of compilation errors depending on which type I specify.

Does anyone know a good way to address this problem?

Many Thanks.
hyperion is offline   Reply With Quote

Old   December 8, 2010, 02:34
Default
  #2
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Daniel,

you can try this:
Code:
    volScalarField cv
    (
            IOobject
            (
                    "cv",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar("zero",dimVolume,0.0)
    );
    
    cv.internalField() = mesh.V();
Martin
Fine and mm.abdollahzadeh like this.
MartinB is offline   Reply With Quote

Old   December 8, 2010, 03:58
Default
  #3
Member
 
Daniel
Join Date: Jul 2010
Location: California
Posts: 39
Rep Power: 7
hyperion is on a distinguished road
Hi Martin - That did the trick. Thank you very much.
hyperion is offline   Reply With Quote

Old   December 8, 2010, 04:34
Default
  #4
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Hi,
i have the same problem. I want to know the volume of each cell and let the program prints them somewhere. With this trick where openFoam prints the volume?

Thank you
Andrea_85 is offline   Reply With Quote

Old   December 8, 2010, 05:04
Default
  #5
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Andrea,

you can loop over the cells volumes and print them with
Code:
	forAll(mesh.V(),celli)
	{
	  Info << mesh.V()[celli] << nl;
	}
Martin
MartinB is offline   Reply With Quote

Old   December 8, 2010, 05:11
Default
  #6
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Hi Martin
thanks for the reply. Where i have to put these lines?

Andrea
Andrea_85 is offline   Reply With Quote

Old   December 8, 2010, 05:22
Default
  #7
Member
 
Daniel
Join Date: Jul 2010
Location: California
Posts: 39
Rep Power: 7
hyperion is on a distinguished road
Hi Andrea - Just in case you want it I have attached the source file that I use to get cell Volumes (Thanks to Martin). I put the file in a directory I created called writeCellVolumes that is in the same directory as writeCellCentres. I also copied over the the other make files (Make/files and Make/options) from writeCellCentres then update those files with the correct EXE = path names then wmake it. The application is executed within a case directory. Good Luck
Attached Files
File Type: c writeCellVolumes.C (2.6 KB, 236 views)
David* and Nucleophobe like this.
hyperion is offline   Reply With Quote

Old   December 8, 2010, 05:25
Default
  #8
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Andrea,

it depends ;-)

Without knowing your solver and the operation you want to perform with cell volumes I can only guess:
- you might want to create a volScalarField that will contain the results
- you loop over the mesh as described above and calculate something important, which is then stored in your new volScalarField
- after the loop you let your volScalarField be written out to the current time directory

This action can be done after the solver is finished completely (saves time), or during each writing of time step, or even in a separated utility that can be run independently of the solver.
You have to make sure that all fields you need for the computation are valid, since some solver use "{ ... }" to end field's life early and reduce memory usage.

May be you can post or send your solver and your intended calculation...

Martin
MartinB is offline   Reply With Quote

Old   December 8, 2010, 05:51
Default
  #9
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Hi Martin

I just started using OpenFoam as part of my phd and so i am a very new user. At the moment i would like to simulate multiphase flows using interFoam, in simple geometry created with blockMesh. The geometry at the moment is a sphere so the cell volumes are not equal.

I just want the openFoam print volumes of each cell in a file, like "points" or "owner" in the polyMesh directory, named for example "Volumes".

Otherwise that openfoam prints after the simulation these volumes during post-processing, like "writeCellCentres" (something like writeCellVolumes).

I dont know which is the best way to do that.
Thanks

Andrea_85 is offline   Reply With Quote

Old   December 8, 2010, 06:57
Default
  #10
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Andrea,

if you don't use a dynamic mesh (i.e. cell volumes stay constant) you can write the values immediately before "Info<< "\nStarting time loop\n" << endl;" in interFoam.

To write the cell volumes to "constant/polyMesh":
Code:
        volScalarField cv
    (
            IOobject
            (
                    "Volumes",
	      runTime.constant(),
	      polyMesh::meshSubDir,
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar("zero",dimVolume,0.0)
    );
    
    cv.internalField() = mesh.V();
    cv.write();
Here are some more useful hints for handling files:
http://openfoamwiki.net/index.php/In...IOobject_class

Martin
MartinB is offline   Reply With Quote

Old   December 9, 2010, 04:23
Default
  #11
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Thanks a lot for your help.

After adding this line i guess that you have to recompile the program...or not?


andrea
Andrea_85 is offline   Reply With Quote

Old   December 9, 2010, 04:39
Default
  #12
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Yes, you must recompile. Simply call wmake in the interFoam directory to make a quick shot.

By the way it might be a good idea to make these kind of changes in your user application folder. Copy and rename interFoam to your favorite project name and do changes there... here is a good explanation how to do it:
http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam

Martin
MartinB is offline   Reply With Quote

Old   December 9, 2010, 05:18
Default
  #13
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Great...it runs...thanks a lot again!

andrea
Andrea_85 is offline   Reply With Quote

Old   January 12, 2011, 05:40
Default
  #14
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Hi,

i got another problem. Now I'm trying to import an external geometry using gmsh (so without using BlockMesh). I would like to do the same thing and print the volume of each cell in the constant/polymesh subdirectory. (while importing the geometry for example). Is possible? Where i have to put those line now?

Thanks

Andrea
Andrea_85 is offline   Reply With Quote

Old   January 12, 2011, 11:24
Default
  #15
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Andrea,
you could append the lines from post #10 to each mesh converting tool, of course... for gmshToFoam.C it would be directly before the line "Info<< "End\n" << endl;".

However it would make more sense to create a separate tool like Daniels "writeCellVolumes.C" from post #7. Replace the location where to write in the IOobject part.

In your every day work you call the mesh converter for gmsh, followed by a writeCellVolumes call.

Martin
MartinB is offline   Reply With Quote

Old   January 13, 2011, 11:52
Default
  #16
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Thanks Martin

Your advice were very helpful!!

Andrea
Andrea_85 is offline   Reply With Quote

Old   January 27, 2011, 10:58
Default
  #17
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
I continue this discussion because the new problem is similar to the old.
I am using the sample Dictionary to get the value of the color function alpha1 on a particular patch (wall for my problem). Is possible (similarly to obtain the volume of each cell) to let OF prints not only the value of the color function but also the area of the surface (in m^2) for each cell on the wall patch??

Thanks

Andrea
Andrea_85 is offline   Reply With Quote

Old   January 28, 2011, 08:27
Default
  #18
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Andrea,

you might have a look at this thread:
Improve simpleFoam convergence

In post #29 there is a piece of code that calculates the wall shear stress and the surface area for each cell surface on a given patch. Replace the wall shear stress by your alpha1 and you have your code... installation hints are given in this thread, too.

A parallel version is given in post #44.

I recommend, that you try to implement a uniprocessor version that you can call like the writeCellVolumes tool discussed earlier. The patch name of the wall you want to examine can be passed as a command line argument or via dictionary.

Martin
MartinB is offline   Reply With Quote

Old   January 28, 2011, 09:29
Default
  #19
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Thanks a lot Martin for all your help!!

I only need the values of cell's area on the patch and i can get alpha1 on the same patch directly from the sample.
So i changed your code (for parallel run, because i really need parallel) in this way. If I have not made mistakes should calculate the cell area on the patch "obstacle " and write it in area.dat.
Please check if there are any errors (sure there will be...)

I commented out all lines that I did not need..


#include "RASModel.H"
// write matlab file
{
Info << "Writing Area of the cells on the patch!" << endl;

word wallPatchName = "obstacle";
Info << "Searching for patch " << wallPatchName << endl;

// Check, if the patch name exists:
label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
if (patchID < 0)
{
Info << "No patch found with name \"" << wallPatchName << "\"!"
<< endl;
}
// else
// {
// calculating tau

//# include "createFields2.H"
// mesh.readUpdate();
// volSymmTensorField Reff(RASModel->devReff());
// volVectorField tau
// (
// IOobject
// (
// "tau",
// runTime.timeName(),
// mesh,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE
// ),
// mesh,
// dimensionedVector
// (
// "tau",
// Reff.dimensions(),
// vector::zero
// )
// );

// forAll(tau.boundaryField(), patchi)
// {
// tau.boundaryField()[patchi] =
// (
// -mesh.Sf().boundaryField()[patchi]
// /mesh.magSf().boundaryField()[patchi]
// ) & Reff.boundaryField()[patchi];
// }



// Pointer to the patch we want to have shear stress for:
const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
label nFaces = cPatch.size();

// pointField faceCenterCoords(nFaces, vector(0,0,0));
pointField faceAreas(nFaces, vector(0,0,0));
// pointField wallShearStress(nFaces, vector(0,0,0));

label runFaces = 0;

// collecting the interesting data
for (int i=0; i<nFaces; i++)
{
// faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i];
faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
// wallShearStress[runFaces] = tau.boundaryField()[patchID][i];
runFaces++;
}

// for parallel
reduce(nFaces, sumOp<label>());
Info << "Number of global faces: " << nFaces << endl;


// make a cute filename
fileName casePath = cwd();//casePath.name();
fileName caseName = casePath.name();
fileName myFile = "area.dat";
Info << "Data is written to file " << myFile << endl;

// create a new file
OFstream *myStream = NULL;
if (Pstream::master())
{
myStream = new OFstream(myFile);
}

// slaves send data to master, master writes data to file
if (Pstream:arRun())
{
if (Pstream::myProcNo() != Pstream::masterNo())
{
OPstream toMaster(Pstream::blocking, Pstream::masterNo());
// toMaster << faceCenterCoords;
toMaster << faceAreas;
// toMaster << wallShearStress;
}
else // master part
{
// first write own data
for (int i=0; i<faceCenterCoords.size(); i++)
{
*myStream << i << ";"
// << faceCenterCoords[i].component(0) << "; "
// << faceCenterCoords[i].component(1) << "; "
// << faceCenterCoords[i].component(2) << "; "

// << wallShearStress[i].component(0) << "; "
// << wallShearStress[i].component(1) << "; "
// << wallShearStress[i].component(2) << "; "

// << mag(wallShearStress[i]) << endl;

<< faceAreas[i] << "; "

}

// label run = faceCenterCoords.size();

// then receive slave data and write
// pointField bufferFaceCenterCoords;
pointField bufferFaceAreas;
// pointField bufferWallShearStress;


for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
{
IPstream fromSlave(Pstream::blocking, slave);
// fromSlave >> bufferFaceCenterCoords;
fromSlave >> bufferFaceAreas;
// fromSlave >> bufferWallShearStress;

for (int i=0; i<bufferFaceCenterCoords.size(); i++)
{
*myStream << run << ";"
// << bufferFaceCenterCoords[i].component(0) << "; "
// << bufferFaceCenterCoords[i].component(1) << "; "
// << bufferFaceCenterCoords[i].component(2) << "; "
//
// << bufferWallShearStress[i].component(0) << "; "
// << bufferWallShearStress[i].component(1) << "; "
// << bufferWallShearStress[i].component(2) << "; "

// << mag(bufferWallShearStress[i]) << endl;

<< bufferFaceAreas << "; "
run++;
}
}
}
}
else // single processor only
{
for (int i=0; i<nFaces; i++)
{
*myStream << i << "; "
// coordinates on face centers:
// << mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
// << mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
// << mesh.Cf().boundaryField()[patchID][i].component(2) << "; "

// wall shear stress as a vector:
// << tau.boundaryField()[patchID][i].component(0) << "; "
// << tau.boundaryField()[patchID][i].component(1) << "; "
// << tau.boundaryField()[patchID][i].component(2) << "; "

// magnitude of wall shear stress:
// << mag(tau.boundaryField()[patchID][i]) << endl;

<< mesh.Sf().boundaryField()[patchID][i] << "; "
}
}

}
}



Thanks
andrea
Andrea_85 is offline   Reply With Quote

Old   January 28, 2011, 10:34
Default
  #20
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Andrea,

some minor problems removed:
Code:
// write matlab file
{
    Info << "Writing Area of the cells on the patch!" << endl;

    word wallPatchName = "obstacle";
    Info << "Searching for patch " << wallPatchName << endl;

    // Check, if the patch name exists:
    label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
    if (patchID < 0)
    {
        Info << "No patch found with name \"" << wallPatchName << "\"!"
        << endl;
    }

    // Pointer to the patch we want to have shear stress for:
    const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
    label nFaces = cPatch.size();

    pointField faceAreas(nFaces, vector(0,0,0));

    label runFaces = 0;

    // collecting the interesting data
    for (int i=0; i<nFaces; i++)
    {
        faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
        runFaces++;
    }

    // for parallel
    reduce(nFaces, sumOp<label>());
    Info << "Number of global faces: " << nFaces << endl;

    // make a cute filename
    fileName casePath = cwd();//casePath.name();
    fileName caseName = casePath.name();
    fileName myFile = "area.dat";
    Info << "Data is written to file " << myFile << endl;

    // create a new file
    OFstream *myStream = NULL;
    if (Pstream::master())
    {
        myStream = new OFstream(myFile);
    }

    // slaves send data to master, master writes data to file
    if (Pstream::parRun())
    {
        if (Pstream::myProcNo() != Pstream::masterNo())
        {
            OPstream toMaster(Pstream::blocking, Pstream::masterNo());
            toMaster << faceAreas;
        }
        else // master part
        {
            // first write own data
            for (int i=0; i<faceAreas.size(); i++)
            {
                *myStream << i << ";"
                << mag(faceAreas[i]) << endl;

            }
        }
        label run = faceAreas.size();

        // then receive slave data and write
        pointField bufferFaceAreas;

        for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
        {
            IPstream fromSlave(Pstream::blocking, slave);
            fromSlave >> bufferFaceAreas;

            for (int i=0; i<bufferFaceAreas.size(); i++)
            {
                *myStream << run << "; "
                << mag(bufferFaceAreas[i]) << endl;
                run++;
            }
        }
    }
    else // single processor only
    {
        for (int i=0; i<nFaces; i++)
        {
            *myStream << i << "; "
            << mag(mesh.Sf().boundaryField()[patchID][i]) << endl;
        }
    }

}
Single processor version seems to work fine, but parallel version does not terminate correctly... you may have a look, I'm a little bit short of time right now...

Martin
MartinB is offline   Reply With Quote

Reply

Tags
cell volume, mesh.v, volume

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Cells with t below lower limit Purushothama CD-adapco 2 May 31, 2010 21:58
mesh file for flow over a circular cylinder Ardalan Main CFD Forum 6 April 17, 2010 23:40
how to write the value of every cell in a volume? Ralf Schmidt FLUENT 0 January 9, 2006 06:18
Trimmed cell and embedded refinement mesh conversion issues michele OpenFOAM Other Meshers: ICEM, Star, Ansys, Pointwise, GridPro, Ansa, ... 2 July 15, 2005 04:15
Warning 097- AB CD-adapco 6 November 15, 2004 05:41


All times are GMT -4. The time now is 10:26.