mapFields.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2018-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 mapFields
29
30Group
31 grpPreProcessingUtilities
32
33Description
34 Maps volume fields from one mesh to another, reading and
35 interpolating all fields present in the time directory of both cases.
36
37 Parallel and non-parallel cases are handled without the need to reconstruct
38 them first.
39
40\*---------------------------------------------------------------------------*/
41
42#include "fvCFD.H"
43#include "meshToMesh0.H"
44#include "processorFvPatch.H"
45#include "MapMeshes.H"
46#include "decompositionModel.H"
47
48// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
49
50int readNumProcs
51(
52 const argList& args,
53 const word& optionName,
54 const Time& runTime
55)
56{
57 return decompositionMethod::nDomains
58 (
59 IOdictionary
60 (
61 IOobject::selectIO
62 (
63 IOobject
64 (
65 decompositionModel::canonicalName,
66 runTime.system(),
67 runTime,
68 IOobject::MUST_READ,
69 IOobject::NO_WRITE,
70 false // do not register
71 ),
72 args.getOrDefault<fileName>(optionName, "")
73 )
74 )
75 );
76}
77
78
79void mapConsistentMesh
80(
81 const fvMesh& meshSource,
82 const fvMesh& meshTarget,
83 const meshToMesh0::order& mapOrder,
84 const bool subtract
85)
86{
87 if (subtract)
88 {
89 MapConsistentMesh<minusEqOp>
90 (
91 meshSource,
92 meshTarget,
93 mapOrder
94 );
95 }
96 else
97 {
98 MapConsistentMesh<eqOp>
99 (
100 meshSource,
101 meshTarget,
102 mapOrder
103 );
104 }
105}
106
107
108void mapSubMesh
109(
110 const fvMesh& meshSource,
111 const fvMesh& meshTarget,
112 const HashTable<word>& patchMap,
113 const wordList& cuttingPatches,
114 const meshToMesh0::order& mapOrder,
115 const bool subtract
116)
117{
118 if (subtract)
119 {
120 MapSubMesh<minusEqOp>
121 (
122 meshSource,
123 meshTarget,
124 patchMap,
125 cuttingPatches,
126 mapOrder
127 );
128 }
129 else
130 {
131 MapSubMesh<eqOp>
132 (
133 meshSource,
134 meshTarget,
135 patchMap,
136 cuttingPatches,
137 mapOrder
138 );
139 }
140}
141
142
143void mapConsistentSubMesh
144(
145 const fvMesh& meshSource,
146 const fvMesh& meshTarget,
147 const meshToMesh0::order& mapOrder,
148 const bool subtract
149)
150{
151 if (subtract)
152 {
153 MapConsistentSubMesh<minusEqOp>
154 (
155 meshSource,
156 meshTarget,
157 mapOrder
158 );
159 }
160 else
161 {
162 MapConsistentSubMesh<eqOp>
163 (
164 meshSource,
165 meshTarget,
166 mapOrder
167 );
168 }
169}
170
171
172wordList addProcessorPatches
173(
174 const fvMesh& meshTarget,
175 const wordList& cuttingPatches
176)
177{
178 // Add the processor patches to the cutting list
179 HashTable<label> cuttingPatchTable;
180 forAll(cuttingPatches, i)
181 {
182 cuttingPatchTable.insert(cuttingPatches[i], i);
183 }
184
185 forAll(meshTarget.boundary(), patchi)
186 {
187 if (isA<processorFvPatch>(meshTarget.boundary()[patchi]))
188 {
189 if
190 (
191 !cuttingPatchTable.found
192 (
193 meshTarget.boundaryMesh()[patchi].name()
194 )
195 )
196 {
197 cuttingPatchTable.insert
198 (
199 meshTarget.boundaryMesh()[patchi].name(),
200 -1
201 );
202 }
203 }
204 }
205
206 return cuttingPatchTable.toc();
207}
208
209
210// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
211
212int main(int argc, char *argv[])
213{
214 argList::addNote
215 (
216 "Map volume fields from one mesh to another"
217 );
218 argList::noParallel();
219 argList::addArgument("sourceCase");
220
221 argList::addOption
222 (
223 "sourceTime",
224 "scalar|'latestTime'",
225 "Specify the source time"
226 );
227 argList::addOption
228 (
229 "sourceRegion",
230 "word",
231 "Specify the source region"
232 );
233 argList::addOption
234 (
235 "targetRegion",
236 "word",
237 "Specify the target region"
238 );
239 argList::addBoolOption
240 (
241 "parallelSource",
242 "The source is decomposed"
243 );
244 argList::addBoolOption
245 (
246 "parallelTarget",
247 "The target is decomposed"
248 );
249 argList::addBoolOption
250 (
251 "consistent",
252 "Source and target geometry and boundary conditions identical"
253 );
254 argList::addOption
255 (
256 "mapMethod",
257 "word",
258 "Specify the mapping method"
259 );
260 argList::addBoolOption
261 (
262 "subtract",
263 "Subtract mapped source from target"
264 );
265 argList::addOption
266 (
267 "sourceDecomposeParDict",
268 "file",
269 "Read decomposePar dictionary from specified location"
270 );
271 argList::addOption
272 (
273 "targetDecomposeParDict",
274 "file",
275 "Read decomposePar dictionary from specified location"
276 );
277
278
279 argList args(argc, argv);
280 if (!args.check())
281 {
283 }
284 #include "foamDlOpenLibs.H"
285
286 fileName rootDirTarget(args.rootPath());
287 fileName caseDirTarget(args.globalCaseName());
288
289 const auto casePath = args.get<fileName>(1);
290 const fileName rootDirSource = casePath.path().toAbsolute();
291 const fileName caseDirSource = casePath.name();
292
293 Info<< "Source: " << rootDirSource << ' ' << caseDirSource;
294 word sourceRegion(polyMesh::defaultRegion);
295 if (args.readIfPresent("sourceRegion", sourceRegion))
296 {
297 Info<< " (region " << sourceRegion << ')';
298 }
299 Info<< endl;
300
301 Info<< "Target: " << rootDirTarget << ' ' << caseDirTarget;
302 word targetRegion(polyMesh::defaultRegion);
303 if (args.readIfPresent("targetRegion", targetRegion))
304 {
305 Info<< " (region " << targetRegion << ')';
306 }
307 Info<< endl;
308
309 const bool parallelSource = args.found("parallelSource");
310 const bool parallelTarget = args.found("parallelTarget");
311 const bool consistent = args.found("consistent");
312
313 meshToMesh0::order mapOrder = meshToMesh0::INTERPOLATE;
314 if (args.found("mapMethod"))
315 {
316 const word mapMethod(args["mapMethod"]);
317 if (mapMethod == "mapNearest")
318 {
319 mapOrder = meshToMesh0::MAP;
320 }
321 else if (mapMethod == "interpolate")
322 {
323 mapOrder = meshToMesh0::INTERPOLATE;
324 }
325 else if (mapMethod == "cellPointInterpolate")
326 {
327 mapOrder = meshToMesh0::CELL_POINT_INTERPOLATE;
328 }
329 else
330 {
332 << "Unknown mapMethod " << mapMethod << ". Valid options are: "
333 << "mapNearest, interpolate and cellPointInterpolate"
334 << exit(FatalError);
335 }
336
337 Info<< "Mapping method: " << mapMethod << endl;
338 }
339
340 const bool subtract = args.found("subtract");
341 if (subtract)
342 {
343 Info<< "Subtracting mapped source field from target" << endl;
344 }
345
346
347 #include "createTimes.H"
348
349 HashTable<word> patchMap;
350 wordList cuttingPatches;
351
352 if (!consistent)
353 {
354 IOdictionary mapFieldsDict
355 (
356 IOobject
357 (
358 "mapFieldsDict",
359 runTimeTarget.system(),
361 IOobject::MUST_READ_IF_MODIFIED,
362 IOobject::NO_WRITE,
363 false
364 )
365 );
366
367 mapFieldsDict.readEntry("patchMap", patchMap);
368 mapFieldsDict.readEntry("cuttingPatches", cuttingPatches);
369 }
370
371 if (parallelSource && !parallelTarget)
372 {
373 const int nProcs = readNumProcs
374 (
375 args,
376 "sourceDecomposeParDict",
378 );
379
380 Info<< "Create target mesh\n" << endl;
381
382 fvMesh meshTarget
383 (
384 IOobject
385 (
386 targetRegion,
387 runTimeTarget.timeName(),
389 )
390 );
391
392 Info<< "Target mesh size: " << meshTarget.nCells() << endl;
393
394 for (int proci=0; proci<nProcs; proci++)
395 {
396 Info<< nl << "Source processor " << proci << endl;
397
398 Time runTimeSource
399 (
400 Time::controlDictName,
401 rootDirSource,
402 caseDirSource/("processor" + Foam::name(proci))
403 );
404
405 #include "setTimeIndex.H"
406
407 fvMesh meshSource
408 (
409 IOobject
410 (
411 sourceRegion,
412 runTimeSource.timeName(),
414 )
415 );
416
417 Info<< "mesh size: " << meshSource.nCells() << endl;
418
419 if (consistent)
420 {
421 mapConsistentSubMesh
422 (
423 meshSource,
424 meshTarget,
425 mapOrder,
426 subtract
427 );
428 }
429 else
430 {
431 mapSubMesh
432 (
433 meshSource,
434 meshTarget,
435 patchMap,
436 cuttingPatches,
437 mapOrder,
438 subtract
439 );
440 }
441 }
442 }
443 else if (!parallelSource && parallelTarget)
444 {
445 const int nProcs = readNumProcs
446 (
447 args,
448 "targetDecomposeParDict",
450 );
451
452
453 Info<< "Create source mesh\n" << endl;
454
455 #include "setTimeIndex.H"
456
457 fvMesh meshSource
458 (
459 IOobject
460 (
461 sourceRegion,
462 runTimeSource.timeName(),
464 )
465 );
466
467 Info<< "Source mesh size: " << meshSource.nCells() << endl;
468
469 for (int proci=0; proci<nProcs; proci++)
470 {
471 Info<< nl << "Target processor " << proci << endl;
472
473 Time runTimeTarget
474 (
475 Time::controlDictName,
476 rootDirTarget,
477 caseDirTarget/("processor" + Foam::name(proci))
478 );
479
480 fvMesh meshTarget
481 (
482 IOobject
483 (
484 targetRegion,
485 runTimeTarget.timeName(),
487 )
488 );
489
490 Info<< "mesh size: " << meshTarget.nCells() << endl;
491
492 if (consistent)
493 {
494 mapConsistentSubMesh
495 (
496 meshSource,
497 meshTarget,
498 mapOrder,
499 subtract
500 );
501 }
502 else
503 {
504 mapSubMesh
505 (
506 meshSource,
507 meshTarget,
508 patchMap,
509 addProcessorPatches(meshTarget, cuttingPatches),
510 mapOrder,
511 subtract
512 );
513 }
514 }
515 }
516 else if (parallelSource && parallelTarget)
517 {
518 const int nProcsSource = readNumProcs
519 (
520 args,
521 "sourceDecomposeParDict",
523 );
524 const int nProcsTarget = readNumProcs
525 (
526 args,
527 "targetDecomposeParDict",
529 );
530
531 List<boundBox> bbsTarget(nProcsTarget);
532 List<bool> bbsTargetSet(nProcsTarget, false);
533
534 for (int procISource=0; procISource<nProcsSource; procISource++)
535 {
536 Info<< nl << "Source processor " << procISource << endl;
537
538 Time runTimeSource
539 (
540 Time::controlDictName,
541 rootDirSource,
542 caseDirSource/("processor" + Foam::name(procISource))
543 );
544
545 #include "setTimeIndex.H"
546
547 fvMesh meshSource
548 (
549 IOobject
550 (
551 sourceRegion,
552 runTimeSource.timeName(),
554 )
555 );
556
557 Info<< "mesh size: " << meshSource.nCells() << endl;
558
559 boundBox bbSource(meshSource.bounds());
560
561 for (int procITarget=0; procITarget<nProcsTarget; procITarget++)
562 {
563 if
564 (
565 !bbsTargetSet[procITarget]
566 || (
567 bbsTargetSet[procITarget]
568 && bbsTarget[procITarget].overlaps(bbSource)
569 )
570 )
571 {
572 Info<< nl << "Target processor " << procITarget << endl;
573
574 Time runTimeTarget
575 (
576 Time::controlDictName,
577 rootDirTarget,
578 caseDirTarget/("processor" + Foam::name(procITarget))
579 );
580
581 fvMesh meshTarget
582 (
583 IOobject
584 (
585 targetRegion,
586 runTimeTarget.timeName(),
588 )
589 );
590
591 Info<< "mesh size: " << meshTarget.nCells() << endl;
592
593 bbsTarget[procITarget] = meshTarget.bounds();
594 bbsTargetSet[procITarget] = true;
595
596 if (bbsTarget[procITarget].overlaps(bbSource))
597 {
598 if (consistent)
599 {
600 mapConsistentSubMesh
601 (
602 meshSource,
603 meshTarget,
604 mapOrder,
605 subtract
606 );
607 }
608 else
609 {
610 mapSubMesh
611 (
612 meshSource,
613 meshTarget,
614 patchMap,
615 addProcessorPatches(meshTarget, cuttingPatches),
616 mapOrder,
617 subtract
618 );
619 }
620 }
621 }
622 }
623 }
624 }
625 else
626 {
627 #include "setTimeIndex.H"
628
629 Info<< "Create meshes\n" << endl;
630
631 fvMesh meshSource
632 (
633 IOobject
634 (
635 sourceRegion,
636 runTimeSource.timeName(),
638 )
639 );
640
641 fvMesh meshTarget
642 (
643 IOobject
644 (
645 targetRegion,
646 runTimeTarget.timeName(),
648 )
649 );
650
651 Info<< "Source mesh size: " << meshSource.nCells() << tab
652 << "Target mesh size: " << meshTarget.nCells() << nl << endl;
653
654 if (consistent)
655 {
656 mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
657 }
658 else
659 {
660 mapSubMesh
661 (
662 meshSource,
663 meshTarget,
664 patchMap,
665 cuttingPatches,
666 mapOrder,
667 subtract
668 );
669 }
670 }
671
672 Info<< "\nEnd\n" << endl;
673
674 return 0;
675}
676
677
678// ************************************************************************* //
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
const fileName & rootPath() const noexcept
Return root path.
Definition: argListI.H:63
bool check(bool checkArgs=argList::argsMandatory(), bool checkOpts=true) const
Definition: argList.C:1913
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
const fileName & globalCaseName() const noexcept
Return global case name.
Definition: argListI.H:75
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
void exit(const int errNo=1)
Exit : can be called for any error to exit program.
Definition: error.C:331
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
List< word > wordList
A List of words.
Definition: fileName.H:63
void subtract(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr char tab
The tab '\t' character(0x09)
Definition: Ostream.H:52
Time runTimeSource(Time::controlDictName, argsSrc)
Time runTimeTarget(Time::controlDictName, args)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333