
[Sponsors] 
October 30, 2009, 03:35 
How to compute cellZone volume

#1 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 15 
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 

October 30, 2009, 03:53 

#2 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Copenhagen, Denmark
Posts: 1,884
Rep Power: 34 
Hi Vincent
I have never worked with cellZones, however, I do expect that addressing these will give you a list of cellnumbers 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 

October 30, 2009, 04:30 

#3 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 15 
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 

October 30, 2009, 05:01 

#4 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 15 
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//OpenFOAM1.6/lib/linux64GccDPOpt/libOpenFOAM.so" #1 Foam::sigSegv::sigSegvHandler(int) in "/CFDBB/vincent/OpenFOAM//OpenFOAM1.6/lib/linux64GccDPOpt/libOpenFOAM.so" #2 __restore_rt at sigaction.c:0 #3 Foam::cellZone::zoneMesh() const in "/CFDBB/vincent/OpenFOAM//OpenFOAM1.6/lib/linux64GccDPOpt/libOpenFOAM.so" #4 main in "/CFDBB/vincent/OpenFOAM/OpenFOAM1.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/OpenFOAM1.6/applications/bin/linux64GccDPOpt/adSimpleFoam" Segmentation fault Any idea on what is going on? Thanks! Last edited by vinz; October 30, 2009 at 05:23. 

October 30, 2009, 10:03 

#5 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 15 
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 

November 1, 2009, 12:21 

#6 
Senior Member
Francesco Del Citto
Join Date: Mar 2009
Location: Zürich Area, Switzerland
Posts: 236
Rep Power: 15 
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 

November 2, 2009, 01:26 

#7 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 15 
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 02:49. 

July 13, 2010, 10:50 

#8 
New Member
Robert
Join Date: Mar 2010
Posts: 16
Rep Power: 13 
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. 

July 15, 2010, 02:32 

#9  
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,294
Rep Power: 31 
Quote:
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. 

July 30, 2010, 16:31 

#10 
New Member
Join Date: Jun 2010
Posts: 14
Rep Power: 13 
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 

July 31, 2010, 02:46 

#11  
Senior Member

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

August 2, 2010, 17:38 

#12  
New Member
Join Date: Jun 2010
Posts: 14
Rep Power: 13 
Quote:
Everything works perfectly! Thanks again for your help. Maxime 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
[blockMesh] BlockMesh FOAM warning  gaottino  OpenFOAM Meshing & Mesh Conversion  7  July 19, 2010 14:11 
On the damBreak4phaseFine cases  paean  OpenFOAM Running, Solving & CFD  0  November 14, 2008 21:14 
How to compute liquid volume of the NOT whole domain  snak  OpenFOAM PostProcessing  4  June 25, 2008 07:34 
fluent add additional zones for the mesh file  SSL  FLUENT  2  January 26, 2008 11:55 
[blockMesh] Axisymmetrical mesh  Rasmus Gjesing (Gjesing)  OpenFOAM Meshing & Mesh Conversion  10  April 2, 2007 14:00 