CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Pre-Processing (http://www.cfd-online.com/Forums/openfoam-pre-processing/)
-   -   assign values to selected cells (http://www.cfd-online.com/Forums/openfoam-pre-processing/103799-assign-values-selected-cells.html)

zxj160 June 27, 2012 05:29

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?

niklas June 27, 2012 06:19

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

what you're after is probably setFields and the setFieldsDict

zxj160 June 27, 2012 06:40

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.

niklas June 27, 2012 06:53

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]

zxj160 June 29, 2012 12:17

Quote:

Originally Posted by niklas (Post 368527)
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]

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?


All times are GMT -4. The time now is 15:39.