[hyperion]

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.

[MartinB]

Hi Daniel,

you can try this:

    volScalarField cv
    (
            IOobject
            (
                    "cv",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar("zero",dimVolume,0.0)
    );
    
    cv.internalField() = mesh.V();

Martin

[hyperion]

Hi Martin - That did the trick. Thank you very much.

[Andrea_85]

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

[MartinB]

Hi Andrea,

you can loop over the cells volumes and print them with

	forAll(mesh.V(),celli)
	{
	  Info << mesh.V()[celli] << nl;
	}

Martin

[Andrea_85]

Hi Martin
thanks for the reply. Where i have to put these lines?

Andrea

[hyperion]

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

[MartinB]

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

[Andrea_85]

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

[MartinB]

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

        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

[Andrea_85]

Thanks a lot for your help.

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


andrea

[MartinB]

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

[Andrea_85]

Great...it runs...thanks a lot again!

andrea

[Andrea_85]

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

[MartinB]

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

[Andrea_85]

Thanks Martin

Your advice were very helpful!!

Andrea

[Andrea_85]

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

[MartinB]

Hi Andrea,

you might have a look at this thread:
http://www.cfd-online.com/Forums/ope...ergence-2.html

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

[Andrea_85]

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

[MartinB]

Hi Andrea,

some minor problems removed:

// 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