CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Field value inside a boundary condition class (http://www.cfd-online.com/Forums/openfoam-programming-development/74980-field-value-inside-boundary-condition-class.html)

pcaron April 13, 2010 15:12

Field value inside a boundary condition class
 
Hello Forum,

I'm using OF-1.5-dev in openSUSE 11.2. I want to implement a boundary condition class which depends on a field value at a point.

I used the parabolicVelocity class as start point. I modified it to fix the boundary value according to different parameters and it works fine.

The next step is to read a field value at a point and use it to modify the boundary value.

I read several post, but I can't understand how to know a field value at a point. These are the posts I read

http://www.cfd-online.com/Forums/ope...ary-point.html
http://openfoamwiki.net/index.php/Sn...value_at_point

Can somebody help me? Thanks in advance

Pablo

su_junwei April 13, 2010 22:58

Hi Pablo
try the following code segment

point p(x,y,z) ;
label cellNo=mesh_.findCell(p);
//using the cell center value or make an interpolation
type value=Variable[cellNO];
If you want to do an interpolation, please see the following thread
http://www.cfd-online.com/Forums/ope...tml#post254038

regards,Junwei

pcaron April 14, 2010 09:52

Hi Junwei

thanks for your reply!

I read your idea before. The problem is that I want to do this evaluation inside a BC class, so I don't have any reference to mesh_ and Variable.

So, I think the improved question could be: How should I reference a Variable (like p or U) and mesh_ inside a BC class?

BTW: These are my first steps in OF programming.

Regards

Pablo

CedricVH April 14, 2010 10:04

You have the access member function const objectRegistry& db() const. Accessing for example U is done by the code:

Code:

const volVectorField& U = db().lookupObject<volVectorField>("U");
For accessing the mesh: when creating a boundary condition, you have access to the patch. You can thus access the mesh with:

Code:

const fvMesh& mesh = patch().boundaryMesh().mesh();

mchurchf April 14, 2010 10:48

Boundary Condition Using Information from Opposite Face
 
Hello,

I have a related question. I am using OpenFOAM-1.6, and I want to set up a boundary condition on velocity that requires knowledge of the velocity at the face opposite the boundary face. Is there a straightforward way to do this from within a derived fvPatchFields class? Or will I have to access the entire velocity volume field from within that derived fvPatchFields class and use something like "oppositeCellFace"? If I need to do the latter, can someone give some guidance, please?

Thank you,

Matt

pcaron April 14, 2010 16:23

Hi Cedric and Junwei

I want to thank you. I could implement my new boundary condition and is working now! :D

BTW, Cedric, As I told you these are my first steps. I found very difficult to understand and/or find functions and classes in OF Oxygen help. Is there another place to find a more user friendly help?

Thanks for your help!

Pablo

Quote:

Originally Posted by CedricVH (Post 254591)
You have the access member function const objectRegistry& db() const. Accessing for example U is done by the code:

Code:

const volVectorField& U = db().lookupObject<volVectorField>("U");
For accessing the mesh: when creating a boundary condition, you have access to the patch. You can thus access the mesh with:

Code:

const fvMesh& mesh = patch().boundaryMesh().mesh();


su_junwei April 14, 2010 22:48

Quote:

Originally Posted by mchurchf (Post 254599)
Hello,

I have a related question. I am using OpenFOAM-1.6, and I want to set up a boundary condition on velocity that requires knowledge of the velocity at the face opposite the boundary face. Is there a straightforward way to do this from within a derived fvPatchFields class?
Thank you,
Matt

Hi Matt

Is your "opposite boundary face" a boundary patch in your case? If it is, You can search the patch by its index or patch name.

label patchId=patch().boundaryMesh().lookupPatchID(patch Name);
and using the following code

const vectorField =U.boundaryField()[patchId];
.... // do your work

U can be obtained using credricVH's method mentioned above.(lookup using ObjectRegistry)

regards, Junwei

mchurchf April 16, 2010 09:38

Accessing internal field and mesh from within boundary condition
 
Thank you all for your help,

I was able to implement my boundary condition as well. I used the following code to access the internal velocity field and the mesh two layers of cells inward from the boundary:

Code:

//  Set up access to the internal velocity field and mesh
    const volVectorField& U =
        db().objectRegistry::lookupObject<volVectorField>(UName_);
    const fvMesh& mesh = patch().boundaryMesh().mesh();

    forAll(patch(), facei)
    {
    // get global cell indices for cells adjacent to patch   
        label celli = patch().faceCells()[facei];

    // get global face indices for faces opposite patch face
    label oppFacei = mesh.cells()[celli].opposingFaceLabel(facei+patch().patch().start(),mesh.faces());

    // get coordinates of center of cell adjacent to patch (patch cells)
    vector cellCentreO = mesh.cellCentres()[mesh.owner()[oppFacei]];

    // get coordinates of center of cell on the side opposite the patch of
    // the patch cell
    vector cellCentreN = mesh.cellCentres()[mesh.neighbour()[oppFacei]];

    // get coordinates of center of face opposite the patch boundary face
    vector faceCentre  = mesh.faceCentres()[oppFacei];

    // get coordinates of center of patch boundary face;
    vector patchFaceCentre = mesh.faceCentres()[facei+patch().patch().start()];
    }


jorkolino July 27, 2010 06:04

How can I access patch() function? I try to include it in my solver files with a compile error "patch is not defined in this scope". My wish is to inject a particle in each boundary cell of a patch. I have the code for particle generation, I just can not find the cell index for which to generate it.
I can not find the function lookupPatchID in the doxygen documentation, either. Though it may still work, I can not test it without first fixing the patch() error message.

JasonWang3 October 1, 2015 12:00

Quote:

Originally Posted by mchurchf (Post 254946)
Thank you all for your help,

I was able to implement my boundary condition as well. I used the following code to access the internal velocity field and the mesh two layers of cells inward from the boundary:

Code:

//  Set up access to the internal velocity field and mesh
    const volVectorField& U =
        db().objectRegistry::lookupObject<volVectorField>(UName_);
    const fvMesh& mesh = patch().boundaryMesh().mesh();

    forAll(patch(), facei)
    {
    // get global cell indices for cells adjacent to patch   
        label celli = patch().faceCells()[facei];

    // get global face indices for faces opposite patch face
    label oppFacei = mesh.cells()[celli].opposingFaceLabel(facei+patch().patch().start(),mesh.faces());

    // get coordinates of center of cell adjacent to patch (patch cells)
    vector cellCentreO = mesh.cellCentres()[mesh.owner()[oppFacei]];

    // get coordinates of center of cell on the side opposite the patch of
    // the patch cell
    vector cellCentreN = mesh.cellCentres()[mesh.neighbour()[oppFacei]];

    // get coordinates of center of face opposite the patch boundary face
    vector faceCentre  = mesh.faceCentres()[oppFacei];

    // get coordinates of center of patch boundary face;
    vector patchFaceCentre = mesh.faceCentres()[facei+patch().patch().start()];
    }


Hi

I follow your code, and write a code to test if it could work or not:
Quote:

#include "fvCFD.H"
#include "wallFvPatch.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "setRootCase.H"
#include "createTime.H"

Info<< "Create mesh, no clear-out\n" << endl;
fvMesh mesh
(
IOobject
(
fvMesh::defaultRegion,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
);

const fvPatchList& patches = mesh.boundary();

forAll(mesh.boundary(), patchi)
{
const fvPatch& patch = patches[patchi];
if (isA<wallFvPatch>(patch)) // is wall or not//
{
const fvMesh& mesh = patch.boundaryMesh().mesh();

forAll(patch, facei)
{
// get global cell indices for cells adjacent to patch
label celli = patch.faceCells()[facei];

// get global face indices for faces opposite patch face
label oppFacei = mesh.cells()[celli].opposingFaceLabel(facei+patch.patch.start(),mesh. faces());

Info<< facei << endl;
Info<< celli << endl;
Info<< oppFacei << endl;
Info<< ' ' << endl;
}

}

}

}
But I can not compile this code, some errors have happens in the opposiongFaceLabel. Do you have any suggestions about the code?

JasonWang3 October 1, 2015 12:44

I figure it out, just delete one patch from opposingFaceLabel(facei+patch.patch.start(),mesh. faces());

Sorry to trouble you.

Jung hoo September 11, 2016 07:14

Hi, I'm trying to create my own boundary condition.

I'm modifying the Member Function from parabolicVelocityFvPatchVectorField.C(i.e, I chose parabolicVelocityFvPatchVectorField.C as the base)

When I tried to read field values like U or p on the boundary member function, compiling was successfully done, but I get the error message as belows when I run "setFields" on an example case based on interFoam.

Quote:

--> FOAM Warning :
From function dlLibraryTableOpen(const fileName& functionLibName)
in file db/dlLibraryTable/dlLibraryTable.C at line 85
could not load /home/lee/foam/lee-3.1/lib/linux64GccDPOpt/libmyfunctionfinal2.so: undefined symbol: _ZTVN4Foam15twoPhaseMixtureE
And then, when I run application(interFoam), an error message is shown :

Quote:

Reading field T

--> FOAM FATAL ERROR:
request for volVectorField U from objectRegistry region0 failed
available objects of type volVectorField are
1
(
T
)

From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/lee/foam/foam-extend-3.1/src/foam/lnInclude/objectRegistryTemplates.C at line 139.
FOAM aborting
Aborted (core dumped)
I used these commans for reading the field values(U, phi, alpha1)

Quote:

const volVectorField& U = db().lookupObject<volVectorField>("U");
const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi");
const volScalarField& alpha1 = db().lookupObject<volScalarField>("alpha1");
I'm beginner in programming, so I guess I'm missing a simple thing...

this is a part of Member Function:

Quote:

void cbVelocityFvPatchVectorField::updateCoeffs()
{
if (updated())
{
return;
}
const fvMesh& mesh = dimensionedInternalField().mesh();
//const fvMesh& mesh = patch().boundaryMesh().mesh();
const word& patchName = this->patch().name();
const label patchID = mesh.boundaryMesh().findPatchID(patchName);
const volVectorField& U = db().lookupObject<volVectorField>("U");
const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi");
const volScalarField& alpha1 = db().lookupObject<volScalarField>(alphaName());
//volVectorField U("U", U);
//surfaceScalarField phi("phi", phi);
//volScalarField alpha1("alpha1", alpha1);
twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);
volSymmTensorField Reff(turbulence->devReff());

scalar rho = 1000;
scalar nu = 1e-06;
scalar mu = nu*rho;
scalar g = 9.81;
scalar s = 1.4;
scalar d = 0.005;
scalar shieldsc = 0.05;
scalar kappa = 0.8;
scalar c0 = 3;
//scalar Alpha1 = 0.8*mathematicalConstantPi;
scalar Alpha1 = 2;
scalarField wallShearStress = rho*mag((-mesh.Sf().boundaryField()[patchID]/mesh.magSf().boundaryField()[patchID]) & Reff.boundaryField()[patchID]);
scalarField shields = (wallShearStress)/(rho*(s-1)*g*d);
scalarField PEF = (pow((1+pow4((((1/6)*mathematicalConstantPi*mu)/(shields-shieldsc)))), -(1/4)));
scalarField linearC = sqrt(((sqr(kappa)*sqr(Alpha1))*(shields-shieldsc-((1/6)*mathematicalConstantPi*mu*PEF)))/(0.013*s*shields));
scalarField cb = c0/pow3((1+(1/linearC)));
vectorField n = patch().nf();
vectorFieldOperator=(-n*cb);
Any comments will help me a lot!
Thank you in advance!!


All times are GMT -4. The time now is 09:51.