CFD Online Logo CFD Online URL
Home > Forums > OpenFOAM Programming & Development

fvc::div for surfaceVectorField

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By Zeppo

LinkBack Thread Tools Display Modes
Old   March 16, 2011, 05:18
Default fvc::div for surfaceVectorField
Artem Shaklein
Join Date: Feb 2010
Location: Russia, Izhevsk
Posts: 43
Rep Power: 9
ARTem is on a distinguished road
Greetings !

I want to calculate divergence of velocity.
Firstly, I have volVectorField U. After some manipulations with U, I obtain surfaceScalarField Usurface. Next, i try to calculate divergence: fvc::div(Usurface), but this action leads to unexpected result.
For example, when i compute div of mesh.Sf() (dimension [m^2]), i recieve volVectorField (dimension [m^-1]) instead of volScalarField (dimension [m^1]).
But computation grad of mag(mesg.Sf()) (surfaceScalarField, dimension [m^2]) leads to volVectorField with dimension [m^1]. This is right,

There are my questions:
-Why fvc::grad(surfaceScalarField) gives correct result, but fvc::div(surfaceVectorField) - incorrect? May be i make mistake in calculations?
-Can i calculate divergence of surfaceVectorField using fvc::div?

Thank you!
Best regards!

Last edited by ARTem; March 16, 2011 at 11:10.
ARTem is offline   Reply With Quote

Old   January 9, 2017, 08:21
New Member
Join Date: Jun 2016
Posts: 9
Rep Power: 2
toboto is on a distinguished road
I've exactly the same problem. When I interpolate velocity field to the surface I got a surfaceVectorField. When I try to calculate the divergence of this field I got a volVectorField while I expect a volScalarField.

According to the programmers' manual:
The fvc::div function can take as its argument either a surface<Type>Field, in which case
φf is specified directly, or a vol<Type>Field which is interpolated to the face by central differencing
I found in the source code of fvcDiv.H that if the fvc::div is applied to a volField it will return a field of one rank lower
    GeometricField <typename innerProduct<vector, Type>::type, fvPatchField, volMesh> > div
         const GeometricField<Type, fvPatchField, volMesh>&,
         const word& name
While if fvc::div is applied to a scalarField it will return a field of the same rank!!!
tmp<GeometricField<Type, fvPatchField, volMesh>> div
const GeometricField<Type, fvsPatchField, surfaceMesh>&

Any explanation is highly appreciated.
Best regards.

Last edited by toboto; January 9, 2017 at 09:52. Reason: Add additional information
toboto is offline   Reply With Quote

Old   January 14, 2017, 17:47
Senior Member
Join Date: Dec 2009
Posts: 150
Rep Power: 10
Zeppo is on a distinguished road
If you check source code you'll see there're 14 overloaded functions fvc:div() for several combinations of input parameters. This set of functions can be broken down into three groups with one "primary" function within each group being called by other functions in the same group. Those primary functions are declared as:

template<class Type >
tmp< GeometricField< typename innerProduct< vector, Type >::type, fvPatchField, volMesh > > div
    const GeometricField< Type, fvPatchField, volMesh > & vf,
    const word &name
template<class Type >
tmp< GeometricField< Type, fvPatchField, volMesh > > div
    const surfaceScalarField &flux,
    const GeometricField< Type, fvPatchField, volMesh > &vf,
    const word &name
template<class Type >
tmp< GeometricField< Type, fvPatchField, volMesh > > div
    const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf
The first function on the list is the most general out of three:

75 template<class Type>
76 tmp
77 <
78     GeometricField
79     <
80         typename innerProduct<vector, Type>::type, fvPatchField, volMesh
81     >
82 >
83 div
84 (
85     const GeometricField<Type, fvPatchField, volMesh>& vf,
86     const word& name
87 )
88 {
89     return fv::divScheme<Type>::New
90     (
91         vf.mesh(), vf.mesh().divScheme(name)
92     )().fvcDiv(vf);
93 }
It selects (through the Runtime Selection Mechanism) one of the available div schemes with name "name" and calls its member function fvcDiv. Currently the only one available div scheme in OpenFOAM is Gauss scheme (see Divergence theorem). Its fvcDiv is:

42 template<class Type>
43 tmp
44 <
45     GeometricField
46     <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
47 >
48 gaussDivScheme<Type>::fvcDiv
49 (
50     const GeometricField<Type, fvPatchField, volMesh>& vf
51 )
52 {
53     tmp
54     <
55         GeometricField
56         <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
57     > tDiv
58     (
59         fvc::surfaceIntegrate
60         (
61             this->mesh_.Sf() & this->tinterpScheme_().interpolate(vf)
62         )
63     );
65     tDiv().rename("div(" + + ')');
67     return tDiv;
68 }
Firstly, it interpolates your volume field (defined in cell centres) onto a surface field (defined in cell face centres) :

Secondly, it makes a dot (inner) product (&) of the mesh face vectors and the interpolated field. It is that operation which reduces the rank! :

this->mesh_.Sf() & this->tinterpScheme_().interpolate(vf)
This dot product can be thought of as an operation of getting surface flux of the field through the cell faces.

Lastly, it loops over the cells and for each cell it sums the fluxes through the cell's faces. This total sum of fluxes for each cell is then divided by a cell's volume

42 template<class Type>
43 void surfaceIntegrate
44 (
45     Field<Type>& ivf,
46     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
47 )
48 {
49     const fvMesh& mesh = ssf.mesh();
51     const labelUList& owner = mesh.owner();
52     const labelUList& neighbour = mesh.neighbour();
54     const Field<Type>& issf = ssf;
56     forAll(owner, facei)
57     {
58         ivf[owner[facei]] += issf[facei];
59         ivf[neighbour[facei]] -= issf[facei];
60     }
62     forAll(mesh.boundary(), patchi)
63     {
64         const labelUList& pFaceCells =
65             mesh.boundary()[patchi].faceCells();
67         const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
69         forAll(mesh.boundary()[patchi], facei)
70         {
71             ivf[pFaceCells[facei]] += pssf[facei];
72         }
73     }
75     ivf /= mesh.Vsc();
76 }
The type of an object returned from fvcDiv is

45     GeometricField
46     <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
It is a GeometricField with elements of type innerProduct<vector, Type>::type :

88 template<class arg1, class arg2>
89 class innerProduct
90 {
91 public:
93     typedef typename typeOfRank
94     <
95         typename pTraits<arg1>::cmptType,
96         int(pTraits<arg1>::rank) + int(pTraits<arg2>::rank) - 2
97     >::type type;
98 };
innerProduct is a class template. type nested in innerProduct defines the type of an object which inner product (basically, dot operation) of two objects (of classes arg1 and arg2) results in. typeOfRank is an example of a well known design pattern called "Trait". It can be instantiated with a specific value parameter rank:

44 template<class Cmpt, int rank>
45 class typeOfRank
46 {};
You just provide it with a rank and it tells you of what type an object of this rank is. OpenFOAM defines three partial specialisations of class template typeOfRank with template parameters (ranks) 0, 1, 2 for types scalar, vector and tensor, accordingly

41 template<class Cmpt>
42 class typeOfRank<Cmpt, 0>
43 {
44 public:
46     typedef Cmpt type;
47 };
129 template<class Cmpt>
130 class typeOfRank<Cmpt, 1>
131 {
132 public:
134     typedef Vector<Cmpt> type;
135 };
178 template<class Cmpt>
179 class typeOfRank<Cmpt, 2>
180 {
181 public:
183     typedef Tensor<Cmpt> type;
184 };
pTraits is ... well another example of design pattern "Trait". If you specialise it for a specific type it can tell you the rank of any object of that type. Of course, in OpenFOMA it is allready specialised for scalar, vector and tensor.

The other two "primary" fvc:div functions mentioned in the very beginning of this very lengthy post can be dissected in the similar way.

If something is unclear please feel free to ask.
qjh888 likes this.

Last edited by Zeppo; January 15, 2017 at 06:21.
Zeppo is offline   Reply With Quote


Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
fvc::div() strange behaviour ivan_cozza OpenFOAM Running, Solving & CFD 2 February 6, 2010 07:09
help needed with fvc::div() johanna OpenFOAM Programming & Development 2 August 31, 2009 07:12
What type return fvcdiv su_junwei OpenFOAM Running, Solving & CFD 6 October 13, 2008 07:09

All times are GMT -4. The time now is 01:46.