setFields.C 12.7 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
9
10
11
12
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

13
14
15
16
    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.
17
18
19
20
21
22
23

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

26
27
28
29
30
31
Application
    setFields

Group
    grpPreProcessingUtilities

32
Description
33
    Set values on a selected set of cells/patch-faces via a dictionary.
34
35
36
37
38
39
40
41

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

#include "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "topoSetSource.H"
#include "cellSet.H"
42
#include "faceSet.H"
43
44
45
46
#include "volFields.H"

using namespace Foam;

Andrew Heather's avatar
Andrew Heather committed
47
template<class Type>
48
bool setCellFieldType
49
(
Andrew Heather's avatar
Andrew Heather committed
50
    const word& fieldTypeDesc,
51
52
53
54
55
    const fvMesh& mesh,
    const labelList& selectedCells,
    Istream& fieldValueStream
)
{
Andrew Heather's avatar
Andrew Heather committed
56
57
58
59
60
61
62
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

    if (fieldTypeDesc != fieldType::typeName + "Value")
    {
        return false;
    }

63
64
    word fieldName(fieldValueStream);

65
    // Check the current time directory
66
67
68
69
70
71
72
73
    IOobject fieldHeader
    (
        fieldName,
        mesh.time().timeName(),
        mesh,
        IOobject::MUST_READ
    );

74
    // Check the "constant" directory
75
    if (!fieldHeader.typeHeaderOk<fieldType>(true))
76
77
78
79
80
81
82
83
84
85
    {
        fieldHeader = IOobject
        (
            fieldName,
            mesh.time().constant(),
            mesh,
            IOobject::MUST_READ
        );
    }

86
    // Check field exists
87
    if (fieldHeader.typeHeaderOk<fieldType>(true))
88
    {
89
90
        Info<< "    Setting internal values of "
            << fieldHeader.headerClassName()
91
92
            << " " << fieldName << endl;

93
        fieldType field(fieldHeader, mesh, false);
94

Andrew Heather's avatar
Andrew Heather committed
95
        const Type& value = pTraits<Type>(fieldValueStream);
96
97
98

        if (selectedCells.size() == field.size())
        {
99
            field.primitiveFieldRef() = value;
100
101
102
103
104
105
106
107
108
        }
        else
        {
            forAll(selectedCells, celli)
            {
                field[selectedCells[celli]] = value;
            }
        }

109
        typename GeometricField<Type, fvPatchField, volMesh>::
110
            Boundary& fieldBf = field.boundaryFieldRef();
111

112
113
        forAll(field.boundaryField(), patchi)
        {
114
            fieldBf[patchi] = fieldBf[patchi].patchInternalField();
115
116
        }

117
118
        if (!field.write())
        {
119
120
            FatalErrorInFunction
              << "Failed writing field " << fieldName << endl;
121
        }
122
123
124
    }
    else
    {
125
126
        WarningInFunction
          << "Field " << fieldName << " not found" << endl;
127
128
129

        // Consume value
        (void)pTraits<Type>(fieldValueStream);
130
    }
Andrew Heather's avatar
Andrew Heather committed
131
132

    return true;
133
134
135
}


136
class setCellField
137
138
139
140
{

public:

141
    setCellField()
142
143
    {}

144
    autoPtr<setCellField> clone() const
145
    {
146
        return autoPtr<setCellField>::New();
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
    }

    class iNew
    {
        const fvMesh& mesh_;
        const labelList& selectedCells_;

    public:

        iNew(const fvMesh& mesh, const labelList& selectedCells)
        :
            mesh_(mesh),
            selectedCells_(selectedCells)
        {}

162
163
164
165
166
167
        iNew(const fvMesh& mesh, labelList&& selectedCells)
        :
            mesh_(mesh),
            selectedCells_(std::move(selectedCells))
        {}

168
        autoPtr<setCellField> operator()(Istream& fieldValues) const
169
170
171
        {
            word fieldType(fieldValues);

Andrew Heather's avatar
Andrew Heather committed
172
173
174
            if
            (
               !(
175
                    setCellFieldType<scalar>
Andrew Heather's avatar
Andrew Heather committed
176
                        (fieldType, mesh_, selectedCells_, fieldValues)
177
                 || setCellFieldType<vector>
Andrew Heather's avatar
Andrew Heather committed
178
                        (fieldType, mesh_, selectedCells_, fieldValues)
179
                 || setCellFieldType<sphericalTensor>
Andrew Heather's avatar
Andrew Heather committed
180
                        (fieldType, mesh_, selectedCells_, fieldValues)
181
                 || setCellFieldType<symmTensor>
Andrew Heather's avatar
Andrew Heather committed
182
                        (fieldType, mesh_, selectedCells_, fieldValues)
183
                 || setCellFieldType<tensor>
Andrew Heather's avatar
Andrew Heather committed
184
185
186
                        (fieldType, mesh_, selectedCells_, fieldValues)
                )
            )
187
            {
188
                WarningInFunction
189
190
191
192
                    << "field type " << fieldType << " not currently supported"
                    << endl;
            }

193
            return autoPtr<setCellField>::New();
194
195
196
197
198
        }
    };
};


