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

How to compute cellZone volume

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

Like Tree5Likes
  • 2 Post By vinz
  • 2 Post By fra76
  • 1 Post By su_junwei

Reply
 
LinkBack Thread Tools Display Modes
Old   October 30, 2009, 04:35
Default How to compute cellZone volume
  #1
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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
vinz is offline   Reply With Quote

Old   October 30, 2009, 04:53
Default
  #2
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
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
ngj is offline   Reply With Quote

Old   October 30, 2009, 05:30
Default
  #3
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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 is offline   Reply With Quote

Old   October 30, 2009, 06:01
Default
  #4
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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:rintStack(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!

Last edited by vinz; October 30, 2009 at 06:23.
vinz is offline   Reply With Quote

Old   October 30, 2009, 11:03
Default
  #5
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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
vinz is offline   Reply With Quote

Old   November 1, 2009, 13:21
Default
  #6
Senior Member
 
Francesco Del Citto
Join Date: Mar 2009
Location: Zürich Area, Switzerland
Posts: 215
Rep Power: 9
fra76 is on a distinguished road
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

fra76 is offline   Reply With Quote

Old   November 2, 2009, 02:26
Default
  #7
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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

Last edited by vinz; November 2, 2009 at 03:49.
vinz is offline   Reply With Quote

Old   July 13, 2010, 10:50
Default
  #8
New Member
 
Robert
Join Date: Mar 2010
Posts: 16
Rep Power: 7
rob3rt is on a distinguished road
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.
rob3rt is offline   Reply With Quote

Old   July 15, 2010, 02:32
Default
  #9
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
Quote:
Originally Posted by rob3rt View Post
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.
olesen is offline   Reply With Quote

Old   July 30, 2010, 16:31
Default
  #10
New Member
 
Join Date: Jun 2010
Posts: 14
Rep Power: 7
maximeg is on a distinguished road
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
maximeg is offline   Reply With Quote

Old   July 31, 2010, 02:46
Default
  #11
Senior Member
 
su_junwei's Avatar
 
su junwei
Join Date: Mar 2009
Location: Xi'an China
Posts: 151
Rep Power: 10
su_junwei is on a distinguished road
Send a message via MSN to su_junwei
Quote:
Originally Posted by maximeg View Post
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
Nucleophobe likes this.
su_junwei is offline   Reply With Quote

Old   August 2, 2010, 17:38
Default
  #12
New Member
 
Join Date: Jun 2010
Posts: 14
Rep Power: 7
maximeg is on a distinguished road
Quote:
Originally Posted by su_junwei View Post
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
maximeg is offline   Reply With Quote

Reply

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
BlockMesh FOAM warning gaottino OpenFOAM Native Meshers: blockMesh 7 July 19, 2010 14:11
On the damBreak4phaseFine cases paean OpenFOAM Running, Solving & CFD 0 November 14, 2008 22:14
How to compute liquid volume of the NOT whole domain snak OpenFOAM Post-Processing 4 June 25, 2008 07:34
fluent add additional zones for the mesh file SSL FLUENT 2 January 26, 2008 12:55
Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Native Meshers: blockMesh 10 April 2, 2007 14:00


All times are GMT -4. The time now is 04:28.