FreeFOAM The Cross-Platform CFD Toolkit
surfaceSplitNonManifolds.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 Application
25  surfaceSplitNonManifolds
26 
27 Description
28  Takes multiply connected surface and tries to split surface at
29  multiply connected edges by duplicating points.
30 
31  Introduces concept of
32  - borderEdge. Edge with 4 faces connected to it.
33  - borderPoint. Point connected to exactly 2 borderEdges.
34  - borderLine. Connected list of borderEdges.
35 
36  By duplicating borderPoints this will split 'borderLines'. As a
37  preprocessing step it can detect borderEdges without any borderPoints
38  and explicitly split these triangles.
39 
40  The problems in this algorithm are:
41  - determining which two (of the four) faces form a surface. Done by walking
42  face-edge-face while keeping and edge or point on the borderEdge
43  borderPoint.
44  - determining the outwards pointing normal to be used to slightly offset the
45  duplicated point.
46 
47  Uses sortedEdgeFaces quite a bit.
48 
49  Is tested on simple borderLines resulting from extracting a surface
50  from a hex mesh. Will quite possibly go wrong on more complicated border
51  lines (i.e. ones forming a loop).
52 
53  Dumps surface every so often since might take a long time to complete.
54 
55 Usage
56 
57  - surfaceSplitNonManifolds [OPTIONS] <Foam surface file> <surface output file>
58 
59  @param <Foam surface file> \n
60  @todo Detailed description of argument.
61 
62  @param <surface output file> \n
63  @todo Detailed description of argument.
64 
65  @param -debug \n
66  Generate debugging output.
67 
68  @param -case <dir>\n
69  Case directory.
70 
71  @param -help \n
72  Display help message.
73 
74  @param -doc \n
75  Display Doxygen API documentation page for this application.
76 
77  @param -srcDoc \n
78  Display Doxygen source documentation page for this application.
79 
80 \*---------------------------------------------------------------------------*/
81 
82 #include <OpenFOAM/argList.H>
83 #include <triSurface/triSurface.H>
84 #include <OpenFOAM/OFstream.H>
85 #include <OpenFOAM/ListOps.H>
87 
88 using namespace Foam;
89 
90 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
91 
92 void writeOBJ(Ostream& os, const pointField& pts)
93 {
94  forAll(pts, i)
95  {
96  const point& pt = pts[i];
97 
98  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
99  }
100 }
101 
102 
103 void dumpPoints(const triSurface& surf, const labelList& borderPoint)
104 {
105  fileName fName("borderPoints.obj");
106 
107  Info<< "Dumping borderPoints as Lightwave .obj file to " << fName
108  << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
109 
110  OFstream os(fName);
111 
112  forAll(borderPoint, pointI)
113  {
114  if (borderPoint[pointI] != -1)
115  {
116  const point& pt = surf.localPoints()[pointI];
117 
118  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
119  }
120  }
121 }
122 
123 
124 void dumpEdges(const triSurface& surf, const boolList& borderEdge)
125 {
126  fileName fName("borderEdges.obj");
127 
128  Info<< "Dumping borderEdges as Lightwave .obj file to " << fName
129  << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
130 
131  OFstream os(fName);
132 
133  writeOBJ(os, surf.localPoints());
134 
135  forAll(borderEdge, edgeI)
136  {
137  if (borderEdge[edgeI])
138  {
139  const edge& e = surf.edges()[edgeI];
140 
141  os << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
142  }
143  }
144 }
145 
146 
147 void dumpFaces
148 (
149  const fileName& fName,
150  const triSurface& surf,
151  const Map<label>& connectedFaces
152 )
153 {
154  Info<< "Dumping connectedFaces as Lightwave .obj file to " << fName
155  << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
156 
157  OFstream os(fName);
158 
159  for
160  (
161  Map<label>::const_iterator iter = connectedFaces.begin();
162  iter != connectedFaces.end();
163  ++iter
164  )
165  {
166  const labelledTri& f = surf.localFaces()[iter.key()];
167 
168  point ctr(f.centre(surf.localPoints()));
169 
170  os << "v " << ctr.x() << ' ' << ctr.y() << ' ' << ctr.z() << endl;
171  }
172 }
173 
174 
175 void testSortedEdgeFaces(const triSurface& surf)
176 {
177  const labelListList& edgeFaces = surf.edgeFaces();
178  const labelListList& sortedEdgeFaces = surf.sortedEdgeFaces();
179 
180  forAll(edgeFaces, edgeI)
181  {
182  const labelList& myFaces = edgeFaces[edgeI];
183  const labelList& sortMyFaces = sortedEdgeFaces[edgeI];
184 
185  forAll(myFaces, i)
186  {
187  if (findIndex(sortMyFaces, myFaces[i]) == -1)
188  {
189  FatalErrorIn("testSortedEdgeFaces(..)") << abort(FatalError);
190  }
191  }
192  forAll(sortMyFaces, i)
193  {
194  if (findIndex(myFaces, sortMyFaces[i]) == -1)
195  {
196  FatalErrorIn("testSortedEdgeFaces(..)") << abort(FatalError);
197  }
198  }
199  }
200 }
201 
202 
203 // Mark off all border edges. Return number of border edges.
204 label markBorderEdges
205 (
206  const bool debug,
207  const triSurface& surf,
208  boolList& borderEdge
209 )
210 {
211  label nBorderEdges = 0;
212 
213  const labelListList& edgeFaces = surf.edgeFaces();
214 
215  forAll(edgeFaces, edgeI)
216  {
217  if (edgeFaces[edgeI].size() == 4)
218  {
219  borderEdge[edgeI] = true;
220 
221  nBorderEdges++;
222  }
223  }
224 
225  if (debug)
226  {
227  dumpEdges(surf, borderEdge);
228  }
229 
230  return nBorderEdges;
231 }
232 
233 
234 // Mark off all border points. Return number of border points. Border points
235 // marked by setting value to newly introduced point.
236 label markBorderPoints
237 (
238  const bool debug,
239  const triSurface& surf,
240  const boolList& borderEdge,
241  labelList& borderPoint
242 )
243 {
244  label nPoints = surf.nPoints();
245 
246  const labelListList& pointEdges = surf.pointEdges();
247 
248  forAll(pointEdges, pointI)
249  {
250  const labelList& pEdges = pointEdges[pointI];
251 
252  label nBorderEdges = 0;
253 
254  forAll(pEdges, i)
255  {
256  if (borderEdge[pEdges[i]])
257  {
258  nBorderEdges++;
259  }
260  }
261 
262  if (nBorderEdges == 2 && borderPoint[pointI] == -1)
263  {
264  borderPoint[pointI] = nPoints++;
265  }
266  }
267 
268  label nBorderPoints = nPoints - surf.nPoints();
269 
270  if (debug)
271  {
272  dumpPoints(surf, borderPoint);
273  }
274 
275  return nBorderPoints;
276 }
277 
278 
279 // Get minumum length of edges connected to pointI
280 // Serves to get some local length scale.
281 scalar minEdgeLen(const triSurface& surf, const label pointI)
282 {
283  const pointField& points = surf.localPoints();
284 
285  const labelList& pEdges = surf.pointEdges()[pointI];
286 
287  scalar minLen = GREAT;
288 
289  forAll(pEdges, i)
290  {
291  label edgeI = pEdges[i];
292 
293  scalar len = surf.edges()[edgeI].mag(points);
294 
295  if (len < minLen)
296  {
297  minLen = len;
298  }
299  }
300  return minLen;
301 }
302 
303 
304 // Find edge among edgeLabels with endpoints v0,v1
305 label findEdge
306 (
307  const triSurface& surf,
308  const labelList& edgeLabels,
309  const label v0,
310  const label v1
311 )
312 {
313  forAll(edgeLabels, i)
314  {
315  label edgeI = edgeLabels[i];
316 
317  const edge& e = surf.edges()[edgeI];
318 
319  if
320  (
321  (
322  e.start() == v0
323  && e.end() == v1
324  )
325  || (
326  e.start() == v1
327  && e.end() == v0
328  )
329  )
330  {
331  return edgeI;
332  }
333  }
334 
335 
336  FatalErrorIn("findEdge(..)") << "Cannot find edge with labels " << v0
337  << ' ' << v1 << " in candidates " << edgeLabels
338  << " with vertices:" << UIndirectList<edge>(surf.edges(), edgeLabels)()
339  << abort(FatalError);
340 
341  return -1;
342 }
343 
344 
345 // Get the other edge connected to pointI on faceI.
346 label otherEdge
347 (
348  const triSurface& surf,
349  const label faceI,
350  const label otherEdgeI,
351  const label pointI
352 )
353 {
354  const labelList& fEdges = surf.faceEdges()[faceI];
355 
356  forAll(fEdges, i)
357  {
358  label edgeI = fEdges[i];
359 
360  const edge& e = surf.edges()[edgeI];
361 
362  if
363  (
364  edgeI != otherEdgeI
365  && (
366  e.start() == pointI
367  || e.end() == pointI
368  )
369  )
370  {
371  return edgeI;
372  }
373  }
374 
375  FatalErrorIn("otherEdge(..)") << "Cannot find other edge on face " << faceI
376  << " verts:" << surf.localPoints()[faceI]
377  << " connected to point " << pointI
378  << " faceEdges:" << UIndirectList<edge>(surf.edges(), fEdges)()
379  << abort(FatalError);
380 
381  return -1;
382 }
383 
384 
385 // Starting from startPoint on startEdge on startFace walk along border
386 // and insert faces along the way. Walk keeps always one point or one edge
387 // on the border line.
388 void walkSplitLine
389 (
390  const triSurface& surf,
391  const boolList& borderEdge,
392  const labelList& borderPoint,
393 
394  const label startFaceI,
395  const label startEdgeI, // is border edge
396  const label startPointI, // is border point
397 
398  Map<label>& faceToEdge,
400 )
401 {
402  label faceI = startFaceI;
403  label edgeI = startEdgeI;
404  label pointI = startPointI;
405 
406  do
407  {
408  //
409  // Stick to pointI and walk face-edge-face until back on border edge.
410  //
411 
412  do
413  {
414  // Cross face to next edge.
415  edgeI = otherEdge(surf, faceI, edgeI, pointI);
416 
417  if (borderEdge[edgeI])
418  {
419  if (!faceToEdge.insert(faceI, edgeI))
420  {
421  // Was already visited.
422  return;
423  }
424  else
425  {
426  // First visit to this borderEdge. We're back on borderline.
427  break;
428  }
429  }
430  else if (!faceToPoint.insert(faceI, pointI))
431  {
432  // Was already visited.
433  return;
434  }
435 
436  // Cross edge to other face
437  const labelList& eFaces = surf.edgeFaces()[edgeI];
438 
439  if (eFaces.size() != 2)
440  {
441  FatalErrorIn("walkSplitPoint(..)")
442  << "Can only handle edges with 2 or 4 edges for now."
443  << abort(FatalError);
444  }
445 
446  if (eFaces[0] == faceI)
447  {
448  faceI = eFaces[1];
449  }
450  else if (eFaces[1] == faceI)
451  {
452  faceI = eFaces[0];
453  }
454  else
455  {
456  FatalErrorIn("walkSplitPoint(..)") << abort(FatalError);
457  }
458  }
459  while (true);
460 
461 
462  //
463  // Back on border edge. Cross to other point on edge.
464  //
465 
466  pointI = surf.edges()[edgeI].otherVertex(pointI);
467 
468  if (borderPoint[pointI] == -1)
469  {
470  return;
471  }
472 
473  }
474  while (true);
475 }
476 
477 
478 // Find second face which is from same surface i.e. has points on the
479 // shared edge in reverse order.
480 label sharedFace
481 (
482  const triSurface& surf,
483  const label firstFaceI,
484  const label sharedEdgeI
485 )
486 {
487  // Find ordering of face points in edge.
488 
489  const edge& e = surf.edges()[sharedEdgeI];
490 
491  const labelledTri& f = surf.localFaces()[firstFaceI];
492 
493  label startIndex = findIndex(f, e.start());
494 
495  // points in face in same order as edge
496  bool edgeOrder = (f[f.fcIndex(startIndex)] == e.end());
497 
498  // Get faces using edge in sorted order. (sorted such that walking
499  // around them in anti-clockwise order corresponds to edge vector
500  // acc. to right-hand rule)
501  const labelList& eFaces = surf.sortedEdgeFaces()[sharedEdgeI];
502 
503  // Get position of face in sorted edge faces
504  label faceIndex = findIndex(eFaces, firstFaceI);
505 
506  if (edgeOrder)
507  {
508  // Get face before firstFaceI
509  return eFaces[eFaces.rcIndex(faceIndex)];
510  }
511  else
512  {
513  // Get face after firstFaceI
514  return eFaces[eFaces.fcIndex(faceIndex)];
515  }
516 }
517 
518 
519 // Calculate (inward pointing) normals on edges shared by faces in faceToEdge and
520 // averages them to pointNormals.
521 void calcPointVecs
522 (
523  const triSurface& surf,
524  const Map<label>& faceToEdge,
525  const Map<label>& faceToPoint,
526  vectorField& borderPointVec
527 )
528 {
529  const labelListList& sortedEdgeFaces = surf.sortedEdgeFaces();
530  const edgeList& edges = surf.edges();
531  const pointField& points = surf.localPoints();
532 
533  boolList edgeDone(surf.nEdges(), false);
534 
535  for
536  (
537  Map<label>::const_iterator iter = faceToEdge.begin();
538  iter != faceToEdge.end();
539  ++iter
540  )
541  {
542  label edgeI = iter();
543 
544  if (!edgeDone[edgeI])
545  {
546  edgeDone[edgeI] = true;
547 
548  // Get the two connected faces in sorted order
549  // Note: should have stored this when walking ...
550 
551  label face0I = -1;
552  label face1I = -1;
553 
554  const labelList& eFaces = sortedEdgeFaces[edgeI];
555 
556  forAll(eFaces, i)
557  {
558  label faceI = eFaces[i];
559 
560  if (faceToEdge.found(faceI))
561  {
562  if (face0I == -1)
563  {
564  face0I = faceI;
565  }
566  else if (face1I == -1)
567  {
568  face1I = faceI;
569 
570  break;
571  }
572  }
573  }
574 
575  if (face0I == -1 && face1I == -1)
576  {
577  Info<< "Writing surface to errorSurf.ftr" << endl;
578 
579  surf.write("errorSurf.ftr");
580 
581  FatalErrorIn("calcPointVecs(..)")
582  << "Cannot find two faces using border edge " << edgeI
583  << " verts:" << edges[edgeI]
584  << " eFaces:" << eFaces << endl
585  << "face0I:" << face0I
586  << " face1I:" << face1I << nl
587  << "faceToEdge:" << faceToEdge << nl
588  << "faceToPoint:" << faceToPoint
589  << "Written surface to errorSurf.ftr"
590  << abort(FatalError);
591  }
592 
593  // Now we have edge and the two faces in counter-clockwise order
594  // as seen from edge vector. Calculate normal.
595  const edge& e = edges[edgeI];
596  vector eVec = e.vec(points);
597 
598  // Determine vector as average of the vectors in the two faces.
599  // If there is only one face available use only one vector.
600  vector midVec(vector::zero);
601 
602  if (face0I != -1)
603  {
604  label v0 = triSurfaceTools::oppositeVertex(surf, face0I, edgeI);
605  vector e0 = (points[v0] - points[e.start()]) ^ eVec;
606  e0 /= mag(e0);
607  midVec = e0;
608  }
609 
610  if (face1I != -1)
611  {
612  label v1 = triSurfaceTools::oppositeVertex(surf, face1I, edgeI);
613  vector e1 = (points[e.start()] - points[v1]) ^ eVec;
614  e1 /= mag(e1);
615  midVec += e1;
616  }
617 
618  scalar magMidVec = mag(midVec);
619 
620  if (magMidVec > SMALL)
621  {
622  midVec /= magMidVec;
623 
624  // Average to pointVec
625  borderPointVec[e.start()] += midVec;
626  borderPointVec[e.end()] += midVec;
627  }
628  }
629  }
630 }
631 
632 
633 // Renumbers vertices (of triangles in faceToEdge) of which the pointMap is
634 // not -1.
635 void renumberFaces
636 (
637  const triSurface& surf,
638  const labelList& pointMap,
639  const Map<label>& faceToEdge,
640  List<labelledTri>& newTris
641 )
642 {
643  for
644  (
645  Map<label>::const_iterator iter = faceToEdge.begin();
646  iter != faceToEdge.end();
647  ++iter
648  )
649  {
650  label faceI = iter.key();
651 
652  const labelledTri& f = surf.localFaces()[faceI];
653 
654  forAll(f, fp)
655  {
656  if (pointMap[f[fp]] != -1)
657  {
658  newTris[faceI][fp] = pointMap[f[fp]];
659  }
660  }
661  }
662 }
663 
664 
665 // Split all borderEdges that don't have borderPoint. Return true if split
666 // anything.
667 bool splitBorderEdges
668 (
669  triSurface& surf,
670  const boolList& borderEdge,
671  const labelList& borderPoint
672 )
673 {
674  labelList edgesToBeSplit(surf.nEdges());
675  label nSplit = 0;
676 
677  forAll(borderEdge, edgeI)
678  {
679  if (borderEdge[edgeI])
680  {
681  const edge& e = surf.edges()[edgeI];
682 
683  if (borderPoint[e.start()] == -1 && borderPoint[e.end()] == -1)
684  {
685  // None of the points of the edge is borderPoint. Split edge
686  // to introduce border point.
687  edgesToBeSplit[nSplit++] = edgeI;
688  }
689  }
690  }
691  edgesToBeSplit.setSize(nSplit);
692 
693  if (nSplit > 0)
694  {
695  Info<< "Splitting surface along " << nSplit << " borderEdges that don't"
696  << " neighbour other borderEdges" << nl << endl;
697 
698  surf = triSurfaceTools::greenRefine(surf, edgesToBeSplit);
699 
700  return true;
701  }
702  else
703  {
704  Info<< "No edges to be split" <<nl << endl;
705 
706  return false;
707  }
708 }
709 
710 
711 // Main program:
712 
713 int main(int argc, char *argv[])
714 {
717 
718  argList::validArgs.append("surface file");
719  argList::validArgs.append("output surface file");
720  argList::validOptions.insert("debug", "");
721 
722  argList args(argc, argv);
723 
724  fileName inSurfName(args.additionalArgs()[0]);
725  fileName outSurfName(args.additionalArgs()[1]);
726  bool debug = args.optionFound("debug");
727 
728 
729  Info<< "Reading surface from " << inSurfName << endl;
730 
731  triSurface surf(inSurfName);
732 
733  // Make sure sortedEdgeFaces is calculated correctly
734  testSortedEdgeFaces(surf);
735 
736  // Get all quad connected edges. These are seen as borders when walking.
737  boolList borderEdge(surf.nEdges(), false);
738  markBorderEdges(debug, surf, borderEdge);
739 
740  // Points on two sides connected to borderEdges are called borderPoints and
741  // will be duplicated. borderPoint contains label of newly introduced vertex.
742  labelList borderPoint(surf.nPoints(), -1);
743  markBorderPoints(debug, surf, borderEdge, borderPoint);
744 
745  // Split edges where there would be no borderPoint to duplicate.
746  splitBorderEdges(surf, borderEdge, borderPoint);
747 
748 
749  Info<< "Writing split surface to " << outSurfName << nl << endl;
750  surf.write(outSurfName);
751  Info<< "Finished writing surface to " << outSurfName << nl << endl;
752 
753 
754  // Last iteration values.
755  label nOldBorderEdges = -1;
756  label nOldBorderPoints = -1;
757 
758  label iteration = 0;
759 
760  do
761  {
762  // Redo borderEdge/borderPoint calculation.
763  boolList borderEdge(surf.nEdges(), false);
764 
765  label nBorderEdges = markBorderEdges(debug, surf, borderEdge);
766 
767  if (nBorderEdges == 0)
768  {
769  Info<< "Found no border edges. Exiting." << nl << nl;
770 
771  break;
772  }
773 
774  // Label of newly introduced duplicate.
775  labelList borderPoint(surf.nPoints(), -1);
776 
777  label nBorderPoints =
778  markBorderPoints
779  (
780  debug,
781  surf,
782  borderEdge,
783  borderPoint
784  );
785 
786  if (nBorderPoints == 0)
787  {
788  Info<< "Found no border points. Exiting." << nl << nl;
789 
790  break;
791  }
792 
793  Info<< "Found:\n"
794  << " border edges :" << nBorderEdges << nl
795  << " border points:" << nBorderPoints << nl
796  << endl;
797 
798  if
799  (
800  nBorderPoints == nOldBorderPoints
801  && nBorderEdges == nOldBorderEdges
802  )
803  {
804  Info<< "Stopping since number of border edges and point is same"
805  << " as in previous iteration" << nl << endl;
806 
807  break;
808  }
809 
810  //
811  // Define splitLine as a series of connected borderEdges. Find start
812  // of one (as edge and point on edge)
813  //
814 
815  label startEdgeI = -1;
816  label startPointI = -1;
817 
818  forAll(borderEdge, edgeI)
819  {
820  if (borderEdge[edgeI])
821  {
822  const edge& e = surf.edges()[edgeI];
823 
824  if ((borderPoint[e[0]] != -1) && (borderPoint[e[1]] == -1))
825  {
826  startEdgeI = edgeI;
827  startPointI = e[0];
828 
829  break;
830  }
831  else if ((borderPoint[e[0]] == -1) && (borderPoint[e[1]] != -1))
832  {
833  startEdgeI = edgeI;
834  startPointI = e[1];
835 
836  break;
837  }
838  }
839  }
840 
841  if (startEdgeI == -1)
842  {
843  Info<< "Cannot find starting point of splitLine\n" << endl;
844 
845  break;
846  }
847 
848  // Pick any face using edge to start from.
849  const labelList& eFaces = surf.edgeFaces()[startEdgeI];
850 
851  label firstFaceI = eFaces[0];
852 
853  // Find second face which is from same surface i.e. has outwards
854  // pointing normal as well (actually bit more complex than this)
855  label secondFaceI = sharedFace(surf, firstFaceI, startEdgeI);
856 
857  Info<< "Starting local walk from:" << endl
858  << " edge :" << startEdgeI << endl
859  << " point:" << startPointI << endl
860  << " face0:" << firstFaceI << endl
861  << " face1:" << secondFaceI << endl
862  << endl;
863 
864  // From face on border edge to edge.
865  Map<label> faceToEdge(2*nBorderEdges);
866  // From face connected to border point (but not border edge) to point.
867  Map<label> faceToPoint(2*nBorderPoints);
868 
869  faceToEdge.insert(firstFaceI, startEdgeI);
870 
871  walkSplitLine
872  (
873  surf,
874  borderEdge,
875  borderPoint,
876 
877  firstFaceI,
878  startEdgeI,
879  startPointI,
880 
881  faceToEdge,
882  faceToPoint
883  );
884 
885  faceToEdge.insert(secondFaceI, startEdgeI);
886 
887  walkSplitLine
888  (
889  surf,
890  borderEdge,
891  borderPoint,
892 
893  secondFaceI,
894  startEdgeI,
895  startPointI,
896 
897  faceToEdge,
898  faceToPoint
899  );
900 
901  Info<< "Finished local walk and visited" << nl
902  << " border edges :" << faceToEdge.size() << nl
903  << " border points(but not edges):" << faceToPoint.size() << nl
904  << endl;
905 
906  if (debug)
907  {
908  dumpFaces("faceToEdge.obj", surf, faceToEdge);
909  dumpFaces("faceToPoint.obj", surf, faceToPoint);
910  }
911 
912  //
913  // Create coordinates for borderPoints by duplicating the existing
914  // point and then slightly shifting it inwards. To determine the
915  // inwards direction get the average normal of both connectedFaces on
916  // the edge and then interpolate this to the (border)point.
917  //
918 
919  vectorField borderPointVec(surf.nPoints(), vector(GREAT, GREAT, GREAT));
920 
921  calcPointVecs(surf, faceToEdge, faceToPoint, borderPointVec);
922 
923 
924  // New position. Start off from copy of old points.
925  pointField newPoints(surf.localPoints());
926  newPoints.setSize(newPoints.size() + nBorderPoints);
927 
928  forAll(borderPoint, pointI)
929  {
930  label newPointI = borderPoint[pointI];
931 
932  if (newPointI != -1)
933  {
934  scalar minLen = minEdgeLen(surf, pointI);
935 
936  vector n = borderPointVec[pointI];
937  n /= mag(n);
938 
939  newPoints[newPointI] = newPoints[pointI] + 0.1 * minLen * n;
940  }
941  }
942 
943 
944  //
945  // Renumber all faces in connectedFaces
946  //
947 
948  // Start off from copy of faces.
949  List<labelledTri> newTris(surf.size());
950 
951  forAll(surf, faceI)
952  {
953  newTris[faceI] = surf.localFaces()[faceI];
954 
955  newTris[faceI].region() = surf[faceI].region();
956  }
957 
958  // Renumber all faces in faceToEdge
959  renumberFaces(surf, borderPoint, faceToEdge, newTris);
960  // Renumber all faces in faceToPoint
961  renumberFaces(surf, borderPoint, faceToPoint, newTris);
962 
963 
964  // Check if faces use unmoved points.
965  forAll(newTris, faceI)
966  {
967  const labelledTri& f = newTris[faceI];
968 
969  forAll(f, fp)
970  {
971  const point& pt = newPoints[f[fp]];
972 
973  if (mag(pt) >= GREAT/2)
974  {
975  Info<< "newTri:" << faceI << " verts:" << f
976  << " vert:" << f[fp] << " point:" << pt << endl;
977  }
978  }
979  }
980 
981 
982  surf = triSurface(newTris, surf.patches(), newPoints);
983 
984  if (debug || (iteration != 0 && (iteration % 20) == 0))
985  {
986  Info<< "Writing surface to " << outSurfName << nl << endl;
987 
988  surf.write(outSurfName);
989 
990  Info<< "Finished writing surface to " << outSurfName << nl << endl;
991  }
992 
993  // Save prev iteration values.
994  nOldBorderEdges = nBorderEdges;
995  nOldBorderPoints = nBorderPoints;
996 
997  iteration++;
998  }
999  while (true);
1000 
1001  Info<< "Writing final surface to " << outSurfName << nl << endl;
1002 
1003  surf.write(outSurfName);
1004 
1005  Info << "End\n" << endl;
1006 
1007  return 0;
1008 }
1009 
1010 
1011 // ************************ vim: set sw=4 sts=4 et: ************************ //