extrude2DMesh.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
26\*---------------------------------------------------------------------------*/
27
28#include "extrude2DMesh.H"
29#include "polyMesh.H"
30#include "polyTopoChange.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
36 defineTypeNameAndDebug(extrude2DMesh, 0);
37}
38
39// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40
41void Foam::extrude2DMesh::check2D() const
42{
43 const faceList& faces = mesh_.faces();
44 forAll(faces, facei)
45 {
46 if (faces[facei].size() != 2)
47 {
49 << "Face " << facei << " size " << faces[facei].size()
50 << " is not of size 2: mesh is not a valid two-dimensional "
51 << "mesh" << exit(FatalError);
52 }
53 }
54}
55
56
57//void Foam::extrude2DMesh::findExtrudeDirection()
58//{
59// scalar minRange = GREAT;
60
61// for (direction dir = 0; dir < 3; dir++)
62// {
63// scalarField cmpts(mesh_.points().component(dir));
64
65// scalar range = max(cmpts)-min(cmpts);
66
67// Info<< "Direction:" << dir << " range:" << range << endl;
68
69// if (range < minRange)
70// {
71// minRange = range;
72// extrudeDir_ = dir;
73// }
74// }
75
76// Info<< "Extruding in direction " << extrudeDir_
77// << " with thickness " << thickness_ << nl
78// << endl;
79//}
80
81
82// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
83
85(
86 polyMesh& mesh,
87 const dictionary& dict,
88 const extrudeModel& model
89)
90:
91 mesh_(mesh),
92 dict_(dict),
93 //patchDict_(dict.subDict("patchInfo")),
94 model_(model),
95 modelType_(dict.get<word>("extrudeModel")),
96 patchType_(dict.get<word>("patchType")),
97 frontPatchi_(-1),
98 backPatchi_(-1)
99{
100 check2D();
101}
102
103
104// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
105
107{
108 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
109
110 frontPatchi_ = patches.findPatchID("front");
111 backPatchi_ = patches.findPatchID("back");
112
113 // Add patch.
114 List<polyPatch*> newPatches(patches.size() + 2);
115
116 forAll(patches, patchi)
117 {
118 const polyPatch& pp = patches[patchi];
119
120 newPatches[patchi] =
121 pp.clone
122 (
123 patches,
124 newPatches.size(),
125 pp.size(),
126 pp.start()
127 ).ptr();
128 }
129
130 if (frontPatchi_ == -1)
131 {
132 frontPatchi_ = patches.size();
133
134 newPatches[frontPatchi_] =
136 (
137 patchType_,
138 "front",
139 0,
140 mesh_.nFaces(),
141 frontPatchi_,
142 patches
143 ).ptr();
144
145 Info<< "Adding patch " << newPatches[frontPatchi_]->name()
146 << " at index " << frontPatchi_
147 << " for front faces." << nl << endl;
148 }
149
150 if (backPatchi_ == -1)
151 {
152 backPatchi_ = patches.size() + 1;
153
154 newPatches[backPatchi_] =
156 (
157 patchType_,
158 "back",
159 0,
160 mesh_.nFaces(),
161 backPatchi_,
162 patches
163 ).ptr();
164
165 Info<< "Adding patch " << newPatches[backPatchi_]->name()
166 << " at index " << backPatchi_
167 << " for back faces." << nl << endl;
168 }
169
170 mesh_.removeBoundary();
171 mesh_.addPatches(newPatches);
172}
173
174
176(
177 polyTopoChange& meshMod
178)
179{
180 const label nLayers = model_.nLayers();
181 const pointField& points = mesh_.points();
182 label nFaces = 0;
183
184 for (label layer = 0; layer < nLayers; ++layer)
185 {
186 label offset = layer * mesh_.nCells();
187
188 forAll(mesh_.cells(), celli)
189 {
190 meshMod.addCell
191 (
192 -1, //masterPointID,
193 -1, //masterEdgeID,
194 -1, //masterFaceID,
195 celli + offset, //masterCellID,
196 mesh_.cellZones().whichZone(celli) //zoneID
197 );
198 }
199 }
200
201
202 // Generate points
203 // ~~~~~~~~~~~~~~~
204
205 for (label layer = 0; layer <= nLayers; ++layer)
206 {
207 label offset = layer * points.size();
208
209 forAll(points, pointi)
210 {
211 // Don't need the surface normal for either linearDirection or
212 // wedge. Will need to add to be able to use others.
213 point newPoint = model_
214 (
215 points[pointi],
216 vector(),
217 layer
218 );
219
220 meshMod.addPoint
221 (
222 newPoint,
223 pointi + offset,
224 -1, // zoneID
225 true // inCell
226 );
227 }
228
229 Pout<< "Added " << points.size() << " points to layer "
230 << layer << endl;
231 }
232
233
234 // Generate faces
235 // ~~~~~~~~~~~~~~
236
237 const faceList& faces = mesh_.faces();
238 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
239
240 for (label layer = 0; layer < nLayers; ++layer)
241 {
242 label currentLayerOffset = layer * mesh_.nPoints();
243 label nextLayerOffset = currentLayerOffset + mesh_.nPoints();
244
245 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
246 {
247 label zoneID = mesh_.faceZones().whichZone(facei);
248 bool zoneFlip = false;
249 if (zoneID != -1)
250 {
251 const faceZone& fZone = mesh_.faceZones()[zoneID];
252 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
253 }
254
255 face newFace(4);
256 const face& f = faces[facei];
257 newFace[0] = f[0] + currentLayerOffset;
258 newFace[1] = f[1] + currentLayerOffset;
259 newFace[2] = f[1] + nextLayerOffset;
260 newFace[3] = f[0] + nextLayerOffset;
261
262//{
263// vector n = newFace.normal(pointField(meshMod.points()));
264// label own = mesh_.faceOwner()[facei];
265// const labelList& ownPoints = mesh_.cellPoints()[own];
266// point ownCc = sum(pointField(mesh_.points(), ownPoints))/ownPoints.size();
267// label nei = mesh_.faceNeighbour()[facei];
268// const labelList& neiPoints = mesh_.cellPoints()[nei];
269// point neiCc = sum(pointField(mesh_.points(), neiPoints))/neiPoints.size();
270// vector d = neiCc - ownCc;
271
272// Pout<< "face:" << facei << " at:" << f.centre(mesh_.points()) << endl
273// << " own:" << own << " at:" << ownCc << endl
274// << " nei:" << nei << " at:" << neiCc << endl
275// << " sign:" << (n & d) << endl
276// << endl;
277//}
278
279 label offset = layer * mesh_.nCells();
280
281 meshMod.addFace
282 (
283 newFace,
284 mesh_.faceOwner()[facei] + offset, // own
285 mesh_.faceNeighbour()[facei] + offset, // nei
286 -1, // masterPointID
287 -1, // masterEdgeID
288 nFaces++, // masterFaceID
289 false, // flipFaceFlux
290 -1, // patchID
291 zoneID, // zoneID
292 zoneFlip // zoneFlip
293 );
294
295 if (debug)
296 {
297 Info<< newFace << " "
298 << mesh_.faceOwner()[facei] + offset << " "
299 << mesh_.faceNeighbour()[facei] + offset << " "
300 << nFaces - 1
301 << endl;
302 }
303 }
304 }
305
306 forAll(patches, patchi)
307 {
308 for (label layer=0; layer < nLayers; layer++)
309 {
310 label currentLayerOffset = layer*mesh_.nPoints();
311 label nextLayerOffset = currentLayerOffset + mesh_.nPoints();
312
313 label startFacei = patches[patchi].start();
314 label endFacei = startFacei + patches[patchi].size();
315
316 for (label facei = startFacei; facei < endFacei; facei++)
317 {
318 label zoneID = mesh_.faceZones().whichZone(facei);
319 bool zoneFlip = false;
320 if (zoneID != -1)
321 {
322 const faceZone& fZone = mesh_.faceZones()[zoneID];
323 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
324 }
325
326 face newFace(4);
327 const face& f = faces[facei];
328 newFace[0] = f[0] + currentLayerOffset;
329 newFace[1] = f[1] + currentLayerOffset;
330 newFace[2] = f[1] + nextLayerOffset;
331 newFace[3] = f[0] + nextLayerOffset;
332
333 label offset = layer * mesh_.nCells();
334
335 meshMod.addFace
336 (
337 newFace,
338 mesh_.faceOwner()[facei] + offset, // own
339 -1, // nei
340 -1, // masterPointID
341 -1, // masterEdgeID
342 nFaces++, // masterFaceID
343 false, // flipFaceFlux
344 patchi, // patchID
345 zoneID, // zoneID
346 zoneFlip // zoneFlip
347 );
348
349 if (debug)
350 {
351 Info<< newFace << " "
352 << mesh_.faceOwner()[facei] + offset << " "
353 << nFaces - 1
354 << endl;
355 }
356 }
357 }
358 }
359
360 // Add extra internal faces that need special treatment for owners and
361 // neighbours.
362 forAll(mesh_.cells(), celli)
363 {
364 const cell& cFaces = mesh_.cells()[celli];
365
366 face frontFace(cFaces.size());
367
368 // Make a loop out of faces.
369 label nextFacei = cFaces[0];
370
371 const face& f = faces[nextFacei];
372
373 label nextPointi;
374 if (mesh_.faceOwner()[nextFacei] == celli)
375 {
376 frontFace[0] = f[0];
377 nextPointi = f[1];
378 }
379 else
380 {
381 frontFace[0] = f[1];
382 nextPointi = f[0];
383 }
384
385
386 for (label i = 1; i < frontFace.size(); i++)
387 {
388 frontFace[i] = nextPointi;
389
390 // Find face containing pointi
391 forAll(cFaces, cFacei)
392 {
393 label facei = cFaces[cFacei];
394 if (facei != nextFacei)
395 {
396 const face& f = faces[facei];
397
398 if (f[0] == nextPointi)
399 {
400 nextPointi = f[1];
401 nextFacei = facei;
402 break;
403 }
404 else if (f[1] == nextPointi)
405 {
406 nextPointi = f[0];
407 nextFacei = facei;
408 break;
409 }
410 }
411 }
412 }
413
414 for (label layer = 0; layer < nLayers - 1; ++layer)
415 {
416 // Offset to create front face.
417 forAll(frontFace, fp)
418 {
419 frontFace[fp] += mesh_.nPoints();
420 }
421
422 label offset = layer * mesh_.nCells();
423
424 label nei = -1;
425 if (layer != nLayers - 1)
426 {
427 nei = celli + offset + mesh_.nCells();
428 }
429
430 meshMod.addFace
431 (
432 frontFace,
433 celli + offset, // own
434 nei, // nei
435 -1, // masterPointID
436 -1, // masterEdgeID
437 nFaces++, // masterFaceID
438 false, // flipFaceFlux
439 -1, // patchID
440 -1, // zoneID
441 false // zoneFlip
442 );
443
444 if (debug)
445 {
446 Info<< frontFace << " "
447 << celli + offset << " "
448 << nei << " "
449 << nFaces - 1
450 << endl;
451 }
452 }
453 }
454
455 // Generate front and back faces
456 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
457
458 forAll(mesh_.cells(), celli)
459 {
460 const cell& cFaces = mesh_.cells()[celli];
461
462 face frontFace(cFaces.size());
463
464 // Make a loop out of faces.
465 label nextFacei = cFaces[0];
466
467 const face& f = faces[nextFacei];
468
469 label nextPointi;
470 if (mesh_.faceOwner()[nextFacei] == celli)
471 {
472 frontFace[0] = f[0];
473 nextPointi = f[1];
474 }
475 else
476 {
477 frontFace[0] = f[1];
478 nextPointi = f[0];
479 }
480
481
482 for (label i = 1; i < frontFace.size(); i++)
483 {
484 frontFace[i] = nextPointi;
485
486 // Find face containing pointi
487 forAll(cFaces, cFacei)
488 {
489 label facei = cFaces[cFacei];
490 if (facei != nextFacei)
491 {
492 const face& f = faces[facei];
493
494 if (f[0] == nextPointi)
495 {
496 nextPointi = f[1];
497 nextFacei = facei;
498 break;
499 }
500 else if (f[1] == nextPointi)
501 {
502 nextPointi = f[0];
503 nextFacei = facei;
504 break;
505 }
506 }
507 }
508 }
509
510 // Add back face.
511 meshMod.addFace
512 (
513 frontFace.reverseFace(),
514 celli, // own
515 -1, // nei
516 -1, // masterPointID
517 -1, // masterEdgeID
518 nFaces++, // masterFaceID
519 false, // flipFaceFlux
520 backPatchi_, // patchID
521 -1, // zoneID
522 false // zoneFlip
523 );
524
525 if (debug)
526 {
527 Info<< nl<<frontFace.reverseFace() << " "
528 << celli << " "
529 << nFaces - 1
530 << endl;
531 }
532
533 // Offset to create front face.
534 forAll(frontFace, fp)
535 {
536 frontFace[fp] += mesh_.nPoints()* (nLayers);
537 }
538
539 label offset = (nLayers - 1) * mesh_.nCells();
540
541 meshMod.addFace
542 (
543 frontFace,
544 celli + offset, // own
545 -1, // nei
546 -1, // masterPointID
547 -1, // masterEdgeID
548 nFaces++, // masterFaceID
549 false, // flipFaceFlux
550 frontPatchi_, // patchID
551 -1, // zoneID
552 false // zoneFlip
553 );
554
555 if (debug)
556 {
557 Info<< frontFace << " "
558 << celli + offset << " "
559 << nFaces - 1
560 << endl;
561 }
562 }
563}
564
565
566// ************************************************************************* //
PtrList< T > clone(Args &&... args) const
Make a copy by cloning each of the list elements.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
Given a 2D mesh insert all the topology changes to extrude. Does not work in parallel.
Definition: extrude2DMesh.H:64
void setRefinement(polyTopoChange &)
Play commands into polyTopoChange to extrude mesh.
void addFrontBackPatches()
Add front and back patches.
label findPatchID(const word &patchName, const bool allowNotFound=true) const
Find patch index given a name, return -1 if not found.
label start() const noexcept
The start label of boundary faces in the polyMesh face list.
void addPatches(polyPatchList &plist, const bool validBoundary=true)
Add boundary patches.
Definition: polyMesh.C:975
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
void removeBoundary()
Remove boundary patches.
Definition: polyMeshClear.C:39
label nFaces() const noexcept
Number of mesh faces.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const pointField & points
const labelIOList & zoneID
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
error FatalError
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Vector< scalar > vector
Definition: vector.H:61
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
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
labelList f(nPoints)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333