FreeFOAM The Cross-Platform CFD Toolkit
surfaceCheck.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  surfaceCheck
26 
27 Description
28  Performs various checks on surface.
29 
30 Usage
31 
32  - surfaceCheck [OPTIONS] <Foam surface file>
33 
34  @param <Foam surface file> \n
35  @todo Detailed description of argument.
36 
37  @param -noSelfIntersection \n
38  Check for self-intersection.
39 
40  @param -case <dir>\n
41  Case directory.
42 
43  @param -help \n
44  Display help message.
45 
46  @param -doc \n
47  Display Doxygen API documentation page for this application.
48 
49  @param -srcDoc \n
50  Display Doxygen source documentation page for this application.
51 
52 \*---------------------------------------------------------------------------*/
53 
54 #include <OpenFOAM/triangle.H>
55 #include <triSurface/triSurface.H>
58 #include <OpenFOAM/argList.H>
59 #include <OpenFOAM/OFstream.H>
61 #include <OpenFOAM/SortableList.H>
62 #include <OpenFOAM/PatchTools.H>
63 
64 using namespace Foam;
65 
66 // Does face use valid vertices?
67 bool validTri(const bool verbose, const triSurface& surf, const label faceI)
68 {
69  // Simple check on indices ok.
70 
71  const labelledTri& f = surf[faceI];
72 
73  if
74  (
75  (f[0] < 0) || (f[0] >= surf.points().size())
76  || (f[1] < 0) || (f[1] >= surf.points().size())
77  || (f[2] < 0) || (f[2] >= surf.points().size())
78  )
79  {
80  WarningIn("validTri(const triSurface&, const label)")
81  << "triangle " << faceI << " vertices " << f
82  << " uses point indices outside point range 0.."
83  << surf.points().size()-1 << endl;
84 
85  return false;
86  }
87 
88  if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
89  {
90  WarningIn("validTri(const triSurface&, const label)")
91  << "triangle " << faceI
92  << " uses non-unique vertices " << f
93  << " coords:" << f.points(surf.points())
94  << endl;
95  return false;
96  }
97 
98  // duplicate triangle check
99 
100  const labelList& fFaces = surf.faceFaces()[faceI];
101 
102  // Check if faceNeighbours use same points as this face.
103  // Note: discards normal information - sides of baffle are merged.
104  forAll(fFaces, i)
105  {
106  label nbrFaceI = fFaces[i];
107 
108  if (nbrFaceI <= faceI)
109  {
110  // lower numbered faces already checked
111  continue;
112  }
113 
114  const labelledTri& nbrF = surf[nbrFaceI];
115 
116  if
117  (
118  ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
119  && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
120  && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
121  )
122  {
123  WarningIn("validTri(const triSurface&, const label)")
124  << "triangle " << faceI << " vertices " << f
125  << " has the same vertices as triangle " << nbrFaceI
126  << " vertices " << nbrF
127  << " coords:" << f.points(surf.points())
128  << endl;
129 
130  return false;
131  }
132  }
133  return true;
134 }
135 
136 
137 labelList countBins
138 (
139  const scalar min,
140  const scalar max,
141  const label nBins,
142  const scalarField& vals
143 )
144 {
145  scalar dist = nBins/(max - min);
146 
147  labelList binCount(nBins, 0);
148 
149  forAll(vals, i)
150  {
151  scalar val = vals[i];
152 
153  label index = -1;
154 
155  if (Foam::mag(val - min) < SMALL)
156  {
157  index = 0;
158  }
159  else if (val >= max - SMALL)
160  {
161  index = nBins - 1;
162  }
163  else
164  {
165  index = label((val - min)*dist);
166 
167  if ((index < 0) || (index >= nBins))
168  {
169  WarningIn
170  (
171  "countBins(const scalar, const scalar, const label"
172  ", const scalarField&)"
173  ) << "value " << val << " at index " << i
174  << " outside range " << min << " .. " << max << endl;
175 
176  if (index < 0)
177  {
178  index = 0;
179  }
180  else
181  {
182  index = nBins - 1;
183  }
184  }
185  }
186  binCount[index]++;
187  }
188 
189  return binCount;
190 }
191 
192 
193 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
194 // Main program:
195 
196 int main(int argc, char *argv[])
197 {
199 
201  argList::validArgs.append("surface file");
202  argList::validOptions.insert("checkSelfIntersection", "");
203  argList::validOptions.insert("verbose", "");
204  argList args(argc, argv);
205 
206  bool checkSelfIntersection = args.optionFound("checkSelfIntersection");
207  bool verbose = args.optionFound("verbose");
208 
209  fileName surfFileName(args.additionalArgs()[0]);
210  Pout<< "Reading surface from " << surfFileName << " ..." << nl << endl;
211 
212 
213  // Read
214  // ~~~~
215 
216  triSurface surf(surfFileName);
217 
218 
219  Pout<< "Statistics:" << endl;
220  surf.writeStats(Pout);
221  Pout<< endl;
222 
223 
224  // Region sizes
225  // ~~~~~~~~~~~~
226 
227  {
228  labelList regionSize(surf.patches().size(), 0);
229 
230  forAll(surf, faceI)
231  {
232  label region = surf[faceI].region();
233 
234  if (region < 0 || region >= regionSize.size())
235  {
237  << "Triangle " << faceI << " vertices " << surf[faceI]
238  << " has region " << region << " which is outside the range"
239  << " of regions 0.." << surf.patches().size()-1
240  << endl;
241  }
242  else
243  {
244  regionSize[region]++;
245  }
246  }
247 
248  Pout<< "Region\tSize" << nl
249  << "------\t----" << nl;
250  forAll(surf.patches(), patchI)
251  {
252  Pout<< surf.patches()[patchI].name() << '\t'
253  << regionSize[patchI] << nl;
254  }
255  Pout<< nl << endl;
256  }
257 
258 
259  // Check triangles
260  // ~~~~~~~~~~~~~~~
261 
262  {
263  DynamicList<label> illegalFaces(surf.size()/100 + 1);
264 
265  forAll(surf, faceI)
266  {
267  if (!validTri(verbose, surf, faceI))
268  {
269  illegalFaces.append(faceI);
270  }
271  }
272 
273  if (illegalFaces.size())
274  {
275  Pout<< "Surface has " << illegalFaces.size()
276  << " illegal triangles." << endl;
277 
278  OFstream str("illegalFaces");
279  Pout<< "Dumping conflicting face labels to " << str.name() << endl
280  << "Paste this into the input for surfaceSubset" << endl;
281  str << illegalFaces;
282  }
283  else
284  {
285  Pout<< "Surface has no illegal triangles." << endl;
286  }
287  Pout<< endl;
288  }
289 
290 
291 
292  // Triangle quality
293  // ~~~~~~~~~~~~~~~~
294 
295  {
296  scalarField triQ(surf.size(), 0);
297  forAll(surf, faceI)
298  {
299  const labelledTri& f = surf[faceI];
300 
301  if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
302  {
303  //WarningIn(args.executable())
304  // << "Illegal triangle " << faceI << " vertices " << f
305  // << " coords " << f.points(surf.points()) << endl;
306  }
307  else
308  {
309  triPointRef tri
310  (
311  surf.points()[f[0]],
312  surf.points()[f[1]],
313  surf.points()[f[2]]
314  );
315 
316  vector ba(tri.b() - tri.a());
317  ba /= mag(ba) + VSMALL;
318 
319  vector ca(tri.c() - tri.a());
320  ca /= mag(ca) + VSMALL;
321 
322  if (mag(ba&ca) > 1-1E-3)
323  {
324  triQ[faceI] = SMALL;
325  }
326  else
327  {
328  triQ[faceI] = triPointRef
329  (
330  surf.points()[f[0]],
331  surf.points()[f[1]],
332  surf.points()[f[2]]
333  ).quality();
334  }
335  }
336  }
337 
338  labelList binCount = countBins(0, 1, 20, triQ);
339 
340  Pout<< "Triangle quality (equilateral=1, collapsed=0):"
341  << endl;
342 
343 
344  OSstream& os = Pout;
345  os.width(4);
346 
347  scalar dist = (1.0 - 0.0)/20.0;
348  scalar min = 0;
349  forAll(binCount, binI)
350  {
351  Pout<< " " << min << " .. " << min+dist << " : "
352  << 1.0/surf.size() * binCount[binI]
353  << endl;
354  min += dist;
355  }
356  Pout<< endl;
357 
358  label minIndex = findMin(triQ);
359  label maxIndex = findMax(triQ);
360 
361  Pout<< " min " << triQ[minIndex] << " for triangle " << minIndex
362  << nl
363  << " max " << triQ[maxIndex] << " for triangle " << maxIndex
364  << nl
365  << endl;
366 
367 
368  if (triQ[minIndex] < SMALL)
369  {
370  WarningIn(args.executable()) << "Minimum triangle quality is "
371  << triQ[minIndex] << ". This might give problems in"
372  << " self-intersection testing later on." << endl;
373  }
374 
375  // Dump for subsetting
376  {
377  DynamicList<label> problemFaces(surf.size()/100+1);
378 
379  forAll(triQ, faceI)
380  {
381  if (triQ[faceI] < 1E-11)
382  {
383  problemFaces.append(faceI);
384  }
385  }
386  OFstream str("badFaces");
387 
388  Pout<< "Dumping bad quality faces to " << str.name() << endl
389  << "Paste this into the input for surfaceSubset" << nl
390  << nl << endl;
391 
392  str << problemFaces;
393  }
394  }
395 
396 
397 
398  // Edges
399  // ~~~~~
400  {
401  const edgeList& edges = surf.edges();
402  const pointField& localPoints = surf.localPoints();
403 
404  scalarField edgeMag(edges.size());
405 
406  forAll(edges, edgeI)
407  {
408  edgeMag[edgeI] = edges[edgeI].mag(localPoints);
409  }
410 
411  label minEdgeI = findMin(edgeMag);
412  label maxEdgeI = findMax(edgeMag);
413 
414  const edge& minE = edges[minEdgeI];
415  const edge& maxE = edges[maxEdgeI];
416 
417 
418  Pout<< "Edges:" << nl
419  << " min " << edgeMag[minEdgeI] << " for edge " << minEdgeI
420  << " points " << localPoints[minE[0]] << localPoints[minE[1]]
421  << nl
422  << " max " << edgeMag[maxEdgeI] << " for edge " << maxEdgeI
423  << " points " << localPoints[maxE[0]] << localPoints[maxE[1]]
424  << nl
425  << endl;
426  }
427 
428 
429 
430  // Close points
431  // ~~~~~~~~~~~~
432  {
433  const edgeList& edges = surf.edges();
434  const pointField& localPoints = surf.localPoints();
435 
436  const boundBox bb(localPoints);
437  scalar smallDim = 1E-6 * bb.mag();
438 
439  Pout<< "Checking for points less than 1E-6 of bounding box ("
440  << bb.span() << " meter) apart."
441  << endl;
442 
443  // Sort points
444  SortableList<scalar> sortedMag(mag(localPoints));
445 
446  label nClose = 0;
447 
448  for (label i = 1; i < sortedMag.size(); i++)
449  {
450  label ptI = sortedMag.indices()[i];
451 
452  label prevPtI = sortedMag.indices()[i-1];
453 
454  if (mag(localPoints[ptI] - localPoints[prevPtI]) < smallDim)
455  {
456  // Check if neighbours.
457  const labelList& pEdges = surf.pointEdges()[ptI];
458 
459  label edgeI = -1;
460 
461  forAll(pEdges, i)
462  {
463  const edge& e = edges[pEdges[i]];
464 
465  if (e[0] == prevPtI || e[1] == prevPtI)
466  {
467  // point1 and point0 are connected through edge.
468  edgeI = pEdges[i];
469 
470  break;
471  }
472  }
473 
474  nClose++;
475 
476  if (edgeI == -1)
477  {
478  Pout<< " close unconnected points "
479  << ptI << ' ' << localPoints[ptI]
480  << " and " << prevPtI << ' '
481  << localPoints[prevPtI]
482  << " distance:"
483  << mag(localPoints[ptI] - localPoints[prevPtI])
484  << endl;
485  }
486  else
487  {
488  Pout<< " small edge between points "
489  << ptI << ' ' << localPoints[ptI]
490  << " and " << prevPtI << ' '
491  << localPoints[prevPtI]
492  << " distance:"
493  << mag(localPoints[ptI] - localPoints[prevPtI])
494  << endl;
495  }
496  }
497  }
498 
499  Pout<< "Found " << nClose << " nearby points." << nl
500  << endl;
501  }
502 
503 
504 
505  // Check manifold
506  // ~~~~~~~~~~~~~~
507 
508  DynamicList<label> problemFaces(surf.size()/100 + 1);
509 
510  const labelListList& eFaces = surf.edgeFaces();
511 
512  label nSingleEdges = 0;
513  forAll(eFaces, edgeI)
514  {
515  const labelList& myFaces = eFaces[edgeI];
516 
517  if (myFaces.size() == 1)
518  {
519  problemFaces.append(myFaces[0]);
520 
521  nSingleEdges++;
522  }
523  }
524 
525  label nMultEdges = 0;
526  forAll(eFaces, edgeI)
527  {
528  const labelList& myFaces = eFaces[edgeI];
529 
530  if (myFaces.size() > 2)
531  {
532  forAll(myFaces, myFaceI)
533  {
534  problemFaces.append(myFaces[myFaceI]);
535  }
536 
537  nMultEdges++;
538  }
539  }
540  problemFaces.shrink();
541 
542  if ((nSingleEdges != 0) || (nMultEdges != 0))
543  {
544  Pout<< "Surface is not closed since not all edges connected to "
545  << "two faces:" << endl
546  << " connected to one face : " << nSingleEdges << endl
547  << " connected to >2 faces : " << nMultEdges << endl;
548 
549  Pout<< "Conflicting face labels:" << problemFaces.size() << endl;
550 
551  OFstream str("problemFaces");
552 
553  Pout<< "Dumping conflicting face labels to " << str.name() << endl
554  << "Paste this into the input for surfaceSubset" << endl;
555 
556  str << problemFaces;
557  }
558  else
559  {
560  Pout<< "Surface is closed. All edges connected to two faces." << endl;
561  }
562  Pout<< endl;
563 
564 
565 
566  // Check singly connected domain
567  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
568 
570  label numZones = surf.markZones(boolList(surf.nEdges(), false), faceZone);
571 
572  Pout<< "Number of unconnected parts : " << numZones << endl;
573 
574  if (numZones > 1)
575  {
576  Pout<< "Splitting surface into parts ..." << endl << endl;
577 
578  fileName surfFileNameBase(surfFileName.name());
579 
580  for(label zone = 0; zone < numZones; zone++)
581  {
582  boolList includeMap(surf.size(), false);
583 
584  forAll(faceZone, faceI)
585  {
586  if (faceZone[faceI] == zone)
587  {
588  includeMap[faceI] = true;
589  }
590  }
591 
592  labelList pointMap;
593  labelList faceMap;
594 
595  triSurface subSurf
596  (
597  surf.subsetMesh
598  (
599  includeMap,
600  pointMap,
601  faceMap
602  )
603  );
604 
605  fileName subFileName
606  (
607  surfFileNameBase.lessExt()
608  + "_"
609  + name(zone)
610  + ".ftr"
611  );
612 
613  Pout<< "writing part " << zone << " size " << subSurf.size()
614  << " to " << subFileName << endl;
615 
616  subSurf.write(subFileName);
617  }
618 
619  return 0;
620  }
621 
622 
623 
624  // Check orientation
625  // ~~~~~~~~~~~~~~~~~
626 
627  labelHashSet borderEdge(surf.size()/1000);
628  PatchTools::checkOrientation(surf, false, &borderEdge);
629 
630  //
631  // Colour all faces into zones using borderEdge
632  //
633  labelList normalZone;
634  label numNormalZones = PatchTools::markZones(surf, borderEdge, normalZone);
635 
636  Pout<< endl
637  << "Number of zones (connected area with consistent normal) : "
638  << numNormalZones << endl;
639 
640  if (numNormalZones > 1)
641  {
642  Pout<< "More than one normal orientation." << endl;
643  }
644  Pout<< endl;
645 
646 
647 
648  // Check self-intersection
649  // ~~~~~~~~~~~~~~~~~~~~~~~
650 
651  if (checkSelfIntersection)
652  {
653  Pout<< "Checking self-intersection." << endl;
654 
655  triSurfaceSearch querySurf(surf);
656  surfaceIntersection inter(querySurf);
657 
658  if (inter.cutEdges().empty() && inter.cutPoints().empty())
659  {
660  Pout<< "Surface is not self-intersecting" << endl;
661  }
662  else
663  {
664  Pout<< "Surface is self-intersecting" << endl;
665  Pout<< "Writing edges of intersection to selfInter.obj" << endl;
666 
667  OFstream intStream("selfInter.obj");
668  forAll(inter.cutPoints(), cutPointI)
669  {
670  const point& pt = inter.cutPoints()[cutPointI];
671 
672  intStream << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
673  << endl;
674  }
675  forAll(inter.cutEdges(), cutEdgeI)
676  {
677  const edge& e = inter.cutEdges()[cutEdgeI];
678 
679  intStream << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
680  }
681  }
682  Pout<< endl;
683  }
684 
685 
686  Pout<< "End\n" << endl;
687 
688  return 0;
689 }
690 
691 
692 // ************************ vim: set sw=4 sts=4 et: ************************ //