TriPatchIntersection.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2023-2026 OpenFOAM Foundation
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 "TriPatchIntersection.H"
27 #include "barycentricTensor2D.H"
28 #include "boundSphere.H"
29 #include "cpuTime.H"
30 #include "indexedOctree.H"
31 #include "OFstream.H"
32 #include "treeDataPrimitivePatch.H"
33 #include "triIntersect.H"
34 #include "vtkWritePolyData.H"
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 template<class SrcPatchType, class TgtPatchType>
39 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::checkPatchFace
40 (
41  const label patchFacei,
42  const bool isSrc
43 ) const
44 {
45  #ifdef FULLDEBUG
46 
47  const labelList& patchFaceTris =
48  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
49 
50  forAll(patchFaceTris, patchFaceTrii)
51  {
52  const label trii = patchFaceTris[patchFaceTrii];
53 
54  if (triPoints_[trii] != FixedList<label, 3>({-1, -1, -1}))
55  {
56  forAll(triEdges_[trii], triEdgei)
57  {
58  const label edgei = triEdges_[trii][triEdgei];
59 
60  // Check that there are no duplicate points
61  const edge e = triEdgePoints(trii, triEdgei);
62  if (e[0] == e[1])
63  {
64  FatalErrorInFunction
65  << "Tri #" << trii << "'s tri-edge #" << triEdgei
66  << " (edge #" << triEdges_[trii][triEdgei] << ") has "
67  << "duplicate points " << e << exit(FatalError);
68  }
69 
70  // Check that this edge also references this tri
71  const label edgeTrii = findIndex(edgeTris_[edgei], trii);
72  if (edgeTrii == -1)
73  {
74  FatalErrorInFunction
75  << "Tri #" << trii << " references edge #" << edgei
76  << " but the reverse is not true" << exit(FatalError);
77  }
78 
79  // If there is a connected tri ...
80  const label otherTrii = edgeTris_[edgei][!edgeTrii];
81  if (otherTrii != -1)
82  {
83  // Check that the connected tri also references this edge
84  const label otherTriEdgei =
85  findIndex(triEdges_[otherTrii], edgei);
86  if (otherTriEdgei == -1)
87  {
88  FatalErrorInFunction
89  << "Edge #" << edgei << " references tri #"
90  << otherTrii << " but the reverse is not true"
91  << exit(FatalError);
92  }
93 
94  // Check that the edge is correspondingly oriented in its
95  // two connected tris
96  const edge otherE = triEdgePoints(otherTrii, otherTriEdgei);
97  if (edge::compare(e, otherE) != -1)
98  {
99  FatalErrorInFunction
100  << "Edge #" << edgei << ' ' << edgePoints(edgei)
101  << " is not the same in adjacent tri #" << trii
102  << ' ' << triPoints(trii) << " and tri #"
103  << otherTrii << ' ' << triPoints(otherTrii)
104  << exit(FatalError);
105  }
106  }
107 
108  // Check edge-points and point-edges ... !!!
109 
110  // If this is a front edge ...
111  if (edgeFrontEdges_[edgei] != -1)
112  {
113  // Check that the front arrays are consistent
114  const label frontEdgei = edgeFrontEdges_[edgei];
115  if
116  (
117  frontEdgeEdges_.size() <= frontEdgei
118  || frontEdgeEdges_[frontEdgei] != edgei
119  )
120  {
121  FatalErrorInFunction
122  << "Edge #" << edgei << " is marked as part of the "
123  << "front but is not in the front edge list"
124  << exit(FatalError);
125  }
126 
127  // Check that the edge connects source to target
128  const label trii0 = edgeTris_[edgei][0];
129  const label trii1 = edgeTris_[edgei][1];
130  if
131  (
132  trii0 == -1
133  || trii1 == -1
134  || (triSrcFace_[trii0] == -1) == (triSrcFace_[trii1] == -1)
135  )
136  {
137  FatalErrorInFunction
138  << "Front edge #" << edgei
139  << " does not connect a source-tri to a target-tri"
140  << exit(FatalError);
141  }
142  }
143 
144  // If this is a star edge ... !!!
145  }
146 
147  // If this is a candidate tri ... !!!
148 
149  // If this is a marked tri ... !!!
150  }
151  }
152 
153  #endif
154 }
155 
156 
157 template<class SrcPatchType, class TgtPatchType>
158 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::checkPatchEdge
159 (
160  const label patchEdgei,
161  const bool isSrc
162 ) const
163 {
164  #ifdef FULLDEBUG
165 
166  const labelListList& patchEdgePatchFaces =
167  isSrc ? this->srcPatch_.edgeFaces() : this->tgtPatch_.edgeFaces();
168 
169  const labelList patchFaceis = patchEdgePatchFaces[patchEdgei];
170 
171  forAll(patchFaceis, i)
172  {
173  checkPatchFace(patchFaceis[i], isSrc);
174  }
175 
176  #endif
177 }
178 
179 
180 template<class SrcPatchType, class TgtPatchType>
181 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::checkPatchFaces
182 (
183  const bool isSrc
184 ) const
185 {
186  #ifdef FULLDEBUG
187 
188  const List<DynamicList<label>>& patchFaceTris =
189  isSrc ? srcFaceTris_ : tgtFaceTris_;
190 
191  forAll(patchFaceTris, patchFacei)
192  {
193  checkPatchFace(patchFacei, isSrc);
194  }
195 
196  #endif
197 }
198 
199 
200 template<class SrcPatchType, class TgtPatchType>
201 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::removeEdge
202 (
203  const label edgei
204 )
205 {
206  #ifdef FULLDEBUG
207  if
208  (
209  edgeTris_[edgei] != labelPair(-1, -1)
210  && intersectEdgeFaces_[edgei] != labelPair(-1, -1)
211  )
212  {
213  FatalErrorInFunction
214  << "Attempted to remove edge #" << edgei << " which is still "
215  << "connected to triangles " << edgeTris_[edgei] << " and faces "
216  << intersectEdgeFaces_[edgei]
217  << exit(FatalError);
218  }
219  #endif
220 
221  removedEdges_.append(edgei);
222 
223  if (edgeFrontEdges_[edgei] != -1)
224  {
225  frontEdgeEdges_[edgeFrontEdges_[edgei]] = -1;
226  edgeFrontEdges_[edgei] = -1;
227  }
228 }
229 
230 
231 template<class SrcPatchType, class TgtPatchType>
232 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::removeTri
233 (
234  const label trii
235 )
236 {
237  triPoints_[trii] = triFace(-1, -1, -1);
238 
239  {
240  const bool isSrc = triSrcFace_[trii] != -1;
241 
242  const label patchFacei = isSrc ? triSrcFace_[trii] : triTgtFace_[trii];
243 
244  DynamicList<label>& patchFaceTris =
245  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
246 
247  label patchFaceTrii = patchFaceTris.size() - 1;
248 
249  for (; patchFaceTris[patchFaceTrii] != trii; -- patchFaceTrii);
250 
251  for (; patchFaceTrii < patchFaceTris.size() - 1; ++ patchFaceTrii)
252  {
253  patchFaceTris[patchFaceTrii] = patchFaceTris[patchFaceTrii + 1];
254  }
255 
256  patchFaceTris.resize(patchFaceTris.size() - 1);
257  }
258 
259  triSrcFace_[trii] = -1;
260  triTgtFace_[trii] = -1;
261 
262  if (triCandidateTris_[trii] != -1)
263  {
264  candidateTriTris_[triCandidateTris_[trii]] = -1;
265  triCandidateTris_[trii] = -1;
266  }
267 
268  if (triMarkedTris_[trii] != -1)
269  {
270  markedTriTris_[triMarkedTris_[trii]] = -1;
271  triMarkedTris_[trii] = -1;
272  }
273 
274  forAll(triEdges_[trii], i)
275  {
276  const label edgei = triEdges_[trii][i];
277 
278  edgeTris_[edgei][edgeTris_[edgei][1] == trii] = -1;
279 
280  if
281  (
282  edgeTris_[edgei] == labelPair(-1, -1)
283  && intersectEdgeFaces_[edgei] == labelPair(-1, -1)
284  )
285  {
286  removeEdge(edgei);
287  }
288  }
289 
290  triEdges_[trii] = FixedList<label, 3>({-1, -1, -1});
291 
292  removedTris_.append(trii);
293 }
294 
295 
296 template<class SrcPatchType, class TgtPatchType>
297 Foam::label Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::newTrii()
298 {
299  if (removedTris_.size())
300  {
301  return removedTris_.remove();
302  }
303  else
304  {
305  triPoints_.append(triFace(-1, -1, -1));
306  triEdges_.append(FixedList<label, 3>({-1, -1, -1}));
307  triSrcFace_.append(-1);
308  triTgtFace_.append(-1);
309  triCandidateTris_.append(-1);
310  triMarkedTris_.append(-1);
311  return triPoints_.size() - 1;
312  }
313 }
314 
315 
316 template<class SrcPatchType, class TgtPatchType>
317 Foam::label Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::newEdgei()
318 {
319  if (removedEdges_.size())
320  {
321  return removedEdges_.remove();
322  }
323  else
324  {
325  edgeTris_.append({-1, -1});
326  intersectEdgeFaces_.append({-1, -1});
327  edgeFrontEdges_.append(-1);
328  return edgeTris_.size() - 1;
329  }
330 }
331 
332 
333 template<class SrcPatchType, class TgtPatchType>
334 Foam::label Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::newPointi()
335 {
336  const label pointi = srcPoints_.size();
337 
338  srcPoints_.append(point::uniform(NaN));
339  srcPointNormals_.append(vector::uniform(NaN));
340  tgtPoints_.append(point::uniform(NaN));
341  pointPoints_.append(pointi);
342  this->pointSrcFaces_.append(-1);
343  this->pointTgtFaces_.append(-1);
344  this->pointSrcEdges_.append(-1);
345  this->pointTgtEdges_.append(-1);
346  this->pointSrcPoints_.append(-1);
347  this->pointTgtPoints_.append(-1);
348  return pointi;
349 }
350 
351 
352 template<class SrcPatchType, class TgtPatchType>
353 Foam::label Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triPoint
354 (
355  const label trii,
356  const label triPointi
357 ) const
358 {
359  return pointPoints_[triPoints_[trii][triPointi]];
360 }
361 
362 
363 template<class SrcPatchType, class TgtPatchType>
364 Foam::triFace Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triPoints
365 (
366  const label trii
367 ) const
368 {
369  triFace result;
370  forAll(triPoints_[trii], triPointi)
371  {
372  result[triPointi] = triPoint(trii, triPointi);
373  }
374  return result;
375 }
376 
377 
378 template<class SrcPatchType, class TgtPatchType>
379 template<class Type>
380 Foam::FixedList<Type, 3>
381 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triPointValues
382 (
383  const label trii,
384  const UList<Type> values
385 ) const
386 {
387  FixedList<Type, 3> result;
388  forAll(result, triPointi)
389  {
390  result[triPointi] = values[triPoint(trii, triPointi)];
391  }
392  return result;
393 }
394 
395 
396 template<class SrcPatchType, class TgtPatchType>
397 Foam::FixedList<bool, 3>
398 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triOwns
399 (
400  const label trii
401 ) const
402 {
403  FixedList<bool, 3> result;
404  forAll(triEdges_[trii], triEdgei)
405  {
406  result[triEdgei] = edgeTris_[triEdges_[trii][triEdgei]][0] == trii;
407  }
408  return result;
409 }
410 
411 
412 template<class SrcPatchType, class TgtPatchType>
413 Foam::FixedList<Foam::label, 3>
414 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triOtherTriPoints
415 (
416  const label trii,
417  const label otherTrii
418 ) const
419 {
420  FixedList<label, 3> result({-1, -1, -1});
421  forAll(triPoints_[trii], triPointi)
422  {
423  forAll(triPoints_[otherTrii], otherTriPointi)
424  {
425  if
426  (
427  triPoint(trii, triPointi)
428  == triPoint(otherTrii, otherTriPointi)
429  )
430  {
431  result[triPointi] = otherTriPointi;
432  }
433  }
434  }
435  return result;
436 }
437 
438 
439 template<class SrcPatchType, class TgtPatchType>
440 Foam::FixedList<Foam::label, 3>
441 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triOtherTriEdges
442 (
443  const label trii,
444  const label otherTrii
445 ) const
446 {
447  FixedList<label, 3> result({-1, -1, -1});
448  forAll(triEdges_[trii], triEdgei)
449  {
450  forAll(triEdges_[otherTrii], otherTriEdgei)
451  {
452  if
453  (
454  triEdges_[trii][triEdgei]
455  == triEdges_[otherTrii][otherTriEdgei]
456  )
457  {
458  result[triEdgei] = otherTriEdgei;
459  }
460  }
461  }
462  return result;
463 }
464 
465 
466 template<class SrcPatchType, class TgtPatchType>
467 Foam::edge Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::triEdgePoints
468 (
469  const label trii,
470  const label triEdgei
471 ) const
472 {
473  return edge
474  (
475  triPoint(trii, triEdgei),
476  triPoint(trii, (triEdgei + 1) % 3)
477  );
478 }
479 
480 
481 template<class SrcPatchType, class TgtPatchType>
482 Foam::edge Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::edgePoints
483 (
484  const label edgei
485 ) const
486 {
487  const label edgeSidei = edgeTris_[edgei][0] == -1;
488 
489  const label trii = edgeTris_[edgei][edgeSidei];
490  const label triEdgei = findIndex(triEdges_[trii], edgei);
491 
492  const edge e = triEdgePoints(trii, triEdgei);
493 
494  return edgeSidei == 0 ? e : e.reverseEdge();
495 }
496 
497 
498 template<class SrcPatchType, class TgtPatchType>
499 Foam::label
500 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePoint
501 (
502  const label patchFacei,
503  const label patchFacePointi,
504  const bool isSrc
505 ) const
506 {
507  return
508  isSrc
509  ? this->srcPointPoints_
510  [
511  this->srcPatch_.localFaces()[patchFacei][patchFacePointi]
512  ]
513  : this->tgtPointPoints_
514  [
515  this->tgtPatch_.localFaces()[patchFacei][patchFacePointi]
516  ];
517 }
518 
519 
520 template<class SrcPatchType, class TgtPatchType>
521 Foam::triFace
522 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePoints
523 (
524  const label patchFacei,
525  const bool isSrc
526 ) const
527 {
528  triFace result;
529  forAll(result, patchFacePointi)
530  {
531  result[patchFacePointi] =
532  patchFacePoint(patchFacei, patchFacePointi, isSrc);
533  }
534  return result;
535 }
536 
537 
538 template<class SrcPatchType, class TgtPatchType>
539 template<class Type>
540 Foam::FixedList<Type, 3>
541 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePointValues
542 (
543  const label patchFacei,
544  const bool isSrc,
545  const UList<Type>& values
546 ) const
547 {
548  FixedList<vector, 3> result;
549  forAll(result, patchFacePointi)
550  {
551  result[patchFacePointi] =
552  values[patchFacePoint(patchFacei, patchFacePointi, isSrc)];
553  }
554  return result;
555 }
556 
557 
558 template<class SrcPatchType, class TgtPatchType>
559 Foam::FixedList<bool, 3>
560 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFaceOwns
561 (
562  const label patchFacei,
563  const bool isSrc
564 ) const
565 {
566  const UList<triFace>& localFaces =
567  isSrc ? this->srcPatch_.localFaces() : this->tgtPatch_.localFaces();
568  const labelListList& faceEdges =
569  isSrc ? this->srcPatch_.faceEdges() : this->tgtPatch_.faceEdges();
570  const edgeList& localEdges =
571  isSrc ? this->srcPatch_.edges() : this->tgtPatch_.edges();
572  FixedList<bool, 3> result;
573  forAll(result, i)
574  {
575  const edge& e = localEdges[faceEdges[patchFacei][i]];
576  const edge fe = localFaces[patchFacei].faceEdge(i);
577  result[i] = edge::compare(e, fe) > 0;
578  }
579  return result;
580 }
581 
582 
583 template<class SrcPatchType, class TgtPatchType>
584 Foam::FixedList<Foam::label, 3>
585 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::
586 patchFaceOtherPatchFacePoints
587 (
588  const label patchFacei,
589  const label otherPatchFacei,
590  const bool isSrc
591 ) const
592 {
593  const triFace pointis = patchFacePoints(patchFacei, isSrc);
594  const triFace otherPointis = patchFacePoints(otherPatchFacei, !isSrc);
595  FixedList<label, 3> result({-1, -1, -1});
596  forAll(pointis, i)
597  {
598  forAll(otherPointis, otheri)
599  {
600  if (pointis[i] == otherPointis[otheri])
601  {
602  result[i] = otheri;
603  }
604  }
605  }
606  return result;
607 }
608 
609 
610 template<class SrcPatchType, class TgtPatchType>
611 Foam::label
612 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePatchPoint
613 (
614  const label patchFacei,
615  const label patchFacePointi,
616  const bool isSrc
617 ) const
618 {
619  return
620  isSrc
621  ? this->srcPatch_.localFaces()[patchFacei][patchFacePointi]
622  : this->tgtPatch_.localFaces()[patchFacei][patchFacePointi];
623 }
624 
625 
626 template<class SrcPatchType, class TgtPatchType>
627 Foam::triFace
628 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePatchPoints
629 (
630  const label patchFacei,
631  const bool isSrc
632 ) const
633 {
634  triFace result;
635  forAll(result, i)
636  {
637  result[i] = patchFacePatchPoint(patchFacei, i, isSrc);
638  }
639  return result;
640 }
641 
642 
643 template<class SrcPatchType, class TgtPatchType>
644 Foam::label
645 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePatchEdge
646 (
647  const label patchFacei,
648  const label patchFaceEdgei,
649  const bool isSrc
650 ) const
651 {
652  return
653  isSrc
654  ? this->srcPatch_.faceEdges()[patchFacei][patchFaceEdgei]
655  : this->tgtPatch_.faceEdges()[patchFacei][patchFaceEdgei];
656 }
657 
658 
659 template<class SrcPatchType, class TgtPatchType>
660 Foam::triFace
661 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::patchFacePatchEdges
662 (
663  const label patchFacei,
664  const bool isSrc
665 ) const
666 {
667  triFace result;
668  forAll(result, i)
669  {
670  result[i] = patchFacePatchEdge(patchFacei, i, isSrc);
671  }
672  return result;
673 }
674 
675 
676 template<class SrcPatchType, class TgtPatchType>
677 Foam::label
678 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::edgePatchEdge
679 (
680  const label edgei,
681  const bool isSrc
682 ) const
683 {
684  const labelListList& patchPointPatchEdges =
685  isSrc ? this->srcPatch_.pointEdges() : this->tgtPatch_.pointEdges();
686 
687  const labelList& pointPatchPoints =
688  isSrc ? this->pointSrcPoints_ : this->pointTgtPoints_;
689  const labelList& pointPatchEdges =
690  isSrc ? this->pointSrcEdges_ : this->pointTgtEdges_;
691 
692  auto nPointPatchEdges = [&](const label pointi)
693  {
694  const label patchPointi = pointPatchPoints[pointi];
695 
696  if (patchPointi == -1)
697  {
698  return pointPatchEdges[pointi] == -1 ? label(0) : label(1);
699  }
700  else
701  {
702  return patchPointPatchEdges[patchPointi].size();
703  }
704  };
705 
706  auto pointPatchEdge = [&]
707  (
708  const label pointi,
709  const label i
710  )
711  {
712  const label patchPointi = pointPatchPoints[pointi];
713 
714  if (patchPointi == -1)
715  {
716  return pointPatchEdges[pointi];
717  }
718  else
719  {
720  return patchPointPatchEdges[patchPointi][i];
721  }
722  };
723 
724  const edge& e = edgePoints(edgei);
725  const label pointi0 = e.first(), pointi1 = e.last();
726 
727  label patchEdgei = -1;
728 
729  for (label i0 = 0; i0 < nPointPatchEdges(pointi0); ++ i0)
730  {
731  for (label i1 = 0; i1 < nPointPatchEdges(pointi1); ++ i1)
732  {
733  const label patchEdgei0 = pointPatchEdge(pointi0, i0);
734  const label patchEdgei1 = pointPatchEdge(pointi1, i1);
735 
736  if (patchEdgei0 == patchEdgei1)
737  {
738  if (patchEdgei != -1 && patchEdgei != patchEdgei0)
739  {
740  FatalErrorInFunction
741  << "Edge #" << edgei << " (points " << e
742  << ") is associated with two or more different "
743  << (isSrc ? "source" : "target")
744  << " patch edges. This should not be possible."
745  << exit(FatalError);
746  }
747 
748  patchEdgei = patchEdgei0;
749  }
750  }
751  }
752 
753  return patchEdgei;
754 }
755 
756 
757 template<class SrcPatchType, class TgtPatchType>
758 Foam::labelPair
759 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::edgePatchEdges
760 (
761  const label edgei
762 ) const
763 {
764  return labelPair(edgePatchEdge(edgei, true), edgePatchEdge(edgei, false));
765 }
766 
767 
768 template<class SrcPatchType, class TgtPatchType>
769 Foam::label Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::addTri
770 (
771  const triFace& pointis,
772  const FixedList<label, 3>& edgeis,
773  const label patchFacei,
774  const bool isSrc
775 )
776 {
777  // Check that edge points match the existing mesh
778  #ifdef FULLDEBUG
779  forAll(edgeis, triEdgei)
780  {
781  const label edgei = edgeis[triEdgei];
782 
783  if (edgei == -1) continue;
784 
785  const label otherTrii = edgeTris_[edgei][edgeTris_[edgei][0] == -1];
786 
787  if (otherTrii == -1) continue;
788 
789  const label otherTriEdgei = findIndex(triEdges_[otherTrii], edgei);
790 
791  const edge eThis = pointis.faceEdge(triEdgei);
792  const edge eOther = triEdgePoints(otherTrii, otherTriEdgei);
793 
794  if (edge::compare(eThis, eOther) != -1)
795  {
796  FatalErrorInFunction
797  << "Edge #" << edgei << ' ' << edgePoints(edgei)
798  << " is not the same in the new tri " << pointis << ' '
799  << edgeis << " and the existing adjacent tri "
800  << triPoints(otherTrii) << ' ' << triEdges_[otherTrii]
801  << exit(FatalError);
802  }
803  }
804  #endif
805 
806  // Get the new triangle label
807  const label trii = newTrii();
808 
809  // Set the points
810  triPoints_[trii] = pointis;
811 
812  // Set the edges, creating new ones where specified, and completing the
813  // edge-tri associations
814  triEdges_[trii] = edgeis;
815  forAll(triEdges_[trii], triEdgei)
816  {
817  label& edgei = triEdges_[trii][triEdgei];
818 
819  if (edgei == -1)
820  {
821  edgei = newEdgei();
822  }
823 
824  edgeTris_[edgei][edgeTris_[edgei][0] != -1] = trii;
825  }
826 
827  // Set the patch face association
828  if (isSrc)
829  {
830  triSrcFace_[trii] = patchFacei;
831  srcFaceTris_[patchFacei].append(trii);
832  }
833  else
834  {
835  triTgtFace_[trii] = patchFacei;
836  tgtFaceTris_[patchFacei].append(trii);
837  }
838 
839  return trii;
840 }
841 
842 
843 template<class SrcPatchType, class TgtPatchType>
844 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::flipEdge
845 (
846  const label edgei
847 )
848 {
849  // Return if this edge does not have two adjacent triangles
850  forAll(edgeTris_[edgei], edgeTrii)
851  {
852  const label trii = edgeTris_[edgei][edgeTrii];
853 
854  if (trii == -1) return;
855  }
856 
857  // Check that the flip doesn't break patch face associations
858  #ifdef FULLDEBUG
859  if (edgePatchEdges(edgei) != labelPair(-1, -1))
860  {
861  FatalErrorInFunction
862  << "Flipping an original edge"
863  << exit(FatalError);
864  }
865  #endif
866 
867  // Store stuff that will disappear when the existing triangles are modified
868  label pointi0 = -1;
869  label pointi1 = -1;
870  labelPair pointiOpps(-1, -1);
871  labelPair edgei0s(-1, -1), trii0s(-1, -1);
872  labelPair edgei1s(-1, -1), trii1s(-1, -1);
873  forAll(edgeTris_[edgei], edgeTrii)
874  {
875  const label trii = edgeTris_[edgei][edgeTrii];
876  const label triEdgei = findIndex(triEdges_[trii], edgei);
877 
878  pointi0 = triPoint(trii, (triEdgei + edgeTrii) % 3);
879  pointi1 = triPoint(trii, (triEdgei + !edgeTrii) % 3);
880  pointiOpps[edgeTrii] = triPoint(trii, (triEdgei + 2) % 3);
881 
882  edgei0s[edgeTrii] = triEdges_[trii][(triEdgei + !edgeTrii + 1) % 3];
883  trii0s[edgeTrii] =
884  edgeTris_
885  [edgei0s[edgeTrii]]
886  [edgeTris_[edgei0s[edgeTrii]][0] == trii];
887  edgei1s[edgeTrii] = triEdges_[trii][(triEdgei + edgeTrii + 1) % 3];
888  trii1s[edgeTrii] =
889  edgeTris_
890  [edgei1s[edgeTrii]]
891  [edgeTris_[edgei1s[edgeTrii]][0] == trii];
892  }
893 
894  // Modify the triangles
895  {
896  const label trii0 = edgeTris_[edgei][0], trii1 = edgeTris_[edgei][1];
897 
898  triPoints_[trii0] = {pointi0, pointiOpps[1], pointiOpps[0]};
899  triEdges_[trii0] = {edgei0s[1], edgei, edgei0s[0]};
900  edgeTris_[edgei0s[1]][edgeTris_[edgei0s[1]][1] == trii1] = trii0;
901 
902  triPoints_[trii1] = {pointi1, pointiOpps[0], pointiOpps[1]};
903  triEdges_[trii1] = {edgei1s[0], edgei, edgei1s[1]};
904  edgeTris_[edgei1s[0]][edgeTris_[edgei1s[0]][1] == trii0] = trii1;
905  }
906 }
907 
908 
909 template<class SrcPatchType, class TgtPatchType>
910 Foam::scalar
911 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::circumDistSqr
912 (
913  const label trii,
914  const label pointi
915 ) const
916 {
917  // !!! It may be preferable to do this in the plane of the patch face,
918  // rather than the current triangle, but that would mean implementing
919  // custom circumcircle methods.
920 
921  const pointField& points =
922  triSrcFace_[trii] != -1 ? srcPoints_ : tgtPoints_;
923 
924  const triPointRef t = triPoints(trii).tri(points);
925 
926  const Tuple2<point, scalar> crSqr = t.circumCircleSqr();
927 
928  return
929  crSqr.second() < 0
930  ? - vGreat
931  : magSqr(points[pointi] - crSqr.first()) - crSqr.second();
932 }
933 
934 
935 template<class SrcPatchType, class TgtPatchType>
936 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::insertPoints
937 (
938  const label insertionTriOrEdgei,
939  const bool isTri,
940  const UList<label>& pointis,
941  UList<label>& insertionEdgeis,
942  const UList<label>& fixedEdgeis
943 )
944 {
945  // Clear the insertion edge list
946  insertionEdgeis = -1;
947 
948  // If inserting into a tri then initialise the candidates with the given
949  // insertion tri
950  if (isTri)
951  {
952  candidateTriTris_.append(insertionTriOrEdgei);
953  triCandidateTris_[insertionTriOrEdgei] = 0;
954  }
955 
956  // If inserting into an edge, get the insertion edge label
957  if (!isTri)
958  {
959  insertionEdgeis[0] = insertionTriOrEdgei;
960  }
961 
962  // Check that we are inserting unambiguously into the source or target side
963  #ifdef FULLDEBUG
964  if (!isTri)
965  {
966  const label trii0 = edgeTris_[insertionTriOrEdgei][0];
967  const label trii1 = edgeTris_[insertionTriOrEdgei][1];
968 
969  if
970  (
971  trii0 != -1
972  && trii1 != -1
973  && (triSrcFace_[trii0] == -1) != (triSrcFace_[trii1] == -1)
974  )
975  {
976  FatalErrorInFunction
977  << "Attempted insertion into front edge #"
978  << insertionTriOrEdgei << exit(FatalError);
979  }
980  }
981  #endif
982 
983  // Determine the patch associations for the insertion tri or edge
984  const label adjacentTri =
985  isTri
986  ? insertionTriOrEdgei
987  : edgeTris_[insertionTriOrEdgei][edgeTris_[insertionTriOrEdgei][0] == -1];
988  const bool isSrc = triSrcFace_[adjacentTri] != -1;
989  const label patchEdgei =
990  isTri
991  ? -1
992  : edgePatchEdge(insertionTriOrEdgei, isSrc);
993  const label patchFacei =
994  patchEdgei != -1
995  ? -1
996  : (isSrc ? triSrcFace_[adjacentTri] : triTgtFace_[adjacentTri]);
997 
998  // Get connectivity arrays that will be modified during insertion
999  List<DynamicList<label>>& patchEdgePoints =
1000  isSrc ? this->srcEdgePoints_ : this->tgtEdgePoints_;
1001  DynamicList<label>& pointPatchEdges =
1002  isSrc ? this->pointSrcEdges_ : this->pointTgtEdges_;
1003  DynamicList<label>& pointPatchFaces =
1004  isSrc ? this->pointSrcFaces_ : this->pointTgtFaces_;
1005 
1006  // Determine whether the insertion star is closed or not
1007  const bool isClosedStar =
1008  isTri || findIndex(edgeTris_[insertionEdgeis[0]], -1) == -1;
1009 
1010  // Loop the points
1011  forAll(pointis, pointii)
1012  {
1013  const label pointi = pointis[pointii];
1014 
1015  // Initial star tri or edge. If a tri, this means searching the
1016  // candidates for the one that best contains the point. If an edge then
1017  // just use the next insertion edge.
1018  label initialTriOrEdgei = -1;
1019  if (isTri)
1020  {
1021  label insertionCandidateTrii = -1;
1022  scalar minDistSqr = vGreat;
1023  forAll(candidateTriTris_, candidateTrii)
1024  {
1025  const label trii = candidateTriTris_[candidateTrii];
1026 
1027  if (trii == -1) continue;
1028  const scalar distSqr = circumDistSqr(trii, pointi);
1029  if (distSqr < minDistSqr)
1030  {
1031  insertionCandidateTrii = candidateTrii;
1032  minDistSqr = distSqr;
1033  }
1034  if (minDistSqr < 0)
1035  {
1036  break;
1037  }
1038  }
1039  initialTriOrEdgei = candidateTriTris_[insertionCandidateTrii];
1040  }
1041  else
1042  {
1043  initialTriOrEdgei = insertionEdgeis[pointii];
1044  }
1045 
1046  // Populate the star
1047  star::context starContext = star_.populate
1048  (
1049  initialTriOrEdgei,
1050  isTri,
1051  [&](const label edgei, const label trii)
1052  {
1053  // If this edge is part of a patch edge, or if it connects to a
1054  // triangle in the other patch, or if it is fixed, then it
1055  // should not be removed, so it should not be crossed
1056  if
1057  (
1058  edgePatchEdges(edgei) != labelPair(-1, -1)
1059  || (isSrc && triTgtFace_[trii] != -1)
1060  || (!isSrc && triSrcFace_[trii] != -1)
1061  || findIndex(fixedEdgeis, edgei) != -1
1062  )
1063  {
1064  return false;
1065  }
1066 
1067  // If the destination tri borders any tris that are already in
1068  // the star (other than the tri from which we walked) then do
1069  // not cross into it. Adding this tri would leave a point
1070  // within the star, which is not allowed.
1071  const label triEdgei = findIndex(triEdges_[trii], edgei);
1072  const label edgej0 = triEdges_[trii][(triEdgei + 1) % 3];
1073  const label edgej1 = triEdges_[trii][(triEdgei + 2) % 3];
1074  const label trij0 =
1075  edgeTris_[edgej0][edgeTris_[edgej0][0] == trii];
1076  const label trij1 =
1077  edgeTris_[edgej1][edgeTris_[edgej1][0] == trii];
1078  if
1079  (
1080  (trij0 != -1 && star_.faceStarFaces()[trij0] != -1)
1081  || (trij1 != -1 && star_.faceStarFaces()[trij1] != -1)
1082  || (!isTri && findIndex(insertionEdgeis, edgej0) != -1)
1083  || (!isTri && findIndex(insertionEdgeis, edgej1) != -1)
1084  )
1085  {
1086  return false;
1087  }
1088 
1089  // Otherwise, cross into the triangle if it's circumcircle
1090  // contains the insertion point
1091  return circumDistSqr(trii, pointi) < 0;
1092  },
1093  triEdges_,
1094  edgeTris_
1095  );
1096 
1097  // Get the end points of the current insertion edge
1098  const edge insertionEdge =
1099  isTri ? edge(-1, -1) : edgePoints(insertionEdgeis[pointii]);
1100 
1101  // Walk around the star polygon disconnecting the old triangles and
1102  // adding in the new ones
1103  label newTrii0 = -1, newTrii00 = -1;
1104  forAllStarEdges(star_, i, starEdgei, edgei)
1105  {
1106  const bool isFirst = i == 0;
1107  const bool isLast = i == star_.starEdgeEdges().size() - 1;
1108 
1109  const bool edgeSidei =
1110  edgeTris_[edgei][0] == -1
1111  || star_.faceStarFaces()[edgeTris_[edgei][0]] == -1;
1112 
1113  const label trii = edgeTris_[edgei][edgeSidei];
1114  const label triEdgei = findIndex(triEdges_[trii], edgei);
1115 
1116  // Disconnect the existing triangle from the star
1117  const label tempEdgei = newEdgei();
1118  triEdges_[trii][triEdgei] = tempEdgei;
1119  edgeTris_[tempEdgei][0] = trii;
1120  edgeTris_[edgei][edgeSidei] = -1;
1121 
1122  // Add the new triangle
1123  const label newTrii =
1124  addTri
1125  (
1126  {
1127  pointi,
1128  triPoint(trii, triEdgei),
1129  triPoint(trii, (triEdgei + 1) % 3),
1130  },
1131  {
1132  !isFirst ? triEdges_[newTrii0][2] : -1,
1133  edgei,
1134  isClosedStar && isLast ? triEdges_[newTrii00][0] : -1
1135  },
1136  isSrc ? triSrcFace_[trii] : triTgtFace_[trii],
1137  isSrc
1138  );
1139  newTrii0 = newTrii;
1140  newTrii00 = isFirst ? newTrii0 : newTrii00;
1141 
1142  // Add the new triangle into the list of candidates
1143  if (isTri)
1144  {
1145  candidateTriTris_.append(newTrii);
1146  triCandidateTris_[newTrii] = candidateTriTris_.size() - 1;
1147  }
1148 
1149  // Get the previous insertion edge
1150  if (triPoint(newTrii, 2) == insertionEdge[0])
1151  {
1152  insertionEdgeis[pointii] = triEdges_[newTrii][2];
1153  }
1154 
1155  // Get the next insertion edge
1156  if (!isTri && triPoint(newTrii, 1) == insertionEdge[1])
1157  {
1158  insertionEdgeis[pointii + 1] = triEdges_[newTrii][0];
1159  }
1160  }
1161 
1162  // Create the edge or face association for the new point
1163  if (patchEdgei != -1)
1164  {
1165  patchEdgePoints[patchEdgei].append(-1);
1166  label patchEdgePointi = patchEdgePoints[patchEdgei].size() - 1;
1167  for (; patchEdgePointi >= 0; -- patchEdgePointi)
1168  {
1169  label& pointi = patchEdgePoints[patchEdgei][patchEdgePointi];
1170  pointi = patchEdgePoints[patchEdgei][patchEdgePointi - 1];
1171  if
1172  (
1173  pointPoints_[pointi] == insertionEdge[0]
1174  || pointPoints_[pointi] == insertionEdge[1]
1175  )
1176  {
1177  break;
1178  }
1179  }
1180  -- patchEdgePointi;
1181  patchEdgePoints[patchEdgei][patchEdgePointi] = pointi;
1182 
1183  pointPatchEdges[pointi] = patchEdgei;
1184  }
1185  else
1186  {
1187  pointPatchFaces[pointi] = patchFacei;
1188  }
1189 
1190  // Check
1191  if (patchEdgei != -1)
1192  {
1193  checkPatchEdge(patchEdgei, isSrc);
1194  }
1195  else
1196  {
1197  checkPatchFace(patchFacei, isSrc);
1198  }
1199 
1200  // If the insertion edge is not the same orientation as the original
1201  // insertion edge then reverse it
1202  if (!isTri)
1203  {
1204  if (edgePoints(insertionEdgeis[pointii + 1])[1] != insertionEdge[1])
1205  {
1206  Swap
1207  (
1208  edgeTris_[insertionEdgeis[pointii + 1]][0],
1209  edgeTris_[insertionEdgeis[pointii + 1]][1]
1210  );
1211  }
1212  }
1213 
1214  // Remove the star triangles
1215  forAllStarFaces(star_, starTrii, trii)
1216  {
1217  removeTri(trii);
1218  }
1219 
1220  // Check
1221  if (patchEdgei != -1)
1222  {
1223  checkPatchEdge(patchEdgei, isSrc);
1224  }
1225  else
1226  {
1227  checkPatchFace(patchFacei, isSrc);
1228  }
1229  }
1230 
1231  // Clear the candidate tris
1232  if (isTri)
1233  {
1234  forAll(candidateTriTris_, candidateTrii)
1235  {
1236  const label trii = candidateTriTris_[candidateTrii];
1237  if (trii != -1)
1238  {
1239  triCandidateTris_[trii] = -1;
1240  }
1241  }
1242  candidateTriTris_.clear();
1243  }
1244 }
1245 
1246 
1247 template<class SrcPatchType, class TgtPatchType>
1248 bool Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::pointCanIntersect
1249 (
1250  const label pointi
1251 ) const
1252 {
1253  return
1254  (
1255  this->pointSrcPoints_[pointi] != -1
1256  && this->pointTgtPoints_[pointi] == -1
1257  && this->pointTgtEdges_[pointi] == -1
1258  && this->pointTgtFaces_[pointi] == -1
1259  )
1260  || (
1261  this->pointTgtPoints_[pointi] != -1
1262  && this->pointSrcPoints_[pointi] == -1
1263  && this->pointSrcEdges_[pointi] == -1
1264  && this->pointSrcFaces_[pointi] == -1
1265  );
1266 }
1267 
1268 
1269 template<class SrcPatchType, class TgtPatchType>
1270 bool Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::edgeCanIntersect
1271 (
1272  const label edgei
1273 ) const
1274 {
1275  return
1276  (
1277  edgePatchEdge(edgei, true) != -1
1278  && edgePatchEdge(edgei, false) == -1
1279  && (
1280  edgeTris_[edgei][0] == -1
1281  || triTgtFace_[edgeTris_[edgei][0]] == -1
1282  )
1283  && (
1284  edgeTris_[edgei][1] == -1
1285  || triTgtFace_[edgeTris_[edgei][1]] == -1
1286  )
1287  )
1288  || (
1289  edgePatchEdge(edgei, false) != -1
1290  && edgePatchEdge(edgei, true) == -1
1291  && (
1292  edgeTris_[edgei][0] == -1
1293  || triSrcFace_[edgeTris_[edgei][0]] == -1
1294  )
1295  && (
1296  edgeTris_[edgei][1] == -1
1297  || triSrcFace_[edgeTris_[edgei][1]] == -1
1298  )
1299  );
1300 }
1301 
1302 
1303 template<class SrcPatchType, class TgtPatchType>
1304 void
1305 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::snapPatchFaceTris
1306 (
1307  const label srcFacei,
1308  const label tgtFacei,
1309  const scalar snapTol
1310 )
1311 {
1312  static const FixedList<label, 3> triNoSnap({-1, -1, -1});
1313 
1314  const vector tgtNormal =
1315  triPointRef(tgtPoints_, patchFacePoints(tgtFacei, false)).normal();
1316 
1317  // Determine what source points snap to
1318  FixedList<barycentric2D, 3> srcFaceTgtTs(barycentric2D::uniform(-vGreat));
1319  FixedList<label, 3> srcFaceSnapTgtFacePoint(triNoSnap);
1320  FixedList<label, 3> srcFaceSnapTgtFaceEdge(triNoSnap);
1321  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1322  {
1323  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1324 
1325  if (!pointCanIntersect(srcPointi)) continue;
1326 
1327  if ((srcPointNormals_[srcPointi] & tgtNormal) < 0)
1328  {
1329  srcFaceTgtTs[srcFacePointi] =
1330  triIntersect::srcPointTgtTriIntersection
1331  (
1332  srcPoints_[srcPointi],
1333  srcPointNormals_[srcPointi],
1334  patchFacePointValues(tgtFacei, false, tgtPoints_)
1335  );
1336 
1337  for (label tgtFacePointi = 0; tgtFacePointi < 3; ++ tgtFacePointi)
1338  {
1339  const label tgtPointi =
1340  patchFacePoint(tgtFacei, tgtFacePointi, false);
1341 
1342  const label tgtFacePointi0 = (tgtFacePointi + 2) % 3;
1343  const label tgtFacePointi1 = (tgtFacePointi + 1) % 3;
1344 
1345  if
1346  (
1347  pointCanIntersect(tgtPointi)
1348  && mag(srcFaceTgtTs[srcFacePointi][tgtFacePointi0]) < snapTol
1349  && mag(srcFaceTgtTs[srcFacePointi][tgtFacePointi1]) < snapTol
1350  )
1351  {
1352  srcFaceSnapTgtFacePoint[srcFacePointi] =
1353  srcFaceSnapTgtFacePoint[srcFacePointi] == -1
1354  ? tgtFacePointi
1355  : -2;
1356 
1357  srcFaceTgtTs[srcFacePointi][tgtFacePointi] = 1;
1358  srcFaceTgtTs[srcFacePointi][tgtFacePointi0] = 0;
1359  srcFaceTgtTs[srcFacePointi][tgtFacePointi1] = 0;
1360  }
1361  }
1362 
1363  srcFaceSnapTgtFacePoint[srcFacePointi] =
1364  max(srcFaceSnapTgtFacePoint[srcFacePointi], -1);
1365 
1366  if (srcFaceSnapTgtFacePoint[srcFacePointi] != -1) continue;
1367 
1368  for (label tgtFaceEdgei = 0; tgtFaceEdgei < 3; ++ tgtFaceEdgei)
1369  {
1370  const label tgtFacePointi0 = tgtFaceEdgei;
1371  const label tgtFacePointi1 = (tgtFaceEdgei + 1) % 3;
1372  const label tgtFacePointiOpp = (tgtFaceEdgei + 2) % 3;
1373 
1374  if
1375  (
1376  srcFaceTgtTs[srcFacePointi][tgtFacePointi0] >= snapTol
1377  && srcFaceTgtTs[srcFacePointi][tgtFacePointi1] >= snapTol
1378  && mag(srcFaceTgtTs[srcFacePointi][tgtFacePointiOpp]) < snapTol
1379  )
1380  {
1381  srcFaceSnapTgtFaceEdge[srcFacePointi] =
1382  srcFaceSnapTgtFaceEdge[srcFacePointi] == -1
1383  ? tgtFaceEdgei
1384  : -2;
1385 
1386  const scalar eps =
1387  srcFaceTgtTs[srcFacePointi][tgtFacePointiOpp];
1388  srcFaceTgtTs[srcFacePointi][tgtFacePointiOpp] = 0;
1389  srcFaceTgtTs[srcFacePointi][tgtFacePointi0] /= 1 - eps;
1390  srcFaceTgtTs[srcFacePointi][tgtFacePointi1] /= 1 - eps;
1391  }
1392  }
1393 
1394  srcFaceSnapTgtFaceEdge[srcFacePointi] =
1395  max(srcFaceSnapTgtFaceEdge[srcFacePointi], -1);
1396  }
1397  }
1398 
1399  // Determine what target points snap to
1400  FixedList<barycentric2D, 3> tgtFaceSrcTs(barycentric2D::uniform(-vGreat));
1401  FixedList<label, 3> tgtFaceSnapSrcFacePoint(triNoSnap);
1402  FixedList<label, 3> tgtFaceSnapSrcFaceEdge(triNoSnap);
1403  forAll(tgtFaceSnapSrcFacePoint, tgtFacePointi)
1404  {
1405  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1406 
1407  if (!pointCanIntersect(tgtPointi)) continue;
1408 
1409  tgtFaceSrcTs[tgtFacePointi] =
1410  triIntersect::srcTriTgtPointIntersection
1411  (
1412  patchFacePointValues(srcFacei, true, srcPoints_),
1413  patchFacePointValues(srcFacei, true, srcPointNormals_),
1414  tgtPoints_[tgtPointi]
1415  );
1416 
1417  const vector srcNormal =
1418  triIntersect::srcTriInterpolate
1419  (
1420  tgtFaceSrcTs[tgtFacePointi],
1421  patchFacePointValues(srcFacei, true, srcPointNormals_)
1422  );
1423 
1424  if ((srcNormal & tgtNormal) < 0)
1425  {
1426  for (label srcFacePointi = 0; srcFacePointi < 3; ++ srcFacePointi)
1427  {
1428  const label srcPointi =
1429  patchFacePoint(srcFacei, srcFacePointi, true);
1430 
1431  const label srcFacePointi0 = (srcFacePointi + 2) % 3;
1432  const label srcFacePointi1 = (srcFacePointi + 1) % 3;
1433 
1434  if
1435  (
1436  pointCanIntersect(srcPointi)
1437  && mag(tgtFaceSrcTs[tgtFacePointi][srcFacePointi0]) < snapTol
1438  && mag(tgtFaceSrcTs[tgtFacePointi][srcFacePointi1]) < snapTol
1439  )
1440  {
1441  tgtFaceSnapSrcFacePoint[tgtFacePointi] =
1442  tgtFaceSnapSrcFacePoint[tgtFacePointi] == -1
1443  ? srcFacePointi
1444  : -2;
1445 
1446  tgtFaceSrcTs[tgtFacePointi][srcFacePointi] = 1;
1447  tgtFaceSrcTs[tgtFacePointi][srcFacePointi0] = 0;
1448  tgtFaceSrcTs[tgtFacePointi][srcFacePointi1] = 0;
1449  }
1450  }
1451 
1452  tgtFaceSnapSrcFacePoint[tgtFacePointi] =
1453  max(tgtFaceSnapSrcFacePoint[tgtFacePointi], -1);
1454 
1455  if (tgtFaceSnapSrcFacePoint[tgtFacePointi] != -1) continue;
1456 
1457  for (label srcFaceEdgei = 0; srcFaceEdgei < 3; ++ srcFaceEdgei)
1458  {
1459  const label srcFacePointi0 = srcFaceEdgei;
1460  const label srcFacePointi1 = (srcFaceEdgei + 1) % 3;
1461  const label srcFacePointiOpp = (srcFaceEdgei + 2) % 3;
1462 
1463  if
1464  (
1465  tgtFaceSrcTs[tgtFacePointi][srcFacePointi0] >= snapTol
1466  && tgtFaceSrcTs[tgtFacePointi][srcFacePointi1] >= snapTol
1467  && mag(tgtFaceSrcTs[tgtFacePointi][srcFacePointiOpp]) < snapTol
1468  )
1469  {
1470  tgtFaceSnapSrcFaceEdge[tgtFacePointi] =
1471  tgtFaceSnapSrcFaceEdge[tgtFacePointi] == -1
1472  ? srcFaceEdgei
1473  : -2;
1474 
1475  const scalar eps =
1476  tgtFaceSrcTs[tgtFacePointi][srcFacePointiOpp];
1477  tgtFaceSrcTs[tgtFacePointi][srcFacePointiOpp] = 0;
1478  tgtFaceSrcTs[tgtFacePointi][srcFacePointi0] /= 1 - eps;
1479  tgtFaceSrcTs[tgtFacePointi][srcFacePointi1] /= 1 - eps;
1480  }
1481  }
1482 
1483  tgtFaceSnapSrcFaceEdge[tgtFacePointi] =
1484  max(tgtFaceSnapSrcFaceEdge[tgtFacePointi], -1);
1485  }
1486  }
1487 
1488  // Return if there is nothing to do
1489  if
1490  (
1491  srcFaceSnapTgtFacePoint == triNoSnap
1492  && srcFaceSnapTgtFaceEdge == triNoSnap
1493  && tgtFaceSnapSrcFacePoint == triNoSnap
1494  && tgtFaceSnapSrcFaceEdge == triNoSnap
1495  )
1496  {
1497  return;
1498  }
1499 
1500  // Make point-point snaps symmetric
1501  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1502  {
1503  const label tgtFacePointi = srcFaceSnapTgtFacePoint[srcFacePointi];
1504 
1505  if (tgtFacePointi != -1)
1506  {
1507  tgtFaceSnapSrcFacePoint[tgtFacePointi] = srcFacePointi;
1508  }
1509  }
1510  forAll(tgtFaceSnapSrcFacePoint, tgtFacePointi)
1511  {
1512  const label srcFacePointi = tgtFaceSnapSrcFacePoint[tgtFacePointi];
1513 
1514  if (srcFacePointi != -1)
1515  {
1516  srcFaceSnapTgtFacePoint[srcFacePointi] = tgtFacePointi;
1517  }
1518  }
1519 
1520  // Do point-point motion
1521  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1522  {
1523  const label tgtFacePointi = srcFaceSnapTgtFacePoint[srcFacePointi];
1524 
1525  if (tgtFacePointi == -1) continue;
1526 
1527  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1528  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1529 
1530  if (!pointCanIntersect(srcPointi)) continue;
1531  if (!pointCanIntersect(tgtPointi)) continue;
1532 
1533  srcPoints_[tgtPointi] = srcPoints_[srcPointi];
1534  srcPointNormals_[tgtPointi] = srcPointNormals_[srcPointi];
1535  tgtPoints_[srcPointi] = tgtPoints_[tgtPointi];
1536 
1537  const point& srcP = srcPoints_[srcPointi];
1538  const vector& srcN = srcPointNormals_[srcPointi];
1539  const point& tgtP = tgtPoints_[tgtPointi];
1540  const vector d = (tensor::I - sqr(srcN)) & (tgtP - srcP);
1541  srcPoints_[srcPointi] += d/2;
1542  tgtPoints_[tgtPointi] -= d/2;
1543  }
1544 
1545  // Do source-point-target-edge motion
1546  forAll(srcFaceSnapTgtFaceEdge, srcFacePointi)
1547  {
1548  const label tgtFaceEdgei = srcFaceSnapTgtFaceEdge[srcFacePointi];
1549 
1550  if (tgtFaceEdgei == -1) continue;
1551 
1552  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1553 
1554  if (!pointCanIntersect(srcPointi)) continue;
1555 
1556  tgtPoints_[srcPointi] =
1557  triIntersect::tgtTriInterpolate
1558  (
1559  srcFaceTgtTs[srcFacePointi],
1560  patchFacePointValues(tgtFacei, false, tgtPoints_)
1561  );
1562 
1563  const point& srcP = srcPoints_[srcPointi];
1564  const vector& srcN = srcPointNormals_[srcPointi];
1565  const point& tgtP = tgtPoints_[srcPointi];
1566  const vector d = (tensor::I - sqr(srcN)) & (tgtP - srcP);
1567  srcPoints_[srcPointi] += d;
1568  }
1569 
1570  // Do target-point-source-edge motion
1571  forAll(tgtFaceSnapSrcFaceEdge, tgtFacePointi)
1572  {
1573  const label srcFaceEdgei = tgtFaceSnapSrcFaceEdge[tgtFacePointi];
1574 
1575  if (srcFaceEdgei == -1) continue;
1576 
1577  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1578 
1579  if (!pointCanIntersect(tgtPointi)) continue;
1580 
1581  srcPoints_[tgtPointi] =
1582  triIntersect::srcTriInterpolate
1583  (
1584  tgtFaceSrcTs[tgtFacePointi],
1585  patchFacePointValues(srcFacei, true, srcPoints_)
1586  );
1587  srcPointNormals_[tgtPointi] =
1588  triIntersect::srcTriInterpolate
1589  (
1590  tgtFaceSrcTs[tgtFacePointi],
1591  patchFacePointValues(srcFacei, true, srcPointNormals_)
1592  );
1593 
1594  const point& srcP = srcPoints_[tgtPointi];
1595  const vector& srcN = srcPointNormals_[tgtPointi];
1596  const point& tgtP = tgtPoints_[tgtPointi];
1597  const vector d = (tensor::I - sqr(srcN)) & (tgtP - srcP);
1598  tgtPoints_[tgtPointi] -= d;
1599  }
1600 
1601  // Do point-point snapping
1602  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1603  {
1604  const label tgtFacePointi = srcFaceSnapTgtFacePoint[srcFacePointi];
1605 
1606  if (tgtFacePointi == -1) continue;
1607 
1608  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1609  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1610 
1611  if (!pointCanIntersect(srcPointi)) continue;
1612  if (!pointCanIntersect(tgtPointi)) continue;
1613 
1614  srcPoints_[tgtPointi] = srcPoints_[srcPointi];
1615  srcPointNormals_[tgtPointi] = srcPointNormals_[srcPointi];
1616  tgtPoints_[srcPointi] = tgtPoints_[tgtPointi];
1617 
1618  pointPoints_[tgtPointi] = srcPointi;
1619 
1620  this->pointTgtEdges_[srcPointi] = this->pointTgtEdges_[tgtPointi];
1621  this->pointSrcEdges_[tgtPointi] = -1;
1622  this->pointTgtEdges_[tgtPointi] = -1;
1623 
1624  this->pointTgtPoints_[srcPointi] = this->pointTgtPoints_[tgtPointi];
1625  this->tgtPointPoints_[this->pointTgtPoints_[tgtPointi]] = srcPointi;
1626  this->pointSrcPoints_[tgtPointi] = -1;
1627  this->pointTgtPoints_[tgtPointi] = -1;
1628  }
1629 
1630  // Check
1631  checkPatchFace(srcFacei, true);
1632  checkPatchFace(tgtFacei, false);
1633 
1634  // Insertion workspace
1635  labelList insertPointis(1), insertEdgeis(2, label(-1)), fixedEdgeis(0);
1636 
1637  // Do source-point-target-edge snapping
1638  forAll(srcFaceSnapTgtFaceEdge, srcFacePointi)
1639  {
1640  const label tgtFaceEdgei = srcFaceSnapTgtFaceEdge[srcFacePointi];
1641 
1642  if (tgtFaceEdgei == -1) continue;
1643 
1644  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1645 
1646  if (!pointCanIntersect(srcPointi)) continue;
1647 
1648  const label tgtPatchEdgei =
1649  patchFacePatchEdge(tgtFacei, tgtFaceEdgei, false);
1650 
1651  // Loop all target face tris and find the best edge into which to
1652  // insert the source point
1653  label insertTgtEdgei = -1;
1654  scalar insertDistSqr = vGreat;
1655  forAll(tgtFaceTris_[tgtFacei], tgtFaceTrii)
1656  {
1657  const label tgtTrii = tgtFaceTris_[tgtFacei][tgtFaceTrii];
1658 
1659  forAll(triEdges_[tgtTrii], tgtTriEdgei)
1660  {
1661  const label tgtEdgei = triEdges_[tgtTrii][tgtTriEdgei];
1662 
1663  if
1664  (
1665  edgeCanIntersect(tgtEdgei)
1666  && tgtPatchEdgei == edgePatchEdge(tgtEdgei, false)
1667  )
1668  {
1669  const scalar distSqr = circumDistSqr(tgtTrii, srcPointi);
1670 
1671  if (distSqr < insertDistSqr)
1672  {
1673  insertDistSqr = distSqr;
1674  insertTgtEdgei = tgtEdgei;
1675  }
1676  }
1677  }
1678  }
1679 
1680  if (insertTgtEdgei != -1)
1681  {
1682  insertPointis[0] = srcPointi;
1683  insertPoints
1684  (
1685  insertTgtEdgei,
1686  false,
1687  insertPointis,
1688  insertEdgeis,
1689  fixedEdgeis
1690  );
1691  }
1692  }
1693 
1694  // Do target-point-source-edge snapping
1695  forAll(tgtFaceSnapSrcFaceEdge, tgtFacePointi)
1696  {
1697  const label srcFaceEdgei = tgtFaceSnapSrcFaceEdge[tgtFacePointi];
1698 
1699  if (srcFaceEdgei == -1) continue;
1700 
1701  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1702 
1703  if (!pointCanIntersect(tgtPointi)) continue;
1704 
1705  const label srcPatchEdgei =
1706  patchFacePatchEdge(srcFacei, srcFaceEdgei, true);
1707 
1708  srcPoints_[tgtPointi] =
1709  triIntersect::srcTriInterpolate
1710  (
1711  tgtFaceSrcTs[tgtFacePointi],
1712  patchFacePointValues(srcFacei, true, srcPoints_)
1713  );
1714  srcPointNormals_[tgtPointi] =
1715  triIntersect::srcTriInterpolate
1716  (
1717  tgtFaceSrcTs[tgtFacePointi],
1718  patchFacePointValues(srcFacei, true, srcPointNormals_)
1719  );
1720 
1721  // Loop all target face tris and find the best edge into which to
1722  // insert the source point
1723  label insertSrcEdgei = -1;
1724  scalar insertDistSqr = vGreat;
1725  forAll(srcFaceTris_[srcFacei], srcFaceTrii)
1726  {
1727  const label srcTrii = srcFaceTris_[srcFacei][srcFaceTrii];
1728 
1729  forAll(triEdges_[srcTrii], srcTriEdgei)
1730  {
1731  const label srcEdgei = triEdges_[srcTrii][srcTriEdgei];
1732 
1733  if
1734  (
1735  edgeCanIntersect(srcEdgei)
1736  && srcPatchEdgei == edgePatchEdge(srcEdgei, true)
1737  )
1738  {
1739  const scalar distSqr = circumDistSqr(srcTrii, tgtPointi);
1740 
1741  if (distSqr < insertDistSqr)
1742  {
1743  insertDistSqr = distSqr;
1744  insertSrcEdgei = srcEdgei;
1745  }
1746  }
1747  }
1748  }
1749 
1750  if (insertSrcEdgei != -1)
1751  {
1752  insertPointis[0] = tgtPointi;
1753  insertPoints
1754  (
1755  insertSrcEdgei,
1756  false,
1757  insertPointis,
1758  insertEdgeis,
1759  fixedEdgeis
1760  );
1761  }
1762  }
1763 
1764  // Check
1765  checkPatchFace(srcFacei, true);
1766  checkPatchFace(tgtFacei, false);
1767 }
1768 
1769 
1770 template<class SrcPatchType, class TgtPatchType>
1771 bool Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::intersectTris
1772 (
1773  const label srcTrii,
1774  const label tgtTrii
1775 )
1776 {
1777  // Get the source and target faces that these tris are associated with
1778  const label srcFacei = triSrcFace_[srcTrii];
1779  const label tgtFacei = triTgtFace_[tgtTrii];
1780  if (srcFacei == -1 || tgtFacei == -1)
1781  {
1782  FatalErrorInFunction
1783  << "Tri-intersections must be between a tri associated with the "
1784  << "source patch and one associated with the target patch"
1785  << exit(FatalError);
1786  }
1787 
1788  // Do intersection
1789  DynamicList<point> ictSrcPoints;
1790  DynamicList<vector> ictSrcPointNormals;
1791  DynamicList<point> ictTgtPoints;
1792  DynamicList<triIntersect::location> ictPointLocations;
1793  triIntersect::intersectTris
1794  (
1795  triPointValues(srcTrii, srcPoints_),
1796  triPointValues(srcTrii, srcPointNormals_),
1797  triOwns(srcTrii),
1798  triOtherTriPoints(srcTrii, tgtTrii),
1799  triPointValues(tgtTrii, tgtPoints_),
1800  triOwns(tgtTrii),
1801  triOtherTriPoints(tgtTrii, srcTrii),
1802  ictSrcPoints,
1803  ictSrcPointNormals,
1804  ictTgtPoints,
1805  ictPointLocations,
1806  this->debug > 3,
1807  this->debug > 4
1808  ? "srcTrii##tgtTrii=" + name(srcTrii) + name(tgtTrii)
1809  : ""
1810  );
1811 
1812  // Return false if no intersection
1813  if (!ictPointLocations.size())
1814  {
1815  return false;
1816  }
1817 
1818  // Loop over the intersection points looking for elements that are already
1819  // intersected with the opposing surface. If any are found, then this
1820  // intersection is considered invalid.
1821  forAll(ictPointLocations, ictPointi)
1822  {
1823  const triIntersect::location& l = ictPointLocations[ictPointi];
1824 
1825  if (l.isSrcPoint() && !l.isTgtPoint())
1826  {
1827  const label pointi = triPoint(srcTrii, l.srcPointi());
1828 
1829  if
1830  (
1831  this->pointTgtPoints_[pointi] != -1
1832  || this->pointTgtEdges_[pointi] != -1
1833  || this->pointTgtFaces_[pointi] != -1
1834  )
1835  {
1836  return false;
1837  }
1838  }
1839  if (!l.isSrcPoint() && l.isTgtPoint())
1840  {
1841  const label pointi = triPoint(tgtTrii, l.tgtPointi());
1842 
1843  if
1844  (
1845  this->pointSrcPoints_[pointi] != -1
1846  || this->pointSrcEdges_[pointi] != -1
1847  || this->pointSrcFaces_[pointi] != -1
1848  )
1849  {
1850  return false;
1851  }
1852  }
1853  if (l.isIntersection())
1854  {
1855  const label srcEdgei = triEdges_[srcTrii][l.srcEdgei()];
1856  const label tgtEdgei = triEdges_[tgtTrii][l.tgtEdgei()];
1857 
1858  const labelPair srcPatchEdges = edgePatchEdges(srcEdgei);
1859  const labelPair tgtPatchEdges = edgePatchEdges(tgtEdgei);
1860 
1861  if
1862  (
1863  (srcPatchEdges[0] != -1 && tgtPatchEdges[1] != -1)
1864  && (srcPatchEdges[1] != -1 || tgtPatchEdges[0] != -1)
1865  )
1866  {
1867  return false;
1868  }
1869 
1870  const label srcTrij =
1871  edgeTris_[srcEdgei][edgeTris_[srcEdgei][0] == srcTrii];
1872  const label tgtTrij =
1873  edgeTris_[tgtEdgei][edgeTris_[tgtEdgei][0] == tgtTrii];
1874 
1875  if
1876  (
1877  (srcTrij != -1 && triTgtFace_[srcTrij] != -1)
1878  || (tgtTrij != -1 && triSrcFace_[tgtTrij] != -1)
1879  )
1880  {
1881  return false;
1882  }
1883  }
1884  }
1885 
1886  // Flags controlling what operations to do
1887  bool insertPointsIntoTri = false;
1888  bool insertPointsIntoEdges = false;
1889 
1890  // Loop over the intersection points creating a map from intersection point
1891  // indices to point indices. Overwrite any intersected source or target
1892  // points, and create points generated by any new intersections. Point-tri
1893  // and point-edge associations are generated later during insertion.
1894  labelList ictPointiToPointi(ictPointLocations.size(), -1);
1895  forAll(ictPointLocations, ictPointi)
1896  {
1897  const triIntersect::location& l = ictPointLocations[ictPointi];
1898 
1899  if (l.isSrcPoint() && !l.isTgtPoint())
1900  {
1901  const label pointi = triPoint(srcTrii, l.srcPointi());
1902 
1903  ictPointiToPointi[ictPointi] = pointi;
1904 
1905  // Store the projected location on the target side
1906  tgtPoints_[pointi] = ictTgtPoints[ictPointi];
1907 
1908  insertPointsIntoTri = true;
1909  }
1910  if (!l.isSrcPoint() && l.isTgtPoint())
1911  {
1912  const label pointi = triPoint(tgtTrii, l.tgtPointi());
1913 
1914  ictPointiToPointi[ictPointi] = pointi;
1915 
1916  // Store the projected location on the source side
1917  srcPoints_[pointi] = ictSrcPoints[ictPointi];
1918  srcPointNormals_[pointi] = ictSrcPointNormals[ictPointi];
1919 
1920  insertPointsIntoTri = true;
1921  }
1922  if (l.isIntersection())
1923  {
1924  const label srcPatchEdgei =
1925  edgePatchEdge(triEdges_[srcTrii][l.srcEdgei()], true);
1926  const label tgtPatchEdgei =
1927  edgePatchEdge(triEdges_[tgtTrii][l.tgtEdgei()], false);
1928 
1929  if (srcPatchEdgei != -1 && tgtPatchEdgei != -1)
1930  {
1931  ictPointiToPointi[ictPointi] = srcPoints_.size();
1932 
1933  const label pointi = newPointi();
1934 
1935  srcPoints_[pointi] = ictSrcPoints[ictPointi];
1936  srcPointNormals_[pointi] = ictSrcPointNormals[ictPointi];
1937  tgtPoints_[pointi] = ictTgtPoints[ictPointi];
1938 
1939  insertPointsIntoEdges = true;
1940  }
1941  }
1942  }
1943 
1944  // Store the source and target tri edges and edge-ownerships
1945  const FixedList<label, 3> srcTriEdges = triEdges_[srcTrii];
1946  const FixedList<bool, 3> srcTriOwns = triOwns(srcTrii);
1947  const FixedList<label, 3> tgtTriEdges = triEdges_[tgtTrii];
1948  const FixedList<bool, 3> tgtTriOwns = triOwns(tgtTrii);
1949 
1950  // Edges which are to be fixed during insertion. These constraints are not
1951  // needed. The only edges that get inserted into are those that associate
1952  // with patch edges, and these are constrained anyway. Constraining all
1953  // edges of the tris leads to poor triangle qualities, as it limits the
1954  // that re-triangulation insertPoints can do. For testing, however, this
1955  // can be quite useful, hence the option to make the system over-constrain
1956  // the problem.
1957  DynamicList<label> fixedSrcEdgeis, fixedTgtEdgeis;
1958  #ifdef OVERCONSTRAIN
1959  forAll(srcTriEdges, srcTriEdgei)
1960  {
1961  fixedSrcEdgeis.append(srcTriEdges[srcTriEdgei]);
1962  }
1963  forAll(tgtTriEdges, tgtTriEdgei)
1964  {
1965  fixedTgtEdgeis.append(tgtTriEdges[tgtTriEdgei]);
1966  }
1967  #endif
1968 
1969  // Insert points into the triangles
1970  if (insertPointsIntoTri)
1971  {
1972  DynamicList<label> insertSrcPointis(3), insertTgtPointis(3);
1973  DynamicList<label> insertEdgeis;
1974  forAll(ictPointLocations, ictPointi)
1975  {
1976  const label pointi = ictPointiToPointi[ictPointi];
1977 
1978  if (pointi == -1) continue;
1979 
1980  const triIntersect::location& l = ictPointLocations[ictPointi];
1981 
1982  if (l.isSrcPoint() && !l.isTgtPoint())
1983  {
1984  insertSrcPointis.append(pointi);
1985  }
1986  if (!l.isSrcPoint() && l.isTgtPoint())
1987  {
1988  insertTgtPointis.append(pointi);
1989  }
1990  }
1991 
1992  if (insertSrcPointis.size())
1993  {
1994  insertPoints
1995  (
1996  tgtTrii,
1997  true,
1998  insertSrcPointis,
1999  insertEdgeis,
2000  fixedTgtEdgeis
2001  );
2002  }
2003 
2004  if (insertTgtPointis.size())
2005  {
2006  insertPoints
2007  (
2008  srcTrii,
2009  true,
2010  insertTgtPointis,
2011  insertEdgeis,
2012  fixedSrcEdgeis
2013  );
2014  }
2015 
2016  // Check
2017  checkPatchFace(srcFacei, true);
2018  checkPatchFace(tgtFacei, false);
2019  }
2020 
2021  // Insert points into the edges
2022  if (insertPointsIntoEdges)
2023  {
2024  DynamicList<label> insertPointis(2), insertEdgeis(3);
2025  forAll(ictPointLocations, ictPointi0)
2026  {
2027  const label ictPointi1 = ictPointLocations.fcIndex(ictPointi0);
2028 
2029  const label pointi0 = ictPointiToPointi[ictPointi0];
2030  const label pointi1 = ictPointiToPointi[ictPointi1];
2031 
2032  const triIntersect::location& l0 = ictPointLocations[ictPointi0];
2033  const triIntersect::location& l1 = ictPointLocations[ictPointi1];
2034 
2035  insertPointis.clear();
2036  insertEdgeis.clear();
2037  insertEdgeis.append(-1);
2038 
2039  if
2040  (
2041  (
2042  l0.isIntersection()
2043  && l1.isSrcPoint()
2044  && l0.srcEdgei() != (l1.srcPointi() + 1) % 3
2045  )
2046  || (
2047  l0.isSrcPoint()
2048  && l1.isIntersection()
2049  && (l0.srcPointi() + 1) % 3 != l1.srcEdgei()
2050  )
2051  || (
2052  l0.isIntersection()
2053  && l1.isIntersection()
2054  && l0.srcEdgei() == l1.srcEdgei()
2055  )
2056  )
2057  {
2058  const label srcTriEdgei =
2059  l0.isIntersection() ? l0.srcEdgei() : l1.srcEdgei();
2060 
2061  if (!l0.isSrcPoint() && pointi0 != -1)
2062  {
2063  insertPointis.append(pointi0);
2064  insertEdgeis.append(-1);
2065  }
2066  if (!l1.isSrcPoint() && pointi1 != -1)
2067  {
2068  insertPointis.append(pointi1);
2069  insertEdgeis.append(-1);
2070  }
2071  if (!srcTriOwns[srcTriEdgei])
2072  {
2073  inplaceReverseList(insertPointis);
2074  }
2075 
2076  if (insertPointis.size())
2077  {
2078  insertPoints
2079  (
2080  srcTriEdges[srcTriEdgei],
2081  false,
2082  insertPointis,
2083  insertEdgeis,
2084  fixedSrcEdgeis
2085  );
2086  }
2087 
2088  #ifdef OVERCONSTRAIN
2089  fixedSrcEdgeis[srcTriEdgei] = insertEdgeis.remove();
2090  fixedSrcEdgeis.append(insertEdgeis);
2091  #endif
2092  }
2093 
2094  if
2095  (
2096  (
2097  l0.isIntersection()
2098  && l1.isTgtPoint()
2099  && l0.tgtEdgei() != (l1.tgtPointi() + 1) % 3
2100  )
2101  || (
2102  l0.isTgtPoint()
2103  && l1.isIntersection()
2104  && (l0.tgtPointi() + 1) % 3 != l1.tgtEdgei()
2105  )
2106 
2107  || (
2108  l0.isIntersection()
2109  && l1.isIntersection()
2110  && l0.tgtEdgei() == l1.tgtEdgei()
2111  )
2112  )
2113  {
2114  const label tgtTriEdgei =
2115  l0.isIntersection() ? l0.tgtEdgei() : l1.tgtEdgei();
2116 
2117  if (!l0.isTgtPoint() && pointi0 != -1)
2118  {
2119  insertPointis.append(pointi0);
2120  insertEdgeis.append(-1);
2121  }
2122  if (!l1.isTgtPoint() && pointi1 != -1)
2123  {
2124  insertPointis.append(pointi1);
2125  insertEdgeis.append(-1);
2126  }
2127  if (tgtTriOwns[tgtTriEdgei])
2128  {
2129  inplaceReverseList(insertPointis);
2130  }
2131 
2132  if (insertPointis.size())
2133  {
2134  insertPoints
2135  (
2136  tgtTriEdges[tgtTriEdgei],
2137  false,
2138  insertPointis,
2139  insertEdgeis,
2140  fixedTgtEdgeis
2141  );
2142  }
2143 
2144  #ifdef OVERCONSTRAIN
2145  fixedTgtEdgeis[tgtTriEdgei] = insertEdgeis.remove();
2146  fixedTgtEdgeis.append(insertEdgeis);
2147  #endif
2148  }
2149  }
2150 
2151  // Check
2152  checkPatchFace(srcFacei, true);
2153  checkPatchFace(tgtFacei, false);
2154  }
2155 
2156  return true;
2157 }
2158 
2159 
2160 template<class SrcPatchType, class TgtPatchType>
2161 void
2162 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::intersectPatchFaceTris
2163 (
2164  const label srcFacei,
2165  const label tgtFacei
2166 )
2167 {
2168  // Get the source face tris and make sure none are marked
2169  const DynamicList<label>& srcFaceTris = srcFaceTris_[srcFacei];
2170  forAll(srcFaceTris, srcFaceTrii)
2171  {
2172  const label srcTrii = srcFaceTris[srcFaceTrii];
2173 
2174  if (triMarkedTris_[srcTrii] != -1)
2175  {
2176  markedTriTris_[triMarkedTris_[srcTrii]] = -1;
2177  triMarkedTris_[srcTrii] = -1;
2178  }
2179  }
2180 
2181  // Loop the source face tris
2182  label srcFaceTrii = 0;
2183  while (srcFaceTrii < srcFaceTris.size())
2184  {
2185  // Mark this source tri
2186  const label srcTrii = srcFaceTris[srcFaceTrii];
2187  markedTriTris_.append(srcTrii);
2188  triMarkedTris_[srcTrii] = markedTriTris_.size() - 1;
2189 
2190  // Get the target face tris and make sure none are marked
2191  const DynamicList<label>& tgtFaceTris = tgtFaceTris_[tgtFacei];
2192  forAll(tgtFaceTris, tgtFaceTrii)
2193  {
2194  const label tgtTrii = tgtFaceTris[tgtFaceTrii];
2195 
2196  if (triMarkedTris_[tgtTrii] != -1)
2197  {
2198  markedTriTris_[triMarkedTris_[tgtTrii]] = -1;
2199  triMarkedTris_[tgtTrii] = -1;
2200  }
2201  }
2202 
2203  // Loop the target face tris
2204  label tgtFaceTrii = 0;
2205  while (tgtFaceTrii < tgtFaceTris.size())
2206  {
2207  // Mark this target tri
2208  const label tgtTrii = tgtFaceTris[tgtFaceTrii];
2209  markedTriTris_.append(tgtTrii);
2210  triMarkedTris_[tgtTrii] = markedTriTris_.size() - 1;
2211 
2212  /*
2213  // If the target tri has not unmarked then snap
2214  if (triMarkedTris_[tgtFaceTris[tgtFaceTrii]] != -1)
2215  {
2216  snapTris(srcTrii, tgtTrii);
2217  }
2218 
2219  // If the source tri has unmarked then break
2220  if (triMarkedTris_[srcFaceTris[srcFaceTrii]] == -1) break;
2221  */
2222 
2223  // If the target tri has not unmarked then intersect
2224  if (triMarkedTris_[tgtFaceTris[tgtFaceTrii]] != -1)
2225  {
2226  intersectTris(srcTrii, tgtTrii);
2227  }
2228 
2229  // If the source tri has unmarked then break
2230  if (triMarkedTris_[srcFaceTris[srcFaceTrii]] == -1) break;
2231 
2232  // Iterate backwards to find the next marked target tri
2233  while
2234  (
2235  tgtFaceTrii >= 0
2236  && triMarkedTris_[tgtFaceTris[tgtFaceTrii]] == -1
2237  )
2238  {
2239  -- tgtFaceTrii;
2240  }
2241 
2242  // Increment to the next unmarked target tri
2243  ++ tgtFaceTrii;
2244  }
2245 
2246  // Iterate backwards to find the next marked source tri
2247  while
2248  (
2249  srcFaceTrii >= 0
2250  && triMarkedTris_[srcFaceTris[srcFaceTrii]] == -1
2251  )
2252  {
2253  -- srcFaceTrii;
2254  }
2255 
2256  // Increment to the next unmarked source tri
2257  ++ srcFaceTrii;
2258  }
2259 
2260  // Clear the marked tris
2261  forAll(markedTriTris_, markedTrii)
2262  {
2263  const label trii = markedTriTris_[markedTrii];
2264  if (trii != -1)
2265  {
2266  triMarkedTris_[trii] = -1;
2267  }
2268  }
2269  markedTriTris_.clear();
2270 }
2271 
2272 
2273 template<class SrcPatchType, class TgtPatchType>
2274 bool
2275 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::conformPatchFaceTris
2276 (
2277  const label patchFacei,
2278  const label otherPatchFacei,
2279  const bool isSrc
2280 )
2281 {
2282  const labelList& patchFaceTris =
2283  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
2284  const labelList& otherPatchFaceTris =
2285  isSrc ? tgtFaceTris_[otherPatchFacei] : srcFaceTris_[otherPatchFacei];
2286 
2287  const labelList& triOtherPatchFaces =
2288  isSrc ? triTgtFace_ : triSrcFace_;
2289 
2290  const pointField& points = isSrc ? srcPoints_ : tgtPoints_;
2291 
2292  const vector& patchFaceN =
2293  isSrc
2294  ? this->srcPatch_.faceNormals()[patchFacei]
2295  : this->tgtPatch_.faceNormals()[patchFacei];
2296 
2297  auto pointIntersectsPatchFace = [&]
2298  (
2299  const label pointi,
2300  const label patchFacei,
2301  const bool isSrc
2302  )
2303  {
2304  const labelList& pointPatchPoints =
2305  isSrc ? this->pointSrcPoints_ : this->pointTgtPoints_;
2306  const labelList& pointPatchEdges =
2307  isSrc ? this->pointSrcEdges_ : this->pointTgtEdges_;
2308  const labelList& pointPatchFaces =
2309  isSrc ? this->pointSrcFaces_ : this->pointTgtFaces_;
2310 
2311  const triFace patchFacePatchPoints
2312  (
2313  this->patchFacePatchPoints(patchFacei, isSrc)
2314  );
2315  const triFace patchFacePatchEdges
2316  (
2317  this->patchFacePatchEdges(patchFacei, isSrc)
2318  );
2319 
2320  return
2321  findIndex
2322  (
2323  patchFacePatchPoints,
2324  pointPatchPoints[pointi]
2325  ) != -1
2326  || findIndex
2327  (
2328  patchFacePatchEdges,
2329  pointPatchEdges[pointi]
2330  ) != -1
2331  || pointPatchFaces[pointi] == patchFacei;
2332  };
2333 
2334  auto edgeIntersectsPatchFace = [&]
2335  (
2336  const label edgei,
2337  const label patchFacei,
2338  const bool isSrc
2339  )
2340  {
2341  const edge e = edgePoints(edgei);
2342  return
2343  pointIntersectsPatchFace(e[0], patchFacei, isSrc)
2344  && pointIntersectsPatchFace(e[1], patchFacei, isSrc);
2345  };
2346 
2347  // Get all the edges that need conforming to
2348  DynamicList<label> conformEdgeis;
2349  forAll(otherPatchFaceTris, otherPatchFaceTrii)
2350  {
2351  const label trii = otherPatchFaceTris[otherPatchFaceTrii];
2352 
2353  forAll(triEdges_[trii], triEdgei)
2354  {
2355  const label edgei = triEdges_[trii][triEdgei];
2356 
2357  const label trij =
2358  edgeTris_[edgei][edgeTris_[edgei][0] == trii];
2359 
2360  const label otherPatchFacej =
2361  trij == -1 ? -1 : triOtherPatchFaces[trij];
2362 
2363  if
2364  (
2365  otherPatchFacei != otherPatchFacej
2366  && edgeIntersectsPatchFace(edgei, patchFacei, isSrc)
2367  )
2368  {
2369  conformEdgeis.append(edgei);
2370  }
2371  }
2372  }
2373 
2374  // Remove edges already present in the triangulation by shuffling up
2375  auto isConformed = [&](const label edgei)
2376  {
2377  const edge e = edgePoints(edgei);
2378 
2379  forAll(patchFaceTris, patchFaceTrii)
2380  {
2381  const label trii = patchFaceTris[patchFaceTrii];
2382 
2383  forAll(triEdges_[trii], triEdgei)
2384  {
2385  const label edgei = triEdges_[trii][triEdgei];
2386 
2387  if (edge::compare(e, edgePoints(edgei)) != 0)
2388  {
2389  return true;
2390  }
2391  }
2392  }
2393 
2394  return false;
2395  };
2396  {
2397  label conformi = 0;
2398  forAll(conformEdgeis, conformj)
2399  {
2400  if (!isConformed(conformEdgeis[conformj]))
2401  {
2402  conformEdgeis[conformi] = conformEdgeis[conformj];
2403  ++ conformi;
2404  }
2405  }
2406  conformEdgeis.resize(conformi);
2407  }
2408 
2409  // Quit if there is nothing to conform to
2410  if (conformEdgeis.empty()) return true;
2411 
2412  // Get the next patch face tri index connected to a given point
2413  auto nextConnectedPatchFaceTri = [&]
2414  (
2415  const label pointi,
2416  const label patchFaceTrii0 = 0
2417  )
2418  {
2419  for
2420  (
2421  label patchFaceTrii = patchFaceTrii0;
2422  patchFaceTrii < patchFaceTris.size();
2423  ++ patchFaceTrii
2424  )
2425  {
2426  const label trii = patchFaceTris[patchFaceTrii];
2427 
2428  if (findIndex(triPoints(trii), pointi) != -1)
2429  {
2430  return patchFaceTrii;
2431  }
2432  }
2433 
2434  return label(-1);
2435  };
2436 
2437  // Track from a point on a triangle towards a given point. Stop at an edge
2438  // and set the index of that edge and update the local coordinates.
2439  auto trackToEdge = [&]
2440  (
2441  const label trii,
2442  label& triEdgei,
2443  barycentric2D& y,
2444  const label pointi1
2445  )
2446  {
2447  const triPointRef t = triPoints(trii).tri(points);
2448 
2449  const barycentricTensor2D A(t.a(), t.b(), t.c());
2450 
2451  const point p = A & y;
2452  const vector dp = points[pointi1] - p;
2453 
2454  const vector ab = t.b() - t.a();
2455  const vector ac = t.c() - t.a();
2456  const vector bc = t.c() - t.b();
2457  const scalar detA = (ab ^ ac) & patchFaceN;
2458  const barycentricTensor2D T
2459  (
2460  patchFaceN ^ bc,
2461  ac ^ patchFaceN,
2462  patchFaceN ^ ab
2463  );
2464 
2465  const barycentric2D TDp = dp & T;
2466 
2467  label iH = -1;
2468  scalar lambdaByDetAH =
2469  !std::isnormal(detA) || detA < 0 ? vGreat : 1/detA;
2470 
2471  forAll(TDp, i)
2472  {
2473  if (TDp[i] < - detA*small)
2474  {
2475  const scalar lambdaByDetA = - y[i]/TDp[i];
2476 
2477  if (0 <= lambdaByDetA && lambdaByDetA < lambdaByDetAH)
2478  {
2479  iH = i;
2480  lambdaByDetAH = lambdaByDetA;
2481  }
2482  }
2483  }
2484 
2485  y += lambdaByDetAH*TDp;
2486  forAll(y, i)
2487  {
2488  y.replace(i, i == iH ? 0 : max(0, y[i]));
2489  }
2490  if (iH == -1)
2491  {
2492  y /= cmptSum(y);
2493  }
2494 
2495  triEdgei = iH == -1 ? -1 : (iH + 1) % 3;
2496  };
2497 
2498  // Cross the edge. Return true if the edge can be crossed, and false
2499  // otherwise. Set the new triangle, new edge index, and transform the local
2500  // coordinates appropriately.
2501  auto crossEdge = [&](label& trii, label& triEdgei, barycentric2D& y)
2502  {
2503  if (triEdgei == -1) return false;
2504 
2505  const label edgei = triEdges_[trii][triEdgei];
2506 
2507  if (edgePatchEdges(edgei) != labelPair(-1, -1)) return false;
2508 
2509  const label trij = edgeTris_[edgei][edgeTris_[edgei][0] == trii];
2510 
2511  if (trij == -1) return false;
2512 
2513  if (triSrcFace_[trii] != triSrcFace_[trij]) return false;
2514  if (triTgtFace_[trii] != triTgtFace_[trij]) return false;
2515 
2516  const label triEdgej = findIndex(triEdges_[trij], edgei);
2517 
2518  auto inplaceRotate = [](barycentric2D& y, label n)
2519  {
2520  n = n % 3;
2521 
2522  if (n == 1)
2523  {
2524  scalar temp = y.a();
2525  y.a() = y.b();
2526  y.b() = y.c();
2527  y.c() = temp;
2528  }
2529 
2530  if (n == 2)
2531  {
2532  scalar temp = y.c();
2533  y.c() = y.b();
2534  y.b() = y.a();
2535  y.a() = temp;
2536  }
2537  };
2538 
2539  inplaceRotate(y, triEdgei - 1 + 3);
2540  Swap(y.b(), y.c());
2541  inplaceRotate(y, 1 - triEdgej + 3);
2542 
2543  trii = trij;
2544  triEdgei = triEdgej;
2545 
2546  return true;
2547  };
2548 
2549  // Assume successful
2550  bool success = true;
2551 
2552  // Conform each identified edge in turn
2553  forAll(conformEdgeis, conformi)
2554  {
2555  const label edgei = conformEdgeis[conformi];
2556 
2557  // Skip if earlier operations have inadvertently resulted in this
2558  // edge being conformed to
2559  if (isConformed(edgei)) continue;
2560 
2561  // Get the edge and the points
2562  const edge e = edgePoints(edgei);
2563  const label pointi0 = e[0], pointi1 = e[1];
2564 
2565  // Find triangles which contains the first and last points
2566  const label patchFaceTrii0 = nextConnectedPatchFaceTri(pointi0);
2567  const label patchFaceTrii1 = nextConnectedPatchFaceTri(pointi1);
2568 
2569  // Conformation to this edge is not possible if either end point is not
2570  // present in the patch face triangulation
2571  if (patchFaceTrii0 == -1 || patchFaceTrii1 == -1) continue;
2572 
2573  // Track from point zero until in a triangle which contains point one.
2574  // Build a route of tris and edges along the track.
2575  DynamicList<label> routeTriis;
2576  DynamicList<label> routeEdgeis;
2577  label routePatchFaceTrii0 = -1;
2578  while (true)
2579  {
2580  routeTriis.clear();
2581  routeEdgeis.clear();
2582 
2583  routePatchFaceTrii0 =
2584  nextConnectedPatchFaceTri(pointi0, routePatchFaceTrii0 + 1);
2585 
2586  if (routePatchFaceTrii0 == -1) break;
2587 
2588  label trii = patchFaceTris[routePatchFaceTrii0];
2589  label triEdgei = -1;
2590  barycentric2D y(0, 0, 0);
2591  y[findIndex(triPoints(trii), pointi0)] = 1;
2592 
2593  forAll(patchFaceTris, iter)
2594  {
2595  if (findIndex(triPoints(trii), pointi0) != -1)
2596  {
2597  routeTriis.clear();
2598  routeEdgeis.clear();
2599  }
2600 
2601  routeTriis.append(trii);
2602 
2603  if (findIndex(triPoints(trii), pointi1) != -1) break;
2604 
2605  trackToEdge(trii, triEdgei, y, pointi1);
2606 
2607  if (triEdgei == -1) break;
2608 
2609  routeEdgeis.append(triEdges_[trii][triEdgei]);
2610 
2611  if (!crossEdge(trii, triEdgei, y)) break;
2612  }
2613 
2614  if (findIndex(triPoints(trii), pointi1) != -1) break;
2615  }
2616 
2617  // Conform the triangulation to the route
2618  if (routeEdgeis.size() == 0)
2619  {
2620  // No route was found. Don't do anything. This edge will not be
2621  // conformed to and the intersection will disconnect here.
2622  if (this->debug > 1)
2623  {
2624  WarningInFunction
2625  << indent << "Failed to route edge " << e
2626  << " through the triangulation of "
2627  << (isSrc ? "source" : "target") << " face #"
2628  << patchFacei << endl;
2629  writePatchFace(patchFacei, isSrc);
2630  }
2631  success = false;
2632  }
2633  else if (routeEdgeis.size() == 1)
2634  {
2635  // The route only spans a single edge. Flip it to conform.
2636  flipEdge(routeEdgeis.first());
2637  }
2638  else
2639  {
2640  // The route spans multiple edges. Remove all triangles in the
2641  // route and then split the resulting polygon along the conform
2642  // edge and triangulate both sides independently.
2643 
2644  // !!! Check that removal of the tris does not result in any points
2645  // becoming disconnected from the mesh
2646  // ...
2647 
2648  // Form the polygons on either side of the conform edge
2649  DynamicList<label> leftPolyPointis, rightPolyPointis;
2650  DynamicList<label> leftPolyEdgeis, rightPolyEdgeis;
2651  DynamicList<label> leftPolyTriis, rightPolyTriis;
2652  {
2653  // First point
2654  leftPolyPointis.append(pointi0);
2655  rightPolyPointis.append(pointi0);
2656 
2657  // First triangle
2658  {
2659  const label trii = routeTriis.first();
2660  const label edgei1 = routeEdgeis.first();
2661  const label triEdgei1 = findIndex(triEdges_[trii], edgei1);
2662  const label triEdgeiLeft = (triEdgei1 + 2) % 3;
2663  const label triEdgeiRight = (triEdgei1 + 1) % 3;
2664  leftPolyEdgeis.append(triEdges_[trii][triEdgeiLeft]);
2665  rightPolyEdgeis.append(triEdges_[trii][triEdgeiRight]);
2666  leftPolyTriis.append(trii);
2667  rightPolyTriis.append(trii);
2668  leftPolyPointis.append(triPoint(trii, triEdgei1));
2669  rightPolyPointis.append(triPoint(trii, triEdgeiRight));
2670  }
2671 
2672  // Intermediate triangles
2673  for
2674  (
2675  label routeTrii = 1;
2676  routeTrii < routeTriis.size() - 1;
2677  ++ routeTrii
2678  )
2679  {
2680  const label trii = routeTriis[routeTrii];
2681  const label edgei0 = routeEdgeis[routeTrii - 1];
2682  const label edgei1 = routeEdgeis[routeTrii];
2683  const label triEdgei0 = findIndex(triEdges_[trii], edgei0);
2684  const label triEdgei1 = findIndex(triEdges_[trii], edgei1);
2685  if ((triEdgei0 + 2) % 3 == triEdgei1)
2686  {
2687  const label triEdgeiLeft = (triEdgei1 + 2) % 3;
2688  leftPolyEdgeis.append(triEdges_[trii][triEdgeiLeft]);
2689  leftPolyTriis.append(trii);
2690  leftPolyPointis.append(triPoint(trii, triEdgei1));
2691  }
2692  else
2693  {
2694  const label triEdgeiRight = (triEdgei1 + 1) % 3;
2695  rightPolyEdgeis.append(triEdges_[trii][triEdgeiRight]);
2696  rightPolyTriis.append(trii);
2697  rightPolyPointis.append(triPoint(trii, triEdgeiRight));
2698  }
2699  }
2700 
2701  // Last triangle
2702  {
2703  const label trii = routeTriis.last();
2704  const label edgei0 = routeEdgeis.last();
2705  const label triEdgei0 = findIndex(triEdges_[trii], edgei0);
2706  const label triEdgeiLeft = (triEdgei0 + 1) % 3;
2707  const label triEdgeiRight = (triEdgei0 + 2) % 3;
2708  leftPolyEdgeis.append(triEdges_[trii][triEdgeiLeft]);
2709  rightPolyEdgeis.append(triEdges_[trii][triEdgeiRight]);
2710  leftPolyTriis.append(trii);
2711  rightPolyTriis.append(trii);
2712  leftPolyPointis.append(pointi1);
2713  rightPolyPointis.append(pointi1);
2714  }
2715  }
2716 
2717  // Disconnect and remove existing tris
2718  auto disconnect = [&]
2719  (
2720  const labelUList& edgeis,
2721  const labelList& triis
2722  )
2723  {
2724  forAll(edgeis, i)
2725  {
2726  const label edgei = edgeis[i];
2727  const label trii = triis[i];
2728 
2729  const label edgeTrii = edgeTris_[edgei][0] != trii;
2730  const label triEdgei = findIndex(triEdges_[trii], edgei);
2731 
2732  const label tempEdgei = newEdgei();
2733  triEdges_[trii][triEdgei] = tempEdgei;
2734  edgeTris_[tempEdgei][0] = trii;
2735  edgeTris_[edgei][edgeTrii] = -1;
2736  }
2737  };
2738  disconnect(leftPolyEdgeis, leftPolyTriis);
2739  disconnect(rightPolyEdgeis, rightPolyTriis);
2740  forAll(routeTriis, routeTrii)
2741  {
2742  removeTri(routeTriis[routeTrii]);
2743  }
2744 
2745  // Create the new conformed edge
2746  const label middleEdgei = newEdgei();
2747  leftPolyEdgeis.append(middleEdgei);
2748  rightPolyEdgeis.append(middleEdgei);
2749 
2750  // Reverse the right polygon
2751  SubList<label> subRightPolyPointis
2752  (
2753  rightPolyPointis,
2754  rightPolyPointis.size() - 1,
2755  1
2756  );
2757  inplaceReverseList(subRightPolyPointis);
2758  inplaceReverseList(rightPolyEdgeis);
2759  inplaceReverseList(rightPolyTriis);
2760 
2761  // Insert the polygons into the triangulation
2762  auto insert = [&]
2763  (
2764  DynamicList<label>& polyPointis,
2765  DynamicList<label>& polyEdgeis
2766  )
2767  {
2768  // Generate new edges
2769  polyEdgeis.resize(2*polyPointis.size() - 3, -1);
2770  for (label i = polyPointis.size(); i < polyEdgeis.size(); ++ i)
2771  {
2772  polyEdgeis[i] = newEdgei();
2773  }
2774 
2775  // Triangulate
2776  polygonTriangulate_.triangulate
2777  (
2778  UIndirectList<point>(points, polyPointis)
2779  );
2780 
2781  // Insert the triangulation
2782  forAll(polygonTriangulate_.triPoints(), trii)
2783  {
2784  addTri
2785  (
2786  polygonTriangulate_.triPoints(trii, polyPointis),
2787  polygonTriangulate_.triEdges(trii, polyEdgeis),
2788  patchFacei,
2789  isSrc
2790  );
2791  }
2792  };
2793  insert(leftPolyPointis, leftPolyEdgeis);
2794  insert(rightPolyPointis, rightPolyEdgeis);
2795  }
2796  }
2797 
2798  // Check
2799  checkPatchFace(patchFacei, isSrc);
2800 
2801  return success;
2802 }
2803 
2804 
2805 template<class SrcPatchType, class TgtPatchType>
2806 bool
2807 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::conformPatchFaceTris
2808 (
2809  const label srcFacei,
2810  const label tgtFacei
2811 )
2812 {
2813  return
2814  conformPatchFaceTris(srcFacei, tgtFacei, true)
2815  && conformPatchFaceTris(tgtFacei, srcFacei, false);
2816 }
2817 
2818 
2819 template<class SrcPatchType, class TgtPatchType>
2820 bool
2821 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::combinePatchFaceTris
2822 (
2823  const label srcFacei,
2824  const label tgtFacei
2825 )
2826 {
2827  // Function to determine whether or not a given triangle is fully
2828  // intersected with the opposite side
2829  auto triIsIntersected = [&]
2830  (
2831  const label trii,
2832  const label otherPatchFacei,
2833  const triFace& otherPatchFacePatchPoints,
2834  const triFace& otherPatchFacePatchEdges,
2835  const bool isSrc
2836  )
2837  {
2838  const labelList& pointOtherPatchPoints =
2839  isSrc ? this->pointTgtPoints_ : this->pointSrcPoints_;
2840  const labelList& pointOtherPatchEdges =
2841  isSrc ? this->pointTgtEdges_ : this->pointSrcEdges_;
2842  const labelList& pointOtherPatchFaces =
2843  isSrc ? this->pointTgtFaces_ : this->pointSrcFaces_;
2844 
2845  forAll(triPoints_[trii], triPointi)
2846  {
2847  const label pointi = triPoint(trii, triPointi);
2848 
2849  if
2850  (
2851  // If the point does not intersect an other patch point ...
2852  findIndex
2853  (
2854  otherPatchFacePatchPoints,
2855  pointOtherPatchPoints[pointi]
2856  ) == -1
2857  // ... and does not intersect an other patch edge ...
2858  && findIndex
2859  (
2860  otherPatchFacePatchEdges,
2861  pointOtherPatchEdges[pointi]
2862  ) == -1
2863  // ... and not intersect the other patch face ...
2864  && pointOtherPatchFaces[pointi] != otherPatchFacei
2865  )
2866  {
2867  // ... then this triangle does not entirely intersect the
2868  // other patch face.
2869  return false;
2870  }
2871  }
2872 
2873  return true;
2874  };
2875 
2876  // Get the polygons which comprise all of the intersected triangles
2877  auto getIntersectionPolygon = [&]
2878  (
2879  const label patchFacei,
2880  const label otherPatchFacei,
2881  const bool isSrc,
2882  DynamicList<label>& triis,
2883  DynamicList<label>& edgeis
2884  )
2885  {
2886  const labelList& patchFaceTris =
2887  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
2888  const labelList& triPatchFaces =
2889  isSrc ? triSrcFace_ : triTgtFace_;
2890 
2891  const triFace otherPatchFacePatchPoints
2892  (
2893  patchFacePatchPoints(otherPatchFacei, !isSrc)
2894  );
2895  const triFace otherPatchFacePatchEdges
2896  (
2897  patchFacePatchEdges(otherPatchFacei, !isSrc)
2898  );
2899 
2900  // Get an initial intersected tri
2901  label trii0 = -1;
2902  forAll(patchFaceTris, patchFaceTrii)
2903  {
2904  const label trii = patchFaceTris[patchFaceTrii];
2905  if
2906  (
2907  triIsIntersected
2908  (
2909  trii,
2910  otherPatchFacei,
2911  otherPatchFacePatchPoints,
2912  otherPatchFacePatchEdges,
2913  isSrc
2914  )
2915  )
2916  {
2917  trii0 = trii;
2918  break;
2919  }
2920  }
2921  if (trii0 == -1) return;
2922 
2923  // Populate the star
2924  star::context starContext = star_.populate
2925  (
2926  trii0,
2927  true,
2928  [&](const label edgei, const label trii)
2929  {
2930  return
2931  triPatchFaces[trii] == patchFacei
2932  && triIsIntersected
2933  (
2934  trii,
2935  otherPatchFacei,
2936  otherPatchFacePatchPoints,
2937  otherPatchFacePatchEdges,
2938  isSrc
2939  );
2940  },
2941  triEdges_,
2942  edgeTris_
2943  );
2944 
2945  // Add triangles to the list
2946  forAllStarFaces(star_, starTrii, trii)
2947  {
2948  triis.append(trii);
2949  }
2950 
2951  // Walk around the star polygon, adding edges to the list
2952  forAllStarEdges(star_, i, starEdgei, edgei)
2953  {
2954  edgeis.append(edgei);
2955  }
2956  };
2957  DynamicList<label> srcPolyTris, srcPolyEdges;
2958  DynamicList<label> tgtPolyTris, tgtPolyEdges;
2959  getIntersectionPolygon
2960  (
2961  srcFacei,
2962  tgtFacei,
2963  true,
2964  srcPolyTris,
2965  srcPolyEdges
2966  );
2967  getIntersectionPolygon
2968  (
2969  tgtFacei,
2970  srcFacei,
2971  false,
2972  tgtPolyTris,
2973  tgtPolyEdges
2974  );
2975 
2976  // Return if there is nothing to do
2977  if (srcPolyEdges.size() == 0 && tgtPolyEdges.size() == 0)
2978  {
2979  return true;
2980  }
2981 
2982  // Align the source and target polygon edges
2983  if (srcPolyEdges.size() && tgtPolyEdges.size())
2984  {
2985  inplaceReverseList(tgtPolyEdges);
2986 
2987  label tgtPolyEdgei0 = 0;
2988  forAll(tgtPolyEdges, tgtPolyEdgei)
2989  {
2990  const edge srcE = edgePoints(srcPolyEdges[0]);
2991  const edge tgtE = edgePoints(tgtPolyEdges[tgtPolyEdgei]);
2992 
2993  if (edge::compare(srcE, tgtE) != 0)
2994  {
2995  tgtPolyEdgei0 = tgtPolyEdgei;
2996  break;
2997  }
2998  }
2999 
3000  inplaceRotateList
3001  (
3002  static_cast<labelList&>(tgtPolyEdges),
3003  - tgtPolyEdgei0
3004  );
3005  }
3006  bool aligned;
3007  if (srcPolyEdges.size() == tgtPolyEdges.size())
3008  {
3009  aligned = true;
3010  forAll(srcPolyEdges, polyEdgei)
3011  {
3012  const edge srcE = edgePoints(srcPolyEdges[polyEdgei]);
3013  const edge tgtE = edgePoints(tgtPolyEdges[polyEdgei]);
3014 
3015  if (edge::compare(srcE, tgtE) == 0)
3016  {
3017  aligned = false;
3018  break;
3019  }
3020  }
3021  }
3022  else
3023  {
3024  aligned = false;
3025  }
3026  if (!aligned)
3027  {
3028  if (this->debug > 1)
3029  {
3030  WarningInFunction
3031  << indent << "Failed to combine intersected parts of source "
3032  << "face #" << srcFacei << " and target face #" << tgtFacei
3033  << endl;
3034  writePatchFace(srcFacei, true);
3035  writePatchFace(tgtFacei, false);
3036  }
3037  return false;
3038  }
3039 
3040  // Mark triangles that belong to either polygon as candidates
3041  forAll(srcPolyTris, srcPolyTrii)
3042  {
3043  const label srcTrii = srcPolyTris[srcPolyTrii];
3044  markedTriTris_.append(srcTrii);
3045  triMarkedTris_[srcTrii] = markedTriTris_.size() - 1;
3046  }
3047  forAll(tgtPolyTris, tgtPolyTrii)
3048  {
3049  const label tgtTrii = tgtPolyTris[tgtPolyTrii];
3050  markedTriTris_.append(tgtTrii);
3051  triMarkedTris_[tgtTrii] = markedTriTris_.size() - 1;
3052  }
3053 
3054  // Change the aligned edges so that they connect opposite sides, thereby
3055  // disconnecting the intersected triangles from the rest of their surfaces.
3056  // Tris in the polygon are connected to the source edge and non-star tris
3057  // are connected to the target edge. If there are no non-star tris
3058  // connected to a pair of aligned edges then the target edge will not be
3059  // connected to any triangles as a result of this operation. However, we do
3060  // not remove the target edge in this case because we need it later to
3061  // attach the new intersected face to.
3062  forAll(srcPolyEdges, polyEdgei)
3063  {
3064  const label srcEdgei = srcPolyEdges[polyEdgei];
3065  const label tgtEdgei = tgtPolyEdges[polyEdgei];
3066 
3067  if (srcEdgei == tgtEdgei) continue;
3068 
3069  const label srcEdgeTrii =
3070  edgeTris_[srcEdgei][0] == -1
3071  || triMarkedTris_[edgeTris_[srcEdgei][0]] == -1
3072  || triSrcFace_[edgeTris_[srcEdgei][0]] == -1;
3073  const label tgtEdgeTrii =
3074  edgeTris_[tgtEdgei][0] == -1
3075  || triMarkedTris_[edgeTris_[tgtEdgei][0]] == -1
3076  || triTgtFace_[edgeTris_[tgtEdgei][0]] == -1;
3077 
3078  const label tgtTrii = edgeTris_[tgtEdgei][tgtEdgeTrii];
3079  const label srcTrij = edgeTris_[srcEdgei][!srcEdgeTrii];
3080  const label tgtTrij = edgeTris_[tgtEdgei][!tgtEdgeTrii];
3081 
3082  edgeTris_[srcEdgei][!srcEdgeTrii] = tgtTrii;
3083  edgeTris_[tgtEdgei][tgtEdgeTrii] = srcTrij;
3084 
3085  const label tgtTriEdgei = findIndex(triEdges_[tgtTrii], tgtEdgei);
3086  triEdges_[tgtTrii][tgtTriEdgei] = srcEdgei;
3087 
3088  if (srcTrij != -1)
3089  {
3090  const label srcTriEdgej = findIndex(triEdges_[srcTrij], srcEdgei);
3091  triEdges_[srcTrij][srcTriEdgej] = tgtEdgei;
3092  }
3093 
3094  if (srcTrij != -1 && tgtTrij != -1)
3095  {
3096  edgeFrontEdges_[tgtEdgei] = frontEdgeEdges_.size();
3097  frontEdgeEdges_.append(tgtEdgei);
3098  }
3099  }
3100 
3101  // Circle the edges to create the intersection face
3102  const label facei = this->faces().size();
3103  this->faces_.append(face(srcPolyEdges.size()));
3104  faceEdges_.append(labelList(srcPolyEdges.size()));
3105  this->srcFaceFaces_[srcFacei].append(facei);
3106  this->tgtFaceFaces_[tgtFacei].append(facei);
3107  this->faceSrcFaces_.append(srcFacei);
3108  this->faceTgtFaces_.append(tgtFacei);
3109  forAll(srcPolyEdges, polyEdgei)
3110  {
3111  const label srcEdgei = srcPolyEdges[polyEdgei];
3112  const label srcEdgeTrii =
3113  edgeTris_[srcEdgei][0] == -1
3114  || triMarkedTris_[edgeTris_[srcEdgei][0]] == -1
3115  || triSrcFace_[edgeTris_[srcEdgei][0]] == -1;
3116 
3117  const label srcTrii = edgeTris_[srcEdgei][srcEdgeTrii];
3118  const label srcTriEdgei = findIndex(triEdges_[srcTrii], srcEdgei);
3119 
3120  this->faces_.last()[polyEdgei] =
3121  triEdgePoints(srcTrii, srcTriEdgei).start();
3122 
3123  const label tgtEdgei = tgtPolyEdges[polyEdgei];
3124 
3125  faceEdges_.last()[polyEdgei] = tgtEdgei;
3126  intersectEdgeFaces_[tgtEdgei][intersectEdgeFaces_[tgtEdgei][0] != -1] =
3127  facei;
3128  }
3129 
3130  // Clear the marked tris
3131  forAll(markedTriTris_, candidateTrii)
3132  {
3133  const label trii = markedTriTris_[candidateTrii];
3134  if (trii != -1)
3135  {
3136  triMarkedTris_[trii] = -1;
3137  }
3138  }
3139  markedTriTris_.clear();
3140 
3141  // Check
3142  checkPatchFace(srcFacei, true);
3143  checkPatchFace(tgtFacei, false);
3144 
3145  // Remove the intersection triangles
3146  forAll(srcPolyTris, srcPolyTrii)
3147  {
3148  removeTri(srcPolyTris[srcPolyTrii]);
3149  }
3150  forAll(tgtPolyTris, tgtPolyTrii)
3151  {
3152  removeTri(tgtPolyTris[tgtPolyTrii]);
3153  }
3154 
3155  // Check
3156  checkPatchFace(srcFacei, true);
3157  checkPatchFace(tgtFacei, false);
3158 
3159  // Update the propagation front
3160  forAll(tgtPolyEdges, polyEdgei)
3161  {
3162  const label edgei = tgtPolyEdges[polyEdgei];
3163 
3164  const label trii0 = edgeTris_[edgei][0];
3165  const label trii1 = edgeTris_[edgei][1];
3166 
3167  const bool isFront =
3168  trii0 != -1
3169  && trii1 != -1
3170  && (triSrcFace_[trii0] == -1) != (triSrcFace_[trii1] == -1);
3171 
3172  if (isFront && edgeFrontEdges_[edgei] == -1)
3173  {
3174  edgeFrontEdges_[edgei] = frontEdgeEdges_.size();
3175  frontEdgeEdges_.append(edgei);
3176  }
3177 
3178  if (!isFront && edgeFrontEdges_[edgei] != -1)
3179  {
3180  frontEdgeEdges_[edgeFrontEdges_[edgei]] = -1;
3181  edgeFrontEdges_[edgei] = -1;
3182  }
3183  }
3184 
3185  // Check
3186  checkPatchFace(srcFacei, true);
3187  checkPatchFace(tgtFacei, false);
3188 
3189  return true;
3190 }
3191 
3192 
3193 template<class SrcPatchType, class TgtPatchType>
3194 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::initialise
3195 (
3196  const vectorField& srcPointNormals
3197 )
3198 {
3199  // Clear the base class data ...
3200 
3201  this->points_.clear();
3202 
3203  this->srcPointPoints_ = -1;
3204  this->tgtPointPoints_ = -1;
3205  this->pointSrcPoints_.clear();
3206  this->pointTgtPoints_.clear();
3207 
3208  forAll(this->srcEdgePoints_, srcEdgei)
3209  {
3210  this->srcEdgePoints_[srcEdgei].clear();
3211  }
3212  forAll(this->tgtEdgePoints_, tgtEdgei)
3213  {
3214  this->tgtEdgePoints_[tgtEdgei].clear();
3215  }
3216  this->pointSrcEdges_.clear();
3217  this->pointTgtEdges_.clear();
3218 
3219  this->pointSrcFaces_.clear();
3220  this->pointTgtFaces_.clear();
3221 
3222  this->faces_.clear();
3223 
3224  forAll(this->srcFaceFaces_, srcFacei)
3225  {
3226  this->srcFaceFaces_[srcFacei].clear();
3227  }
3228  forAll(this->tgtFaceFaces_, tgtFacei)
3229  {
3230  this->tgtFaceFaces_[tgtFacei].clear();
3231  }
3232  this->faceSrcFaces_.clear();
3233  this->faceTgtFaces_.clear();
3234 
3235  // Initialise with the source and target patch data ...
3236 
3237  // Points
3238  const label nPoints =
3239  this->srcPatch_.nPoints() + this->tgtPatch_.nPoints();
3240  srcPoints_.resize(nPoints, point::uniform(NaN));
3241  srcPointNormals_.resize(nPoints, vector::uniform(NaN));
3242  tgtPoints_.resize(nPoints, point::uniform(NaN));
3243  pointPoints_.resize(nPoints, -1);
3244  this->pointSrcPoints_.resize(nPoints, -1);
3245  this->pointTgtPoints_.resize(nPoints, -1);
3246  this->pointSrcEdges_.resize(nPoints, -1);
3247  this->pointTgtEdges_.resize(nPoints, -1);
3248  this->pointSrcFaces_.resize(nPoints, -1);
3249  this->pointTgtFaces_.resize(nPoints, -1);
3250  forAll(this->srcPatch_.localPoints(), srcPointi)
3251  {
3252  const label pointi = srcPointi;
3253 
3254  srcPoints_[pointi] = this->srcPatch_.localPoints()[srcPointi];
3255  tgtPoints_[pointi] = this->srcPatch_.localPoints()[srcPointi];
3256  srcPointNormals_[pointi] = srcPointNormals[srcPointi];
3257  pointPoints_[pointi] = pointi;
3258  this->srcPointPoints_[srcPointi] = pointi;
3259  this->pointSrcPoints_[pointi] = srcPointi;
3260  }
3261  forAll(this->tgtPatch_.localPoints(), tgtPointi)
3262  {
3263  const label pointi = this->srcPatch_.nPoints() + tgtPointi;
3264 
3265  srcPoints_[pointi] = this->tgtPatch_.localPoints()[tgtPointi];
3266  tgtPoints_[pointi] = this->tgtPatch_.localPoints()[tgtPointi];
3267  pointPoints_[pointi] = pointi;
3268  this->tgtPointPoints_[tgtPointi] = pointi;
3269  this->pointTgtPoints_[pointi] = tgtPointi;
3270  }
3271 
3272  // Edges
3273  const label nEdges = this->srcPatch_.nEdges() + this->tgtPatch_.nEdges();
3274  edgeTris_.resize(nEdges, labelPair(-1, -1));
3275  intersectEdgeFaces_.resize(nEdges, labelPair(-1, -1));
3276  forAll(this->srcPatch_.faceEdges(), srcFacei)
3277  {
3278  const label trii = srcFacei;
3279 
3280  forAll(this->srcPatch_.faceEdges()[srcFacei], srcFaceEdgei)
3281  {
3282  const label srcEdgei =
3283  this->srcPatch_.faceEdges()[srcFacei][srcFaceEdgei];
3284  const label edgei = srcEdgei;
3285 
3286  const edge& e = this->srcPatch_.edges()[srcEdgei];
3287  const edge fe =
3288  this->srcPatch_.localFaces()[srcFacei].faceEdge(srcFaceEdgei);
3289 
3290  edgeTris_[edgei][edge::compare(e, fe) < 0] = trii;
3291  }
3292  }
3293  forAll(this->srcPatch_.edges(), srcEdgei)
3294  {
3295  const edge& e = this->srcPatch_.edges()[srcEdgei];
3296 
3297  this->srcEdgePoints_[srcEdgei].append(this->srcPointPoints_[e[0]]);
3298  this->srcEdgePoints_[srcEdgei].append(this->srcPointPoints_[e[1]]);
3299  }
3300  forAll(this->tgtPatch_.faceEdges(), tgtFacei)
3301  {
3302  const label trii = this->srcPatch_.size() + tgtFacei;
3303 
3304  forAll(this->tgtPatch_.faceEdges()[tgtFacei], tgtFaceEdgei)
3305  {
3306  const label tgtEdgei =
3307  this->tgtPatch_.faceEdges()[tgtFacei][tgtFaceEdgei];
3308  const label edgei = this->srcPatch_.nEdges() + tgtEdgei;
3309 
3310  const edge& e = this->tgtPatch_.edges()[tgtEdgei];
3311  const edge fe =
3312  this->tgtPatch_.localFaces()[tgtFacei].faceEdge(tgtFaceEdgei);
3313 
3314  edgeTris_[edgei][edge::compare(e, fe) < 0] = trii;
3315  }
3316  }
3317  forAll(this->tgtPatch_.edges(), tgtEdgei)
3318  {
3319  const edge& e = this->tgtPatch_.edges()[tgtEdgei];
3320 
3321  this->tgtEdgePoints_[tgtEdgei].append(this->tgtPointPoints_[e[0]]);
3322  this->tgtEdgePoints_[tgtEdgei].append(this->tgtPointPoints_[e[1]]);
3323  }
3324 
3325  // Tris
3326  const label nTris = this->srcPatch_.size() + this->tgtPatch_.size();
3327  triPoints_.resize(nTris);
3328  triEdges_.resize(nTris);
3329  triSrcFace_.resize(nTris, -1);
3330  triTgtFace_.resize(nTris, -1);
3331  forAll(this->srcPatch_.localFaces(), srcFacei)
3332  {
3333  const label trii = srcFacei;
3334 
3335  forAll(this->srcPatch_.localFaces()[srcFacei], i)
3336  {
3337  triPoints_[trii][i] = this->srcPatch_.localFaces()[srcFacei][i];
3338  triEdges_[trii][i] = this->srcPatch_.faceEdges()[srcFacei][i];
3339  }
3340  triSrcFace_[trii] = srcFacei;
3341  srcFaceTris_[srcFacei].resize(1, trii);
3342  }
3343  forAll(this->tgtPatch_.localFaces(), tgtFacei)
3344  {
3345  const label trii = this->srcPatch_.size() + tgtFacei;
3346 
3347  forAll(this->tgtPatch_.localFaces()[tgtFacei], i)
3348  {
3349  triPoints_[trii][i] =
3350  this->srcPatch_.nPoints()
3351  + this->tgtPatch_.localFaces()[tgtFacei][i];
3352  triEdges_[trii][i] =
3353  this->srcPatch_.nEdges()
3354  + this->tgtPatch_.faceEdges()[tgtFacei][i];
3355  }
3356  triTgtFace_[trii] = tgtFacei;
3357  tgtFaceTris_[tgtFacei].resize(1, trii);
3358  }
3359 
3360  // Removal
3361  removedEdges_.clear();
3362  removedTris_.clear();
3363 
3364  // Front propagation
3365  frontEdgeEdges_.clear();
3366  edgeFrontEdges_ = DynamicList<label>(edgeTris_.size(), -1);
3367 
3368  // Insertion
3369  candidateTriTris_.clear();
3370  triCandidateTris_ = DynamicList<label>(triPoints_.size(), -1);
3371 
3372  // Marking
3373  markedTriTris_.clear();
3374  triMarkedTris_ = DynamicList<label>(triPoints_.size(), -1);
3375 
3376  checkPatchFaces(true);
3377  checkPatchFaces(false);
3378 }
3379 
3380 
3381 template<class SrcPatchType, class TgtPatchType>
3382 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::clean()
3383 {
3384  // Resolve point-points
3385  inplaceRenumber(pointPoints_, triPoints_);
3386  inplaceRenumber(pointPoints_, this->srcPointPoints_);
3387  inplaceRenumber(pointPoints_, this->tgtPointPoints_);
3388  inplaceRenumber(pointPoints_, this->srcEdgePoints_);
3389  inplaceRenumber(pointPoints_, this->tgtEdgePoints_);
3390 
3391  // Remove points by shuffling up
3392  labelList oldPointNewPoints(pointPoints_.size(), -1);
3393  label pointi = 0;
3394  forAll(pointPoints_, pointj)
3395  {
3396  if (pointPoints_[pointj] == pointj)
3397  {
3398  oldPointNewPoints[pointj] = pointi;
3399  srcPoints_[pointi] = srcPoints_[pointj];
3400  srcPointNormals_[pointi] = srcPointNormals_[pointj];
3401  tgtPoints_[pointi] = tgtPoints_[pointj];
3402  pointPoints_[pointi] = pointi;
3403  this->pointSrcPoints_[pointi] = this->pointSrcPoints_[pointj];
3404  this->pointTgtPoints_[pointi] = this->pointTgtPoints_[pointj];
3405  this->pointSrcEdges_[pointi] = this->pointSrcEdges_[pointj];
3406  this->pointTgtEdges_[pointi] = this->pointTgtEdges_[pointj];
3407  this->pointSrcFaces_[pointi] = this->pointSrcFaces_[pointj];
3408  this->pointTgtFaces_[pointi] = this->pointTgtFaces_[pointj];
3409  ++ pointi;
3410  }
3411  }
3412  srcPoints_.resize(pointi);
3413  srcPointNormals_.resize(pointi);
3414  tgtPoints_.resize(pointi);
3415  pointPoints_.resize(pointi);
3416  this->pointSrcPoints_.resize(pointi);
3417  this->pointTgtPoints_.resize(pointi);
3418  this->pointSrcEdges_.resize(pointi);
3419  this->pointTgtEdges_.resize(pointi);
3420  this->pointSrcFaces_.resize(pointi);
3421  this->pointTgtFaces_.resize(pointi);
3422 
3423  // Remove edges by shuffling up
3424  labelList oldEdgeNewEdges(edgeTris_.size(), -1);
3425  label edgei = 0;
3426  forAll(edgeTris_, edgej)
3427  {
3428  if
3429  (
3430  edgeTris_[edgej] != labelPair(-1, -1)
3431  || intersectEdgeFaces_[edgej] != labelPair(-1, -1)
3432  )
3433  {
3434  oldEdgeNewEdges[edgej] = edgei;
3435  edgeTris_[edgei] = edgeTris_[edgej];
3436  intersectEdgeFaces_[edgei] = intersectEdgeFaces_[edgej];
3437  edgeFrontEdges_[edgei] = edgeFrontEdges_[edgej];
3438  ++ edgei;
3439  }
3440  }
3441  edgeTris_.resize(edgei);
3442  intersectEdgeFaces_.resize(edgei);
3443  edgeFrontEdges_.resize(edgei);
3444 
3445  // Remove tris by shuffling up
3446  labelList oldTriNewTris(triPoints_.size(), -1);
3447  label trii = 0;
3448  forAll(triPoints_, trij)
3449  {
3450  if (triPoints_[trij] != FixedList<label, 3>({-1, -1, -1}))
3451  {
3452  oldTriNewTris[trij] = trii;
3453  triPoints_[trii] = triPoints_[trij];
3454  triEdges_[trii] = triEdges_[trij];
3455  triSrcFace_[trii] = triSrcFace_[trij];
3456  triTgtFace_[trii] = triTgtFace_[trij];
3457  ++ trii;
3458  }
3459  }
3460  triPoints_.resize(trii);
3461  triEdges_.resize(trii);
3462  triSrcFace_.resize(trii);
3463  triTgtFace_.resize(trii);
3464 
3465  // Map
3466  inplaceRenumber(oldTriNewTris, edgeTris_);
3467  inplaceRenumber(oldPointNewPoints, triPoints_);
3468  inplaceRenumber(oldEdgeNewEdges, triEdges_);
3469  inplaceRenumber(oldTriNewTris, srcFaceTris_);
3470  inplaceRenumber(oldTriNewTris, tgtFaceTris_);
3471  inplaceRenumber(oldEdgeNewEdges, frontEdgeEdges_);
3472  inplaceRenumber(oldPointNewPoints, this->faces_);
3473  inplaceRenumber(oldEdgeNewEdges, faceEdges_);
3474  inplaceRenumber(oldPointNewPoints, this->srcPointPoints_);
3475  inplaceRenumber(oldPointNewPoints, this->tgtPointPoints_);
3476  inplaceRenumber(oldPointNewPoints, this->srcEdgePoints_);
3477  inplaceRenumber(oldPointNewPoints, this->tgtEdgePoints_);
3478 
3479  // Removal
3480  removedEdges_.clear();
3481  removedTris_.clear();
3482 
3483  // Insertion
3484  candidateTriTris_.clear();
3485  triCandidateTris_ = DynamicList<label>(triPoints_.size(), -1);
3486 
3487  // Marking
3488  markedTriTris_.clear();
3489  triMarkedTris_ = DynamicList<label>(triPoints_.size(), -1);
3490 
3491  checkPatchFaces(true);
3492  checkPatchFaces(false);
3493 }
3494 
3495 
3496 template<class SrcPatchType, class TgtPatchType>
3497 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::finalise()
3498 {
3499  clean();
3500 
3501  // Add the remaining triangles as uncoupled faces ...
3502 
3503  const label nFaces = this->faces_.size();
3504 
3505  this->faces_.resize(nFaces + triPoints_.size());
3506  faceEdges_.resize(nFaces + triPoints_.size());
3507  this->faceSrcFaces_.resize(nFaces + triPoints_.size());
3508  this->faceTgtFaces_.resize(nFaces + triPoints_.size());
3509 
3510  forAll(triPoints_, trii)
3511  {
3512  const label facei = nFaces + trii;
3513 
3514  if (triSrcFace_[trii] != -1)
3515  {
3516  this->faces_[facei] = triPoints(trii);
3517  faceEdges_[facei] = labelList(triEdges_[trii]);
3518  this->srcFaceFaces_[triSrcFace_[trii]].append(facei);
3519  this->faceSrcFaces_[facei] = triSrcFace_[trii];
3520  this->faceTgtFaces_[facei] = -1;
3521  }
3522  else
3523  {
3524  this->faces_[facei] = triPoints(trii).reverseFace();
3525  faceEdges_[facei] = labelList(reverseList(triEdges_[trii]));
3526  this->tgtFaceFaces_[triTgtFace_[trii]].append(facei);
3527  this->faceSrcFaces_[facei] = -1;
3528  this->faceTgtFaces_[facei] = triTgtFace_[trii];
3529  }
3530  }
3531 
3532  nonIntersectEdgeFaces_.resize(edgeTris_.size());
3533  forAll(edgeTris_, edgei)
3534  {
3535  forAll(edgeTris_[edgei], edgeTrii)
3536  {
3537  nonIntersectEdgeFaces_[edgei][edgeTrii] =
3538  edgeTris_[edgei][edgeTrii] + nFaces;
3539  }
3540  }
3541 
3542  checkPatchFaces(true);
3543  checkPatchFaces(false);
3544 }
3545 
3546 
3547 template<class SrcPatchType, class TgtPatchType>
3548 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::unFinalise()
3549 {
3550  clean();
3551 
3552  // Remove all uncoupled faces and store them as triangles ...
3553 
3554  label nFaces = 0;
3555  while
3556  (
3557  nFaces < this->faces_.size()
3558  && this->faceSrcFaces_[nFaces] != -1
3559  && this->faceTgtFaces_[nFaces] != -1
3560  )
3561  {
3562  ++ nFaces;
3563  }
3564 
3565  this->faces_.resize(nFaces);
3566  faceEdges_.resize(nFaces);
3567  this->faceSrcFaces_.resize(nFaces);
3568  this->faceTgtFaces_.resize(nFaces);
3569 
3570  forAll(triPoints_, trii)
3571  {
3572  DynamicList<label>& faceis =
3573  triSrcFace_[trii] != -1
3574  ? this->srcFaceFaces_[triSrcFace_[trii]]
3575  : this->tgtFaceFaces_[triTgtFace_[trii]];
3576 
3577  label n = 0;
3578  while (n < faceis.size() && faceis[n] < nFaces)
3579  {
3580  ++ n;
3581  }
3582 
3583  faceis.resize(n);
3584  }
3585 
3586  nonIntersectEdgeFaces_.clear();
3587 
3588  checkPatchFaces(true);
3589  checkPatchFaces(false);
3590 }
3591 
3592 
3593 template<class SrcPatchType, class TgtPatchType>
3594 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::write()
3595 {
3596  if (this->debug > 2)
3597  {
3598  finalise();
3599 
3600  // Use base class to write the patch
3601  this->report(name(writei_));
3602 
3603  unFinalise();
3604 
3605  // Write the edge front
3606  const fileName frontFileName =
3607  type() + "_front_" + name(writei_) + ".vtk";
3608  Info<< indent << "Writing front to " << frontFileName << endl;
3609 
3610  DynamicList<point> writePoints(frontEdgeEdges_.size()*2);
3611  DynamicList<labelPair> writeLines(frontEdgeEdges_.size()*2);
3612  forAll(frontEdgeEdges_, frontEdgei)
3613  {
3614  const label edgei = frontEdgeEdges_[frontEdgei];
3615 
3616  if (edgei == -1) continue;
3617 
3618  const edge e = edgePoints(edgei);
3619  writePoints.append(tgtPoints_[e.start()]);
3620  writePoints.append(tgtPoints_[e.end()]);
3621 
3622  const label i = writePoints.size() - 2;
3623  writeLines.append(labelPair(i, i + 1));
3624  }
3625 
3626  vtkWritePolyData::write
3627  (
3628  frontFileName,
3629  "front",
3630  false,
3631  writePoints,
3632  labelList(),
3633  writeLines,
3634  faceList()
3635  );
3636 
3637  writei_ ++;
3638  }
3639 }
3640 
3641 
3642 template<class SrcPatchType, class TgtPatchType>
3643 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::writePatchFace
3644 (
3645  const label patchFacei,
3646  const bool isSrc
3647 ) const
3648 {
3649  OFstream os
3650  (
3651  word(isSrc ? "src" : "tgt") + "Face_" + name(patchFacei) + ".obj"
3652  );
3653 
3654  OFstream tos
3655  (
3656  word(isSrc ? "src" : "tgt") + "FaceTris_" + name(patchFacei) + ".obj"
3657  );
3658 
3659  Info<< indent << "Writing patch face to " << os.name()
3660  << " and patch face triangulation to " << tos.name()
3661  << incrIndent << endl;
3662 
3663  forAll(patchFacePatchPoints(patchFacei, isSrc), patchFacePatchPointi)
3664  {
3665  const label patchPointi =
3666  patchFacePatchPoints(patchFacei, isSrc)[patchFacePatchPointi];
3667  const point& p =
3668  isSrc
3669  ? this->srcPatch_.localPoints()[patchPointi]
3670  : this->tgtPatch_.localPoints()[patchPointi];
3671  os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
3672  }
3673  os << "f 1 2 3" << nl;
3674 
3675  const labelList& patchFaceTris =
3676  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
3677 
3678  forAll(patchFaceTris, patchFaceTrii)
3679  {
3680  const label trii = patchFaceTris[patchFaceTrii];
3681 
3682  Info<< indent << "tri #" << trii << " points=" << triPoints(trii)
3683  << " edges=" << triEdges_[trii] << endl;
3684 
3685  forAll(triPoints_[trii], triPointi)
3686  {
3687  const label pointi = triPoint(trii, triPointi);
3688  const point& p = isSrc ? srcPoints_[pointi] : tgtPoints_[pointi];
3689  tos << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
3690  }
3691  tos << "f";
3692  forAll(triPoints_[trii], triPointi)
3693  {
3694  tos << " " << 1 + 3*patchFaceTrii + triPointi;
3695  }
3696  tos << nl;
3697  }
3698 
3699  Info<< decrIndent;
3700 }
3701 
3702 
3703 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
3704 
3705 template<class SrcPatchType, class TgtPatchType>
3706 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::TriPatchIntersection
3707 (
3708  const SrcPatchType& srcPatch,
3709  const TgtPatchType& tgtPatch,
3710  const scalar snapTol
3711 )
3712 :
3713  TriPatchIntersection(srcPatch, srcPatch.pointNormals(), tgtPatch, snapTol)
3714 {}
3715 
3716 
3717 template<class SrcPatchType, class TgtPatchType>
3718 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::TriPatchIntersection
3719 (
3720  const SrcPatchType& srcPatch,
3721  const vectorField& srcPointNormals,
3722  const TgtPatchType& tgtPatch,
3723  const scalar snapTol
3724 )
3725 :
3726  PatchIntersection<SrcPatchType, TgtPatchType>(srcPatch, tgtPatch),
3727 
3728  srcPoints_(),
3729  srcPointNormals_(),
3730  tgtPoints_(this->points_),
3731  pointPoints_(),
3732 
3733  edgeTris_(),
3734  intersectEdgeFaces_(),
3735  nonIntersectEdgeFaces_(),
3736 
3737  triPoints_(),
3738  triEdges_(),
3739  triSrcFace_(),
3740  triTgtFace_(),
3741  srcFaceTris_(srcPatch.size()),
3742  tgtFaceTris_(tgtPatch.size()),
3743 
3744  faceEdges_(),
3745 
3746  removedTris_(),
3747  removedEdges_(),
3748 
3749  frontEdgeEdges_(),
3750  edgeFrontEdges_(),
3751 
3752  candidateTriTris_(),
3753  triCandidateTris_(),
3754 
3755  markedTriTris_(),
3756  triMarkedTris_(),
3757 
3758  polygonTriangulate_(),
3759 
3760  star_(),
3761 
3762  writei_(0)
3763 {
3764  cpuTime time;
3765 
3766  Info<< indent << type() << ": Intersecting "
3767  << this->srcPatch_.size() << " source tri faces and "
3768  << this->tgtPatch_.size() << " target tri faces" << incrIndent << endl;
3769 
3770  if (this->debug)
3771  {
3772  Info<< indent << "Writing tri patches" << incrIndent << endl;
3773  const fileName srcFileName = type() + "_srcPatch.vtk";
3774  Info<< indent << "Writing patch to " << srcFileName << endl;
3775  vtkWritePolyData::write
3776  (
3777  srcFileName,
3778  "source",
3779  false,
3780  this->srcPatch_.localPoints(),
3781  labelList(),
3782  labelListList(),
3783  this->srcPatch_.localFaces(),
3784  "normals",
3785  true,
3786  srcPointNormals
3787  );
3788  const fileName tgtFileName = type() + "_tgtPatch.vtk";
3789  Info<< indent << "Writing patch to " << tgtFileName << endl;
3790  vtkWritePolyData::write
3791  (
3792  tgtFileName,
3793  "target",
3794  false,
3795  this->tgtPatch_.localPoints(),
3796  labelList(),
3797  labelListList(),
3798  this->tgtPatch_.localFaces()
3799  );
3800  Info<< decrIndent;
3801  }
3802 
3803  // Create bound spheres for patch faces for proximity testing
3804  List<boundSphere> srcFaceSpheres(this->srcPatch_.size());
3805  List<boundSphere> tgtFaceSpheres(this->tgtPatch_.size());
3806  forAll(this->srcPatch_, srcFacei)
3807  {
3808  const triFace& srcFace = this->srcPatch_.localFaces()[srcFacei];
3809  srcFaceSpheres[srcFacei] =
3810  boundSphere::trivial
3811  (
3812  this->srcPatch_.localPoints(),
3813  {srcFace[0], srcFace[1], srcFace[2], -1},
3814  3
3815  );
3816  srcFaceSpheres[srcFacei].inflate(snapTol);
3817  }
3818  forAll(this->tgtPatch_, tgtFacei)
3819  {
3820  const triFace& tgtFace = this->tgtPatch_.localFaces()[tgtFacei];
3821  tgtFaceSpheres[tgtFacei] =
3822  boundSphere::trivial
3823  (
3824  this->tgtPatch_.localPoints(),
3825  {tgtFace[0], tgtFace[1], tgtFace[2], -1},
3826  3
3827  );
3828  tgtFaceSpheres[tgtFacei].inflate(snapTol);
3829  }
3830 
3831  // Construct table to store what faces have been snapped
3832  HashSet<labelPair, labelPair::Hash<>> srcFaceTgtFaceSnaps
3833  (
3834  12*(this->srcPatch_.size() + this->tgtPatch_.size())
3835  );
3836 
3837  // Count the number of successes and failures
3838  label nIntersections = 0, nIntersectionFailures = 0;
3839 
3840  // Function for intersecting two patch faces
3841  auto intersect = [&](const label srcFacei, const label tgtFacei)
3842  {
3843  // Single snapping stage
3844  /*
3845  snapPatchFaceTris(srcFacei, tgtFacei, snapTol);
3846  */
3847 
3848  // Propagate out and snap everything in advance of the intersections
3849  const List<triFace>& srcLocalFaces = this->srcPatch_.localFaces();
3850  const List<triFace>& tgtLocalFaces = this->tgtPatch_.localFaces();
3851  const labelListList& srcPointFaces = this->srcPatch_.pointFaces();
3852  const labelListList& tgtPointFaces = this->tgtPatch_.pointFaces();
3853  forAll(srcLocalFaces[srcFacei], srcFacePointi)
3854  {
3855  const label srcPointi =
3856  srcLocalFaces[srcFacei][srcFacePointi];
3857 
3858  forAll(srcPointFaces[srcPointi], srcPointFacei)
3859  {
3860  const label srcFacej =
3861  srcPointFaces[srcPointi][srcPointFacei];
3862 
3863  forAll(tgtLocalFaces[tgtFacei], tgtFacePointi)
3864  {
3865  const label tgtPointi =
3866  tgtLocalFaces[tgtFacei][tgtFacePointi];
3867 
3868  forAll(tgtPointFaces[tgtPointi], tgtPointFacei)
3869  {
3870  const label tgtFacej =
3871  tgtPointFaces[tgtPointi][tgtPointFacei];
3872 
3873  if
3874  (
3875  boundSphere::overlap
3876  (
3877  srcFaceSpheres[srcFacej],
3878  tgtFaceSpheres[tgtFacej]
3879  )
3880  && !srcFaceTgtFaceSnaps.found({srcFacej, tgtFacej})
3881  )
3882  {
3883  srcFaceTgtFaceSnaps.insert({srcFacej, tgtFacej});
3884  snapPatchFaceTris(srcFacej, tgtFacej, snapTol);
3885  }
3886  }
3887  }
3888  }
3889  }
3890 
3891  intersectPatchFaceTris(srcFacei, tgtFacei);
3892 
3893  write();
3894 
3895  const bool conformFailure = !conformPatchFaceTris(srcFacei, tgtFacei);
3896 
3897  const bool combineFailure = !combinePatchFaceTris(srcFacei, tgtFacei);
3898 
3899  nIntersections ++;
3900  nIntersectionFailures += conformFailure || combineFailure;
3901 
3902  write();
3903  };
3904 
3905  // Target search tree. Used to find an initial source and target face to
3906  // intersect. Once the first intersection has been done the rest follow via
3907  // a front-propagation algorithm, without reference to this tree.
3908  indexedOctree<treeDataPrimitivePatch<TgtPatchType>> tgtTree
3909  (
3910  treeDataPrimitivePatch<TgtPatchType>
3911  (
3912  false,
3913  this->tgtPatch_,
3914  indexedOctree<treeDataPrimitivePatch<TgtPatchType>>::perturbTol()
3915  ),
3916  treeBoundBox(this->tgtPatch_.localPoints()).extend(1e-4),
3917  8,
3918  10,
3919  3
3920  );
3921 
3922  // Populate local data from the source and target patches
3923  initialise(srcPointNormals);
3924  write();
3925 
3926  // Loop the source points, looking for ones that have not been intersected
3927  forAll(this->points_, pointi)
3928  {
3929  // Get the next source patch point
3930  const label srcPointi = this->pointSrcPoints_[pointi];
3931 
3932  // Continue if this point is already intersected
3933  if
3934  (
3935  srcPointi == -1
3936  || this->pointTgtPoints_[pointi] != -1
3937  || this->pointTgtEdges_[pointi] != -1
3938  || this->pointTgtFaces_[pointi] != -1
3939  ) continue;
3940 
3941  // Get a length scale for this point by averaging the lengths of the
3942  // connected edges
3943  const point& srcP = srcPoints_[pointi];
3944  scalar srcL = 0;
3945  forAll(this->srcPatch_.pointEdges()[srcPointi], srcPointEdgei)
3946  {
3947  const label srcEdgei =
3948  this->srcPatch_.pointEdges()[srcPointi][srcPointEdgei];
3949  srcL +=
3950  this->srcPatch_.edges()[srcEdgei].mag
3951  (
3952  this->srcPatch_.localPoints()
3953  );
3954  }
3955  srcL /= this->srcPatch_.pointEdges()[srcPointi].size();
3956 
3957  // Find all the target triangles that are within the length scale of
3958  // the source point
3959  const labelList tgtFaceis = tgtTree.findSphere(srcP, sqr(srcL));
3960 
3961  // Continue if nothing was found
3962  if (tgtFaceis.empty()) continue;
3963 
3964  // Loop all the potential source/target face intersections until an
3965  // edge front is generated
3966  forAll(this->srcPatch_.pointFaces()[srcPointi], srcPointFacei)
3967  {
3968  const label srcFacei =
3969  this->srcPatch_.pointFaces()[srcPointi][srcPointFacei];
3970 
3971  forAll(tgtFaceis, i)
3972  {
3973  intersect(srcFacei, tgtFaceis[i]);
3974 
3975  if (frontEdgeEdges_.size()) break;
3976  }
3977 
3978  if (frontEdgeEdges_.size()) break;
3979  }
3980 
3981  // Propagate until the edge front is empty
3982  while (frontEdgeEdges_.size())
3983  {
3984  const label edgei = frontEdgeEdges_.remove();
3985 
3986  if (edgei == -1) continue;
3987 
3988  edgeFrontEdges_[edgei] = -1;
3989 
3990  label srcTrii = edgeTris_[edgei][0];
3991  label tgtTrii = edgeTris_[edgei][1];
3992 
3993  if (triSrcFace_[srcTrii] == -1)
3994  {
3995  Swap(srcTrii, tgtTrii);
3996  }
3997 
3998  intersect(triSrcFace_[srcTrii], triTgtFace_[tgtTrii]);
3999  }
4000  }
4001 
4002  // Populate data in the base class which was not generated as part of the
4003  // intersection process
4004  finalise();
4005 
4006  this->report();
4007 
4008  // Warn about any failures
4009  if (nIntersectionFailures)
4010  {
4011  Info<< indent << "*** Topology could not be generated in "
4012  << nIntersectionFailures << "/" << nIntersections << " cases"
4013  << endl << indent << " The intersection may be incomplete"
4014  << endl;
4015  }
4016 
4017  Info<< indent << this->faces_.size() << " faces generated in "
4018  << time.cpuTimeIncrement() << 's' << endl;
4019 
4020  Info<< decrIndent;
4021 }
4022 
4023 
4024 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
4025 
4026 template<class SrcPatchType, class TgtPatchType>
4027 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::
4028 ~TriPatchIntersection()
4029 {}
4030 
4031 
4032 // ************************************************************************* //
success
bool success
Definition: LISASMDCalcMethod1.H:16
y
scalar y
Definition: LISASMDCalcMethod1.H:14
n
label n
Definition: TABSMDCalcMethod2.H:31
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:449
Foam::DynamicList::remove
T remove()
Remove and return the top element.
Definition: DynamicListI.H:351
Foam::HashSet
A HashTable with keys but without contents.
Definition: HashSet.H:62
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:109
Foam::HashTable::found
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:138
Foam::IOstream::name
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.H:297
Foam::List::size
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Foam::PatchIntersection
Base class for patch intersections. Provides storage and access to the intersection points and faces ...
Definition: PatchIntersection.H:57
Foam::PatchIntersection::srcPatch_
const SrcPatchType & srcPatch_
Reference to the source patch.
Definition: PatchIntersection.H:63
Foam::PatchIntersection::tgtPatch_
const TgtPatchType & tgtPatch_
Reference to the target patch.
Definition: PatchIntersection.H:66
Foam::TriPatchIntersection
Patch intersection based on triangular faces. Intersects and combines two triangulated patches increm...
Definition: TriPatchIntersection.H:58
Foam::TriPatchIntersection::~TriPatchIntersection
virtual ~TriPatchIntersection()
Destructor.
Definition: TriPatchIntersection.C:4028
Foam::TriPatchIntersection::TriPatchIntersection
TriPatchIntersection(const SrcPatchType &srcPatch, const TgtPatchType &tgtPatch, const scalar snapTol)
Construct from a source and a target patch.
Definition: TriPatchIntersection.C:3707
Foam::TriPatchIntersection::faceEdges
const DynamicList< labelList > & faceEdges() const
...
Definition: TriPatchIntersection.H:450
Foam::UList::empty
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
Foam::boundSphere::overlap
static bool overlap(const boundSphere &a, const boundSphere &b)
Return whether two spheres overlap.
Definition: boundSphereI.H:183
Foam::boundSphere::trivial
static boundSphere trivial(const PointList &ps, const FixedList< label, 4 > &pis, const label nPs)
Return the sphere bounding the given set of up to four points.
Foam::cpuTime
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTime.H:55
Foam::cpuTime::cpuTimeIncrement
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:60
Foam::edge
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
Foam::fileName
A class for handling file names.
Definition: fileName.H:82
Foam::fileName::name
word name() const
Return file name (part beyond last /)
Definition: fileName.C:195
Foam::indexedOctree
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:72
Foam::patchIntersection::pointTgtFaces_
DynamicList< label > pointTgtFaces_
The intersection points' target faces, or -1 if the point.
Definition: patchIntersection.H:108
Foam::patchIntersection::points_
DynamicField< point > points_
The intersection points.
Definition: patchIntersection.H:66
Foam::patchIntersection::pointTgtPoints_
DynamicList< label > pointTgtPoints_
The intersection points' corresponding target points, or -1.
Definition: patchIntersection.H:80
Foam::patchIntersection::report
void report(const word &writeSuffix=word::null)
Report properties of the intersection process.
Definition: patchIntersection.C:42
Foam::patchIntersection::pointSrcPoints_
DynamicList< label > pointSrcPoints_
The intersection points' corresponding source points, or -1.
Definition: patchIntersection.H:76
Foam::patchIntersection::faces_
DynamicList< face > faces_
The intersection faces.
Definition: patchIntersection.H:114
Foam::patchIntersection::pointTgtEdges_
DynamicList< label > pointTgtEdges_
The intersection points' target edges, or -1 if the point.
Definition: patchIntersection.H:94
List
Motion of the mesh specified as a list of pointMeshMovers.
Foam::treeBoundBox
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:90
Foam::treeBoundBox::extend
treeBoundBox extend(const scalar s) const
Return asymmetrically extended bounding box, with guaranteed.
Definition: treeBoundBoxI.H:321
Foam::treeDataPrimitivePatch
Encapsulation of data needed to search on PrimitivePatches.
Definition: treeDataPrimitivePatch.H:64
Foam::triFace
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:71
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
vectorField
volVectorField vectorField(fieldObject, mesh)
points
const pointField & points
Definition: gmvOutputHeader.H:1
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:250
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Foam::dimensions::time
const dimensionSet time
Foam::saturationModels::A
static const coefficient A("A", dimPressure, 611.21)
Foam::triIntersect::srcPointTgtTriIntersection
barycentric2D srcPointTgtTriIntersection(const point &srcP, const vector &srcN, const FixedList< point, 3 > &tgtPs)
Calculate the intersection of a projected source point with a target.
Definition: triIntersect.C:1786
Foam::triIntersect::srcTriTgtPointIntersection
barycentric2D srcTriTgtPointIntersection(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const point &tgtP)
Calculate the intersection of a target point with a source triangle's.
Definition: triIntersect.C:1538
Foam::triIntersect::tgtTriInterpolate
Type tgtTriInterpolate(const barycentric2D &y, const FixedList< Type, 3 > &tgtPsis)
Use the coordinates obtained from srcPointTgtTriIntersection to interpolate.
Definition: triIntersectTemplates.C:87
Foam::triIntersect::srcTriInterpolate
Type srcTriInterpolate(const barycentric2D &y, const FixedList< Type, 3 > &tgtPsis)
Use the coordinates obtained from srcTriTgtPointIntersection to interpolate.
Definition: triIntersectTemplates.C:54
Foam::triIntersect::intersectTris
void intersectTris(const FixedList< point, 3 > &srcPs, const FixedList< vector, 3 > &srcNs, const FixedList< bool, 3 > &srcOwns, const FixedList< label, 3 > &srcTgtPis, const FixedList< point, 3 > &tgtPs, const FixedList< bool, 3 > &tgtOwns, const FixedList< label, 3 > &tgtSrcPis, DynamicList< point > &srcPoints, DynamicList< vector > &srcPointNormals, DynamicList< point > &tgtPoints, DynamicList< location > &pointLocations, const bool debug, const word &writePrefix=word::null)
Construct the intersection of a source triangle's projected volume and a.
Definition: triIntersect.C:1849
Foam::vtkWriteOps::insert
void insert(const scalar, DynamicList< floatScalar > &)
Append scalar to given DynamicList.
Foam::vtkWriteOps::write
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
Foam::vtkWritePolyData::write
void write(const fileName &file, const word &title, const bool binary, const PointField &points, const VertexList &vertices, const LineList &lines, const FaceList &faces, const wordList &fieldNames, const boolList &fieldIsPointValues, const UPtrList< const Field< label >> &fieldLabelValues #define FieldTypeValuesConstArg(Type, nullArg))
Write VTK polygonal data to a file. Takes a PtrList of fields of labels and.
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Foam::labelPair
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
Foam::decrIndent
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:272
Foam::labelList
List< label > labelList
A List of labels.
Definition: labelList.H:56
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:288
Foam::reverseList
ListType reverseList(const ListType &list)
Reverse a list. First element becomes last element etc.
Definition: ListOpsTemplates.C:780
Foam::inplaceReverseList
void inplaceReverseList(ListType &list)
Inplace reversal of a list using Swap.
Definition: ListOpsTemplates.C:797
Foam::Info
messageStream Info
Foam::inplaceRotateList
void inplaceRotateList(ListType< DataType > &list, label n)
Inplace reversal of a list using the Reversal Block Swapping algorithm.
Definition: ListOpsTemplates.C:837
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Foam::I
static const Identity< scalar > I
Definition: Identity.H:93
Foam::barycentric2D
Barycentric2D< scalar > barycentric2D
A scalar version of the templated Barycentric2D.
Definition: barycentric2D.H:45
Foam::incrIndent
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:265
Foam::point
vector point
Point is a vector.
Definition: point.H:41
Foam::sqr
tmp< DimensionedField< typename outerProduct< Type, Type >::type, GeoMesh, Field >> sqr(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
Definition: DimensionedFieldFunctions.C:134
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Foam::barycentricTensor2D
BarycentricTensor2D< scalar > barycentricTensor2D
A scalar version of the templated BarycentricTensor2D.
Definition: barycentricTensor2D.H:46
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
Definition: ListOpsTemplates.C:73
Foam::name
word name(const LagrangianState state)
Return a string representation of a Lagrangian state enumeration.
Definition: LagrangianState.C:30
Foam::triPointRef
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
Foam::mag
tmp< DimensionedField< scalar, GeoMesh, Field > > mag(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
Definition: DimensionedFieldFunctions.C:238
Foam::findIndex
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Foam::FatalError
error FatalError
Foam::Swap
void Swap(T &a, T &b)
Definition: Swap.H:43
Foam::cmptSum
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:507
Foam::magSqr
tmp< DimensionedField< scalar, GeoMesh, Field > > magSqr(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
Definition: DimensionedFieldFunctions.C:193
Foam::faceList
List< face > faceList
Definition: faceListFwd.H:41
Foam::labelUList
UList< label > labelUList
Definition: UList.H:65
Foam::indent
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:243
Foam::T
void T(GeometricField< Type, GeoMesh, PrimitiveField1 > &gf, const GeometricField< Type, GeoMesh, PrimitiveField2 > &gf1)
Definition: GeometricFieldFunctions.C:81
Foam::nl
static const char nl
Definition: Ostream.H:297
Foam::max
dimensioned< Type > max(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
Definition: DimensionedFieldFunctions.C:373
Foam::type
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
newPointi
label newPointi
Definition: readKivaGrid.H:495
triFace
face triFace(3)
p
volScalarField & p
Definition: createFieldRefs.H:4
forAllStarFaces
#define forAllStarFaces(star, starFacei, facei)
Definition: star.H:219
forAllStarEdges
#define forAllStarEdges(star, i, starEdgei, edgei)
Definition: star.H:240
triIntersect.H
Functions with which to perform an intersection of a pair of triangles; the source and target....