foamToTetDualMesh.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-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 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
26Application
27 foamToTetDualMesh
28
29Group
30 grpPostProcessingUtilities
31
32Description
33 Convert polyMesh results to tetDualMesh.
34
35\*---------------------------------------------------------------------------*/
36
37#include "argList.H"
38#include "fvMesh.H"
39#include "volFields.H"
40#include "pointFields.H"
41#include "Time.H"
42#include "IOobjectList.H"
43
44using namespace Foam;
45
46// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47
48template<class ReadGeoField, class MappedGeoField>
49void ReadAndMapFields
50(
51 const fvMesh& mesh,
52 const IOobjectList& objects,
53 const fvMesh& tetDualMesh,
54 const labelList& map,
55 const typename MappedGeoField::value_type& nullValue,
57)
58{
59 typedef typename MappedGeoField::value_type Type;
60
61 // Search list of objects for wanted type
62 IOobjectList fieldObjects(objects.lookupClass(ReadGeoField::typeName));
63
64 tetFields.setSize(fieldObjects.size());
65
66 label i = 0;
67 forAllConstIters(fieldObjects, iter)
68 {
69 Info<< "Converting " << ReadGeoField::typeName << ' ' << iter.key()
70 << endl;
71
72 ReadGeoField readField(*iter(), mesh);
73
74 tetFields.set
75 (
76 i,
77 new MappedGeoField
78 (
80 (
81 readField.name(),
82 readField.instance(),
83 readField.local(),
84 tetDualMesh,
85 IOobject::NO_READ,
86 IOobject::AUTO_WRITE,
87 readField.registerObject()
88 ),
89 pointMesh::New(tetDualMesh),
90 dimensioned<Type>(readField.dimensions(), Zero)
91 )
92 );
93
94 Field<Type>& fld = tetFields[i].primitiveFieldRef();
95
96 // Map from read field. Set unmapped entries to nullValue.
97 fld.setSize(map.size(), nullValue);
98 forAll(map, pointi)
99 {
100 label index = map[pointi];
101
102 if (index > 0)
103 {
104 label celli = index-1;
105 fld[pointi] = readField[celli];
106 }
107 else if (index < 0)
108 {
109 label facei = -index-1;
110 label bFacei = facei - mesh.nInternalFaces();
111 if (bFacei >= 0)
112 {
113 label patchi = mesh.boundaryMesh().patchID()[bFacei];
114 label localFacei = mesh.boundaryMesh()[patchi].whichFace
115 (
116 facei
117 );
118 fld[pointi] = readField.boundaryField()[patchi][localFacei];
119 }
120 //else
121 //{
122 // FatalErrorInFunction
123 // << "Face " << facei << " from index " << index
124 // << " is not a boundary face." << abort(FatalError);
125 //}
126
127 }
128 //else
129 //{
130 // WarningInFunction
131 // << "Point " << pointi << " at "
132 // << tetDualMesh.points()[pointi]
133 // << " has no dual correspondence." << endl;
134 //}
135 }
136 tetFields[i].correctBoundaryConditions();
137
138 i++;
139 }
140}
141
142
143int main(int argc, char *argv[])
144{
145 argList::addNote
146 (
147 "Convert polyMesh results to tetDualMesh"
148 );
149
150 #include "addOverwriteOption.H"
151 #include "addTimeOptions.H"
152
153 #include "setRootCase.H"
154 #include "createTime.H"
155 // Get times list
156 instantList Times = runTime.times();
157 #include "checkTimeOptions.H"
158 runTime.setTime(Times[startTime], startTime);
159
160 // Read the mesh
161 #include "createNamedMesh.H"
162
163 // Read the tetDualMesh
164 Info<< "Create tetDualMesh for time = "
165 << runTime.timeName() << nl << endl;
166
167 fvMesh tetDualMesh
168 (
170 (
171 "tetDualMesh",
172 runTime.timeName(),
173 runTime,
174 IOobject::MUST_READ
175 )
176 );
177 // From tet vertices to poly cells/faces
178 const labelIOList pointDualAddressing
179 (
181 (
182 "pointDualAddressing",
183 tetDualMesh.facesInstance(),
184 tetDualMesh.meshSubDir,
185 tetDualMesh,
186 IOobject::MUST_READ
187 )
188 );
189
190 if (pointDualAddressing.size() != tetDualMesh.nPoints())
191 {
193 << "Size " << pointDualAddressing.size()
194 << " of addressing map " << pointDualAddressing.objectPath()
195 << " differs from number of points in mesh "
196 << tetDualMesh.nPoints()
197 << exit(FatalError);
198 }
199
200
201 // Some stats on addressing
202 label nCells = 0;
203 label nPatchFaces = 0;
204 label nUnmapped = 0;
205 forAll(pointDualAddressing, pointi)
206 {
207 label index = pointDualAddressing[pointi];
208
209 if (index > 0)
210 {
211 nCells++;
212 }
213 else if (index == 0)
214 {
215 nUnmapped++;
216 }
217 else
218 {
219 label facei = -index-1;
220 if (facei < mesh.nInternalFaces())
221 {
223 << "Face " << facei << " from index " << index
224 << " is not a boundary face."
225 << " nInternalFaces:" << mesh.nInternalFaces()
226 << exit(FatalError);
227 }
228 else
229 {
230 nPatchFaces++;
231 }
232 }
233 }
234
235 reduce(nCells, sumOp<label>());
236 reduce(nPatchFaces, sumOp<label>());
237 reduce(nUnmapped, sumOp<label>());
238 Info<< "tetDualMesh points : " << tetDualMesh.nPoints()
239 << " of which mapped to" << nl
240 << " cells : " << nCells << nl
241 << " patch faces : " << nPatchFaces << nl
242 << " not mapped : " << nUnmapped << nl
243 << endl;
244
245
246 // Read objects in time directory
247 IOobjectList objects(mesh, runTime.timeName());
248
249 // Read vol fields, interpolate onto tet points
251 ReadAndMapFields<volScalarField, pointScalarField>
252 (
253 mesh,
254 objects,
255 tetDualMesh,
256 pointDualAddressing,
257 Zero, // nullValue
258 psFlds
259 );
260
262 ReadAndMapFields<volVectorField, pointVectorField>
263 (
264 mesh,
265 objects,
266 tetDualMesh,
267 pointDualAddressing,
268 Zero, // nullValue
269 pvFlds
270 );
271
273 ReadAndMapFields<volSphericalTensorField, pointSphericalTensorField>
274 (
275 mesh,
276 objects,
277 tetDualMesh,
278 pointDualAddressing,
279 Zero, // nullValue
280 pstFlds
281 );
282
284 ReadAndMapFields<volSymmTensorField, pointSymmTensorField>
285 (
286 mesh,
287 objects,
288 tetDualMesh,
289 pointDualAddressing,
290 Zero, // nullValue
291 psymmtFlds
292 );
293
295 ReadAndMapFields<volTensorField, pointTensorField>
296 (
297 mesh,
298 objects,
299 tetDualMesh,
300 pointDualAddressing,
301 Zero, // nullValue
302 ptFlds
303 );
304
305 tetDualMesh.objectRegistry::write();
306
307 Info<< "End\n" << endl;
308
309 return 0;
310}
311
312
313// ************************************************************************* //
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
reduce(hasMovingMesh, orOp< bool >())
Foam::label startTime
Generic templated field type.
Definition: Field.H:82
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
IOobjectList lookupClass(const char *clsName) const
The list of IOobjects with the given headerClassName.
Definition: IOobjectList.C:313
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
void setSize(const label newLen)
Same as resize()
Definition: PtrList.H:151
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Generic dimensioned Type class.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:324
label nPoints() const noexcept
Number of mesh points.
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
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
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278