patchInjectionBase.C 8.22 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    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.

    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
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

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

#include "patchInjectionBase.H"
#include "polyMesh.H"
#include "SubField.H"
#include "cachedRandom.H"
#include "triPointRef.H"
31
32
#include "volFields.H"
#include "polyMeshTetDecomposition.H"
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::patchInjectionBase::patchInjectionBase
(
    const polyMesh& mesh,
    const word& patchName
)
:
    patchName_(patchName),
    patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
    patchArea_(0.0),
    patchNormal_(),
    cellOwners_(),
    triFace_(),
    triToFace_(),
    triCumulativeMagSf_(),
    sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
{
    if (patchId_ < 0)
    {
54
55
        FatalErrorInFunction
            << "Requested patch " << patchName_ << " not found" << nl
56
57
            << "Available patches are: " << mesh.boundaryMesh().names() << nl
            << exit(FatalError);
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
    }

    updateMesh(mesh);
}


Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
:
    patchName_(pib.patchName_),
    patchId_(pib.patchId_),
    patchArea_(pib.patchArea_),
    patchNormal_(pib.patchNormal_),
    cellOwners_(pib.cellOwners_),
    triFace_(pib.triFace_),
    triToFace_(pib.triToFace_),
    triCumulativeMagSf_(pib.triCumulativeMagSf_),
    sumTriMagSf_(pib.sumTriMagSf_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::patchInjectionBase::~patchInjectionBase()
{}


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

void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
{
    // Set/cache the injector cells
    const polyPatch& patch = mesh.boundaryMesh()[patchId_];
    const pointField& points = patch.points();

    cellOwners_ = patch.faceCells();

94
    // Triangulate the patch faces and create addressing
95
96
97
98
99
    DynamicList<label> triToFace(2*patch.size());
    DynamicList<scalar> triMagSf(2*patch.size());
    DynamicList<face> triFace(2*patch.size());
    DynamicList<face> tris(5);

100
    // Set zero value at the start of the tri area list
101
102
    triMagSf.append(0.0);

103
    forAll(patch, facei)
104
    {
105
        const face& f = patch[facei];
106
107
108
109
110
111

        tris.clear();
        f.triangles(points, tris);

        forAll(tris, i)
        {
112
            triToFace.append(facei);
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
            triFace.append(tris[i]);
            triMagSf.append(tris[i].mag(points));
        }
    }

    forAll(sumTriMagSf_, i)
    {
        sumTriMagSf_[i] = 0.0;
    }

    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);

    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
    Pstream::listCombineScatter(sumTriMagSf_);

    for (label i = 1; i < triMagSf.size(); i++)
    {
        triMagSf[i] += triMagSf[i-1];
    }

133
    // Transfer to persistent storage
134
135
136
137
    triFace_.transfer(triFace);
    triToFace_.transfer(triToFace);
    triCumulativeMagSf_.transfer(triMagSf);

138
    // Convert sumTriMagSf_ into cumulative sum of areas per proc
139
140
    for (label i = 1; i < sumTriMagSf_.size(); i++)
    {
141
        sumTriMagSf_[i] += sumTriMagSf_[i-1];
142
143
144
145
146
147
148
149
150
151
152
    }

    const scalarField magSf(mag(patch.faceAreas()));
    patchArea_ = sum(magSf);
    patchNormal_ = patch.faceAreas()/magSf;
    reduce(patchArea_, sumOp<scalar>());
}


void Foam::patchInjectionBase::setPositionAndCell
(
153
    const fvMesh& mesh,
154
155
156
    cachedRandom& rnd,
    vector& position,
    label& cellOwner,
157
    label& tetFacei,
Henry Weller's avatar
Henry Weller committed
158
    label& tetPti
159
160
)
{
161
    scalar areaFraction = rnd.globalPosition(scalar(0), patchArea_);
162
163
164
165

    if (cellOwners_.size() > 0)
    {
        // Determine which processor to inject from
166
        label proci = 0;
167
168
        forAllReverse(sumTriMagSf_, i)
        {
169
            if (areaFraction >= sumTriMagSf_[i])
170
            {
171
                proci = i;
172
173
174
175
                break;
            }
        }

176
        if (Pstream::myProcNo() == proci)
177
        {
178
            // Find corresponding decomposed face triangle
179
            label trii = 0;
180
            scalar offset = sumTriMagSf_[proci];
181
182
183
184
            forAllReverse(triCumulativeMagSf_, i)
            {
                if (areaFraction > triCumulativeMagSf_[i] + offset)
                {
185
                    trii = i;
186
187
188
189
                    break;
                }
            }

190
            // Set cellOwner
191
            label facei = triToFace_[trii];
192
            cellOwner = cellOwners_[facei];
193

194
            // Find random point in triangle
195
196
            const polyPatch& patch = mesh.boundaryMesh()[patchId_];
            const pointField& points = patch.points();
197
            const face& tf = triFace_[trii];
198
199
200
            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
            const point pf(tri.randomPoint(rnd));

201
            // Position perturbed away from face (into domain)
Henry's avatar
Henry committed
202
            const scalar a = rnd.position(scalar(0.1), scalar(0.5));
203
            const vector& pc = mesh.cellCentres()[cellOwner];
204
205
            const vector d =
                mag((pf - pc) & patchNormal_[facei])*patchNormal_[facei];
206
207
208

            position = pf - a*d;

209
            // Try to find tetFacei and tetPti in the current position
210
211
            mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);

212
            // tetFacei and tetPti not found, check if the cell has changed
213
214
            if (tetFacei == -1 ||tetPti == -1)
            {
215
216
217
218
219
220
221
222
223
                mesh.findCellFacePt(position, cellOwner, tetFacei, tetPti);
            }

            // Both searches failed, choose a random position within
            // the original cell
            if (tetFacei == -1 ||tetPti == -1)
            {
                // Reset cellOwner
                cellOwner = cellOwners_[facei];
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
                const scalarField& V = mesh.V();

                // Construct cell tet indices
                const List<tetIndices> cellTetIs =
                    polyMeshTetDecomposition::cellTetIndices(mesh, cellOwner);

                // Construct cell tet volume fractions
                scalarList cTetVFrac(cellTetIs.size(), 0.0);
                for (label teti=1; teti<cellTetIs.size()-1; teti++)
                {
                    cTetVFrac[teti] =
                        cTetVFrac[teti-1]
                      + cellTetIs[teti].tet(mesh).mag()/V[cellOwner];
                }
                cTetVFrac.last() = 1;

                // Set new particle position
                const scalar volFrac = rnd.sample01<scalar>();
                label teti = 0;
                forAll(cTetVFrac, vfI)
                {
                    if (cTetVFrac[vfI] > volFrac)
                    {
                        teti = vfI;
                        break;
                    }
                }
                position = cellTetIs[teti].tet(mesh).randomPoint(rnd);
                tetFacei = cellTetIs[teti].face();
                tetPti = cellTetIs[teti].tetPt();
            }
255
256
257
258
        }
        else
        {
            cellOwner = -1;
259
            tetFacei = -1;
Henry Weller's avatar
Henry Weller committed
260
            tetPti = -1;
261

262
            // Dummy position
263
264
265
266
267
268
            position = pTraits<vector>::max;
        }
    }
    else
    {
        cellOwner = -1;
269
        tetFacei = -1;
Henry Weller's avatar
Henry Weller committed
270
        tetPti = -1;
271

272
        // Dummy position
273
274
275
276
277
278
        position = pTraits<vector>::max;
    }
}


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