[vinz]

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]

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]

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]

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!

[vinz]

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]

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]

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]

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]

Originally Posted by rob3rt

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:

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]

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]

Originally Posted by maximeg

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]

Originally Posted by su_junwei

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