CFD Online Discussion Forums - assign values to selected cells

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Pre-Processing (https://www.cfd-online.com/Forums/openfoam-pre-processing/)

- - assign values to selected cells (https://www.cfd-online.com/Forums/openfoam-pre-processing/103799-assign-values-selected-cells.html)

assign values to selected cells

Does anyone have the experience to set values to selected cells.

How to check the x y z coordinates of the grid points and store their ID in an array? With the ID, I may then set the values of them by direct assignment?

check the Allrun-script in the tutorial
tutorials/multiphase/interFoam/ras/damBreak

what you're after is probably setFields and the setFieldsDict

Quote:

Originally Posted by niklas (Post 368517)

check the Allrun-script in the tutorial
tutorials/multiphase/interFoam/ras/damBreak

what you're after is probably setFields and the setFieldsDict

Thanks. Is it possible to store a list of cells and defined value for them? I also want to the stored list unchanged (NO_WRITE) when calculated each time step.

sounds like you should use topoSet and setsToZones to create a cell-selection then.

and you can use the code option in functionObjects to do whatever you want with that selection

here's my controlDict for doing some averaging over certain cell-selections (and to calculate total pressure drop) that should get you started

Code:



functions

{



    cuttingAverage

    {

        type               coded;

        functionObjectLibs ( "libutilityFunctionObjects.so" );

        redirectType       average;

        //outputControl      timeStep;

        code

        #{



             vector normal(1.0, 0.0, 0.0);

             const volScalarField& T = mesh().lookupObject<volScalarField>("T");

             const labelList& ids = mesh().cellZones().findIndices("docInlet");

             scalar Tav = 0.0;

             scalar Tdev = 0.0;

             scalar vol = 0.0;

             forAll(ids, i)

             {

                 const labelList& cells = mesh().cellZones()[ids[i]];

                 forAll(cells, j)

                 {

                     label celli = cells[j];

                     scalar V = mesh().V()[celli];

                     vol += V;

                     Tav += T[celli]*V;

                 }

             }

             reduce(vol, sumOp<scalar>());

             reduce(Tav, sumOp<scalar>());

             

             Tav /= vol + VSMALL;

             forAll(ids, i)

             {

                 const labelList& cells = mesh().cellZones()[ids[i]];

                 forAll(cells, j)

                 {

                     label celli = cells[j];

                     scalar V = mesh().V()[celli];

                     Tdev += sqr(T[celli]-Tav)*V;

                 }

             }

             reduce(Tdev, sumOp<scalar>());

             

             Tdev /= vol + VSMALL;

             Tdev = sqrt(Tdev);

             

             Info << "T avg = " << Tav << ", dev = " << Tdev << endl;

             Info << "T min/max = " << min(T).value() << ", " << max(T).value() << endl;

             const volScalarField& p = mesh().lookupObject<volScalarField>("p");

             const volScalarField& rho = mesh().lookupObject<volScalarField>("rho");             

             const volVectorField& U = mesh().lookupObject<volVectorField>("U");

             volScalarField pTot = p + 0.5*rho*(U&U);



             label inletID = mesh().boundaryMesh().findPatchID("STL_inlet");

             label outletID = mesh().boundaryMesh().findPatchID("STL_outlet");



             scalar areaIn = sum(mag(mesh().boundaryMesh()[inletID].faceAreas()));

             reduce(areaIn, sumOp<scalar>());



             scalar areaOut = sum(mag(mesh().boundaryMesh()[outletID].faceAreas()));

             reduce(areaOut, sumOp<scalar>());



             scalar pTotIn = sum(pTot.boundaryField()[inletID]*mag(mesh().boundaryMesh()[inletID].faceAreas()));

             reduce(pTotIn, sumOp<scalar>());



             scalar pTotOut = sum(pTot.boundaryField()[outletID]*mag(mesh().boundaryMesh()[outletID].faceAreas()));

             reduce(pTotOut, sumOp<scalar>());

                  

             Info << "pressure drop: " << pTotIn/areaIn - pTotOut/areaOut 

                  << " in: " << pTotIn/areaIn 

                  << " out: " << pTotOut/areaOut << endl;



             vector Uav = vector::zero;



             forAll(ids, i)

             {

                 const labelList& cells = mesh().cellZones()[ids[i]];

                 forAll(cells, j)

                 {

                     label celli = cells[j];

                     scalar V = mesh().V()[celli];

                     Uav += U[celli]*V;

                 }

             }

             reduce(Uav, sumOp<vector>());

             Uav /= vol + VSMALL;



             scalar Uavn = Uav & normal;

             scalar Udevn = 0.0;



             forAll(ids, i)

             {

                 const labelList& cells = mesh().cellZones()[ids[i]];

                 forAll(cells, j)

                 {

                     label celli = cells[j];

                     scalar V = mesh().V()[celli];

                     scalar Un = U[celli] & normal;

                     Udevn += (Un - Uavn)*(Un - Uavn)*V;

                 }

             }

             reduce(Udevn, sumOp<scalar>());

             Udevn /= vol*Uavn + VSMALL;

             Udevn = sqrt(mag(Udevn));

             Info << "Udevn = " << Udevn << endl;

        #};

    }



}

and this is my topoSetDict

Code:

actions

(

    {

        name            doc;

        action          new ;

        type            cellSet;

        source          cylinderToCell;

        sourceInfo

        {

            p1 (  0.0 0.0 0.553);

            p2 (  0.0 0.0 0.656);

            radius  0.17;

        }

    }



    {

        name            docInlet;

        action          new ;

        type            cellSet;

        source          cylinderToCell;

        sourceInfo

        {

            p1 (  0.0 0.0 0.553);

            p2 (  0.0 0.0 0.557);

            radius  0.17;

        }

    }

);

[/code]

Quote:

Originally Posted by niklas (Post 368527)

Code:



functions

{



    cuttingAverage

    {

        type               coded;

    functionObjectLibs ( "libutilityFunctionObjects.so" );

    redirectType       average;

    //outputControl      timeStep;

        code

    #{



         vector normal(1.0, 0.0, 0.0);

         const volScalarField& T = mesh().lookupObject<volScalarField>("T");

         const labelList& ids = mesh().cellZones().findIndices("docInlet");

         scalar Tav = 0.0;

         scalar Tdev = 0.0;

         scalar vol = 0.0;

         forAll(ids, i)

         {

         const labelList& cells = mesh().cellZones()[ids[i]];

         forAll(cells, j)

         {

             label celli = cells[j];

             scalar V = mesh().V()[celli];

             vol += V;

             Tav += T[celli]*V;

         }

         }

         reduce(vol, sumOp<scalar>());

         reduce(Tav, sumOp<scalar>());

         

         Tav /= vol + VSMALL;

         forAll(ids, i)

         {

         const labelList& cells = mesh().cellZones()[ids[i]];

         forAll(cells, j)

         {

             label celli = cells[j];

             scalar V = mesh().V()[celli];

             Tdev += sqr(T[celli]-Tav)*V;

         }

         }

         reduce(Tdev, sumOp<scalar>());

         

         Tdev /= vol + VSMALL;

         Tdev = sqrt(Tdev);

         

         Info << "T avg = " << Tav << ", dev = " << Tdev << endl;

         Info << "T min/max = " << min(T).value() << ", " << max(T).value() << endl;

             const volScalarField& p = mesh().lookupObject<volScalarField>("p");

             const volScalarField& rho = mesh().lookupObject<volScalarField>("rho");         

         const volVectorField& U = mesh().lookupObject<volVectorField>("U");

         volScalarField pTot = p + 0.5*rho*(U&U);



         label inletID = mesh().boundaryMesh().findPatchID("STL_inlet");

         label outletID = mesh().boundaryMesh().findPatchID("STL_outlet");



             scalar areaIn = sum(mag(mesh().boundaryMesh()[inletID].faceAreas()));

         reduce(areaIn, sumOp<scalar>());



             scalar areaOut = sum(mag(mesh().boundaryMesh()[outletID].faceAreas()));

         reduce(areaOut, sumOp<scalar>());



         scalar pTotIn = sum(pTot.boundaryField()[inletID]*mag(mesh().boundaryMesh()[inletID].faceAreas()));

         reduce(pTotIn, sumOp<scalar>());



         scalar pTotOut = sum(pTot.boundaryField()[outletID]*mag(mesh().boundaryMesh()[outletID].faceAreas()));

         reduce(pTotOut, sumOp<scalar>());

              

         Info << "pressure drop: " << pTotIn/areaIn - pTotOut/areaOut 

          << " in: " << pTotIn/areaIn 

          << " out: " << pTotOut/areaOut << endl;



         vector Uav = vector::zero;



         forAll(ids, i)

         {

         const labelList& cells = mesh().cellZones()[ids[i]];

         forAll(cells, j)

         {

             label celli = cells[j];

             scalar V = mesh().V()[celli];

             Uav += U[celli]*V;

         }

         }

         reduce(Uav, sumOp<vector>());

         Uav /= vol + VSMALL;



         scalar Uavn = Uav & normal;

         scalar Udevn = 0.0;



         forAll(ids, i)

         {

         const labelList& cells = mesh().cellZones()[ids[i]];

         forAll(cells, j)

         {

             label celli = cells[j];

             scalar V = mesh().V()[celli];

             scalar Un = U[celli] & normal;

             Udevn += (Un - Uavn)*(Un - Uavn)*V;

         }

         }

         reduce(Udevn, sumOp<scalar>());

         Udevn /= vol*Uavn + VSMALL;

         Udevn = sqrt(mag(Udevn));

         Info << "Udevn = " << Udevn << endl;

    #};

    }



}

and this is my topoSetDict

Code:

actions

(

    {

        name            doc;

        action          new ;

        type            cellSet;

        source          cylinderToCell;

        sourceInfo

    {

            p1 (  0.0 0.0 0.553);

            p2 (  0.0 0.0 0.656);

            radius  0.17;

    }

    }



    {

        name            docInlet;

        action          new ;

        type            cellSet;

        source          cylinderToCell;

        sourceInfo

    {

            p1 (  0.0 0.0 0.553);

            p2 (  0.0 0.0 0.557);

            radius  0.17;

    }

    }

);

[/code]

Many thanks. Actually, my purpose is to select the first layers near the inlet and outlet and then give values to them. Since the inlet and outlet are cyclic patches and they do not accept other type of BCs. Do you have this kind of experience?