199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
template<class Type>
bool setFaceFieldType
(
    const word& fieldTypeDesc,
    const fvMesh& mesh,
    const labelList& selectedFaces,
    Istream& fieldValueStream
)
{
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

    if (fieldTypeDesc != fieldType::typeName + "Value")
    {
        return false;
    }

    word fieldName(fieldValueStream);

217
    // Check the current time directory
218
219
220
221
222
223
224
225
    IOobject fieldHeader
    (
        fieldName,
        mesh.time().timeName(),
        mesh,
        IOobject::MUST_READ
    );

226
    // Check the "constant" directory
227
    if (!fieldHeader.typeHeaderOk<fieldType>(true))
228
229
230
231
232
233
234
235
236
237
    {
        fieldHeader = IOobject
        (
            fieldName,
            mesh.time().constant(),
            mesh,
            IOobject::MUST_READ
        );
    }

238
    // Check field exists
239
    if (fieldHeader.typeHeaderOk<fieldType>(true))
240
241
242
243
244
245
246
247
248
249
    {
        Info<< "    Setting patchField values of "
            << fieldHeader.headerClassName()
            << " " << fieldName << endl;

        fieldType field(fieldHeader, mesh);

        const Type& value = pTraits<Type>(fieldValueStream);

        // Create flat list of selected faces and their value.
250
        Field<Type> allBoundaryValues(mesh.nBoundaryFaces());
251
252
253
254
255
256
257
258
        forAll(field.boundaryField(), patchi)
        {
            SubField<Type>
            (
                allBoundaryValues,
                field.boundaryField()[patchi].size(),
                field.boundaryField()[patchi].patch().start()
              - mesh.nInternalFaces()
259
            ) = field.boundaryField()[patchi];
260
261
262
        }

        // Override
263
        bool hasWarned = false;
264
265
266
267
268
        labelList nChanged
        (
            returnReduce(field.boundaryField().size(), maxOp<label>()),
            0
        );
269
270
271
272
273
        forAll(selectedFaces, i)
        {
            label facei = selectedFaces[i];
            if (mesh.isInternalFace(facei))
            {
274
275
276
                if (!hasWarned)
                {
                    hasWarned = true;
277
                    WarningInFunction
278
279
280
                        << "Ignoring internal face " << facei
                        << ". Suppressing further warnings." << endl;
                }
281
282
283
            }
            else
            {
284
285
286
                label bFacei = facei-mesh.nInternalFaces();
                allBoundaryValues[bFacei] = value;
                nChanged[mesh.boundaryMesh().patchID()[bFacei]]++;
287
288
289
290
291
292
            }
        }

        Pstream::listCombineGather(nChanged, plusEqOp<label>());
        Pstream::listCombineScatter(nChanged);

293
        typename GeometricField<Type, fvPatchField, volMesh>::
294
            Boundary& fieldBf = field.boundaryFieldRef();
295

296
297
298
299
300
301
302
303
        // Reassign.
        forAll(field.boundaryField(), patchi)
        {
            if (nChanged[patchi] > 0)
            {
                Info<< "    On patch "
                    << field.boundaryField()[patchi].patch().name()
                    << " set " << nChanged[patchi] << " values" << endl;
304
                fieldBf[patchi] == SubField<Type>
305
306
                (
                    allBoundaryValues,
307
308
                    fieldBf[patchi].size(),
                    fieldBf[patchi].patch().start()
309
310
311
312
313
                  - mesh.nInternalFaces()
                );
            }
        }

314
315
        if (!field.write())
        {
316
317
            FatalErrorInFunction
                << "Failed writing field " << field.name() << exit(FatalError);
318
        }
319
320
321
    }
    else
    {
322
323
        WarningInFunction
          << "Field " << fieldName << " not found" << endl;
324
325
326

        // Consume value
        (void)pTraits<Type>(fieldValueStream);
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
    }

    return true;
}


