CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   How to compute cellZone volume (http://www.cfd-online.com/Forums/openfoam-programming-development/69661-how-compute-cellzone-volume.html)

vinz October 30, 2009 04:35

How to compute cellZone volume
 
Dear OpenFoamers,

I am currently trying to develop a solver derived from simpleFoam. In this one I add a source term S which is zero almost everywhere except on a cellZone. In my equation I would need to divide S by the volume of this cellZone.
Does anybody know how I could compute this volume in the code. For example, I know there is something like mesh.V() in order to get the volume of one cell. How should I proceed to get the volume of my whole cellZone.
Any help would be appreciated.

Regards,

Vincent

ngj October 30, 2009 04:53

Hi Vincent

I have never worked with cellZones, however, I do expect that addressing these will give you a list of cell-numbers referring to the mesh. Then loop over these indices and sum the volumes.

E.g.

scalar vol(0);
forAll(cellZonesIndices, cI)
{
vol += mesh.V()[cellZonesIndices[cI]];
}

Best regards

Niels

vinz October 30, 2009 05:30

Hi Niels,

And thanks a lot for your fast reply.
I am going to try your solution which seems rather logical.
However, one more question, what do you mean by "adressing" the cellZone? and how do you get the cellIndices?
I guess it is a newbie question but I never went deep into the openFoam code, so I don't know how to do it.
If you have already some hints, that would be nice!

Regars,

Vincent

vinz October 30, 2009 06:01

I tried this in my code since the name of my zone in cellZones file is "boundaryZone"

word boundaryZone;
label cellZoneID = mesh.cellZones().findZoneID(boundaryZone);
const cellZone& zone = mesh.cellZones()[cellZoneID];
const cellZoneMesh& zoneMesh = zone.zoneMesh();
const labelList& cellsZone = zoneMesh[cellZoneID];


scalar cellZoneVol(0);
forAll(cellsZone, cI)
{
cellZoneVol += mesh.V()[cellsZone[cI]];
}

It compiles but crashes straight at the first iteration with this message:

#0 Foam::error::printStack(Foam::Ostream&) in "/CFDBB/vincent/OpenFOAM//OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::sigSegv::sigSegvHandler(int) in "/CFDBB/vincent/OpenFOAM//OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 __restore_rt at sigaction.c:0
#3 Foam::cellZone::zoneMesh() const in "/CFDBB/vincent/OpenFOAM//OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#4 main in "/CFDBB/vincent/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/adSimpleFoam"
#5 __libc_start_main in "/lib64/libc.so.6"
#6 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/CFDBB/vincent/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/adSimpleFoam"
Segmentation fault

Any idea on what is going on?
Thanks!

vinz October 30, 2009 11:03

Ok I reply to myself, here is the solution which works now:

const label cellZoneID = mesh.cellZones().findZoneID("boundaryZone");
const cellZone& zone = mesh.cellZones()[cellZoneID];
const cellZoneMesh& zoneMesh = zone.zoneMesh();
const labelList& cellsZone = zoneMesh[cellZoneID]; //list of all radiatorZone cell ID's

scalar cellZoneVol(0);
forAll(cellsZone, cI)
{
cellZoneVol += mesh.V()[cellsZone[cI]];
}

Info<< "cellZoneVol = " << cellZoneVol << nl << endl;

tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
);

UEqn().relax();

eqnResidual = solve
(
UEqn() == -fvc::grad(p) - force/cellZoneVol
).initialResidual();


However this solution only works in serial mode. If I decompose my case and try to run it in parallel it crashes since apparently the cellZoneVol is equal to 0 at some point??
Does someone know how to run this in parallel?

Regards,

Vincent

fra76 November 1, 2009 12:21

