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