FreeFOAM The Cross-Platform CFD Toolkit
triSurfaceMesh.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 "triSurfaceMesh.H"
27 #include <OpenFOAM/Random.H>
29 #include <OpenFOAM/EdgeMap.H>
31 #include <OpenFOAM/Time.H>
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 defineTypeNameAndDebug(triSurfaceMesh, 0);
40 addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
41 
42 }
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
48 //Foam::word Foam::triSurfaceMesh::findRawInstance
49 //(
50 // const Time& runTime,
51 // const fileName& dir,
52 // const word& name
53 //)
54 //{
55 // // Check current time first
56 // if (isFile(runTime.path()/runTime.timeName()/dir/name))
57 // {
58 // return runTime.timeName();
59 // }
60 // instantList ts = runTime.times();
61 // label instanceI;
62 //
63 // for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
64 // {
65 // if (ts[instanceI].value() <= runTime.timeOutputValue())
66 // {
67 // break;
68 // }
69 // }
70 //
71 // // continue searching from here
72 // for (; instanceI >= 0; --instanceI)
73 // {
74 // if (isFile(runTime.path()/ts[instanceI].name()/dir/name))
75 // {
76 // return ts[instanceI].name();
77 // }
78 // }
79 //
80 //
81 // // not in any of the time directories, try constant
82 //
83 // // Note. This needs to be a hard-coded constant, rather than the
84 // // constant function of the time, because the latter points to
85 // // the case constant directory in parallel cases
86 //
87 // if (isFile(runTime.path()/runTime.constant()/dir/name))
88 // {
89 // return runTime.constant();
90 // }
91 //
92 // FatalErrorIn
93 // (
94 // "searchableSurfaces::findRawInstance"
95 // "(const Time&, const fileName&, const word&)"
96 // ) << "Cannot find file \"" << name << "\" in directory "
97 // << runTime.constant()/dir
98 // << exit(FatalError);
99 //
100 // return runTime.constant();
101 //}
102 
103 
104 //- Check file existence
105 const Foam::fileName& Foam::triSurfaceMesh::checkFile
106 (
107  const fileName& fName,
108  const fileName& objectName
109 )
110 {
111  if (fName.empty())
112  {
114  (
115  "triSurfaceMesh::checkFile(const fileName&, const fileName&)"
116  ) << "Cannot find triSurfaceMesh starting from "
117  << objectName << exit(FatalError);
118  }
119  return fName;
120 }
121 
122 
123 bool Foam::triSurfaceMesh::addFaceToEdge
124 (
125  const edge& e,
126  EdgeMap<label>& facesPerEdge
127 )
128 {
129  EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
130  if (eFnd != facesPerEdge.end())
131  {
132  if (eFnd() == 2)
133  {
134  return false;
135  }
136  eFnd()++;
137  }
138  else
139  {
140  facesPerEdge.insert(e, 1);
141  }
142  return true;
143 }
144 
145 
146 bool Foam::triSurfaceMesh::isSurfaceClosed() const
147 {
148  // Construct pointFaces. Let's hope surface has compact point
149  // numbering ...
150  labelListList pointFaces;
151  invertManyToMany(points().size(), *this, pointFaces);
152 
153  // Loop over all faces surrounding point. Count edges emanating from point.
154  // Every edge should be used by two faces exactly.
155  // To prevent doing work twice per edge only look at edges to higher
156  // point
157  EdgeMap<label> facesPerEdge(100);
158  forAll(pointFaces, pointI)
159  {
160  const labelList& pFaces = pointFaces[pointI];
161 
162  facesPerEdge.clear();
163  forAll(pFaces, i)
164  {
165  const labelledTri& f = triSurface::operator[](pFaces[i]);
166  label fp = findIndex(f, pointI);
167 
168  // Something weird: if I expand the code of addFaceToEdge in both
169  // below instances it gives a segmentation violation on some
170  // surfaces. Compiler (4.3.2) problem?
171 
172 
173  // Forward edge
174  label nextPointI = f[f.fcIndex(fp)];
175 
176  if (nextPointI > pointI)
177  {
178  bool okFace = addFaceToEdge
179  (
180  edge(pointI, nextPointI),
181  facesPerEdge
182  );
183 
184  if (!okFace)
185  {
186  return false;
187  }
188  }
189  // Reverse edge
190  label prevPointI = f[f.rcIndex(fp)];
191 
192  if (prevPointI > pointI)
193  {
194  bool okFace = addFaceToEdge
195  (
196  edge(pointI, prevPointI),
197  facesPerEdge
198  );
199 
200  if (!okFace)
201  {
202  return false;
203  }
204  }
205  }
206 
207  // Check for any edges used only once.
208  forAllConstIter(EdgeMap<label>, facesPerEdge, iter)
209  {
210  if (iter() != 2)
211  {
212  return false;
213  }
214  }
215  }
216 
217  return true;
218 }
219 
220 
221 // Gets all intersections after initial one. Adds smallVec and starts tracking
222 // from there.
223 void Foam::triSurfaceMesh::getNextIntersections
224 (
225  const indexedOctree<treeDataTriSurface>& octree,
226  const point& start,
227  const point& end,
228  const vector& smallVec,
229  DynamicList<pointIndexHit, 1, 1>& hits
230 )
231 {
232  const vector dirVec(end-start);
233  const scalar magSqrDirVec(magSqr(dirVec));
234 
235  // Initial perturbation amount
236  vector perturbVec(smallVec);
237 
238  while (true)
239  {
240  // Start tracking from last hit.
241  point pt = hits[hits.size()-1].hitPoint() + perturbVec;
242 
243  if (((pt-start)&dirVec) > magSqrDirVec)
244  {
245  return;
246  }
247 
248  // See if any intersection between pt and end
249  pointIndexHit inter = octree.findLine(pt, end);
250 
251  if (!inter.hit())
252  {
253  return;
254  }
255 
256  // Check if already found this intersection
257  bool duplicateHit = false;
258  forAllReverse(hits, i)
259  {
260  if (hits[i].index() == inter.index())
261  {
262  duplicateHit = true;
263  break;
264  }
265  }
266 
267 
268  if (duplicateHit)
269  {
270  // Hit same triangle again. Increase perturbVec and try again.
271  perturbVec *= 2;
272  }
273  else
274  {
275  // Proper hit
276  hits.append(inter);
277  // Restore perturbVec
278  perturbVec = smallVec;
279  }
280  }
281 }
282 
283 
285 {
286  // Unfortunately nPoints constructs meshPoints() so do compact version
287  // ourselves
288 
289  const triSurface& s = static_cast<const triSurface&>(*this);
290 
291  PackedBoolList pointIsUsed(points().size());
292 
293  nPoints = 0;
295 
296  forAll(s, triI)
297  {
298  const labelledTri& f = s[triI];
299 
300  forAll(f, fp)
301  {
302  label pointI = f[fp];
303  if (pointIsUsed.set(pointI, 1u))
304  {
305  bb.min() = ::Foam::min(bb.min(), points()[pointI]);
306  bb.max() = ::Foam::max(bb.max(), points()[pointI]);
307  nPoints++;
308  }
309  }
310  }
311 }
312 
313 
314 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
315 
316 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
317 :
318  searchableSurface(io),
320  (
321  IOobject
322  (
323  io.name(),
324  io.instance(),
325  io.local(),
326  io.db(),
327  io.readOpt(),
328  io.writeOpt(),
329  false // searchableSurface already registered under name
330  )
331  ),
332  triSurface(s),
333  tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
334  maxTreeDepth_(10),
335  surfaceClosed_(-1)
336 {}
337 
338 
339 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
340 :
341  // Find instance for triSurfaceMesh
342  searchableSurface(io),
343  //(
344  // IOobject
345  // (
346  // io.name(),
347  // io.time().findInstance(io.local(), word::null),
348  // io.local(),
349  // io.db(),
350  // io.readOpt(),
351  // io.writeOpt(),
352  // io.registerObject()
353  // )
354  //),
355  // Reused found instance in objectRegistry
357  (
358  IOobject
359  (
360  io.name(),
361  static_cast<const searchableSurface&>(*this).instance(),
362  io.local(),
363  io.db(),
364  io.readOpt(),
365  io.writeOpt(),
366  false // searchableSurface already registered under name
367  )
368  ),
369  triSurface
370  (
371  checkFile
372  (
373  searchableSurface::filePath(),
374  searchableSurface::objectPath()
375  )
376  ),
377  tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
378  maxTreeDepth_(10),
379  surfaceClosed_(-1)
380 {}
381 
382 
383 Foam::triSurfaceMesh::triSurfaceMesh
384 (
385  const IOobject& io,
386  const dictionary& dict
387 )
388 :
389  searchableSurface(io),
390  //(
391  // IOobject
392  // (
393  // io.name(),
394  // io.time().findInstance(io.local(), word::null),
395  // io.local(),
396  // io.db(),
397  // io.readOpt(),
398  // io.writeOpt(),
399  // io.registerObject()
400  // )
401  //),
402  // Reused found instance in objectRegistry
404  (
405  IOobject
406  (
407  io.name(),
408  static_cast<const searchableSurface&>(*this).instance(),
409  io.local(),
410  io.db(),
411  io.readOpt(),
412  io.writeOpt(),
413  false // searchableSurface already registered under name
414  )
415  ),
416  triSurface
417  (
418  checkFile
419  (
422  )
423  ),
425  maxTreeDepth_(10),
426  surfaceClosed_(-1)
427 {
428  scalar scaleFactor = 0;
429 
430  // allow rescaling of the surface points
431  // eg, CAD geometries are often done in millimeters
432  if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0)
433  {
434  Info<< searchableSurface::name() << " : using scale " << scaleFactor
435  << endl;
436  triSurface::scalePoints(scaleFactor);
437  }
438 
439 
440  // Have optional non-standard search tolerance for gappy surfaces.
441  if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
442  {
443  Info<< searchableSurface::name() << " : using intersection tolerance "
444  << tolerance_ << endl;
445  }
446 
447 
448  // Have optional non-standard tree-depth to limit storage.
449  if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
450  {
451  Info<< searchableSurface::name() << " : using maximum tree depth "
452  << maxTreeDepth_ << endl;
453  }
454 }
455 
456 
457 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
458 
460 {
461  clearOut();
462 }
463 
464 
466 {
467  tree_.clear();
468  edgeTree_.clear();
470 }
471 
472 
473 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
474 
476 {
477  // Use copy to calculate face centres so they don't get stored
479  (
482  ).faceCentres();
483 }
484 
485 
487 {
488  tree_.clear();
489  edgeTree_.clear();
490  triSurface::movePoints(newPoints);
491 }
492 
493 
496 {
497  if (tree_.empty())
498  {
499  // Calculate bb without constructing local point numbering.
500  treeBoundBox bb;
501  label nPoints;
502  calcBounds(bb, nPoints);
503 
504  if (nPoints != points().size())
505  {
506  WarningIn("triSurfaceMesh::tree() const")
507  << "Surface " << searchableSurface::name()
508  << " does not have compact point numbering."
509  << " Of " << points().size() << " only " << nPoints
510  << " are used. This might give problems in some routines."
511  << endl;
512  }
513 
514 
515  // Random number generator. Bit dodgy since not exactly random ;-)
516  Random rndGen(65431);
517 
518  // Slightly extended bb. Slightly off-centred just so on symmetric
519  // geometry there are less face/edge aligned items.
520  bb = bb.extend(rndGen, 1E-4);
521  bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
522  bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
523 
526 
527  tree_.reset
528  (
530  (
531  treeDataTriSurface(*this),
532  bb,
533  maxTreeDepth_, // maxLevel
534  10, // leafsize
535  3.0 // duplicity
536  )
537  );
538 
540  }
541 
542  return tree_();
543 }
544 
545 
548 {
549  if (edgeTree_.empty())
550  {
551  // Boundary edges
552  labelList bEdges
553  (
554  identity
555  (
556  nEdges()
557  -nInternalEdges()
558  )
559  + nInternalEdges()
560  );
561 
562  treeBoundBox bb;
563  label nPoints;
564  calcBounds(bb, nPoints);
565 
566  // Random number generator. Bit dodgy since not exactly random ;-)
567  Random rndGen(65431);
568 
569  // Slightly extended bb. Slightly off-centred just so on symmetric
570  // geometry there are less face/edge aligned items.
571 
572  bb = bb.extend(rndGen, 1E-4);
573  bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
574  bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
575 
576  edgeTree_.reset
577  (
579  (
581  (
582  false, // cachebb
583  edges(), // edges
584  localPoints(), // points
585  bEdges // selected edges
586  ),
587  bb, // bb
588  maxTreeDepth_, // maxLevel
589  10, // leafsize
590  3.0 // duplicity
591  )
592  );
593  }
594  return edgeTree_();
595 }
596 
597 
599 {
600  if (regions_.empty())
601  {
602  regions_.setSize(patches().size());
603  forAll(regions_, regionI)
604  {
605  regions_[regionI] = patches()[regionI].name();
606  }
607  }
608  return regions_;
609 }
610 
611 
612 // Find out if surface is closed.
614 {
615  if (surfaceClosed_ == -1)
616  {
617  if (isSurfaceClosed())
618  {
619  surfaceClosed_ = 1;
620  }
621  else
622  {
623  surfaceClosed_ = 0;
624  }
625  }
626 
627  return surfaceClosed_ == 1;
628 }
629 
630 
632 (
633  const pointField& samples,
634  const scalarField& nearestDistSqr,
635  List<pointIndexHit>& info
636 ) const
637 {
638  const indexedOctree<treeDataTriSurface>& octree = tree();
639 
640  info.setSize(samples.size());
641 
644 
645  forAll(samples, i)
646  {
647  static_cast<pointIndexHit&>(info[i]) = octree.findNearest
648  (
649  samples[i],
650  nearestDistSqr[i]
651  );
652  }
653 
655 }
656 
657 
659 (
660  const pointField& start,
661  const pointField& end,
662  List<pointIndexHit>& info
663 ) const
664 {
665  const indexedOctree<treeDataTriSurface>& octree = tree();
666 
667  info.setSize(start.size());
668 
671 
672  forAll(start, i)
673  {
674  static_cast<pointIndexHit&>(info[i]) = octree.findLine
675  (
676  start[i],
677  end[i]
678  );
679  }
680 
682 }
683 
684 
686 (
687  const pointField& start,
688  const pointField& end,
689  List<pointIndexHit>& info
690 ) const
691 {
692  const indexedOctree<treeDataTriSurface>& octree = tree();
693 
694  info.setSize(start.size());
695 
698 
699  forAll(start, i)
700  {
701  static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
702  (
703  start[i],
704  end[i]
705  );
706  }
707 
709 }
710 
711 
713 (
714  const pointField& start,
715  const pointField& end,
716  List<List<pointIndexHit> >& info
717 ) const
718 {
719  const indexedOctree<treeDataTriSurface>& octree = tree();
720 
721  info.setSize(start.size());
722 
725 
726  // Work array
728 
729  // Tolerances:
730  // To find all intersections we add a small vector to the last intersection
731  // This is chosen such that
732  // - it is significant (SMALL is smallest representative relative tolerance;
733  // we need something bigger since we're doing calculations)
734  // - if the start-end vector is zero we still progress
735  const vectorField dirVec(end-start);
736  const scalarField magSqrDirVec(magSqr(dirVec));
737  const vectorField smallVec
738  (
740  + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
741  );
742 
743  forAll(start, pointI)
744  {
745  // See if any intersection between pt and end
746  pointIndexHit inter = octree.findLine(start[pointI], end[pointI]);
747 
748  if (inter.hit())
749  {
750  hits.clear();
751  hits.append(inter);
752 
753  getNextIntersections
754  (
755  octree,
756  start[pointI],
757  end[pointI],
758  smallVec[pointI],
759  hits
760  );
761 
762  info[pointI].transfer(hits);
763  }
764  else
765  {
766  info[pointI].clear();
767  }
768  }
769 
771 }
772 
773 
775 (
776  const List<pointIndexHit>& info,
777  labelList& region
778 ) const
779 {
780  region.setSize(info.size());
781  forAll(info, i)
782  {
783  if (info[i].hit())
784  {
785  region[i] = triSurface::operator[](info[i].index()).region();
786  }
787  else
788  {
789  region[i] = -1;
790  }
791  }
792 }
793 
794 
796 (
797  const List<pointIndexHit>& info,
799 ) const
800 {
801  normal.setSize(info.size());
802 
803  forAll(info, i)
804  {
805  if (info[i].hit())
806  {
807  label triI = info[i].index();
808  //- Cached:
809  //normal[i] = faceNormals()[triI];
810 
811  //- Uncached
812  normal[i] = triSurface::operator[](triI).normal(points());
813  normal[i] /= mag(normal[i]) + VSMALL;
814  }
815  else
816  {
817  // Set to what?
818  normal[i] = vector::zero;
819  }
820  }
821 }
822 
823 
825 {
827  (
829  (
830  IOobject
831  (
832  "values",
833  objectRegistry::time().timeName(), // instance
834  "triSurface", // local
835  *this,
838  ),
839  *this,
840  dimless,
841  labelField(values)
842  )
843  );
844 
845  // Store field on triMesh
846  fldPtr.ptr()->store();
847 }
848 
849 
851 (
852  const List<pointIndexHit>& info,
853  labelList& values
854 ) const
855 {
856  if (foundObject<triSurfaceLabelField>("values"))
857  {
858  values.setSize(info.size());
859 
860  const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
861  (
862  "values"
863  );
864 
865  forAll(info, i)
866  {
867  if (info[i].hit())
868  {
869  values[i] = fld[info[i].index()];
870  }
871  }
872  }
873 }
874 
875 
877 (
878  const pointField& points,
879  List<volumeType>& volType
880 ) const
881 {
882  volType.setSize(points.size());
883 
886 
887  forAll(points, pointI)
888  {
889  const point& pt = points[pointI];
890 
891  // - use cached volume type per each tree node
892  // - cheat conversion since same values
893  volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
894  }
895 
897 }
898 
899 
900 //- Write using given format, version and compression
902 (
906 ) const
907 {
909 
910  if (!mkDir(fullPath.path()))
911  {
912  return false;
913  }
914 
915  triSurface::write(fullPath);
916 
917  if (!isFile(fullPath))
918  {
919  return false;
920  }
921 
922  //return objectRegistry::writeObject(fmt, ver, cmp);
923  return true;
924 }
925 
926 
927 // ************************ vim: set sw=4 sts=4 et: ************************ //