Just reduce the value of the volume after computing it, so that all the processors have the whole result (I've added output command to show you the result):

scalar cellZoneVol(0);
forAll(cellsZone, cI)
{
cellZoneVol += mesh.V()[cellsZone[cI]];
}

Pout << "before reduce: cellZoneVol = " << cellZoneVol << endl;
reduce(cellZoneVol,sumOp<scalar>());
Pout << "after reduce: cellZoneVol = " << cellZoneVol << endl;


Hope this helps,
Francesco


vinz November 2, 2009 01:26

Hi everybody,

And Thanks a lot Francesco, that works perfectly now!
I hope I'll get some results that I can post here later on.

Regards,

Vincent

rob3rt July 13, 2010 10:50

Dear Vincent,

How to compute cellZone area? my problem is similar like you had but only this time S = area/volume. I know there is a function like mag(mesh.Sf().boundaryField()[cellZoneID][cI]), but that is for the boundary area i think. My object is a hemisphere.

Please anyone give me a feedback. Very much appreciate it.

Kind Regards,
Robert.

olesen July 15, 2010 02:32

Quote:

Originally Posted by rob3rt (Post 267069)
Dear Vincent,

How to compute cellZone area? my problem is similar like you had but only this time S = area/volume. I know there is a function like mag(mesh.Sf().boundaryField()[cellZoneID][cI]), but that is for the boundary area i think. My object is a hemisphere.

Please anyone give me a feedback. Very much appreciate it.

Here's a rough idea of how you could calculate the outside surface area of a cellZone:


Code:

const cellList& cells  = mesh.cells();
const vectorField& faceAreas = mesh.faceAreas();

List<bool> outsideFaces(faceAreas.size(), false);

forAll(mesh.cellZones(), zoneI)
{
    const labelList& cellLabels = mesh_.cellZones()[zoneI];
    outsideFaces = false;

    // mark all faces that are NOT internal to the cellZone:
    forAll(cellLabels, i)
    {
        const cell& c = cells[cellLabels[i]];
        forAll(c, cFaceI)
        {
            const label faceI = c[cFaceI];

            // xor operation
            // internal faces get marked twice, outside faces get marked once
            if (outsideFaces[faceI])
            {
                outsideFaces[faceI] = false;
            }
            else
            {
                outsideFaces[faceI] = true;
            }
        }
    }

    // now calculate the area
    scalar zoneOutsideArea = 0;
    label  zoneOutsideNFaces = 0;

    forAll(outsideFaces, faceI)
    {
        if (outsideFaces[faceI])
        {
            zoneSurfaceArea += mag(faceAreas[faceI]);
            zoneOutsideNFaces++;
        }
    }

    Info<<"zone:" << zoneI
        << " nFaces:" << zoneOutsideNFaces
        << " area:" << zoneOutsideArea << endl;
}


NOTE: this code has not been tested (at all). You'll need to add some extra operations to ensure that it also works in parallel. Nonetheless, it should help get you going.

maximeg July 30, 2010 16:31

Hi all,
Anybody knows how to access the normal vector of a specific face (facei) of a specific -cubical- cell (celli)?
It should look like something like this:
vector normal_vector=mesh[celli].faceAreas[facei];
Because they are all cubical cells, facei varies from 0 to 5.
Thanks a lot for your help!
Maxime

su_junwei July 31, 2010 02:46

Quote:

Originally Posted by maximeg (Post 269588)
Hi all,
Anybody knows how to access the normal vector of a specific face (facei) of a specific -cubical- cell (celli)?
It should look like something like this:
vector normal_vector=mesh[celli].faceAreas[facei];
Because they are all cubical cells, facei varies from 0 to 5.
Thanks a lot for your help!
Maxime

If your faceI is a globe index or a local index.

For fvMesh
if a globe one, you can use
mesh.Sf()[facei];
a local one
convert it to globe face index
mesh.Sf()[mesh.cells()[cellI][faceI] ]

for polyMesh
you can also use mesh.faceArea() instead of mesh.Sf()

Junwei

maximeg August 2, 2010 17:38

Quote:

Originally Posted by su_junwei (Post 269607)
If your faceI is a globe index or a local index.

For fvMesh
if a globe one, you can use
mesh.Sf()[facei];
a local one
convert it to globe face index
mesh.Sf()[mesh.cells()[cellI][faceI] ]

for polyMesh
you can also use mesh.faceArea() instead of mesh.Sf()

Junwei

Hello Junwei,
Everything works perfectly!
Thanks again for your help.
Maxime


All times are GMT -4. The time now is 21:53.