meshDualiser.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) 2020-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
27\*---------------------------------------------------------------------------*/
28
29#include "meshDualiser.H"
30#include "meshTools.H"
31#include "polyMesh.H"
32#include "polyTopoChange.H"
33#include "mapPolyMesh.H"
34#include "edgeFaceCirculator.H"
35#include "mergePoints.H"
36#include "OFstream.H"
37
38// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39
40namespace Foam
41{
42 defineTypeNameAndDebug(meshDualiser, 0);
43}
44
45
46// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47
48void Foam::meshDualiser::checkPolyTopoChange(const polyTopoChange& meshMod)
49{
50 // Assume no removed points
51 pointField points(meshMod.points().size());
52 forAll(meshMod.points(), i)
53 {
54 points[i] = meshMod.points()[i];
55 }
56
57 labelList oldToNew;
58 label nUnique = Foam::mergePoints
59 (
60 points,
61 1e-6,
62 false,
63 oldToNew
64 );
65
66 if (nUnique < points.size())
67 {
68 labelListList newToOld(invertOneToMany(nUnique, oldToNew));
69
70 forAll(newToOld, newI)
71 {
72 if (newToOld[newI].size() != 1)
73 {
75 << "duplicate verts:" << newToOld[newI]
76 << " coords:"
77 << UIndirectList<point>(points, newToOld[newI])
78 << abort(FatalError);
79 }
80 }
81 }
82}
83
84
85// Dump state so far.
86void Foam::meshDualiser::dumpPolyTopoChange
87(
88 const polyTopoChange& meshMod,
89 const fileName& prefix
90)
91{
92 OFstream str1(prefix + "Faces.obj");
93 OFstream str2(prefix + "Edges.obj");
94
95 Info<< "Dumping current polyTopoChange. Faces to " << str1.name()
96 << " , points and edges to " << str2.name() << endl;
97
98 const DynamicList<point>& points = meshMod.points();
99
100 forAll(points, pointi)
101 {
102 meshTools::writeOBJ(str1, points[pointi]);
103 meshTools::writeOBJ(str2, points[pointi]);
104 }
105
106 const DynamicList<face>& faces = meshMod.faces();
107
108 forAll(faces, facei)
109 {
110 const face& f = faces[facei];
111
112 str1<< 'f';
113 forAll(f, fp)
114 {
115 str1<< ' ' << f[fp]+1;
116 }
117 str1<< nl;
118
119 str2<< 'l';
120 forAll(f, fp)
121 {
122 str2<< ' ' << f[fp]+1;
123 }
124 str2<< ' ' << f[0]+1 << nl;
125 }
126}
127
128
129Foam::label Foam::meshDualiser::findDualCell
130(
131 const label celli,
132 const label pointi
133) const
134{
135 const labelList& dualCells = pointToDualCells_[pointi];
136
137 if (dualCells.size() == 1)
138 {
139 return dualCells[0];
140 }
141 else
142 {
143 label index = mesh_.pointCells()[pointi].find(celli);
144
145 return dualCells[index];
146 }
147}
148
149
150void Foam::meshDualiser::generateDualBoundaryEdges
151(
152 const bitSet& isBoundaryEdge,
153 const label pointi,
154 polyTopoChange& meshMod
155)
156{
157 const labelList& pEdges = mesh_.pointEdges()[pointi];
158
159 forAll(pEdges, pEdgeI)
160 {
161 label edgeI = pEdges[pEdgeI];
162
163 if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.test(edgeI))
164 {
165 const edge& e = mesh_.edges()[edgeI];
166
167 edgeToDualPoint_[edgeI] = meshMod.addPoint
168 (
169 e.centre(mesh_.points()),
170 pointi, // masterPoint
171 -1, // zoneID
172 true // inCell
173 );
174 }
175 }
176}
177
178
179// Return true if point on face has same dual cells on both owner and neighbour
180// sides.
181bool Foam::meshDualiser::sameDualCell
182(
183 const label facei,
184 const label pointi
185) const
186{
187 if (!mesh_.isInternalFace(facei))
188 {
190 << "face:" << facei << " is not internal face."
191 << abort(FatalError);
192 }
193
194 label own = mesh_.faceOwner()[facei];
195 label nei = mesh_.faceNeighbour()[facei];
196
197 return findDualCell(own, pointi) == findDualCell(nei, pointi);
198}
199
200
201Foam::label Foam::meshDualiser::addInternalFace
202(
203 const label masterPointi,
204 const label masterEdgeI,
205 const label masterFacei,
206
207 const bool edgeOrder,
208 const label dualCell0,
209 const label dualCell1,
210 const DynamicList<label>& verts,
211 polyTopoChange& meshMod
212) const
213{
214 face newFace(verts);
215
216 if (edgeOrder != (dualCell0 < dualCell1))
217 {
218 reverse(newFace);
219 }
220
221 if (debug)
222 {
223 pointField facePoints(meshMod.points(), newFace);
224
225 labelList oldToNew;
226 label nUnique = Foam::mergePoints
227 (
228 facePoints,
229 1e-6,
230 false,
231 oldToNew
232 );
233
234 if (nUnique < facePoints.size())
235 {
237 << "verts:" << verts << " newFace:" << newFace
238 << " face points:" << facePoints
239 << abort(FatalError);
240 }
241 }
242
243
244 label zoneID = -1;
245 bool zoneFlip = false;
246 if (masterFacei != -1)
247 {
248 zoneID = mesh_.faceZones().whichZone(masterFacei);
249
250 if (zoneID != -1)
251 {
252 const faceZone& fZone = mesh_.faceZones()[zoneID];
253
254 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
255 }
256 }
257
258 label dualFacei;
259
260 if (dualCell0 < dualCell1)
261 {
262 dualFacei = meshMod.addFace
263 (
264 newFace,
265 dualCell0, // own
266 dualCell1, // nei
267 masterPointi, // masterPointID
268 masterEdgeI, // masterEdgeID
269 masterFacei, // masterFaceID
270 false, // flipFaceFlux
271 -1, // patchID
272 zoneID, // zoneID
273 zoneFlip // zoneFlip
274 );
275
276 //pointField dualPoints(meshMod.points());
277 //const vector n(newFace.unitNormal(dualPoints));
278 //
279 //Pout<< "Generated internal dualFace:" << dualFacei
280 // << " verts:" << newFace
281 // << " points:" << UIndirectList<point>(meshMod.points(), newFace)
282 // << " n:" << n
283 // << " between dualowner:" << dualCell0
284 // << " dualneighbour:" << dualCell1
285 // << endl;
286 }
287 else
288 {
289 dualFacei = meshMod.addFace
290 (
291 newFace,
292 dualCell1, // own
293 dualCell0, // nei
294 masterPointi, // masterPointID
295 masterEdgeI, // masterEdgeID
296 masterFacei, // masterFaceID
297 false, // flipFaceFlux
298 -1, // patchID
299 zoneID, // zoneID
300 zoneFlip // zoneFlip
301 );
302
303 //pointField dualPoints(meshMod.points());
304 //const vector n(newFace.unitNormal(dualPoints));
305 //
306 //Pout<< "Generated internal dualFace:" << dualFacei
307 // << " verts:" << newFace
308 // << " points:" << UIndirectList<point>(meshMod.points(), newFace)
309 // << " n:" << n
310 // << " between dualowner:" << dualCell1
311 // << " dualneighbour:" << dualCell0
312 // << endl;
313 }
314 return dualFacei;
315}
316
317
318Foam::label Foam::meshDualiser::addBoundaryFace
319(
320 const label masterPointi,
321 const label masterEdgeI,
322 const label masterFacei,
323
324 const label dualCelli,
325 const label patchi,
326 const DynamicList<label>& verts,
327 polyTopoChange& meshMod
328) const
329{
330 face newFace(verts);
331
332 label zoneID = -1;
333 bool zoneFlip = false;
334 if (masterFacei != -1)
335 {
336 zoneID = mesh_.faceZones().whichZone(masterFacei);
337
338 if (zoneID != -1)
339 {
340 const faceZone& fZone = mesh_.faceZones()[zoneID];
341
342 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
343 }
344 }
345
346 label dualFacei = meshMod.addFace
347 (
348 newFace,
349 dualCelli, // own
350 -1, // nei
351 masterPointi, // masterPointID
352 masterEdgeI, // masterEdgeID
353 masterFacei, // masterFaceID
354 false, // flipFaceFlux
355 patchi, // patchID
356 zoneID, // zoneID
357 zoneFlip // zoneFlip
358 );
359
360 //pointField dualPoints(meshMod.points());
361 //const vector n(newFace.unitNormal(dualPoints));
362 //
363 //Pout<< "Generated boundary dualFace:" << dualFacei
364 // << " verts:" << newFace
365 // << " points:" << UIndirectList<point>(meshMod.points(), newFace)
366 // << " n:" << n
367 // << " on dualowner:" << dualCelli
368 // << endl;
369 return dualFacei;
370}
371
372
373// Walks around edgeI.
374// splitFace=true : creates multiple faces
375// splitFace=false: creates single face if same dual cells on both sides,
376// multiple faces otherwise.
377void Foam::meshDualiser::createFacesAroundEdge
378(
379 const bool splitFace,
380 const bitSet& isBoundaryEdge,
381 const label edgeI,
382 const label startFacei,
383 polyTopoChange& meshMod,
384 boolList& doneEFaces
385) const
386{
387 const edge& e = mesh_.edges()[edgeI];
388 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
389
391 (
392 mesh_.faces()[startFacei],
393 e[0],
394 e[1]
395 );
396
397 edgeFaceCirculator ie
398 (
399 mesh_,
400 startFacei, // face
401 true, // ownerSide
402 fp, // fp
403 isBoundaryEdge.test(edgeI) // isBoundaryEdge
404 );
405 ie.setCanonical();
406
407 bool edgeOrder = ie.sameOrder(e[0], e[1]);
408 label startFaceLabel = ie.faceLabel();
409
410 //Pout<< "At edge:" << edgeI << " verts:" << e
411 // << " points:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
412 // << " started walking at face:" << ie.faceLabel()
413 // << " verts:" << mesh_.faces()[ie.faceLabel()]
414 // << " edgeOrder:" << edgeOrder
415 // << " in direction of cell:" << ie.cellLabel()
416 // << endl;
417
418 // Walk and collect face.
419 DynamicList<label> verts(100);
420
421 if (edgeToDualPoint_[edgeI] != -1)
422 {
423 verts.append(edgeToDualPoint_[edgeI]);
424 }
425 if (faceToDualPoint_[ie.faceLabel()] != -1)
426 {
427 doneEFaces[eFaces.find(ie.faceLabel())] = true;
428 verts.append(faceToDualPoint_[ie.faceLabel()]);
429 }
430 if (cellToDualPoint_[ie.cellLabel()] != -1)
431 {
432 verts.append(cellToDualPoint_[ie.cellLabel()]);
433 }
434
435 label currentDualCell0 = findDualCell(ie.cellLabel(), e[0]);
436 label currentDualCell1 = findDualCell(ie.cellLabel(), e[1]);
437
438 ++ie;
439
440 while (true)
441 {
442 label facei = ie.faceLabel();
443
444 // Mark face as visited.
445 doneEFaces[eFaces.find(facei)] = true;
446
447 if (faceToDualPoint_[facei] != -1)
448 {
449 verts.append(faceToDualPoint_[facei]);
450 }
451
452 label celli = ie.cellLabel();
453
454 if (celli == -1)
455 {
456 // At ending boundary face. We've stored the face point above
457 // so this is the whole face.
458 break;
459 }
460
461
462 label dualCell0 = findDualCell(celli, e[0]);
463 label dualCell1 = findDualCell(celli, e[1]);
464
465 // Generate face. (always if splitFace=true; only if needed to
466 // separate cells otherwise)
467 if
468 (
469 splitFace
470 || (
471 dualCell0 != currentDualCell0
472 || dualCell1 != currentDualCell1
473 )
474 )
475 {
476 // Close current face.
477 addInternalFace
478 (
479 -1, // masterPointi
480 edgeI, // masterEdgeI
481 -1, // masterFacei
482 edgeOrder,
483 currentDualCell0,
484 currentDualCell1,
485 verts.shrink(),
486 meshMod
487 );
488
489 // Restart
490 currentDualCell0 = dualCell0;
491 currentDualCell1 = dualCell1;
492
493 verts.clear();
494 if (edgeToDualPoint_[edgeI] != -1)
495 {
496 verts.append(edgeToDualPoint_[edgeI]);
497 }
498 if (faceToDualPoint_[facei] != -1)
499 {
500 verts.append(faceToDualPoint_[facei]);
501 }
502 }
503
504 if (cellToDualPoint_[celli] != -1)
505 {
506 verts.append(cellToDualPoint_[celli]);
507 }
508
509 ++ie;
510
511 if (ie == ie.end())
512 {
513 // Back at start face (for internal edge only). See if this needs
514 // adding.
515 if (!isBoundaryEdge.test(edgeI))
516 {
517 label startDual = faceToDualPoint_[startFaceLabel];
518
519 if (startDual != -1)
520 {
521 verts.appendUniq(startDual);
522 }
523 }
524 break;
525 }
526 }
527
528 verts.shrink();
529 addInternalFace
530 (
531 -1, // masterPointi
532 edgeI, // masterEdgeI
533 -1, // masterFacei
534 edgeOrder,
535 currentDualCell0,
536 currentDualCell1,
537 verts,
538 meshMod
539 );
540}
541
542
543// Walks around circumference of facei. Creates single face. Gets given
544// starting (feature) edge to start from. Returns ending edge. (all edges
545// in form of index in faceEdges)
546void Foam::meshDualiser::createFaceFromInternalFace
547(
548 const label facei,
549 label& fp,
550 polyTopoChange& meshMod
551) const
552{
553 const face& f = mesh_.faces()[facei];
554 const labelList& fEdges = mesh_.faceEdges()[facei];
555 label own = mesh_.faceOwner()[facei];
556 label nei = mesh_.faceNeighbour()[facei];
557
558 //Pout<< "createFaceFromInternalFace : At face:" << facei
559 // << " verts:" << f
560 // << " points:" << UIndirectList<point>(mesh_.points(), f)
561 // << " started walking at edge:" << fEdges[fp]
562 // << " verts:" << mesh_.edges()[fEdges[fp]]
563 // << endl;
564
565
566 // Walk and collect face.
567 DynamicList<label> verts(100);
568
569 verts.append(faceToDualPoint_[facei]);
570 verts.append(edgeToDualPoint_[fEdges[fp]]);
571
572 // Step to vertex after edge mid
573 fp = f.fcIndex(fp);
574
575 // Get cells on either side of face at that point
576 label currentDualCell0 = findDualCell(own, f[fp]);
577 label currentDualCell1 = findDualCell(nei, f[fp]);
578
579 forAll(f, i)
580 {
581 // Check vertex
582 if (pointToDualPoint_[f[fp]] != -1)
583 {
584 verts.append(pointToDualPoint_[f[fp]]);
585 }
586
587 // Edge between fp and fp+1
588 label edgeI = fEdges[fp];
589
590 if (edgeToDualPoint_[edgeI] != -1)
591 {
592 verts.append(edgeToDualPoint_[edgeI]);
593 }
594
595 // Next vertex on edge
596 label nextFp = f.fcIndex(fp);
597
598 // Get dual cells on nextFp to check whether face needs closing.
599 label dualCell0 = findDualCell(own, f[nextFp]);
600 label dualCell1 = findDualCell(nei, f[nextFp]);
601
602 if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
603 {
604 // Check: make sure that there is a midpoint on the edge.
605 if (edgeToDualPoint_[edgeI] == -1)
606 {
608 << "face:" << facei << " verts:" << f
609 << " points:" << UIndirectList<point>(mesh_.points(), f)
610 << " no feature edge between " << f[fp]
611 << " and " << f[nextFp] << " although have different"
612 << " dual cells." << endl
613 << "point " << f[fp] << " has dual cells "
614 << currentDualCell0 << " and " << currentDualCell1
615 << " ; point "<< f[nextFp] << " has dual cells "
616 << dualCell0 << " and " << dualCell1
617 << abort(FatalError);
618 }
619
620
621 // Close current face.
622 verts.shrink();
623 addInternalFace
624 (
625 -1, // masterPointi
626 -1, // masterEdgeI
627 facei, // masterFacei
628 true, // edgeOrder,
629 currentDualCell0,
630 currentDualCell1,
631 verts,
632 meshMod
633 );
634 break;
635 }
636
637 fp = nextFp;
638 }
639}
640
641
642// Given a point on a face converts the faces around the point.
643// (pointFaces()). Gets starting face and marks off visited faces in donePFaces.
644void Foam::meshDualiser::createFacesAroundBoundaryPoint
645(
646 const label patchi,
647 const label patchPointi,
648 const label startFacei,
649 polyTopoChange& meshMod,
650 boolList& donePFaces // pFaces visited
651) const
652{
653 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
654 const polyPatch& pp = patches[patchi];
655 const labelList& pFaces = pp.pointFaces()[patchPointi];
656 const labelList& own = mesh_.faceOwner();
657
658 label pointi = pp.meshPoints()[patchPointi];
659
660 if (pointToDualPoint_[pointi] == -1)
661 {
662 // Not a feature point. Loop over all connected
663 // pointFaces.
664
665 // Starting face
666 label facei = startFacei;
667
668 DynamicList<label> verts(4);
669
670 while (true)
671 {
672 label index = pFaces.find(facei-pp.start());
673
674 // Has face been visited already?
675 if (donePFaces[index])
676 {
677 break;
678 }
679 donePFaces[index] = true;
680
681 // Insert face centre
682 verts.append(faceToDualPoint_[facei]);
683
684 label dualCelli = findDualCell(own[facei], pointi);
685
686 // Get the edge before the patchPointi
687 const face& f = mesh_.faces()[facei];
688 label fp = f.find(pointi);
689 label prevFp = f.rcIndex(fp);
690 label edgeI = mesh_.faceEdges()[facei][prevFp];
691
692 if (edgeToDualPoint_[edgeI] != -1)
693 {
694 verts.append(edgeToDualPoint_[edgeI]);
695 }
696
697 // Get next boundary face (whilst staying on edge).
698 edgeFaceCirculator circ
699 (
700 mesh_,
701 facei,
702 true, // ownerSide
703 prevFp, // index of edge in face
704 true // isBoundaryEdge
705 );
706
707 do
708 {
709 ++circ;
710 }
711 while (mesh_.isInternalFace(circ.faceLabel()));
712
713 // Step to next face
714 facei = circ.faceLabel();
715
716 if (facei < pp.start() || facei >= pp.start()+pp.size())
717 {
719 << "Walked from face on patch:" << patchi
720 << " to face:" << facei
721 << " fc:" << mesh_.faceCentres()[facei]
722 << " on patch:" << patches.whichPatch(facei)
723 << abort(FatalError);
724 }
725
726 // Check if different cell.
727 if (dualCelli != findDualCell(own[facei], pointi))
728 {
730 << "Different dual cells but no feature edge"
731 << " inbetween point:" << pointi
732 << " coord:" << mesh_.points()[pointi]
733 << abort(FatalError);
734 }
735 }
736
737 verts.shrink();
738
739 label dualCelli = findDualCell(own[facei], pointi);
740
741 //Bit dodgy: create dualface from the last face (instead of from
742 // the central point). This will also use the original faceZone to
743 // put the new face (which might span multiple original faces) in.
744
745 addBoundaryFace
746 (
747 //pointi, // masterPointi
748 -1, // masterPointi
749 -1, // masterEdgeI
750 facei, // masterFacei
751 dualCelli,
752 patchi,
753 verts,
754 meshMod
755 );
756 }
757 else
758 {
759 label facei = startFacei;
760
761 // Storage for face
762 DynamicList<label> verts(mesh_.faces()[facei].size());
763
764 // Starting point.
765 verts.append(pointToDualPoint_[pointi]);
766
767 // Find edge between pointi and next point on face.
768 const labelList& fEdges = mesh_.faceEdges()[facei];
769 label nextEdgeI = fEdges[mesh_.faces()[facei].find(pointi)];
770 if (edgeToDualPoint_[nextEdgeI] != -1)
771 {
772 verts.append(edgeToDualPoint_[nextEdgeI]);
773 }
774
775 do
776 {
777 label index = pFaces.find(facei-pp.start());
778
779 // Has face been visited already?
780 if (donePFaces[index])
781 {
782 break;
783 }
784 donePFaces[index] = true;
785
786 // Face centre
787 verts.append(faceToDualPoint_[facei]);
788
789 // Find edge before pointi on facei
790 const labelList& fEdges = mesh_.faceEdges()[facei];
791 const face& f = mesh_.faces()[facei];
792 label prevFp = f.rcIndex(f.find(pointi));
793 label edgeI = fEdges[prevFp];
794
795 if (edgeToDualPoint_[edgeI] != -1)
796 {
797 // Feature edge. Close any face so far. Note: uses face to
798 // create dualFace from. Could use pointi instead.
799 verts.append(edgeToDualPoint_[edgeI]);
800 addBoundaryFace
801 (
802 -1, // masterPointi
803 -1, // masterEdgeI
804 facei, // masterFacei
805 findDualCell(own[facei], pointi),
806 patchi,
807 verts.shrink(),
808 meshMod
809 );
810 verts.clear();
811
812 verts.append(pointToDualPoint_[pointi]);
813 verts.append(edgeToDualPoint_[edgeI]);
814 }
815
816 // Cross edgeI to next boundary face
817 edgeFaceCirculator circ
818 (
819 mesh_,
820 facei,
821 true, // ownerSide
822 prevFp, // index of edge in face
823 true // isBoundaryEdge
824 );
825
826 do
827 {
828 ++circ;
829 }
830 while (mesh_.isInternalFace(circ.faceLabel()));
831
832 // Step to next face. Quit if not on same patch.
833 facei = circ.faceLabel();
834 }
835 while
836 (
837 facei != startFacei
838 && facei >= pp.start()
839 && facei < pp.start()+pp.size()
840 );
841
842 if (verts.size() > 2)
843 {
844 // Note: face created from face, not from pointi
845 addBoundaryFace
846 (
847 -1, // masterPointi
848 -1, // masterEdgeI
849 startFacei, // masterFacei
850 findDualCell(own[facei], pointi),
851 patchi,
852 verts.shrink(),
853 meshMod
854 );
855 }
856 }
857}
858
859
860// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
861
863:
864 mesh_(mesh),
865 pointToDualCells_(mesh_.nPoints()),
866 pointToDualPoint_(mesh_.nPoints(), -1),
867 cellToDualPoint_(mesh_.nCells()),
868 faceToDualPoint_(mesh_.nFaces(), -1),
869 edgeToDualPoint_(mesh_.nEdges(), -1)
870{}
871
872
873// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
874
876(
877 const bool splitFace,
878 const labelList& featureFaces,
879 const labelList& featureEdges,
880 const labelList& singleCellFeaturePoints,
881 const labelList& multiCellFeaturePoints,
882 polyTopoChange& meshMod
883)
884{
885 const labelList& own = mesh_.faceOwner();
886 const labelList& nei = mesh_.faceNeighbour();
887 const vectorField& cellCentres = mesh_.cellCentres();
888
889 // Mark boundary edges and points.
890 // (Note: in 1.4.2 we can use the built-in mesh point ordering
891 // facility instead)
892 bitSet isBoundaryEdge(mesh_.nEdges());
893 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
894 {
895 const labelList& fEdges = mesh_.faceEdges()[facei];
896
897 isBoundaryEdge.set(fEdges);
898 }
899
900
901 if (splitFace)
902 {
903 // This is a special mode where whenever we are walking around an edge
904 // every area through a cell becomes a separate dualface. So two
905 // dual cells will probably have more than one dualface between them!
906 // This mode implies that
907 // - all faces have to be feature faces since there has to be a
908 // dualpoint at the face centre.
909 // - all edges have to be feature edges ,,
910 boolList featureFaceSet(mesh_.nFaces(), false);
911 forAll(featureFaces, i)
912 {
913 featureFaceSet[featureFaces[i]] = true;
914 }
915 label facei = featureFaceSet.find(false);
916
917 if (facei != -1)
918 {
920 << "In split-face-mode (splitFace=true) but not all faces"
921 << " marked as feature faces." << endl
922 << "First conflicting face:" << facei
923 << " centre:" << mesh_.faceCentres()[facei]
924 << abort(FatalError);
925 }
926
927 boolList featureEdgeSet(mesh_.nEdges(), false);
928 forAll(featureEdges, i)
929 {
930 featureEdgeSet[featureEdges[i]] = true;
931 }
932 label edgeI = featureEdgeSet.find(false);
933
934 if (edgeI != -1)
935 {
936 const edge& e = mesh_.edges()[edgeI];
938 << "In split-face-mode (splitFace=true) but not all edges"
939 << " marked as feature edges." << endl
940 << "First conflicting edge:" << edgeI
941 << " verts:" << e
942 << " coords:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
943 << abort(FatalError);
944 }
945 }
946 else
947 {
948 // Check that all boundary faces are feature faces.
949
950 boolList featureFaceSet(mesh_.nFaces(), false);
951 forAll(featureFaces, i)
952 {
953 featureFaceSet[featureFaces[i]] = true;
954 }
955 for
956 (
957 label facei = mesh_.nInternalFaces();
958 facei < mesh_.nFaces();
959 facei++
960 )
961 {
962 if (!featureFaceSet[facei])
963 {
965 << "Not all boundary faces marked as feature faces."
966 << endl
967 << "First conflicting face:" << facei
968 << " centre:" << mesh_.faceCentres()[facei]
969 << abort(FatalError);
970 }
971 }
972 }
973
974
975
976
977 // Start creating cells, points, and faces (in that order)
978
979
980 // 1. Mark which cells to create
981 // Mostly every point becomes one cell but sometimes (for feature points)
982 // all cells surrounding a feature point become cells. Also a non-manifold
983 // point can create two cells! So a dual cell is uniquely defined by a
984 // mesh point + cell (as in pointCells index)
985
986 // 2. Mark which face centres to create
987
988 // 3. Internal faces can now consist of
989 // - only cell centres of walk around edge
990 // - cell centres + face centres of walk around edge
991 // - same but now other side is not a single cell
992
993 // 4. Boundary faces (or internal faces between cell zones!) now consist of
994 // - walk around boundary point.
995
996
997
998 autoPtr<OFstream> dualCcStr;
999 if (debug)
1000 {
1001 dualCcStr.reset(new OFstream("dualCc.obj"));
1002 Pout<< "Dumping centres of dual cells to " << dualCcStr().name()
1003 << endl;
1004 }
1005
1006
1007 // Dual cells (from points)
1008 // ~~~~~~~~~~~~~~~~~~~~~~~~
1009
1010 // pointToDualCells_[pointi]
1011 // - single entry : all cells surrounding point all become the same
1012 // cell.
1013 // - multiple entries: in order of pointCells.
1014
1015
1016 // feature points that become single cell
1017 forAll(singleCellFeaturePoints, i)
1018 {
1019 label pointi = singleCellFeaturePoints[i];
1020
1021 pointToDualPoint_[pointi] = meshMod.addPoint
1022 (
1023 mesh_.points()[pointi],
1024 pointi, // masterPoint
1025 mesh_.pointZones().whichZone(pointi), // zoneID
1026 true // inCell
1027 );
1028
1029 // Generate single cell
1030 pointToDualCells_[pointi].setSize(1);
1031 pointToDualCells_[pointi][0] = meshMod.addCell
1032 (
1033 pointi, //masterPointID,
1034 -1, //masterEdgeID,
1035 -1, //masterFaceID,
1036 -1, //masterCellID,
1037 -1 //zoneID
1038 );
1039 if (dualCcStr)
1040 {
1041 meshTools::writeOBJ(*dualCcStr, mesh_.points()[pointi]);
1042 }
1043 }
1044
1045 // feature points that become multiple cells
1046 forAll(multiCellFeaturePoints, i)
1047 {
1048 label pointi = multiCellFeaturePoints[i];
1049
1050 if (pointToDualCells_[pointi].size() > 0)
1051 {
1053 << "Point " << pointi << " at:" << mesh_.points()[pointi]
1054 << " is both in singleCellFeaturePoints"
1055 << " and multiCellFeaturePoints."
1056 << abort(FatalError);
1057 }
1058
1059 pointToDualPoint_[pointi] = meshMod.addPoint
1060 (
1061 mesh_.points()[pointi],
1062 pointi, // masterPoint
1063 mesh_.pointZones().whichZone(pointi), // zoneID
1064 true // inCell
1065 );
1066
1067 // Create dualcell for every cell connected to dual point
1068
1069 const labelList& pCells = mesh_.pointCells()[pointi];
1070
1071 pointToDualCells_[pointi].setSize(pCells.size());
1072
1073 forAll(pCells, pCelli)
1074 {
1075 pointToDualCells_[pointi][pCelli] = meshMod.addCell
1076 (
1077 pointi, //masterPointID
1078 -1, //masterEdgeID
1079 -1, //masterFaceID
1080 -1, //masterCellID
1081 mesh_.cellZones().whichZone(pCells[pCelli]) //zoneID
1082 );
1083 if (dualCcStr)
1084 {
1086 (
1087 *dualCcStr,
1088 0.5*(mesh_.points()[pointi]+cellCentres[pCells[pCelli]])
1089 );
1090 }
1091 }
1092 }
1093 // Normal points
1094 forAll(mesh_.points(), pointi)
1095 {
1096 if (pointToDualCells_[pointi].empty())
1097 {
1098 pointToDualCells_[pointi].setSize(1);
1099 pointToDualCells_[pointi][0] = meshMod.addCell
1100 (
1101 pointi, //masterPointID,
1102 -1, //masterEdgeID,
1103 -1, //masterFaceID,
1104 -1, //masterCellID,
1105 -1 //zoneID
1106 );
1107
1108 if (dualCcStr)
1109 {
1110 meshTools::writeOBJ(*dualCcStr, mesh_.points()[pointi]);
1111 }
1112 }
1113 }
1114
1115
1116 // Dual points (from cell centres, feature faces, feature edges)
1117 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1118
1119 forAll(cellToDualPoint_, celli)
1120 {
1121 cellToDualPoint_[celli] = meshMod.addPoint
1122 (
1123 cellCentres[celli],
1124 mesh_.faces()[mesh_.cells()[celli][0]][0], // masterPoint
1125 -1, // zoneID
1126 true // inCell
1127 );
1128 }
1129
1130 // From face to dual point
1131
1132 forAll(featureFaces, i)
1133 {
1134 label facei = featureFaces[i];
1135
1136 faceToDualPoint_[facei] = meshMod.addPoint
1137 (
1138 mesh_.faceCentres()[facei],
1139 mesh_.faces()[facei][0], // masterPoint
1140 -1, // zoneID
1141 true // inCell
1142 );
1143 }
1144 // Detect whether different dual cells on either side of a face. This
1145 // would necessitate having a dual face built from the face and thus a
1146 // dual point at the face centre.
1147 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1148 {
1149 if (faceToDualPoint_[facei] == -1)
1150 {
1151 const face& f = mesh_.faces()[facei];
1152
1153 forAll(f, fp)
1154 {
1155 label ownDualCell = findDualCell(own[facei], f[fp]);
1156 label neiDualCell = findDualCell(nei[facei], f[fp]);
1157
1158 if (ownDualCell != neiDualCell)
1159 {
1160 faceToDualPoint_[facei] = meshMod.addPoint
1161 (
1162 mesh_.faceCentres()[facei],
1163 f[fp], // masterPoint
1164 -1, // zoneID
1165 true // inCell
1166 );
1167
1168 break;
1169 }
1170 }
1171 }
1172 }
1173
1174 // From edge to dual point
1175
1176 forAll(featureEdges, i)
1177 {
1178 label edgeI = featureEdges[i];
1179
1180 const edge& e = mesh_.edges()[edgeI];
1181
1182 edgeToDualPoint_[edgeI] = meshMod.addPoint
1183 (
1184 e.centre(mesh_.points()),
1185 e[0], // masterPoint
1186 -1, // zoneID
1187 true // inCell
1188 );
1189 }
1190
1191 // Detect whether different dual cells on either side of an edge. This
1192 // would neccesitate having a dual face built perpendicular to the edge
1193 // and thus a dual point at the mid of the edge.
1194 // Note: not really true - the face can be built without the edge centre!
1195 const labelListList& edgeCells = mesh_.edgeCells();
1196
1197 forAll(edgeCells, edgeI)
1198 {
1199 if (edgeToDualPoint_[edgeI] == -1)
1200 {
1201 const edge& e = mesh_.edges()[edgeI];
1202
1203 // We need a point on the edge if not all cells on both sides
1204 // are the same.
1205
1206 const labelList& eCells = mesh_.edgeCells()[edgeI];
1207
1208 label dualE0 = findDualCell(eCells[0], e[0]);
1209 label dualE1 = findDualCell(eCells[0], e[1]);
1210
1211 for (label i = 1; i < eCells.size(); i++)
1212 {
1213 label newDualE0 = findDualCell(eCells[i], e[0]);
1214
1215 if (dualE0 != newDualE0)
1216 {
1217 edgeToDualPoint_[edgeI] = meshMod.addPoint
1218 (
1219 e.centre(mesh_.points()),
1220 e[0], // masterPoint
1221 mesh_.pointZones().whichZone(e[0]), // zoneID
1222 true // inCell
1223 );
1224
1225 break;
1226 }
1227
1228 label newDualE1 = findDualCell(eCells[i], e[1]);
1229
1230 if (dualE1 != newDualE1)
1231 {
1232 edgeToDualPoint_[edgeI] = meshMod.addPoint
1233 (
1234 e.centre(mesh_.points()),
1235 e[1], // masterPoint
1236 mesh_.pointZones().whichZone(e[1]), // zoneID
1237 true // inCell
1238 );
1239
1240 break;
1241 }
1242 }
1243 }
1244 }
1245
1246 // Make sure all boundary edges emanating from feature points are
1247 // feature edges as well.
1248 forAll(singleCellFeaturePoints, i)
1249 {
1250 generateDualBoundaryEdges
1251 (
1252 isBoundaryEdge,
1253 singleCellFeaturePoints[i],
1254 meshMod
1255 );
1256 }
1257 forAll(multiCellFeaturePoints, i)
1258 {
1259 generateDualBoundaryEdges
1260 (
1261 isBoundaryEdge,
1262 multiCellFeaturePoints[i],
1263 meshMod
1264 );
1265 }
1266
1267
1268 // Check for duplicate points
1269 if (debug)
1270 {
1271 dumpPolyTopoChange(meshMod, "generatedPoints_");
1272 checkPolyTopoChange(meshMod);
1273 }
1274
1275
1276 // Now we have all points and cells
1277 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1278 // - pointToDualCells_ : per point a single dualCell or multiple dualCells
1279 // - pointToDualPoint_ : per point -1 or the dual point at the coordinate
1280 // - edgeToDualPoint_ : per edge -1 or the edge centre
1281 // - faceToDualPoint_ : per face -1 or the face centre
1282 // - cellToDualPoint_ : per cell the cell centre
1283 // Now we have to walk all edges and construct faces. Either single face
1284 // per edge or multiple (-if nonmanifold edge -if different dualcells)
1285
1286 const edgeList& edges = mesh_.edges();
1287
1288 forAll(edges, edgeI)
1289 {
1290 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
1291
1292 boolList doneEFaces(eFaces.size(), false);
1293
1294 forAll(eFaces, i)
1295 {
1296 if (!doneEFaces[i])
1297 {
1298 // We found a face that hasn't yet been visited. This might
1299 // happen for non-manifold edges where a single edge can
1300 // become multiple faces.
1301
1302 label startFacei = eFaces[i];
1303
1304 //Pout<< "Walking edge:" << edgeI
1305 // << " points:" << mesh_.points()[e[0]]
1306 // << mesh_.points()[e[1]]
1307 // << " startFace:" << startFacei
1308 // << " at:" << mesh_.faceCentres()[startFacei]
1309 // << endl;
1310
1311 createFacesAroundEdge
1312 (
1313 splitFace,
1314 isBoundaryEdge,
1315 edgeI,
1316 startFacei,
1317 meshMod,
1318 doneEFaces
1319 );
1320 }
1321 }
1322 }
1323
1324 if (debug)
1325 {
1326 dumpPolyTopoChange(meshMod, "generatedFacesFromEdges_");
1327 }
1328
1329 // Create faces from feature faces. These can be internal or external faces.
1330 // - feature face : centre needs to be included.
1331 // - single cells on either side: triangulate
1332 // - multiple cells: create single face between unique cell pair. Only
1333 // create face where cells differ on either side.
1334 // - non-feature face : inbetween cell zones.
1335 forAll(faceToDualPoint_, facei)
1336 {
1337 if (faceToDualPoint_[facei] != -1 && mesh_.isInternalFace(facei))
1338 {
1339 const face& f = mesh_.faces()[facei];
1340 const labelList& fEdges = mesh_.faceEdges()[facei];
1341
1342 // Starting edge
1343 label fp = 0;
1344
1345 do
1346 {
1347 // Find edge that is in dual mesh and where the
1348 // next point (fp+1) has different dual cells on either side.
1349 bool foundStart = false;
1350
1351 do
1352 {
1353 if
1354 (
1355 edgeToDualPoint_[fEdges[fp]] != -1
1356 && !sameDualCell(facei, f.nextLabel(fp))
1357 )
1358 {
1359 foundStart = true;
1360 break;
1361 }
1362 fp = f.fcIndex(fp);
1363 }
1364 while (fp != 0);
1365
1366 if (!foundStart)
1367 {
1368 break;
1369 }
1370
1371 // Walk from edge fp and generate a face.
1372 createFaceFromInternalFace
1373 (
1374 facei,
1375 fp,
1376 meshMod
1377 );
1378 }
1379 while (fp != 0);
1380 }
1381 }
1382
1383 if (debug)
1384 {
1385 dumpPolyTopoChange(meshMod, "generatedFacesFromFeatFaces_");
1386 }
1387
1388
1389 // Create boundary faces. Every boundary point has one or more dualcells.
1390 // These need to be closed.
1391 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1392
1393 forAll(patches, patchi)
1394 {
1395 const polyPatch& pp = patches[patchi];
1396
1397 const labelListList& pointFaces = pp.pointFaces();
1398
1399 forAll(pointFaces, patchPointi)
1400 {
1401 const labelList& pFaces = pointFaces[patchPointi];
1402
1403 boolList donePFaces(pFaces.size(), false);
1404
1405 forAll(pFaces, i)
1406 {
1407 if (!donePFaces[i])
1408 {
1409 // Starting face
1410 label startFacei = pp.start()+pFaces[i];
1411
1412 //Pout<< "Walking around point:" << pointi
1413 // << " coord:" << mesh_.points()[pointi]
1414 // << " on patch:" << patchi
1415 // << " startFace:" << startFacei
1416 // << " at:" << mesh_.faceCentres()[startFacei]
1417 // << endl;
1418
1419 createFacesAroundBoundaryPoint
1420 (
1421 patchi,
1422 patchPointi,
1423 startFacei,
1424 meshMod,
1425 donePFaces // pFaces visited
1426 );
1427 }
1428 }
1429 }
1430 }
1431
1432 if (debug)
1433 {
1434 dumpPolyTopoChange(meshMod, "generatedFacesFromBndFaces_");
1435 }
1436}
1437
1438
1439// ************************************************************************* //
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:330
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
label rcIndex(const label i) const noexcept
Definition: UListI.H:67
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
static label getMinIndex(const face &f, const label v0, const label v1)
Helper: find index in face of edge or -1. Index is such that edge is.
Creates dual of polyMesh. Every point becomes a cell (or multiple cells for feature points),...
Definition: meshDualiser.H:69
void setRefinement(const bool splitFace, const labelList &featureFaces, const labelList &featureEdges, const labelList &singleCellFeaturePoints, const labelList &multiCellFeaturePoints, polyTopoChange &meshMod)
Insert all changes into meshMod to convert the polyMesh into.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
Definition: polyMesh.H:504
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
const pointZoneMesh & pointZones() const noexcept
Return point zone mesh.
Definition: polyMesh.H:492
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
const labelListList & pointCells() const
label nInternalFaces() const noexcept
Number of internal faces.
const vectorField & cellCentres() const
label nFaces() const noexcept
Number of mesh faces.
const labelListList & edgeFaces() const
const labelListList & faceEdges() const
const labelListList & edgeCells() const
const cellList & cells() const
label nEdges() const
Number of mesh edges.
#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 labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
const pointField & points
label nPoints
const labelIOList & zoneID
Geometric merging of points. See below.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Definition: UListI.H:449
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
label mergePoints(const PointList &points, labelList &pointToUnique, labelList &uniquePoints, const scalar mergeTol=SMALL, const bool verbose=false)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< bool > boolList
A List of bools.
Definition: List.H:64
error FatalError
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:114
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333