class setFaceField
{

public:

    setFaceField()
    {}

    autoPtr<setFaceField> clone() const
    {
343
        return autoPtr<setFaceField>::New();
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
    }

    class iNew
    {
        const fvMesh& mesh_;
        const labelList& selectedFaces_;

    public:

        iNew(const fvMesh& mesh, const labelList& selectedFaces)
        :
            mesh_(mesh),
            selectedFaces_(selectedFaces)
        {}

359
360
361
362
363
364
        iNew(const fvMesh& mesh, labelList&& selectedFaces)
        :
            mesh_(mesh),
            selectedFaces_(std::move(selectedFaces))
        {}

365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
        autoPtr<setFaceField> operator()(Istream& fieldValues) const
        {
            word fieldType(fieldValues);

            if
            (
               !(
                    setFaceFieldType<scalar>
                        (fieldType, mesh_, selectedFaces_, fieldValues)
                 || setFaceFieldType<vector>
                        (fieldType, mesh_, selectedFaces_, fieldValues)
                 || setFaceFieldType<sphericalTensor>
                        (fieldType, mesh_, selectedFaces_, fieldValues)
                 || setFaceFieldType<symmTensor>
                        (fieldType, mesh_, selectedFaces_, fieldValues)
                 || setFaceFieldType<tensor>
                        (fieldType, mesh_, selectedFaces_, fieldValues)
                )
            )
            {
385
                WarningInFunction
386
387
388
389
                    << "field type " << fieldType << " not currently supported"
                    << endl;
            }

390
            return autoPtr<setFaceField>::New();
391
392
393
394
395
396
        }
    };
};



397
398
399
400
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
401
402
403
404
405
    argList::addNote
    (
        "Set values on a selected set of cells/patch-faces via a dictionary"
    );

406
    argList::addOption("dict", "file", "Alternative setFieldsDict");
407

andy's avatar
andy committed
408
409
410
411
    #include "addRegionOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createNamedMesh.H"
412

Henry Weller's avatar
Henry Weller committed
413
414
    const word dictName("setFieldsDict");
    #include "setSystemMeshDictionaryIO.H"
415

416
    Info<< "Reading " << dictIO.name() << nl << endl;
Henry Weller's avatar
Henry Weller committed
417
418

    IOdictionary setFieldsDict(dictIO);
419
420
421
422

    if (setFieldsDict.found("defaultFieldValues"))
    {
        Info<< "Setting field default values" << endl;
423
        PtrList<setCellField> defaultFieldValues
424
425
        (
            setFieldsDict.lookup("defaultFieldValues"),
426
            setCellField::iNew(mesh, labelList(mesh.nCells()))
427
428
429
430
431
432
433
434
435
        );
        Info<< endl;
    }


    Info<< "Setting field region values" << endl;

    PtrList<entry> regions(setFieldsDict.lookup("regions"));

436
    forAll(regions, regionI)
437
    {
438
        const entry& region = regions[regionI];
439

440
        autoPtr<topoSetSource> source =
441
442
            topoSetSource::New(region.keyword(), mesh, region.dict());

443
444
445
446
447
448
449
450
        if (source().setType() == topoSetSource::CELLSETSOURCE)
        {
            cellSet selectedCellSet
            (
                mesh,
                "cellSet",
                mesh.nCells()/10+1  // Reasonable size estimate.
            );
451

452
453
454
455
456
            source->applyToSet
            (
                topoSetSource::NEW,
                selectedCellSet
            );
457

458
459
460
            PtrList<setCellField> fieldValues
            (
                region.dict().lookup("fieldValues"),
461
                setCellField::iNew(mesh, selectedCellSet.sortedToc())
462
463
464
465
466
467
468
469
            );
        }
        else if (source().setType() == topoSetSource::FACESETSOURCE)
        {
            faceSet selectedFaceSet
            (
                mesh,
                "faceSet",
470
                mesh.nBoundaryFaces()/10+1
471
472
473
474
475
476
477
478
479
480
481
            );

            source->applyToSet
            (
                topoSetSource::NEW,
                selectedFaceSet
            );

            PtrList<setFaceField> fieldValues
            (
                region.dict().lookup("fieldValues"),
482
                setFaceField::iNew(mesh, selectedFaceSet.sortedToc())
483
484
            );
        }
485
486
    }

487

Mark Olesen's avatar
Mark Olesen committed
488
    Info<< "\nEnd\n" << endl;
489
490
491
492
493
494

    return 0;
}


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