boundary mesh and interpolation questions
 

Dear all,
The boundary mesh handling of the fvMesh and GeometricField classes does not work quite as I expect. I believe the implementation makes sense, but I'll need some support to loop over boundaries faces and get field values on those face centers or vertices.

For example, if I run a modified solidDisplacementFoam on the plateHole tutorial (OF 1.5), I observe that the D.boundaryField() has 140 faces, while mesh.boundary() has 2140 faces. (D is the volVectorField related to deformation).
It looks like D.boundaryField() does not include the 2x1000 faces of patch "frontAndBack" which is defined with "type empty".

1) What is the natural procedure to loop over the values of D on *all* boundary faces (including the "frontAndBack" patch)? Right now I think I need to first interpolate D onto the face centers, then for each face check whether it is touching only one (->boundary face) or two cells (->internal face). I'm sure, OpenFOAM has functions for this purpose. Can you point me to the right functions in the class documentation? I didn't find enough information in the user guide or in the programming guide.

2) How do I get the face label of a Foam::face?
For example, setting
Foam::face fc=D.boundaryField()[patchIndex][faceIndex];
gives me a face fc, which I would like find in mesh.faces(). Without getting the label of fc, I can only compare the geometric coordinates, which is very expensive.

3) Does a function exist which returns the values of a vol(Scalar/Vector)Field on the vertices of the cells? Or do I first have to interpolate onto the face centers and from the face centers onto the vertices?

Thank you
Paul

Hi, Paul

I think, you need to read description of boundaries types in OpenFOAM carefully. Empty - is a PHYSICAL boundary type for excluding dimension (perpendicular to boundary faces) from calculations. For example, if you set all planes (oXoY) to 'empty', then, direction Z will be excluded from calculations - see fied D in tutorial plateHole - 3rd component of displacement vector is about 1e-17. So, it is no need to output values at this boundaries - they are like 'slip' type - if field is scalar, acts as zeroGradient, if field is vector - normal components are zero, tangential are zeroGradient.

Moreover, if field value is uniform across the boundary, then values are not printed, but used spicifier 'uniform' and 'value' to show this. See fixedValue, fixedGardient, zeroGradient (grad(psi) = 0)

to iterate through all cells (values at cell centers of internal field volScalar(or Vector)Field psi)
forAll(psi, celli)
{
psi[celli] = ...
}

to iterate through all faces (values at all face centers at all patches)
forAll(psi.boundaryField(), patchi)
{

forAll(psi.boundaryField()[patchi], facei)
{
psi.boundaryField()[patchi][facei] = ....
}
}

or

forAll(psi.boundaryField(), patchi)
{
fvPatchScalarField &ppsi = psi.boundaryField()[patchi]
forAll(ppsi, facei)
{
ppsi[facei] = ....
}
}

I think, in previous example 'facei' is face label.

For managing point scalar (vector, tensor, etc) fields, use
pointScalar(Vector, Tensor, etc)Field, defined as

typedef GeometricField< scalar, pointPatchField, pointMesh > pointScalarField

defined at
line 45 of file pointFields.H


To volume field vField to point file pField (both of them define on fvMesh mesh), use (for example)

volPointInterpolation::New(mesh).interpolate
(
vField,
pField
);



As you see, the key function here is volPointInterpolation::interpolate(...) member of class volPointInterpolation

hope, this helps you

Dear Matjev,
thank you for the explanations. But some of your suggestions are not helpful, as I had already tried them, and they don't provide a solution to my questions.

In my opinion there is need to output values on these boundaries. A surface plot in paraview shows not only values on the fixedValue, zeroGradient or whatever patches, but also on the empty patches, which is very useful. The displayed value is the value of the adjacent cell, I assume. In a similar fashion I want to evaluate D on all boundaries.

No, in this example facei is the index into the patch psi.boundaryField()[patchi]. mesh.faces()[facei] is a different face than psi.boundaryField()[patchi][facei].

Thanks, this is what I needed.

Paul