BUG: UListIO: byteSize overflowing on really big faceLists
[OpenFOAM-2.0.x.git] / applications / utilities / postProcessing / dataConversion / foamToEnsight / ensightField.C
blobbf33d165d09585160b486e15ed808b1d3d023012
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
11     OpenFOAM is free software: you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by
13     the Free Software Foundation, either version 3 of the License, or
14     (at your option) any later version.
16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19     for more details.
21     You should have received a copy of the GNU General Public License
22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 #include "ensightField.H"
27 #include "fvMesh.H"
28 #include "volFields.H"
29 #include "OFstream.H"
30 #include "IOmanip.H"
31 #include "itoa.H"
32 #include "volPointInterpolation.H"
33 #include "ensightBinaryStream.H"
34 #include "ensightAsciiStream.H"
35 #include "globalIndex.H"
37 using namespace Foam;
39 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
41 template<class Type>
42 Field<Type> map
44     const Field<Type>& vf,
45     const labelList& map1,
46     const labelList& map2
49     Field<Type> mf(map1.size() + map2.size());
51     forAll(map1, i)
52     {
53         mf[i] = vf[map1[i]];
54     }
56     label offset = map1.size();
58     forAll(map2, i)
59     {
60         mf[i + offset] = vf[map2[i]];
61     }
63     return mf;
67 template<class Type>
68 void writeField
70     const char* key,
71     const Field<Type>& vf,
72     ensightStream& ensightFile
75     if (returnReduce(vf.size(), sumOp<label>()) > 0)
76     {
77         if (Pstream::master())
78         {
79             ensightFile.write(key);
81             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
82             {
83                 ensightFile.write(vf.component(cmpt));
85                 for (int slave=1; slave<Pstream::nProcs(); slave++)
86                 {
87                     IPstream fromSlave(Pstream::scheduled, slave);
88                     scalarField slaveData(fromSlave);
89                     ensightFile.write(slaveData);
90                 }
91             }
92         }
93         else
94         {
95             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
96             {
97                 OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
98                 toMaster<< vf.component(cmpt);
99             }
100         }
101     }
105 template<class Type>
106 bool writePatchField
108     const Field<Type>& pf,
109     const label patchi,
110     const label ensightPatchI,
111     const faceSets& boundaryFaceSet,
112     const ensightMesh::nFacePrimitives& nfp,
113     ensightStream& ensightFile
116     if (nfp.nTris || nfp.nQuads || nfp.nPolys)
117     {
118         if (Pstream::master())
119         {
120             ensightFile.writePartHeader(ensightPatchI);
121         }
123         writeField
124         (
125             "tria3",
126             Field<Type>(pf, boundaryFaceSet.tris),
127             ensightFile
128         );
130         writeField
131         (
132             "quad4",
133             Field<Type>(pf, boundaryFaceSet.quads),
134             ensightFile
135         );
137         writeField
138         (
139             "nsided",
140             Field<Type>(pf, boundaryFaceSet.polys),
141             ensightFile
142         );
144         return true;
145     }
146     else
147     {
148         return false;
149     }
153 template<class Type>
154 void writePatchField
156     const word& fieldName,
157     const Field<Type>& pf,
158     const word& patchName,
159     const ensightMesh& eMesh,
160     const fileName& postProcPath,
161     const word& prepend,
162     const label timeIndex,
163     const bool binary,
164     Ostream& ensightCaseFile
167     const Time& runTime = eMesh.mesh().time();
169     const List<faceSets>& boundaryFaceSets = eMesh.boundaryFaceSets();
170     const wordList& allPatchNames = eMesh.allPatchNames();
171     const HashTable<ensightMesh::nFacePrimitives>&
172         nPatchPrims = eMesh.nPatchPrims();
174     label ensightPatchI = eMesh.patchPartOffset();
176     label patchi = -1;
178     forAll(allPatchNames, i)
179     {
180         if (allPatchNames[i] == patchName)
181         {
182             patchi = i;
183             break;
184         }
185         ensightPatchI++;
186     }
189     word pfName = patchName + '.' + fieldName;
191     word timeFile = prepend + itoa(timeIndex);
193     ensightStream* ensightFilePtr = NULL;
194     if (Pstream::master())
195     {
196         if (timeIndex == 0)
197         {
198             ensightCaseFile.setf(ios_base::left);
200             ensightCaseFile
201                 << pTraits<Type>::typeName
202                 << " per element:            1       "
203                 << setw(15) << pfName
204                 << (' ' + prepend + "***." + pfName).c_str()
205                 << nl;
206         }
208         // set the filename of the ensight file
209         fileName ensightFileName(timeFile + "." + pfName);
211         if (binary)
212         {
213             ensightFilePtr = new ensightBinaryStream
214             (
215                 postProcPath/ensightFileName,
216                 runTime
217             );
218         }
219         else
220         {
221             ensightFilePtr = new ensightAsciiStream
222             (
223                 postProcPath/ensightFileName,
224                 runTime
225             );
226         }
227     }
229     ensightStream& ensightFile = *ensightFilePtr;
231     if (Pstream::master())
232     {
233         ensightFile.write(pTraits<Type>::typeName);
234     }
236     if (patchi >= 0)
237     {
238         writePatchField
239         (
240             pf,
241             patchi,
242             ensightPatchI,
243             boundaryFaceSets[patchi],
244             nPatchPrims.find(patchName)(),
245             ensightFile
246         );
247     }
248     else
249     {
250         faceSets nullFaceSets;
252         writePatchField
253         (
254             Field<Type>(),
255             -1,
256             ensightPatchI,
257             nullFaceSets,
258             nPatchPrims.find(patchName)(),
259             ensightFile
260         );
261     }
263     if (Pstream::master())
264     {
265         delete ensightFilePtr;
266     }
270 template<class Type>
271 void ensightField
273     const GeometricField<Type, fvPatchField, volMesh>& vf,
274     const ensightMesh& eMesh,
275     const fileName& postProcPath,
276     const word& prepend,
277     const label timeIndex,
278     const bool binary,
279     Ostream& ensightCaseFile
282     Info<< "Converting field " << vf.name() << endl;
284     word timeFile = prepend + itoa(timeIndex);
286     const fvMesh& mesh = eMesh.mesh();
287     const Time& runTime = mesh.time();
289     const cellSets& meshCellSets = eMesh.meshCellSets();
290     const List<faceSets>& boundaryFaceSets = eMesh.boundaryFaceSets();
291     const wordList& allPatchNames = eMesh.allPatchNames();
292     const wordHashSet& patchNames = eMesh.patchNames();
293     const HashTable<ensightMesh::nFacePrimitives>&
294         nPatchPrims = eMesh.nPatchPrims();
295     const List<faceSets>& faceZoneFaceSets = eMesh.faceZoneFaceSets();
296     const wordHashSet& faceZoneNames = eMesh.faceZoneNames();
297     const HashTable<ensightMesh::nFacePrimitives>&
298         nFaceZonePrims = eMesh.nFaceZonePrims();
300     const labelList& tets = meshCellSets.tets;
301     const labelList& pyrs = meshCellSets.pyrs;
302     const labelList& prisms = meshCellSets.prisms;
303     const labelList& wedges = meshCellSets.wedges;
304     const labelList& hexes = meshCellSets.hexes;
305     const labelList& polys = meshCellSets.polys;
307     ensightStream* ensightFilePtr = NULL;
308     if (Pstream::master())
309     {
310         // set the filename of the ensight file
311         fileName ensightFileName(timeFile + "." + vf.name());
313         if (binary)
314         {
315             ensightFilePtr = new ensightBinaryStream
316             (
317                 postProcPath/ensightFileName,
318                 runTime
319             );
320         }
321         else
322         {
323             ensightFilePtr = new ensightAsciiStream
324             (
325                 postProcPath/ensightFileName,
326                 runTime
327             );
328         }
329     }
331     ensightStream& ensightFile = *ensightFilePtr;
333     if (patchNames.empty())
334     {
335         eMesh.barrier();
337         if (Pstream::master())
338         {
339             if (timeIndex == 0)
340             {
341                 ensightCaseFile.setf(ios_base::left);
343                 ensightCaseFile
344                     << pTraits<Type>::typeName
345                     << " per element:            1       "
346                     << setw(15) << vf.name()
347                     << (' ' + prepend + "***." + vf.name()).c_str()
348                     << nl;
349             }
351             ensightFile.write(pTraits<Type>::typeName);
352             ensightFile.writePartHeader(1);
353         }
355         writeField
356         (
357             "hexa8",
358             map(vf, hexes, wedges),
359             ensightFile
360         );
362         writeField
363         (
364             "penta6",
365             Field<Type>(vf, prisms),
366             ensightFile
367         );
369         writeField
370         (
371             "pyramid5",
372             Field<Type>(vf, pyrs),
373             ensightFile
374         );
376         writeField
377         (
378             "tetra4",
379             Field<Type>(vf, tets),
380             ensightFile
381         );
383         writeField
384         (
385             "nfaced",
386             Field<Type>(vf, polys),
387             ensightFile
388         );
389     }
391     label ensightPatchI = eMesh.patchPartOffset();
393     forAll(allPatchNames, patchi)
394     {
395         const word& patchName = allPatchNames[patchi];
397         eMesh.barrier();
399         if (patchNames.empty() || patchNames.found(patchName))
400         {
401             if
402             (
403                 writePatchField
404                 (
405                     vf.boundaryField()[patchi],
406                     patchi,
407                     ensightPatchI,
408                     boundaryFaceSets[patchi],
409                     nPatchPrims.find(patchName)(),
410                     ensightFile
411                 )
412             )
413             {
414                 ensightPatchI++;
415             }
416         }
417     }
419     // write faceZones, if requested
420     if (faceZoneNames.size())
421     {
422         // Interpolates cell values to faces - needed only when exporting
423         // faceZones...
424         GeometricField<Type, fvsPatchField, surfaceMesh> sf
425         (
426             linearInterpolate(vf)
427         );
429         forAllConstIter(wordHashSet, faceZoneNames, iter)
430         {
431             const word& faceZoneName = iter.key();
433             eMesh.barrier();
435             label zoneID = mesh.faceZones().findZoneID(faceZoneName);
437             const faceZone& fz = mesh.faceZones()[zoneID];
439             // Prepare data to write
440             label nIncluded = 0;
441             forAll(fz, i)
442             {
443                 if (eMesh.faceToBeIncluded(fz[i]))
444                 {
445                     ++nIncluded;
446                 }
447             }
449             Field<Type> values(nIncluded);
451             // Loop on the faceZone and store the needed field values
452             label j = 0;
453             forAll(fz, i)
454             {
455                 label faceI = fz[i];
456                 if (mesh.isInternalFace(faceI))
457                 {
458                     values[j] = sf[faceI];
459                     ++j;
460                 }
461                 else
462                 {
463                     if (eMesh.faceToBeIncluded(faceI))
464                     {
465                         label patchI = mesh.boundaryMesh().whichPatch(faceI);
466                         const polyPatch& pp = mesh.boundaryMesh()[patchI];
467                         label patchFaceI = pp.whichFace(faceI);
468                         Type value = sf.boundaryField()[patchI][patchFaceI];
469                         values[j] = value;
470                         ++j;
471                     }
472                 }
473             }
475             if
476             (
477                 writePatchField
478                 (
479                     values,
480                     zoneID,
481                     ensightPatchI,
482                     faceZoneFaceSets[zoneID],
483                     nFaceZonePrims.find(faceZoneName)(),
484                     ensightFile
485                 )
486             )
487             {
488                 ensightPatchI++;
489             }
490         }
491     }
492     if (Pstream::master())
493     {
494         delete ensightFilePtr;
495     }
499 template<class Type>
500 void ensightPointField
502     const GeometricField<Type, pointPatchField, pointMesh>& pf,
503     const ensightMesh& eMesh,
504     const fileName& postProcPath,
505     const word& prepend,
506     const label timeIndex,
507     const bool binary,
508     Ostream& ensightCaseFile
511     Info<< "Converting field " << pf.name() << endl;
513     word timeFile = prepend + itoa(timeIndex);
515     const fvMesh& mesh = eMesh.mesh();
516     const wordList& allPatchNames = eMesh.allPatchNames();
517     const wordHashSet& patchNames = eMesh.patchNames();
518     const wordHashSet& faceZoneNames = eMesh.faceZoneNames();
521     ensightStream* ensightFilePtr = NULL;
522     if (Pstream::master())
523     {
524         // set the filename of the ensight file
525         fileName ensightFileName(timeFile + "." + pf.name());
527         if (binary)
528         {
529             ensightFilePtr = new ensightBinaryStream
530             (
531                 postProcPath/ensightFileName,
532                 mesh.time()
533             );
534         }
535         else
536         {
537             ensightFilePtr = new ensightAsciiStream
538             (
539                 postProcPath/ensightFileName,
540                 mesh.time()
541             );
542         }
543     }
545     ensightStream& ensightFile = *ensightFilePtr;
547     if (eMesh.patchNames().empty())
548     {
549         eMesh.barrier();
551         if (Pstream::master())
552         {
553             if (timeIndex == 0)
554             {
555                 ensightCaseFile.setf(ios_base::left);
557                 ensightCaseFile
558                     << pTraits<Type>::typeName
559                     << " per node:            1       "
560                     << setw(15) << pf.name()
561                     << (' ' + prepend + "***." + pf.name()).c_str()
562                     << nl;
563             }
565             ensightFile.write(pTraits<Type>::typeName);
566             ensightFile.writePartHeader(1);
567         }
569         writeField
570         (
571             "coordinates",
572             Field<Type>(pf.internalField(), eMesh.uniquePointMap()),
573             ensightFile
574         );
575     }
578     label ensightPatchI = eMesh.patchPartOffset();
580     forAll(allPatchNames, patchi)
581     {
582         const word& patchName = allPatchNames[patchi];
584         eMesh.barrier();
586         if (patchNames.empty() || patchNames.found(patchName))
587         {
588             const fvPatch& p = mesh.boundary()[patchi];
589             if
590             (
591                 returnReduce(p.size(), sumOp<label>())
592               > 0
593             )
594             {
595                 // Renumber the patch points/faces into unique points
596                 labelList pointToGlobal;
597                 labelList uniqueMeshPointLabels;
598                 autoPtr<globalIndex> globalPointsPtr =
599                 mesh.globalData().mergePoints
600                 (
601                     p.patch().meshPoints(),
602                     p.patch().meshPointMap(),
603                     pointToGlobal,
604                     uniqueMeshPointLabels
605                 );
607                 if (Pstream::master())
608                 {
609                     ensightFile.writePartHeader(ensightPatchI);
610                 }
612                 writeField
613                 (
614                     "coordinates",
615                     Field<Type>(pf.internalField(), uniqueMeshPointLabels),
616                     ensightFile
617                 );
619                 ensightPatchI++;
620             }
621         }
622     }
624     // write faceZones, if requested
625     if (faceZoneNames.size())
626     {
627         forAllConstIter(wordHashSet, faceZoneNames, iter)
628         {
629             const word& faceZoneName = iter.key();
631             eMesh.barrier();
633             label zoneID = mesh.faceZones().findZoneID(faceZoneName);
635             const faceZone& fz = mesh.faceZones()[zoneID];
637             if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
638             {
639                 // Renumber the faceZone points/faces into unique points
640                 labelList pointToGlobal;
641                 labelList uniqueMeshPointLabels;
642                 autoPtr<globalIndex> globalPointsPtr =
643                 mesh.globalData().mergePoints
644                 (
645                     fz().meshPoints(),
646                     fz().meshPointMap(),
647                     pointToGlobal,
648                     uniqueMeshPointLabels
649                 );
651                 if (Pstream::master())
652                 {
653                     ensightFile.writePartHeader(ensightPatchI);
654                 }
656                 writeField
657                 (
658                     "coordinates",
659                     Field<Type>
660                     (
661                         pf.internalField(),
662                         uniqueMeshPointLabels
663                     ),
664                     ensightFile
665                 );
667                 ensightPatchI++;
668             }
669         }
670     }
672     if (Pstream::master())
673     {
674         delete ensightFilePtr;
675     }
679 template<class Type>
680 void ensightField
682     const IOobject& fieldObject,
683     const ensightMesh& eMesh,
684     const fileName& postProcPath,
685     const word& prepend,
686     const label timeIndex,
687     const bool binary,
688     const bool nodeValues,
689     Ostream& ensightCaseFile
692     // Read field
693     GeometricField<Type, fvPatchField, volMesh> vf(fieldObject, eMesh.mesh());
695     if (nodeValues)
696     {
697         tmp<GeometricField<Type, pointPatchField, pointMesh> > pfld
698         (
699             volPointInterpolation::New(eMesh.mesh()).interpolate(vf)
700         );
701         pfld().rename(vf.name());
703         ensightPointField<Type>
704         (
705             pfld,
706             eMesh,
707             postProcPath,
708             prepend,
709             timeIndex,
710             binary,
711             ensightCaseFile
712         );
713     }
714     else
715     {
716         ensightField<Type>
717         (
718             vf,
719             eMesh,
720             postProcPath,
721             prepend,
722             timeIndex,
723             binary,
724             ensightCaseFile
725         );
726     }
730 // ************************************************************************* //