extendedCellToFaceStencilTemplates.C 4.73 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
6
7
8
9
10
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
15
16
17
18
19
20
21

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
22
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
23
24
25

\*---------------------------------------------------------------------------*/

mattijs's avatar
mattijs committed
26
#include "extendedCellToFaceStencil.H"
27
28
29
30

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type>
mattijs's avatar
mattijs committed
31
void Foam::extendedCellToFaceStencil::collectData
32
(
33
34
    const mapDistribute& map,
    const labelListList& stencil,
35
36
    const GeometricField<Type, fvPatchField, volMesh>& fld,
    List<List<Type> >& stencilFld
37
)
38
39
{
    // 1. Construct cell data in compact addressing
40
    List<Type> flatFld(map.constructSize(), pTraits<Type>::zero);
41
42
43
44

    // Insert my internal values
    forAll(fld, cellI)
    {
45
        flatFld[cellI] = fld[cellI];
46
47
48
49
50
51
    }
    // Insert my boundary values
    forAll(fld.boundaryField(), patchI)
    {
        const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];

52
53
54
55
56
        label nCompact =
            pfld.patch().start()
           -fld.mesh().nInternalFaces()
           +fld.mesh().nCells();

57
58
        forAll(pfld, i)
        {
59
            flatFld[nCompact++] = pfld[i];
60
61
62
63
        }
    }

    // Do all swapping
64
    map.distribute(flatFld);
65
66

    // 2. Pull to stencil
67
    stencilFld.setSize(stencil.size());
68

69
    forAll(stencil, faceI)
70
    {
71
        const labelList& compactCells = stencil[faceI];
72
73
74
75
76

        stencilFld[faceI].setSize(compactCells.size());

        forAll(compactCells, i)
        {
77
            stencilFld[faceI][i] = flatFld[compactCells[i]];
78
79
80
81
82
83
84
        }
    }
}


template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
mattijs's avatar
mattijs committed
85
Foam::extendedCellToFaceStencil::weightedSum
86
(
87
88
    const mapDistribute& map,
    const labelListList& stencil,
89
90
    const GeometricField<Type, fvPatchField, volMesh>& fld,
    const List<List<scalar> >& stencilWeights
91
)
92
93
94
95
96
{
    const fvMesh& mesh = fld.mesh();

    // Collect internal and boundary values
    List<List<Type> > stencilFld;
97
    collectData(map, stencil, fld, stencilFld);
98
99
100
101
102
103
104
105
106

    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
    (
        new GeometricField<Type, fvsPatchField, surfaceMesh>
        (
            IOobject
            (
                fld.name(),
                mesh.time().timeName(),
mattijs's avatar
mattijs committed
107
108
109
110
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
111
112
113
114
115
116
117
118
119
120
121
122
            ),
            mesh,
            dimensioned<Type>
            (
                fld.name(),
                fld.dimensions(),
                pTraits<Type>::zero
            )
        )
    );
    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsfCorr();

123
    // Internal faces
124
125
126
127
128
129
130
131
132
133
    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
    {
        const List<Type>& stField = stencilFld[faceI];
        const List<scalar>& stWeight = stencilWeights[faceI];

        forAll(stField, i)
        {
            sf[faceI] += stField[i]*stWeight[i];
        }
    }
henry's avatar
henry committed
134

135
136
    // Boundaries. Either constrained or calculated so assign value
    // directly (instead of nicely using operator==)
henry's avatar
henry committed
137
138
    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
        GeometricBoundaryField& bSfCorr = sf.boundaryField();
139

henry's avatar
henry committed
140
141
142
143
    forAll(bSfCorr, patchi)
    {
        fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];

henry's avatar
henry committed
144
        if (pSfCorr.coupled())
henry's avatar
henry committed
145
        {
146
            label faceI = pSfCorr.patch().start();
henry's avatar
henry committed
147

henry's avatar
henry committed
148
            forAll(pSfCorr, i)
henry's avatar
henry committed
149
            {
henry's avatar
henry committed
150
151
                const List<Type>& stField = stencilFld[faceI];
                const List<scalar>& stWeight = stencilWeights[faceI];
152

henry's avatar
henry committed
153
154
155
156
157
158
159
                forAll(stField, j)
                {
                    pSfCorr[i] += stField[j]*stWeight[j];
                }

                faceI++;
            }
henry's avatar
henry committed
160
161
        }
    }
162
163
164
165
166
167

    return tsfCorr;
}


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