OpenFOAM类库介绍（二）Gauss梯度计算

Green-Gauss法梯度计算

计算梯度的核心函数位于/opt/openfoam9/src/finiteVolume/finiteVolume/gradSchemes/gradScheme目录下的gaussGrad.C文件中的gradf函数。根据笔者的理解，添加注释后，将代码列于下方：

 template<class Type>//Type可以是任何类型，如果是scalar那么就是求标量场的梯度，如果是vector那么就是求矢量场梯度（结果是张量场）
 Foam::tmp//函数输出tmp，为了方便C++内存管理，在openFoam中tmp代表任何指定的对象
 <
     Foam::GeometricField//openFoam中最经常使用的场数据，派生自Field，不仅包含场数值，还包含网格数据、单位和边界条件
     <
         typename Foam::outerProduct<Foam::vector, Type>::type,//这里输出场的类型不是指定的，而是推导得到的。使用外积的概念可知，如果对一个标量场求梯度，那么会得到矢量场。如果对矢量场求梯度，那么会得到张量场，以此类推...
         Foam::fvPatchField,
         Foam::volMesh//由fvPatchField和volMesh可知，输出场的数值是存储在网格体心的，而不是面心或者顶点
     >
 >
 Foam::fv::gaussGrad<Type>::gradf
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,//gradf函数的第一个参数：是一个待求梯度的场，注意：由fvsPatchField和surfaceMesh可知，该场的数值是存储在网格面心的，而不是体心或者顶点。
     const word& name//gradf函数的第二个参数：场的名称，这个参数无关紧要
 )
 {
     typedef typename outerProduct<vector, Type>::type GradType;//使用外积推导梯度场的类型，如果对一个标量场求梯度，那么会得到矢量场。如果对矢量场求梯度，那么会得到张量场，以此类推...
 
     const fvMesh& mesh = ssf.mesh();//获取网格
 
     tmp<GeometricField<GradType, fvPatchField, volMesh>> tgGrad//新建一个tmp临时变量，这个变量就是最终返回的梯度场。该场的数值是存储在网格体心的，而不是面心或者顶点
     (
         GeometricField<GradType, fvPatchField, volMesh>::New
         (//GeometricField的构造函数，需要指定名称、网格、量纲、初始值、边界处数值外推方式
             name,
             mesh,
             dimensioned<GradType>
             (
                 "0",
                 ssf.dimensions()/dimLength,
                 Zero
             ),
             extrapolatedCalculatedFvPatchField<GradType>::typeName
         )
     );
     GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad.ref();//将tgGrad(tmp类型）引用赋值给gGrad(GeometricField)。之所以这样做是因为tmp类型仅仅用于储存对象而不能调用对象中的方法（函数），需要把对象的引用拿出来后才能调用其中的方法（函数）
 
     const labelUList& owner = mesh.owner();//owner列表：内部面的所有者单元序号，来自const/polyMesh/owner文件的一部分
     const labelUList& neighbour = mesh.neighbour();//neighbour列表：内部面的邻居单元序号，来自const/polyMesh/neighbour文件
     const vectorField& Sf = mesh.Sf();//一个矢量场（场其实就是列表或数组），存储了内部面的法向量(不包含边界面)，这些法向量的数值是存储在网格面心的，而不是体心或者顶点。法向量的大小是面的面积，方向由owner指向neighbour
 
     Field<GradType>& igGrad = gGrad;//新建一个场，来自gGrad的内部网格（非边界），舍弃了gGrad的边界条件、量纲等信息。该场的数值是存储在网格体心的，而不是面心或者顶点。igGrad中的i表示internal，是内部的含义。
     const Field<Type>& issf = ssf;//新建一个场，来自ssf的内部网格。该场的数值是存储在网格面心的，而不是体心或者顶点
 			//内部面对梯度的贡献
     forAll(owner, facei)//遍历所有内部面(不包括边界面），面的序号是facei
     {
         GradType Sfssf = Sf[facei]*issf[facei];//计算面积与待求梯度场之乘积，Green-Gauss法
 
         igGrad[owner[facei]] += Sfssf;//给第facei个面的owner单元增加Sfssf，因为Sf的方向由owner指向neighbour,朝向owner单元外部
         igGrad[neighbour[facei]] -= Sfssf;//给第facei个面的neighbour单元减少Sfssf，因为Sf的方向由owner指向neighbour,朝向neighbour单元内部
     }
 			//边界面对梯度的贡献
     forAll(mesh.boundary(), patchi)//对每个patch遍历，注：往往边界是由若干个patch组成的
     {
         const labelUList& pFaceCells =
             mesh.boundary()[patchi].faceCells();//一个列表，表示第patchi个patch所对应的单元序号，也就是这个patch上所有面所对应的owner单元序号。
 
         const vectorField& pSf = mesh.Sf().boundaryField()[patchi];//一个向量场，表示第patchi个patch上的面法向量场。很显然，该场的数值是存储在网格面心的，而不是体心或者顶点。
 
         const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];//边界上的ssf场
 
         forAll(mesh.boundary()[patchi], facei)//对第patchi个patch，遍历每一个面
         {
             igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];//给第facei个面的owner单元增加pSf[facei]*pssf[facei]
         }
     }
 
     igGrad /= mesh.V();//igGrad除以单元体积即可得到梯度
 
     gGrad.correctBoundaryConditions();//更新边界条件
 
     return tgGrad;//返回tmp变量
 }

openFoam中计算梯度都需要调用gradf函数，以fvc::grad函数为例，其返回值调用了gradf函数：

 template<class Type>
 tmp
 <
     GeometricField
     <
         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
     >
 >
 grad
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
 )
 {
     return fv::gaussGrad<Type>::gradf(ssf, "grad(" + ssf.name() + ')');
 }

当然还有一系列grad函数的重载，但它们最终都回归到gradf函数