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 |
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 |
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 |
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! |
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 |
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 |
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 |
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. |
Quote:
Code:
const cellList& cells = mesh.cells(); 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. |
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 |
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 |
Quote:
Everything works perfectly! Thanks again for your help. Maxime |
All times are GMT -4. The time now is 10:54. |