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-2024 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> circle = t.circumCircle();
927 
928  return magSqr(points[pointi] - circle.first()) - sqr(circle.second());
929 }
930 
931 
932 template<class SrcPatchType, class TgtPatchType>
933 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::insertPoints
934 (
935  const label insertionTriOrEdgei,
936  const bool isTri,
937  const UList<label>& pointis,
938  UList<label>& insertionEdgeis,
939  const UList<label>& fixedEdgeis
940 )
941 {
942  // Clear the insertion edge list
943  insertionEdgeis = -1;
944 
945  // If inserting into a tri then initialise the candidates with the given
946  // insertion tri
947  if (isTri)
948  {
949  candidateTriTris_.append(insertionTriOrEdgei);
950  triCandidateTris_[insertionTriOrEdgei] = 0;
951  }
952 
953  // If inserting into an edge, get the insertion edge label
954  if (!isTri)
955  {
956  insertionEdgeis[0] = insertionTriOrEdgei;
957  }
958 
959  // Check that we are inserting unambiguously into the source or target side
960  #ifdef FULLDEBUG
961  if (!isTri)
962  {
963  const label trii0 = edgeTris_[insertionTriOrEdgei][0];
964  const label trii1 = edgeTris_[insertionTriOrEdgei][1];
965 
966  if
967  (
968  trii0 != -1
969  && trii1 != -1
970  && (triSrcFace_[trii0] == -1) != (triSrcFace_[trii1] == -1)
971  )
972  {
973  FatalErrorInFunction
974  << "Attempted insertion into front edge #"
975  << insertionTriOrEdgei << exit(FatalError);
976  }
977  }
978  #endif
979 
980  // Determine the patch associations for the insertion tri or edge
981  const label adjacentTri =
982  isTri
983  ? insertionTriOrEdgei
984  : edgeTris_[insertionTriOrEdgei][edgeTris_[insertionTriOrEdgei][0] == -1];
985  const bool isSrc = triSrcFace_[adjacentTri] != -1;
986  const label patchEdgei =
987  isTri
988  ? -1
989  : edgePatchEdge(insertionTriOrEdgei, isSrc);
990  const label patchFacei =
991  patchEdgei != -1
992  ? -1
993  : (isSrc ? triSrcFace_[adjacentTri] : triTgtFace_[adjacentTri]);
994 
995  // Get connectivity arrays that will be modified during insertion
996  List<DynamicList<label>>& patchEdgePoints =
997  isSrc ? this->srcEdgePoints_ : this->tgtEdgePoints_;
998  DynamicList<label>& pointPatchEdges =
999  isSrc ? this->pointSrcEdges_ : this->pointTgtEdges_;
1000  DynamicList<label>& pointPatchFaces =
1001  isSrc ? this->pointSrcFaces_ : this->pointTgtFaces_;
1002 
1003  // Determine whether the insertion star is closed or not
1004  const bool isClosedStar =
1005  isTri || findIndex(edgeTris_[insertionEdgeis[0]], -1) == -1;
1006 
1007  // Loop the points
1008  forAll(pointis, pointii)
1009  {
1010  const label pointi = pointis[pointii];
1011 
1012  // Initial star tri or edge. If a tri, this means searching the
1013  // candidates for the one that best contains the point. If an edge then
1014  // just use the next insertion edge.
1015  label initialTriOrEdgei = -1;
1016  if (isTri)
1017  {
1018  label insertionCandidateTrii = -1;
1019  scalar minDistSqr = vGreat;
1020  forAll(candidateTriTris_, candidateTrii)
1021  {
1022  const label trii = candidateTriTris_[candidateTrii];
1023  if (trii == -1) continue;
1024  const scalar distSqr = circumDistSqr(trii, pointi);
1025  if (distSqr < minDistSqr)
1026  {
1027  insertionCandidateTrii = candidateTrii;
1028  minDistSqr = distSqr;
1029  }
1030  if (minDistSqr < 0)
1031  {
1032  break;
1033  }
1034  }
1035  initialTriOrEdgei = candidateTriTris_[insertionCandidateTrii];
1036  }
1037  else
1038  {
1039  initialTriOrEdgei = insertionEdgeis[pointii];
1040  }
1041 
1042  // Populate the star
1043  star::context starContext = star_.populate
1044  (
1045  initialTriOrEdgei,
1046  isTri,
1047  [&](const label edgei, const label trii)
1048  {
1049  // If this edge is part of a patch edge, or if it connects to a
1050  // triangle in the other patch, or if it is fixed, then it
1051  // should not be removed, so it should not be crossed
1052  if
1053  (
1054  edgePatchEdges(edgei) != labelPair(-1, -1)
1055  || (isSrc && triTgtFace_[trii] != -1)
1056  || (!isSrc && triSrcFace_[trii] != -1)
1057  || findIndex(fixedEdgeis, edgei) != -1
1058  )
1059  {
1060  return false;
1061  }
1062 
1063  // If the destination tri borders any tris that are already in
1064  // the star (other than the tri from which we walked) then do
1065  // not cross into it. Adding this tri would leave a point
1066  // within the star, which is not allowed.
1067  const label triEdgei = findIndex(triEdges_[trii], edgei);
1068  const label edgej0 = triEdges_[trii][(triEdgei + 1) % 3];
1069  const label edgej1 = triEdges_[trii][(triEdgei + 2) % 3];
1070  const label trij0 =
1071  edgeTris_[edgej0][edgeTris_[edgej0][0] == trii];
1072  const label trij1 =
1073  edgeTris_[edgej1][edgeTris_[edgej1][0] == trii];
1074  if
1075  (
1076  (trij0 != -1 && star_.faceStarFaces()[trij0] != -1)
1077  || (trij1 != -1 && star_.faceStarFaces()[trij1] != -1)
1078  || (!isTri && findIndex(insertionEdgeis, edgej0) != -1)
1079  || (!isTri && findIndex(insertionEdgeis, edgej1) != -1)
1080  )
1081  {
1082  return false;
1083  }
1084 
1085  // Otherwise, cross into the triangle if it's circumcircle
1086  // contains the insertion point
1087  return circumDistSqr(trii, pointi) < 0;
1088  },
1089  triEdges_,
1090  edgeTris_
1091  );
1092 
1093  // Get the end points of the current insertion edge
1094  const edge insertionEdge =
1095  isTri ? edge(-1, -1) : edgePoints(insertionEdgeis[pointii]);
1096 
1097  // Walk around the star polygon disconnecting the old triangles and
1098  // adding in the new ones
1099  label newTrii0 = -1, newTrii00 = -1;
1100  forAllStarEdges(star_, i, starEdgei, edgei)
1101  {
1102  const bool isFirst = i == 0;
1103  const bool isLast = i == star_.starEdgeEdges().size() - 1;
1104 
1105  const bool edgeSidei =
1106  edgeTris_[edgei][0] == -1
1107  || star_.faceStarFaces()[edgeTris_[edgei][0]] == -1;
1108 
1109  const label trii = edgeTris_[edgei][edgeSidei];
1110  const label triEdgei = findIndex(triEdges_[trii], edgei);
1111 
1112  // Disconnect the existing triangle from the star
1113  const label tempEdgei = newEdgei();
1114  triEdges_[trii][triEdgei] = tempEdgei;
1115  edgeTris_[tempEdgei][0] = trii;
1116  edgeTris_[edgei][edgeSidei] = -1;
1117 
1118  // Add the new triangle
1119  const label newTrii =
1120  addTri
1121  (
1122  {
1123  pointi,
1124  triPoint(trii, triEdgei),
1125  triPoint(trii, (triEdgei + 1) % 3),
1126  },
1127  {
1128  !isFirst ? triEdges_[newTrii0][2] : -1,
1129  edgei,
1130  isClosedStar && isLast ? triEdges_[newTrii00][0] : -1
1131  },
1132  isSrc ? triSrcFace_[trii] : triTgtFace_[trii],
1133  isSrc
1134  );
1135  newTrii0 = newTrii;
1136  newTrii00 = isFirst ? newTrii0 : newTrii00;
1137 
1138  // Add the new triangle into the list of candidates
1139  if (isTri)
1140  {
1141  candidateTriTris_.append(newTrii);
1142  triCandidateTris_[newTrii] = candidateTriTris_.size() - 1;
1143  }
1144 
1145  // Get the previous insertion edge
1146  if (triPoint(newTrii, 2) == insertionEdge[0])
1147  {
1148  insertionEdgeis[pointii] = triEdges_[newTrii][2];
1149  }
1150 
1151  // Get the next insertion edge
1152  if (!isTri && triPoint(newTrii, 1) == insertionEdge[1])
1153  {
1154  insertionEdgeis[pointii + 1] = triEdges_[newTrii][0];
1155  }
1156  }
1157 
1158  // Create the edge or face association for the new point
1159  if (patchEdgei != -1)
1160  {
1161  patchEdgePoints[patchEdgei].append(-1);
1162  label patchEdgePointi = patchEdgePoints[patchEdgei].size() - 1;
1163  for (; patchEdgePointi >= 0; -- patchEdgePointi)
1164  {
1165  label& pointi = patchEdgePoints[patchEdgei][patchEdgePointi];
1166  pointi = patchEdgePoints[patchEdgei][patchEdgePointi - 1];
1167  if
1168  (
1169  pointPoints_[pointi] == insertionEdge[0]
1170  || pointPoints_[pointi] == insertionEdge[1]
1171  )
1172  {
1173  break;
1174  }
1175  }
1176  -- patchEdgePointi;
1177  patchEdgePoints[patchEdgei][patchEdgePointi] = pointi;
1178 
1179  pointPatchEdges[pointi] = patchEdgei;
1180  }
1181  else
1182  {
1183  pointPatchFaces[pointi] = patchFacei;
1184  }
1185 
1186  // Check
1187  if (patchEdgei != -1)
1188  {
1189  checkPatchEdge(patchEdgei, isSrc);
1190  }
1191  else
1192  {
1193  checkPatchFace(patchFacei, isSrc);
1194  }
1195 
1196  // If the insertion edge is not the same orientation as the original
1197  // insertion edge then reverse it
1198  if (!isTri)
1199  {
1200  if (edgePoints(insertionEdgeis[pointii + 1])[1] != insertionEdge[1])
1201  {
1202  Swap
1203  (
1204  edgeTris_[insertionEdgeis[pointii + 1]][0],
1205  edgeTris_[insertionEdgeis[pointii + 1]][1]
1206  );
1207  }
1208  }
1209 
1210  // Remove the star triangles
1211  forAllStarFaces(star_, starTrii, trii)
1212  {
1213  removeTri(trii);
1214  }
1215 
1216  // Check
1217  if (patchEdgei != -1)
1218  {
1219  checkPatchEdge(patchEdgei, isSrc);
1220  }
1221  else
1222  {
1223  checkPatchFace(patchFacei, isSrc);
1224  }
1225  }
1226 
1227  // Clear the candidate tris
1228  if (isTri)
1229  {
1230  forAll(candidateTriTris_, candidateTrii)
1231  {
1232  const label trii = candidateTriTris_[candidateTrii];
1233  if (trii != -1)
1234  {
1235  triCandidateTris_[trii] = -1;
1236  }
1237  }
1238  candidateTriTris_.clear();
1239  }
1240 }
1241 
1242 
1243 template<class SrcPatchType, class TgtPatchType>
1244 bool Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::pointCanIntersect
1245 (
1246  const label pointi
1247 ) const
1248 {
1249  return
1250  (
1251  this->pointSrcPoints_[pointi] != -1
1252  && this->pointTgtPoints_[pointi] == -1
1253  && this->pointTgtEdges_[pointi] == -1
1254  && this->pointTgtFaces_[pointi] == -1
1255  )
1256  || (
1257  this->pointTgtPoints_[pointi] != -1
1258  && this->pointSrcPoints_[pointi] == -1
1259  && this->pointSrcEdges_[pointi] == -1
1260  && this->pointSrcFaces_[pointi] == -1
1261  );
1262 }
1263 
1264 
1265 template<class SrcPatchType, class TgtPatchType>
1266 bool Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::edgeCanIntersect
1267 (
1268  const label edgei
1269 ) const
1270 {
1271  return
1272  (
1273  edgePatchEdge(edgei, true) != -1
1274  && edgePatchEdge(edgei, false) == -1
1275  && (
1276  edgeTris_[edgei][0] == -1
1277  || triTgtFace_[edgeTris_[edgei][0]] == -1
1278  )
1279  && (
1280  edgeTris_[edgei][1] == -1
1281  || triTgtFace_[edgeTris_[edgei][1]] == -1
1282  )
1283  )
1284  || (
1285  edgePatchEdge(edgei, false) != -1
1286  && edgePatchEdge(edgei, true) == -1
1287  && (
1288  edgeTris_[edgei][0] == -1
1289  || triSrcFace_[edgeTris_[edgei][0]] == -1
1290  )
1291  && (
1292  edgeTris_[edgei][1] == -1
1293  || triSrcFace_[edgeTris_[edgei][1]] == -1
1294  )
1295  );
1296 }
1297 
1298 
1299 template<class SrcPatchType, class TgtPatchType>
1300 void
1301 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::snapPatchFaceTris
1302 (
1303  const label srcFacei,
1304  const label tgtFacei,
1305  const scalar snapTol
1306 )
1307 {
1308  static const FixedList<label, 3> triNoSnap({-1, -1, -1});
1309 
1310  const vector tgtNormal =
1311  triPointRef(tgtPoints_, patchFacePoints(tgtFacei, false)).normal();
1312 
1313  // Determine what source points snap to
1314  FixedList<barycentric2D, 3> srcFaceTgtTs(barycentric2D::uniform(-vGreat));
1315  FixedList<label, 3> srcFaceSnapTgtFacePoint(triNoSnap);
1316  FixedList<label, 3> srcFaceSnapTgtFaceEdge(triNoSnap);
1317  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1318  {
1319  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1320 
1321  if (!pointCanIntersect(srcPointi)) continue;
1322 
1323  if ((srcPointNormals_[srcPointi] & tgtNormal) < 0)
1324  {
1325  srcFaceTgtTs[srcFacePointi] =
1326  triIntersect::srcPointTgtTriIntersection
1327  (
1328  srcPoints_[srcPointi],
1329  srcPointNormals_[srcPointi],
1330  patchFacePointValues(tgtFacei, false, tgtPoints_)
1331  );
1332 
1333  for (label tgtFacePointi = 0; tgtFacePointi < 3; ++ tgtFacePointi)
1334  {
1335  const label tgtPointi =
1336  patchFacePoint(tgtFacei, tgtFacePointi, false);
1337 
1338  const label tgtFacePointi0 = (tgtFacePointi + 2) % 3;
1339  const label tgtFacePointi1 = (tgtFacePointi + 1) % 3;
1340 
1341  if
1342  (
1343  pointCanIntersect(tgtPointi)
1344  && mag(srcFaceTgtTs[srcFacePointi][tgtFacePointi0]) < snapTol
1345  && mag(srcFaceTgtTs[srcFacePointi][tgtFacePointi1]) < snapTol
1346  )
1347  {
1348  srcFaceSnapTgtFacePoint[srcFacePointi] =
1349  srcFaceSnapTgtFacePoint[srcFacePointi] == -1
1350  ? tgtFacePointi
1351  : -2;
1352 
1353  srcFaceTgtTs[srcFacePointi][tgtFacePointi] = 1;
1354  srcFaceTgtTs[srcFacePointi][tgtFacePointi0] = 0;
1355  srcFaceTgtTs[srcFacePointi][tgtFacePointi1] = 0;
1356  }
1357  }
1358 
1359  srcFaceSnapTgtFacePoint[srcFacePointi] =
1360  max(srcFaceSnapTgtFacePoint[srcFacePointi], -1);
1361 
1362  if (srcFaceSnapTgtFacePoint[srcFacePointi] != -1) continue;
1363 
1364  for (label tgtFaceEdgei = 0; tgtFaceEdgei < 3; ++ tgtFaceEdgei)
1365  {
1366  const label tgtFacePointi0 = tgtFaceEdgei;
1367  const label tgtFacePointi1 = (tgtFaceEdgei + 1) % 3;
1368  const label tgtFacePointiOpp = (tgtFaceEdgei + 2) % 3;
1369 
1370  if
1371  (
1372  srcFaceTgtTs[srcFacePointi][tgtFacePointi0] >= snapTol
1373  && srcFaceTgtTs[srcFacePointi][tgtFacePointi1] >= snapTol
1374  && mag(srcFaceTgtTs[srcFacePointi][tgtFacePointiOpp]) < snapTol
1375  )
1376  {
1377  srcFaceSnapTgtFaceEdge[srcFacePointi] =
1378  srcFaceSnapTgtFaceEdge[srcFacePointi] == -1
1379  ? tgtFaceEdgei
1380  : -2;
1381 
1382  const scalar eps =
1383  srcFaceTgtTs[srcFacePointi][tgtFacePointiOpp];
1384  srcFaceTgtTs[srcFacePointi][tgtFacePointiOpp] = 0;
1385  srcFaceTgtTs[srcFacePointi][tgtFacePointi0] /= 1 - eps;
1386  srcFaceTgtTs[srcFacePointi][tgtFacePointi1] /= 1 - eps;
1387  }
1388  }
1389 
1390  srcFaceSnapTgtFaceEdge[srcFacePointi] =
1391  max(srcFaceSnapTgtFaceEdge[srcFacePointi], -1);
1392  }
1393  }
1394 
1395  // Determine what target points snap to
1396  FixedList<barycentric2D, 3> tgtFaceSrcTs(barycentric2D::uniform(-vGreat));
1397  FixedList<label, 3> tgtFaceSnapSrcFacePoint(triNoSnap);
1398  FixedList<label, 3> tgtFaceSnapSrcFaceEdge(triNoSnap);
1399  forAll(tgtFaceSnapSrcFacePoint, tgtFacePointi)
1400  {
1401  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1402 
1403  if (!pointCanIntersect(tgtPointi)) continue;
1404 
1405  tgtFaceSrcTs[tgtFacePointi] =
1406  triIntersect::srcTriTgtPointIntersection
1407  (
1408  patchFacePointValues(srcFacei, true, srcPoints_),
1409  patchFacePointValues(srcFacei, true, srcPointNormals_),
1410  tgtPoints_[tgtPointi]
1411  );
1412 
1413  const vector srcNormal =
1414  triIntersect::srcTriInterpolate
1415  (
1416  tgtFaceSrcTs[tgtFacePointi],
1417  patchFacePointValues(srcFacei, true, srcPointNormals_)
1418  );
1419 
1420  if ((srcNormal & tgtNormal) < 0)
1421  {
1422  for (label srcFacePointi = 0; srcFacePointi < 3; ++ srcFacePointi)
1423  {
1424  const label srcPointi =
1425  patchFacePoint(srcFacei, srcFacePointi, true);
1426 
1427  const label srcFacePointi0 = (srcFacePointi + 2) % 3;
1428  const label srcFacePointi1 = (srcFacePointi + 1) % 3;
1429 
1430  if
1431  (
1432  pointCanIntersect(srcPointi)
1433  && mag(tgtFaceSrcTs[tgtFacePointi][srcFacePointi0]) < snapTol
1434  && mag(tgtFaceSrcTs[tgtFacePointi][srcFacePointi1]) < snapTol
1435  )
1436  {
1437  tgtFaceSnapSrcFacePoint[tgtFacePointi] =
1438  tgtFaceSnapSrcFacePoint[tgtFacePointi] == -1
1439  ? srcFacePointi
1440  : -2;
1441 
1442  tgtFaceSrcTs[tgtFacePointi][srcFacePointi] = 1;
1443  tgtFaceSrcTs[tgtFacePointi][srcFacePointi0] = 0;
1444  tgtFaceSrcTs[tgtFacePointi][srcFacePointi1] = 0;
1445  }
1446  }
1447 
1448  tgtFaceSnapSrcFacePoint[tgtFacePointi] =
1449  max(tgtFaceSnapSrcFacePoint[tgtFacePointi], -1);
1450 
1451  if (tgtFaceSnapSrcFacePoint[tgtFacePointi] != -1) continue;
1452 
1453  for (label srcFaceEdgei = 0; srcFaceEdgei < 3; ++ srcFaceEdgei)
1454  {
1455  const label srcFacePointi0 = srcFaceEdgei;
1456  const label srcFacePointi1 = (srcFaceEdgei + 1) % 3;
1457  const label srcFacePointiOpp = (srcFaceEdgei + 2) % 3;
1458 
1459  if
1460  (
1461  tgtFaceSrcTs[tgtFacePointi][srcFacePointi0] >= snapTol
1462  && tgtFaceSrcTs[tgtFacePointi][srcFacePointi1] >= snapTol
1463  && mag(tgtFaceSrcTs[tgtFacePointi][srcFacePointiOpp]) < snapTol
1464  )
1465  {
1466  tgtFaceSnapSrcFaceEdge[tgtFacePointi] =
1467  tgtFaceSnapSrcFaceEdge[tgtFacePointi] == -1
1468  ? srcFaceEdgei
1469  : -2;
1470 
1471  const scalar eps =
1472  tgtFaceSrcTs[tgtFacePointi][srcFacePointiOpp];
1473  tgtFaceSrcTs[tgtFacePointi][srcFacePointiOpp] = 0;
1474  tgtFaceSrcTs[tgtFacePointi][srcFacePointi0] /= 1 - eps;
1475  tgtFaceSrcTs[tgtFacePointi][srcFacePointi1] /= 1 - eps;
1476  }
1477  }
1478 
1479  tgtFaceSnapSrcFaceEdge[tgtFacePointi] =
1480  max(tgtFaceSnapSrcFaceEdge[tgtFacePointi], -1);
1481  }
1482  }
1483 
1484  // Return if there is nothing to do
1485  if
1486  (
1487  srcFaceSnapTgtFacePoint == triNoSnap
1488  && srcFaceSnapTgtFaceEdge == triNoSnap
1489  && tgtFaceSnapSrcFacePoint == triNoSnap
1490  && tgtFaceSnapSrcFaceEdge == triNoSnap
1491  )
1492  {
1493  return;
1494  }
1495 
1496  // Make point-point snaps symmetric
1497  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1498  {
1499  const label tgtFacePointi = srcFaceSnapTgtFacePoint[srcFacePointi];
1500 
1501  if (tgtFacePointi != -1)
1502  {
1503  tgtFaceSnapSrcFacePoint[tgtFacePointi] = srcFacePointi;
1504  }
1505  }
1506  forAll(tgtFaceSnapSrcFacePoint, tgtFacePointi)
1507  {
1508  const label srcFacePointi = tgtFaceSnapSrcFacePoint[tgtFacePointi];
1509 
1510  if (srcFacePointi != -1)
1511  {
1512  srcFaceSnapTgtFacePoint[srcFacePointi] = tgtFacePointi;
1513  }
1514  }
1515 
1516  // Do point-point motion
1517  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1518  {
1519  const label tgtFacePointi = srcFaceSnapTgtFacePoint[srcFacePointi];
1520 
1521  if (tgtFacePointi == -1) continue;
1522 
1523  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1524  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1525 
1526  if (!pointCanIntersect(srcPointi)) continue;
1527  if (!pointCanIntersect(tgtPointi)) continue;
1528 
1529  srcPoints_[tgtPointi] = srcPoints_[srcPointi];
1530  srcPointNormals_[tgtPointi] = srcPointNormals_[srcPointi];
1531  tgtPoints_[srcPointi] = tgtPoints_[tgtPointi];
1532 
1533  const point& srcP = srcPoints_[srcPointi];
1534  const vector& srcN = srcPointNormals_[srcPointi];
1535  const point& tgtP = tgtPoints_[tgtPointi];
1536  const vector d = (tensor::I - sqr(srcN)) & (tgtP - srcP);
1537  srcPoints_[srcPointi] += d/2;
1538  tgtPoints_[tgtPointi] -= d/2;
1539  }
1540 
1541  // Do source-point-target-edge motion
1542  forAll(srcFaceSnapTgtFaceEdge, srcFacePointi)
1543  {
1544  const label tgtFaceEdgei = srcFaceSnapTgtFaceEdge[srcFacePointi];
1545 
1546  if (tgtFaceEdgei == -1) continue;
1547 
1548  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1549 
1550  if (!pointCanIntersect(srcPointi)) continue;
1551 
1552  tgtPoints_[srcPointi] =
1553  triIntersect::tgtTriInterpolate
1554  (
1555  srcFaceTgtTs[srcFacePointi],
1556  patchFacePointValues(tgtFacei, false, tgtPoints_)
1557  );
1558 
1559  const point& srcP = srcPoints_[srcPointi];
1560  const vector& srcN = srcPointNormals_[srcPointi];
1561  const point& tgtP = tgtPoints_[srcPointi];
1562  const vector d = (tensor::I - sqr(srcN)) & (tgtP - srcP);
1563  srcPoints_[srcPointi] += d;
1564  }
1565 
1566  // Do target-point-source-edge motion
1567  forAll(tgtFaceSnapSrcFaceEdge, tgtFacePointi)
1568  {
1569  const label srcFaceEdgei = tgtFaceSnapSrcFaceEdge[tgtFacePointi];
1570 
1571  if (srcFaceEdgei == -1) continue;
1572 
1573  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1574 
1575  if (!pointCanIntersect(tgtPointi)) continue;
1576 
1577  srcPoints_[tgtPointi] =
1578  triIntersect::srcTriInterpolate
1579  (
1580  tgtFaceSrcTs[tgtFacePointi],
1581  patchFacePointValues(srcFacei, true, srcPoints_)
1582  );
1583  srcPointNormals_[tgtPointi] =
1584  triIntersect::srcTriInterpolate
1585  (
1586  tgtFaceSrcTs[tgtFacePointi],
1587  patchFacePointValues(srcFacei, true, srcPointNormals_)
1588  );
1589 
1590  const point& srcP = srcPoints_[tgtPointi];
1591  const vector& srcN = srcPointNormals_[tgtPointi];
1592  const point& tgtP = tgtPoints_[tgtPointi];
1593  const vector d = (tensor::I - sqr(srcN)) & (tgtP - srcP);
1594  tgtPoints_[tgtPointi] -= d;
1595  }
1596 
1597  // Do point-point snapping
1598  forAll(srcFaceSnapTgtFacePoint, srcFacePointi)
1599  {
1600  const label tgtFacePointi = srcFaceSnapTgtFacePoint[srcFacePointi];
1601 
1602  if (tgtFacePointi == -1) continue;
1603 
1604  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1605  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1606 
1607  if (!pointCanIntersect(srcPointi)) continue;
1608  if (!pointCanIntersect(tgtPointi)) continue;
1609 
1610  srcPoints_[tgtPointi] = srcPoints_[srcPointi];
1611  srcPointNormals_[tgtPointi] = srcPointNormals_[srcPointi];
1612  tgtPoints_[srcPointi] = tgtPoints_[tgtPointi];
1613 
1614  pointPoints_[tgtPointi] = srcPointi;
1615 
1616  this->pointTgtEdges_[srcPointi] = this->pointTgtEdges_[tgtPointi];
1617  this->pointSrcEdges_[tgtPointi] = -1;
1618  this->pointTgtEdges_[tgtPointi] = -1;
1619 
1620  this->pointTgtPoints_[srcPointi] = this->pointTgtPoints_[tgtPointi];
1621  this->tgtPointPoints_[this->pointTgtPoints_[tgtPointi]] = srcPointi;
1622  this->pointSrcPoints_[tgtPointi] = -1;
1623  this->pointTgtPoints_[tgtPointi] = -1;
1624  }
1625 
1626  // Check
1627  checkPatchFace(srcFacei, true);
1628  checkPatchFace(tgtFacei, false);
1629 
1630  // Insertion workspace
1631  labelList insertPointis(1), insertEdgeis(2, label(-1)), fixedEdgeis(0);
1632 
1633  // Do source-point-target-edge snapping
1634  forAll(srcFaceSnapTgtFaceEdge, srcFacePointi)
1635  {
1636  const label tgtFaceEdgei = srcFaceSnapTgtFaceEdge[srcFacePointi];
1637 
1638  if (tgtFaceEdgei == -1) continue;
1639 
1640  const label srcPointi = patchFacePoint(srcFacei, srcFacePointi, true);
1641 
1642  if (!pointCanIntersect(srcPointi)) continue;
1643 
1644  const label tgtPatchEdgei =
1645  patchFacePatchEdge(tgtFacei, tgtFaceEdgei, false);
1646 
1647  // Loop all target face tris and find the best edge into which to
1648  // insert the source point
1649  label insertTgtEdgei = -1;
1650  scalar insertDistSqr = vGreat;
1651  forAll(tgtFaceTris_[tgtFacei], tgtFaceTrii)
1652  {
1653  const label tgtTrii = tgtFaceTris_[tgtFacei][tgtFaceTrii];
1654 
1655  forAll(triEdges_[tgtTrii], tgtTriEdgei)
1656  {
1657  const label tgtEdgei = triEdges_[tgtTrii][tgtTriEdgei];
1658 
1659  if
1660  (
1661  edgeCanIntersect(tgtEdgei)
1662  && tgtPatchEdgei == edgePatchEdge(tgtEdgei, false)
1663  )
1664  {
1665  const scalar distSqr = circumDistSqr(tgtTrii, srcPointi);
1666 
1667  if (distSqr < insertDistSqr)
1668  {
1669  insertDistSqr = distSqr;
1670  insertTgtEdgei = tgtEdgei;
1671  }
1672  }
1673  }
1674  }
1675 
1676  if (insertTgtEdgei != -1)
1677  {
1678  insertPointis[0] = srcPointi;
1679  insertPoints
1680  (
1681  insertTgtEdgei,
1682  false,
1683  insertPointis,
1684  insertEdgeis,
1685  fixedEdgeis
1686  );
1687  }
1688  }
1689 
1690  // Do target-point-source-edge snapping
1691  forAll(tgtFaceSnapSrcFaceEdge, tgtFacePointi)
1692  {
1693  const label srcFaceEdgei = tgtFaceSnapSrcFaceEdge[tgtFacePointi];
1694 
1695  if (srcFaceEdgei == -1) continue;
1696 
1697  const label tgtPointi = patchFacePoint(tgtFacei, tgtFacePointi, false);
1698 
1699  if (!pointCanIntersect(tgtPointi)) continue;
1700 
1701  const label srcPatchEdgei =
1702  patchFacePatchEdge(srcFacei, srcFaceEdgei, true);
1703 
1704  srcPoints_[tgtPointi] =
1705  triIntersect::srcTriInterpolate
1706  (
1707  tgtFaceSrcTs[tgtFacePointi],
1708  patchFacePointValues(srcFacei, true, srcPoints_)
1709  );
1710  srcPointNormals_[tgtPointi] =
1711  triIntersect::srcTriInterpolate
1712  (
1713  tgtFaceSrcTs[tgtFacePointi],
1714  patchFacePointValues(srcFacei, true, srcPointNormals_)
1715  );
1716 
1717  // Loop all target face tris and find the best edge into which to
1718  // insert the source point
1719  label insertSrcEdgei = -1;
1720  scalar insertDistSqr = vGreat;
1721  forAll(srcFaceTris_[srcFacei], srcFaceTrii)
1722  {
1723  const label srcTrii = srcFaceTris_[srcFacei][srcFaceTrii];
1724 
1725  forAll(triEdges_[srcTrii], srcTriEdgei)
1726  {
1727  const label srcEdgei = triEdges_[srcTrii][srcTriEdgei];
1728 
1729  if
1730  (
1731  edgeCanIntersect(srcEdgei)
1732  && srcPatchEdgei == edgePatchEdge(srcEdgei, true)
1733  )
1734  {
1735  const scalar distSqr = circumDistSqr(srcTrii, tgtPointi);
1736 
1737  if (distSqr < insertDistSqr)
1738  {
1739  insertDistSqr = distSqr;
1740  insertSrcEdgei = srcEdgei;
1741  }
1742  }
1743  }
1744  }
1745 
1746  if (insertSrcEdgei != -1)
1747  {
1748  insertPointis[0] = tgtPointi;
1749  insertPoints
1750  (
1751  insertSrcEdgei,
1752  false,
1753  insertPointis,
1754  insertEdgeis,
1755  fixedEdgeis
1756  );
1757  }
1758  }
1759 
1760  // Check
1761  checkPatchFace(srcFacei, true);
1762  checkPatchFace(tgtFacei, false);
1763 }
1764 
1765 
1766 template<class SrcPatchType, class TgtPatchType>
1767 bool Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::intersectTris
1768 (
1769  const label srcTrii,
1770  const label tgtTrii
1771 )
1772 {
1773  // Get the source and target faces that these tris are associated with
1774  const label srcFacei = triSrcFace_[srcTrii];
1775  const label tgtFacei = triTgtFace_[tgtTrii];
1776  if (srcFacei == -1 || tgtFacei == -1)
1777  {
1778  FatalErrorInFunction
1779  << "Tri-intersections must be between a tri associated with the "
1780  << "source patch and one associated with the target patch"
1781  << exit(FatalError);
1782  }
1783 
1784  // Do intersection
1785  DynamicList<point> ictSrcPoints;
1786  DynamicList<vector> ictSrcPointNormals;
1787  DynamicList<point> ictTgtPoints;
1788  DynamicList<triIntersect::location> ictPointLocations;
1789  triIntersect::intersectTris
1790  (
1791  triPointValues(srcTrii, srcPoints_),
1792  triPointValues(srcTrii, srcPointNormals_),
1793  triOwns(srcTrii),
1794  triOtherTriPoints(srcTrii, tgtTrii),
1795  triPointValues(tgtTrii, tgtPoints_),
1796  triOwns(tgtTrii),
1797  triOtherTriPoints(tgtTrii, srcTrii),
1798  ictSrcPoints,
1799  ictSrcPointNormals,
1800  ictTgtPoints,
1801  ictPointLocations,
1802  this->debug > 3,
1803  this->debug > 4
1804  ? "srcTrii##tgtTrii=" + name(srcTrii) + name(tgtTrii)
1805  : ""
1806  );
1807 
1808  // Return false if no intersection
1809  if (!ictPointLocations.size())
1810  {
1811  return false;
1812  }
1813 
1814  // Loop over the intersection points looking for elements that are already
1815  // intersected with the opposing surface. If any are found, then this
1816  // intersection is considered invalid.
1817  forAll(ictPointLocations, ictPointi)
1818  {
1819  const triIntersect::location& l = ictPointLocations[ictPointi];
1820 
1821  if (l.isSrcPoint() && !l.isTgtPoint())
1822  {
1823  const label pointi = triPoint(srcTrii, l.srcPointi());
1824 
1825  if
1826  (
1827  this->pointTgtPoints_[pointi] != -1
1828  || this->pointTgtEdges_[pointi] != -1
1829  || this->pointTgtFaces_[pointi] != -1
1830  )
1831  {
1832  return false;
1833  }
1834  }
1835  if (!l.isSrcPoint() && l.isTgtPoint())
1836  {
1837  const label pointi = triPoint(tgtTrii, l.tgtPointi());
1838 
1839  if
1840  (
1841  this->pointSrcPoints_[pointi] != -1
1842  || this->pointSrcEdges_[pointi] != -1
1843  || this->pointSrcFaces_[pointi] != -1
1844  )
1845  {
1846  return false;
1847  }
1848  }
1849  if (l.isIntersection())
1850  {
1851  const label srcEdgei = triEdges_[srcTrii][l.srcEdgei()];
1852  const label tgtEdgei = triEdges_[tgtTrii][l.tgtEdgei()];
1853 
1854  const labelPair srcPatchEdges = edgePatchEdges(srcEdgei);
1855  const labelPair tgtPatchEdges = edgePatchEdges(tgtEdgei);
1856 
1857  if
1858  (
1859  (srcPatchEdges[0] != -1 && tgtPatchEdges[1] != -1)
1860  && (srcPatchEdges[1] != -1 || tgtPatchEdges[0] != -1)
1861  )
1862  {
1863  return false;
1864  }
1865 
1866  const label srcTrij =
1867  edgeTris_[srcEdgei][edgeTris_[srcEdgei][0] == srcTrii];
1868  const label tgtTrij =
1869  edgeTris_[tgtEdgei][edgeTris_[tgtEdgei][0] == tgtTrii];
1870 
1871  if
1872  (
1873  (srcTrij != -1 && triTgtFace_[srcTrij] != -1)
1874  || (tgtTrij != -1 && triSrcFace_[tgtTrij] != -1)
1875  )
1876  {
1877  return false;
1878  }
1879  }
1880  }
1881 
1882  // Flags controlling what operations to do
1883  bool insertPointsIntoTri = false;
1884  bool insertPointsIntoEdges = false;
1885 
1886  // Loop over the intersection points creating a map from intersection point
1887  // indices to point indices. Overwrite any intersected source or target
1888  // points, and create points generated by any new intersections. Point-tri
1889  // and point-edge associations are generated later during insertion.
1890  labelList ictPointiToPointi(ictPointLocations.size(), -1);
1891  forAll(ictPointLocations, ictPointi)
1892  {
1893  const triIntersect::location& l = ictPointLocations[ictPointi];
1894 
1895  if (l.isSrcPoint() && !l.isTgtPoint())
1896  {
1897  const label pointi = triPoint(srcTrii, l.srcPointi());
1898 
1899  ictPointiToPointi[ictPointi] = pointi;
1900 
1901  // Store the projected location on the target side
1902  tgtPoints_[pointi] = ictTgtPoints[ictPointi];
1903 
1904  insertPointsIntoTri = true;
1905  }
1906  if (!l.isSrcPoint() && l.isTgtPoint())
1907  {
1908  const label pointi = triPoint(tgtTrii, l.tgtPointi());
1909 
1910  ictPointiToPointi[ictPointi] = pointi;
1911 
1912  // Store the projected location on the source side
1913  srcPoints_[pointi] = ictSrcPoints[ictPointi];
1914  srcPointNormals_[pointi] = ictSrcPointNormals[ictPointi];
1915 
1916  insertPointsIntoTri = true;
1917  }
1918  if (l.isIntersection())
1919  {
1920  const label srcPatchEdgei =
1921  edgePatchEdge(triEdges_[srcTrii][l.srcEdgei()], true);
1922  const label tgtPatchEdgei =
1923  edgePatchEdge(triEdges_[tgtTrii][l.tgtEdgei()], false);
1924 
1925  if (srcPatchEdgei != -1 && tgtPatchEdgei != -1)
1926  {
1927  ictPointiToPointi[ictPointi] = srcPoints_.size();
1928 
1929  const label pointi = newPointi();
1930 
1931  srcPoints_[pointi] = ictSrcPoints[ictPointi];
1932  srcPointNormals_[pointi] = ictSrcPointNormals[ictPointi];
1933  tgtPoints_[pointi] = ictTgtPoints[ictPointi];
1934 
1935  insertPointsIntoEdges = true;
1936  }
1937  }
1938  }
1939 
1940  // Store the source and target tri edges and edge-ownerships
1941  const FixedList<label, 3> srcTriEdges = triEdges_[srcTrii];
1942  const FixedList<bool, 3> srcTriOwns = triOwns(srcTrii);
1943  const FixedList<label, 3> tgtTriEdges = triEdges_[tgtTrii];
1944  const FixedList<bool, 3> tgtTriOwns = triOwns(tgtTrii);
1945 
1946  // Edges which are to be fixed during insertion. These constraints are not
1947  // needed. The only edges that get inserted into are those that associate
1948  // with patch edges, and these are constrained anyway. Constraining all
1949  // edges of the tris leads to poor triangle qualities, as it limits the
1950  // that re-triangulation insertPoints can do. For testing, however, this
1951  // can be quite useful, hence the option to make the system over-constrain
1952  // the problem.
1953  DynamicList<label> fixedSrcEdgeis, fixedTgtEdgeis;
1954  #ifdef OVERCONSTRAIN
1955  forAll(srcTriEdges, srcTriEdgei)
1956  {
1957  fixedSrcEdgeis.append(srcTriEdges[srcTriEdgei]);
1958  }
1959  forAll(tgtTriEdges, tgtTriEdgei)
1960  {
1961  fixedTgtEdgeis.append(tgtTriEdges[tgtTriEdgei]);
1962  }
1963  #endif
1964 
1965  // Insert points into the triangles
1966  if (insertPointsIntoTri)
1967  {
1968  DynamicList<label> insertSrcPointis(3), insertTgtPointis(3);
1969  DynamicList<label> insertEdgeis;
1970  forAll(ictPointLocations, ictPointi)
1971  {
1972  const label pointi = ictPointiToPointi[ictPointi];
1973 
1974  if (pointi == -1) continue;
1975 
1976  const triIntersect::location& l = ictPointLocations[ictPointi];
1977 
1978  if (l.isSrcPoint() && !l.isTgtPoint())
1979  {
1980  insertSrcPointis.append(pointi);
1981  }
1982  if (!l.isSrcPoint() && l.isTgtPoint())
1983  {
1984  insertTgtPointis.append(pointi);
1985  }
1986  }
1987 
1988  if (insertSrcPointis.size())
1989  {
1990  insertPoints
1991  (
1992  tgtTrii,
1993  true,
1994  insertSrcPointis,
1995  insertEdgeis,
1996  fixedTgtEdgeis
1997  );
1998  }
1999 
2000  if (insertTgtPointis.size())
2001  {
2002  insertPoints
2003  (
2004  srcTrii,
2005  true,
2006  insertTgtPointis,
2007  insertEdgeis,
2008  fixedSrcEdgeis
2009  );
2010  }
2011 
2012  // Check
2013  checkPatchFace(srcFacei, true);
2014  checkPatchFace(tgtFacei, false);
2015  }
2016 
2017  // Insert points into the edges
2018  if (insertPointsIntoEdges)
2019  {
2020  DynamicList<label> insertPointis(2), insertEdgeis(3);
2021  forAll(ictPointLocations, ictPointi0)
2022  {
2023  const label ictPointi1 = ictPointLocations.fcIndex(ictPointi0);
2024 
2025  const label pointi0 = ictPointiToPointi[ictPointi0];
2026  const label pointi1 = ictPointiToPointi[ictPointi1];
2027 
2028  const triIntersect::location& l0 = ictPointLocations[ictPointi0];
2029  const triIntersect::location& l1 = ictPointLocations[ictPointi1];
2030 
2031  insertPointis.clear();
2032  insertEdgeis.clear();
2033  insertEdgeis.append(-1);
2034 
2035  if
2036  (
2037  (
2038  l0.isIntersection()
2039  && l1.isSrcPoint()
2040  && l0.srcEdgei() != (l1.srcPointi() + 1) % 3
2041  )
2042  || (
2043  l0.isSrcPoint()
2044  && l1.isIntersection()
2045  && (l0.srcPointi() + 1) % 3 != l1.srcEdgei()
2046  )
2047  || (
2048  l0.isIntersection()
2049  && l1.isIntersection()
2050  && l0.srcEdgei() == l1.srcEdgei()
2051  )
2052  )
2053  {
2054  const label srcTriEdgei =
2055  l0.isIntersection() ? l0.srcEdgei() : l1.srcEdgei();
2056 
2057  if (!l0.isSrcPoint() && pointi0 != -1)
2058  {
2059  insertPointis.append(pointi0);
2060  insertEdgeis.append(-1);
2061  }
2062  if (!l1.isSrcPoint() && pointi1 != -1)
2063  {
2064  insertPointis.append(pointi1);
2065  insertEdgeis.append(-1);
2066  }
2067  if (!srcTriOwns[srcTriEdgei])
2068  {
2069  inplaceReverseList(insertPointis);
2070  }
2071 
2072  if (insertPointis.size())
2073  {
2074  insertPoints
2075  (
2076  srcTriEdges[srcTriEdgei],
2077  false,
2078  insertPointis,
2079  insertEdgeis,
2080  fixedSrcEdgeis
2081  );
2082  }
2083 
2084  #ifdef OVERCONSTRAIN
2085  fixedSrcEdgeis[srcTriEdgei] = insertEdgeis.remove();
2086  fixedSrcEdgeis.append(insertEdgeis);
2087  #endif
2088  }
2089 
2090  if
2091  (
2092  (
2093  l0.isIntersection()
2094  && l1.isTgtPoint()
2095  && l0.tgtEdgei() != (l1.tgtPointi() + 1) % 3
2096  )
2097  || (
2098  l0.isTgtPoint()
2099  && l1.isIntersection()
2100  && (l0.tgtPointi() + 1) % 3 != l1.tgtEdgei()
2101  )
2102 
2103  || (
2104  l0.isIntersection()
2105  && l1.isIntersection()
2106  && l0.tgtEdgei() == l1.tgtEdgei()
2107  )
2108  )
2109  {
2110  const label tgtTriEdgei =
2111  l0.isIntersection() ? l0.tgtEdgei() : l1.tgtEdgei();
2112 
2113  if (!l0.isTgtPoint() && pointi0 != -1)
2114  {
2115  insertPointis.append(pointi0);
2116  insertEdgeis.append(-1);
2117  }
2118  if (!l1.isTgtPoint() && pointi1 != -1)
2119  {
2120  insertPointis.append(pointi1);
2121  insertEdgeis.append(-1);
2122  }
2123  if (tgtTriOwns[tgtTriEdgei])
2124  {
2125  inplaceReverseList(insertPointis);
2126  }
2127 
2128  if (insertPointis.size())
2129  {
2130  insertPoints
2131  (
2132  tgtTriEdges[tgtTriEdgei],
2133  false,
2134  insertPointis,
2135  insertEdgeis,
2136  fixedTgtEdgeis
2137  );
2138  }
2139 
2140  #ifdef OVERCONSTRAIN
2141  fixedTgtEdgeis[tgtTriEdgei] = insertEdgeis.remove();
2142  fixedTgtEdgeis.append(insertEdgeis);
2143  #endif
2144  }
2145  }
2146 
2147  // Check
2148  checkPatchFace(srcFacei, true);
2149  checkPatchFace(tgtFacei, false);
2150  }
2151 
2152  return true;
2153 }
2154 
2155 
2156 template<class SrcPatchType, class TgtPatchType>
2157 void
2158 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::intersectPatchFaceTris
2159 (
2160  const label srcFacei,
2161  const label tgtFacei
2162 )
2163 {
2164  // Get the source face tris and make sure none are marked
2165  const DynamicList<label>& srcFaceTris = srcFaceTris_[srcFacei];
2166  forAll(srcFaceTris, srcFaceTrii)
2167  {
2168  const label srcTrii = srcFaceTris[srcFaceTrii];
2169 
2170  if (triMarkedTris_[srcTrii] != -1)
2171  {
2172  markedTriTris_[triMarkedTris_[srcTrii]] = -1;
2173  triMarkedTris_[srcTrii] = -1;
2174  }
2175  }
2176 
2177  // Loop the source face tris
2178  label srcFaceTrii = 0;
2179  while (srcFaceTrii < srcFaceTris.size())
2180  {
2181  // Mark this source tri
2182  const label srcTrii = srcFaceTris[srcFaceTrii];
2183  markedTriTris_.append(srcTrii);
2184  triMarkedTris_[srcTrii] = markedTriTris_.size() - 1;
2185 
2186  // Get the target face tris and make sure none are marked
2187  const DynamicList<label>& tgtFaceTris = tgtFaceTris_[tgtFacei];
2188  forAll(tgtFaceTris, tgtFaceTrii)
2189  {
2190  const label tgtTrii = tgtFaceTris[tgtFaceTrii];
2191 
2192  if (triMarkedTris_[tgtTrii] != -1)
2193  {
2194  markedTriTris_[triMarkedTris_[tgtTrii]] = -1;
2195  triMarkedTris_[tgtTrii] = -1;
2196  }
2197  }
2198 
2199  // Loop the target face tris
2200  label tgtFaceTrii = 0;
2201  while (tgtFaceTrii < tgtFaceTris.size())
2202  {
2203  // Mark this target tri
2204  const label tgtTrii = tgtFaceTris[tgtFaceTrii];
2205  markedTriTris_.append(tgtTrii);
2206  triMarkedTris_[tgtTrii] = markedTriTris_.size() - 1;
2207 
2208  /*
2209  // If the target tri has not unmarked then snap
2210  if (triMarkedTris_[tgtFaceTris[tgtFaceTrii]] != -1)
2211  {
2212  snapTris(srcTrii, tgtTrii);
2213  }
2214 
2215  // If the source tri has unmarked then break
2216  if (triMarkedTris_[srcFaceTris[srcFaceTrii]] == -1) break;
2217  */
2218 
2219  // If the target tri has not unmarked then intersect
2220  if (triMarkedTris_[tgtFaceTris[tgtFaceTrii]] != -1)
2221  {
2222  intersectTris(srcTrii, tgtTrii);
2223  }
2224 
2225  // If the source tri has unmarked then break
2226  if (triMarkedTris_[srcFaceTris[srcFaceTrii]] == -1) break;
2227 
2228  // Iterate backwards to find the next marked target tri
2229  while
2230  (
2231  tgtFaceTrii >= 0
2232  && triMarkedTris_[tgtFaceTris[tgtFaceTrii]] == -1
2233  )
2234  {
2235  -- tgtFaceTrii;
2236  }
2237 
2238  // Increment to the next unmarked target tri
2239  ++ tgtFaceTrii;
2240  }
2241 
2242  // Iterate backwards to find the next marked source tri
2243  while
2244  (
2245  srcFaceTrii >= 0
2246  && triMarkedTris_[srcFaceTris[srcFaceTrii]] == -1
2247  )
2248  {
2249  -- srcFaceTrii;
2250  }
2251 
2252  // Increment to the next unmarked source tri
2253  ++ srcFaceTrii;
2254  }
2255 
2256  // Clear the marked tris
2257  forAll(markedTriTris_, markedTrii)
2258  {
2259  const label trii = markedTriTris_[markedTrii];
2260  if (trii != -1)
2261  {
2262  triMarkedTris_[trii] = -1;
2263  }
2264  }
2265  markedTriTris_.clear();
2266 }
2267 
2268 
2269 template<class SrcPatchType, class TgtPatchType>
2270 bool
2271 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::conformPatchFaceTris
2272 (
2273  const label patchFacei,
2274  const label otherPatchFacei,
2275  const bool isSrc
2276 )
2277 {
2278  const labelList& patchFaceTris =
2279  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
2280  const labelList& otherPatchFaceTris =
2281  isSrc ? tgtFaceTris_[otherPatchFacei] : srcFaceTris_[otherPatchFacei];
2282 
2283  const labelList& triOtherPatchFaces =
2284  isSrc ? triTgtFace_ : triSrcFace_;
2285 
2286  const pointField& points = isSrc ? srcPoints_ : tgtPoints_;
2287 
2288  const vector& patchFaceN =
2289  isSrc
2290  ? this->srcPatch_.faceNormals()[patchFacei]
2291  : this->tgtPatch_.faceNormals()[patchFacei];
2292 
2293  auto pointIntersectsPatchFace = [&]
2294  (
2295  const label pointi,
2296  const label patchFacei,
2297  const bool isSrc
2298  )
2299  {
2300  const labelList& pointPatchPoints =
2301  isSrc ? this->pointSrcPoints_ : this->pointTgtPoints_;
2302  const labelList& pointPatchEdges =
2303  isSrc ? this->pointSrcEdges_ : this->pointTgtEdges_;
2304  const labelList& pointPatchFaces =
2305  isSrc ? this->pointSrcFaces_ : this->pointTgtFaces_;
2306 
2307  const triFace patchFacePatchPoints
2308  (
2309  this->patchFacePatchPoints(patchFacei, isSrc)
2310  );
2311  const triFace patchFacePatchEdges
2312  (
2313  this->patchFacePatchEdges(patchFacei, isSrc)
2314  );
2315 
2316  return
2317  findIndex
2318  (
2319  patchFacePatchPoints,
2320  pointPatchPoints[pointi]
2321  ) != -1
2322  || findIndex
2323  (
2324  patchFacePatchEdges,
2325  pointPatchEdges[pointi]
2326  ) != -1
2327  || pointPatchFaces[pointi] == patchFacei;
2328  };
2329 
2330  auto edgeIntersectsPatchFace = [&]
2331  (
2332  const label edgei,
2333  const label patchFacei,
2334  const bool isSrc
2335  )
2336  {
2337  const edge e = edgePoints(edgei);
2338  return
2339  pointIntersectsPatchFace(e[0], patchFacei, isSrc)
2340  && pointIntersectsPatchFace(e[1], patchFacei, isSrc);
2341  };
2342 
2343  // Get all the edges that need conforming to
2344  DynamicList<label> conformEdgeis;
2345  forAll(otherPatchFaceTris, otherPatchFaceTrii)
2346  {
2347  const label trii = otherPatchFaceTris[otherPatchFaceTrii];
2348 
2349  forAll(triEdges_[trii], triEdgei)
2350  {
2351  const label edgei = triEdges_[trii][triEdgei];
2352 
2353  const label trij =
2354  edgeTris_[edgei][edgeTris_[edgei][0] == trii];
2355 
2356  const label otherPatchFacej =
2357  trij == -1 ? -1 : triOtherPatchFaces[trij];
2358 
2359  if
2360  (
2361  otherPatchFacei != otherPatchFacej
2362  && edgeIntersectsPatchFace(edgei, patchFacei, isSrc)
2363  )
2364  {
2365  conformEdgeis.append(edgei);
2366  }
2367  }
2368  }
2369 
2370  // Remove edges already present in the triangulation by shuffling up
2371  auto isConformed = [&](const label edgei)
2372  {
2373  const edge e = edgePoints(edgei);
2374 
2375  forAll(patchFaceTris, patchFaceTrii)
2376  {
2377  const label trii = patchFaceTris[patchFaceTrii];
2378 
2379  forAll(triEdges_[trii], triEdgei)
2380  {
2381  const label edgei = triEdges_[trii][triEdgei];
2382 
2383  if (edge::compare(e, edgePoints(edgei)) != 0)
2384  {
2385  return true;
2386  }
2387  }
2388  }
2389 
2390  return false;
2391  };
2392  {
2393  label conformi = 0;
2394  forAll(conformEdgeis, conformj)
2395  {
2396  if (!isConformed(conformEdgeis[conformj]))
2397  {
2398  conformEdgeis[conformi] = conformEdgeis[conformj];
2399  ++ conformi;
2400  }
2401  }
2402  conformEdgeis.resize(conformi);
2403  }
2404 
2405  // Quit if there is nothing to conform to
2406  if (conformEdgeis.empty()) return true;
2407 
2408  // Get the next patch face tri index connected to a given point
2409  auto nextConnectedPatchFaceTri = [&]
2410  (
2411  const label pointi,
2412  const label patchFaceTrii0 = 0
2413  )
2414  {
2415  for
2416  (
2417  label patchFaceTrii = patchFaceTrii0;
2418  patchFaceTrii < patchFaceTris.size();
2419  ++ patchFaceTrii
2420  )
2421  {
2422  const label trii = patchFaceTris[patchFaceTrii];
2423 
2424  if (findIndex(triPoints(trii), pointi) != -1)
2425  {
2426  return patchFaceTrii;
2427  }
2428  }
2429 
2430  return label(-1);
2431  };
2432 
2433  // Track from a point on a triangle towards a given point. Stop at an edge
2434  // and set the index of that edge and update the local coordinates.
2435  auto trackToEdge = [&]
2436  (
2437  const label trii,
2438  label& triEdgei,
2439  barycentric2D& y,
2440  const label pointi1
2441  )
2442  {
2443  const triPointRef t = triPoints(trii).tri(points);
2444 
2445  const barycentricTensor2D A(t.a(), t.b(), t.c());
2446 
2447  const point p = A & y;
2448  const vector dp = points[pointi1] - p;
2449 
2450  const vector ab = t.b() - t.a();
2451  const vector ac = t.c() - t.a();
2452  const vector bc = t.c() - t.b();
2453  const scalar detA = (ab ^ ac) & patchFaceN;
2454  const barycentricTensor2D T
2455  (
2456  patchFaceN ^ bc,
2457  ac ^ patchFaceN,
2458  patchFaceN ^ ab
2459  );
2460 
2461  const barycentric2D TDp = dp & T;
2462 
2463  label iH = -1;
2464  scalar lambdaByDetAH =
2465  !std::isnormal(detA) || detA < 0 ? vGreat : 1/detA;
2466 
2467  forAll(TDp, i)
2468  {
2469  if (TDp[i] < - detA*small)
2470  {
2471  const scalar lambdaByDetA = - y[i]/TDp[i];
2472 
2473  if (0 <= lambdaByDetA && lambdaByDetA < lambdaByDetAH)
2474  {
2475  iH = i;
2476  lambdaByDetAH = lambdaByDetA;
2477  }
2478  }
2479  }
2480 
2481  y += lambdaByDetAH*TDp;
2482  forAll(y, i)
2483  {
2484  y.replace(i, i == iH ? 0 : max(0, y[i]));
2485  }
2486  if (iH == -1)
2487  {
2488  y /= cmptSum(y);
2489  }
2490 
2491  triEdgei = iH == -1 ? -1 : (iH + 1) % 3;
2492  };
2493 
2494  // Cross the edge. Return true if the edge can be crossed, and false
2495  // otherwise. Set the new triangle, new edge index, and transform the local
2496  // coordinates appropriately.
2497  auto crossEdge = [&](label& trii, label& triEdgei, barycentric2D& y)
2498  {
2499  if (triEdgei == -1) return false;
2500 
2501  const label edgei = triEdges_[trii][triEdgei];
2502 
2503  if (edgePatchEdges(edgei) != labelPair(-1, -1)) return false;
2504 
2505  const label trij = edgeTris_[edgei][edgeTris_[edgei][0] == trii];
2506 
2507  if (trij == -1) return false;
2508 
2509  if (triSrcFace_[trii] != triSrcFace_[trij]) return false;
2510  if (triTgtFace_[trii] != triTgtFace_[trij]) return false;
2511 
2512  const label triEdgej = findIndex(triEdges_[trij], edgei);
2513 
2514  auto inplaceRotate = [](barycentric2D& y, label n)
2515  {
2516  n = n % 3;
2517 
2518  if (n == 1)
2519  {
2520  scalar temp = y.a();
2521  y.a() = y.b();
2522  y.b() = y.c();
2523  y.c() = temp;
2524  }
2525 
2526  if (n == 2)
2527  {
2528  scalar temp = y.c();
2529  y.c() = y.b();
2530  y.b() = y.a();
2531  y.a() = temp;
2532  }
2533  };
2534 
2535  inplaceRotate(y, triEdgei - 1 + 3);
2536  Swap(y.b(), y.c());
2537  inplaceRotate(y, 1 - triEdgej + 3);
2538 
2539  trii = trij;
2540  triEdgei = triEdgej;
2541 
2542  return true;
2543  };
2544 
2545  // Assume successful
2546  bool success = true;
2547 
2548  // Conform each identified edge in turn
2549  forAll(conformEdgeis, conformi)
2550  {
2551  const label edgei = conformEdgeis[conformi];
2552 
2553  // Skip if earlier operations have inadvertently resulted in this
2554  // edge being conformed to
2555  if (isConformed(edgei)) continue;
2556 
2557  // Get the edge and the points
2558  const edge e = edgePoints(edgei);
2559  const label pointi0 = e[0], pointi1 = e[1];
2560 
2561  // Find triangles which contains the first and last points
2562  const label patchFaceTrii0 = nextConnectedPatchFaceTri(pointi0);
2563  const label patchFaceTrii1 = nextConnectedPatchFaceTri(pointi1);
2564 
2565  // Conformation to this edge is not possible if either end point is not
2566  // present in the patch face triangulation
2567  if (patchFaceTrii0 == -1 || patchFaceTrii1 == -1) continue;
2568 
2569  // Track from point zero until in a triangle which contains point one.
2570  // Build a route of tris and edges along the track.
2571  DynamicList<label> routeTriis;
2572  DynamicList<label> routeEdgeis;
2573  label routePatchFaceTrii0 = -1;
2574  while (true)
2575  {
2576  routeTriis.clear();
2577  routeEdgeis.clear();
2578 
2579  routePatchFaceTrii0 =
2580  nextConnectedPatchFaceTri(pointi0, routePatchFaceTrii0 + 1);
2581 
2582  if (routePatchFaceTrii0 == -1) break;
2583 
2584  label trii = patchFaceTris[routePatchFaceTrii0];
2585  label triEdgei = -1;
2586  barycentric2D y(0, 0, 0);
2587  y[findIndex(triPoints(trii), pointi0)] = 1;
2588 
2589  forAll(patchFaceTris, iter)
2590  {
2591  if (findIndex(triPoints(trii), pointi0) != -1)
2592  {
2593  routeTriis.clear();
2594  routeEdgeis.clear();
2595  }
2596 
2597  routeTriis.append(trii);
2598 
2599  if (findIndex(triPoints(trii), pointi1) != -1) break;
2600 
2601  trackToEdge(trii, triEdgei, y, pointi1);
2602 
2603  if (triEdgei == -1) break;
2604 
2605  routeEdgeis.append(triEdges_[trii][triEdgei]);
2606 
2607  if (!crossEdge(trii, triEdgei, y)) break;
2608  }
2609 
2610  if (findIndex(triPoints(trii), pointi1) != -1) break;
2611  }
2612 
2613  // Conform the triangulation to the route
2614  if (routeEdgeis.size() == 0)
2615  {
2616  // No route was found. Don't do anything. This edge will not be
2617  // conformed to and the intersection will disconnect here.
2618  if (this->debug > 1)
2619  {
2620  WarningInFunction
2621  << indent << "Failed to route edge " << e
2622  << " through the triangulation of "
2623  << (isSrc ? "source" : "target") << " face #"
2624  << patchFacei << endl;
2625  writePatchFace(patchFacei, isSrc);
2626  }
2627  success = false;
2628  }
2629  else if (routeEdgeis.size() == 1)
2630  {
2631  // The route only spans a single edge. Flip it to conform.
2632  flipEdge(routeEdgeis.first());
2633  }
2634  else
2635  {
2636  // The route spans multiple edges. Remove all triangles in the
2637  // route and then split the resulting polygon along the conform
2638  // edge and triangulate both sides independently.
2639 
2640  // !!! Check that removal of the tris does not result in any points
2641  // becoming disconnected from the mesh
2642  // ...
2643 
2644  // Form the polygons on either side of the conform edge
2645  DynamicList<label> leftPolyPointis, rightPolyPointis;
2646  DynamicList<label> leftPolyEdgeis, rightPolyEdgeis;
2647  DynamicList<label> leftPolyTriis, rightPolyTriis;
2648  {
2649  // First point
2650  leftPolyPointis.append(pointi0);
2651  rightPolyPointis.append(pointi0);
2652 
2653  // First triangle
2654  {
2655  const label trii = routeTriis.first();
2656  const label edgei1 = routeEdgeis.first();
2657  const label triEdgei1 = findIndex(triEdges_[trii], edgei1);
2658  const label triEdgeiLeft = (triEdgei1 + 2) % 3;
2659  const label triEdgeiRight = (triEdgei1 + 1) % 3;
2660  leftPolyEdgeis.append(triEdges_[trii][triEdgeiLeft]);
2661  rightPolyEdgeis.append(triEdges_[trii][triEdgeiRight]);
2662  leftPolyTriis.append(trii);
2663  rightPolyTriis.append(trii);
2664  leftPolyPointis.append(triPoint(trii, triEdgei1));
2665  rightPolyPointis.append(triPoint(trii, triEdgeiRight));
2666  }
2667 
2668  // Intermediate triangles
2669  for
2670  (
2671  label routeTrii = 1;
2672  routeTrii < routeTriis.size() - 1;
2673  ++ routeTrii
2674  )
2675  {
2676  const label trii = routeTriis[routeTrii];
2677  const label edgei0 = routeEdgeis[routeTrii - 1];
2678  const label edgei1 = routeEdgeis[routeTrii];
2679  const label triEdgei0 = findIndex(triEdges_[trii], edgei0);
2680  const label triEdgei1 = findIndex(triEdges_[trii], edgei1);
2681  if ((triEdgei0 + 2) % 3 == triEdgei1)
2682  {
2683  const label triEdgeiLeft = (triEdgei1 + 2) % 3;
2684  leftPolyEdgeis.append(triEdges_[trii][triEdgeiLeft]);
2685  leftPolyTriis.append(trii);
2686  leftPolyPointis.append(triPoint(trii, triEdgei1));
2687  }
2688  else
2689  {
2690  const label triEdgeiRight = (triEdgei1 + 1) % 3;
2691  rightPolyEdgeis.append(triEdges_[trii][triEdgeiRight]);
2692  rightPolyTriis.append(trii);
2693  rightPolyPointis.append(triPoint(trii, triEdgeiRight));
2694  }
2695  }
2696 
2697  // Last triangle
2698  {
2699  const label trii = routeTriis.last();
2700  const label edgei0 = routeEdgeis.last();
2701  const label triEdgei0 = findIndex(triEdges_[trii], edgei0);
2702  const label triEdgeiLeft = (triEdgei0 + 1) % 3;
2703  const label triEdgeiRight = (triEdgei0 + 2) % 3;
2704  leftPolyEdgeis.append(triEdges_[trii][triEdgeiLeft]);
2705  rightPolyEdgeis.append(triEdges_[trii][triEdgeiRight]);
2706  leftPolyTriis.append(trii);
2707  rightPolyTriis.append(trii);
2708  leftPolyPointis.append(pointi1);
2709  rightPolyPointis.append(pointi1);
2710  }
2711  }
2712 
2713  // Disconnect and remove existing tris
2714  auto disconnect = [&]
2715  (
2716  const labelUList& edgeis,
2717  const labelList& triis
2718  )
2719  {
2720  forAll(edgeis, i)
2721  {
2722  const label edgei = edgeis[i];
2723  const label trii = triis[i];
2724 
2725  const label edgeTrii = edgeTris_[edgei][0] != trii;
2726  const label triEdgei = findIndex(triEdges_[trii], edgei);
2727 
2728  const label tempEdgei = newEdgei();
2729  triEdges_[trii][triEdgei] = tempEdgei;
2730  edgeTris_[tempEdgei][0] = trii;
2731  edgeTris_[edgei][edgeTrii] = -1;
2732  }
2733  };
2734  disconnect(leftPolyEdgeis, leftPolyTriis);
2735  disconnect(rightPolyEdgeis, rightPolyTriis);
2736  forAll(routeTriis, routeTrii)
2737  {
2738  removeTri(routeTriis[routeTrii]);
2739  }
2740 
2741  // Create the new conformed edge
2742  const label middleEdgei = newEdgei();
2743  leftPolyEdgeis.append(middleEdgei);
2744  rightPolyEdgeis.append(middleEdgei);
2745 
2746  // Reverse the right polygon
2747  SubList<label> subRightPolyPointis
2748  (
2749  rightPolyPointis,
2750  rightPolyPointis.size() - 1,
2751  1
2752  );
2753  inplaceReverseList(subRightPolyPointis);
2754  inplaceReverseList(rightPolyEdgeis);
2755  inplaceReverseList(rightPolyTriis);
2756 
2757  // Insert the polygons into the triangulation
2758  auto insert = [&]
2759  (
2760  DynamicList<label>& polyPointis,
2761  DynamicList<label>& polyEdgeis
2762  )
2763  {
2764  // Generate new edges
2765  polyEdgeis.resize(2*polyPointis.size() - 3, -1);
2766  for (label i = polyPointis.size(); i < polyEdgeis.size(); ++ i)
2767  {
2768  polyEdgeis[i] = newEdgei();
2769  }
2770 
2771  // Triangulate
2772  polygonTriangulate_.triangulate
2773  (
2774  UIndirectList<point>(points, polyPointis)
2775  );
2776 
2777  // Insert the triangulation
2778  forAll(polygonTriangulate_.triPoints(), trii)
2779  {
2780  addTri
2781  (
2782  polygonTriangulate_.triPoints(trii, polyPointis),
2783  polygonTriangulate_.triEdges(trii, polyEdgeis),
2784  patchFacei,
2785  isSrc
2786  );
2787  }
2788  };
2789  insert(leftPolyPointis, leftPolyEdgeis);
2790  insert(rightPolyPointis, rightPolyEdgeis);
2791  }
2792  }
2793 
2794  // Check
2795  checkPatchFace(patchFacei, isSrc);
2796 
2797  return success;
2798 }
2799 
2800 
2801 template<class SrcPatchType, class TgtPatchType>
2802 bool
2803 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::conformPatchFaceTris
2804 (
2805  const label srcFacei,
2806  const label tgtFacei
2807 )
2808 {
2809  return
2810  conformPatchFaceTris(srcFacei, tgtFacei, true)
2811  && conformPatchFaceTris(tgtFacei, srcFacei, false);
2812 }
2813 
2814 
2815 template<class SrcPatchType, class TgtPatchType>
2816 bool
2817 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::combinePatchFaceTris
2818 (
2819  const label srcFacei,
2820  const label tgtFacei
2821 )
2822 {
2823  // Function to determine whether or not a given triangle is fully
2824  // intersected with the opposite side
2825  auto triIsIntersected = [&]
2826  (
2827  const label trii,
2828  const label otherPatchFacei,
2829  const triFace& otherPatchFacePatchPoints,
2830  const triFace& otherPatchFacePatchEdges,
2831  const bool isSrc
2832  )
2833  {
2834  const labelList& pointOtherPatchPoints =
2835  isSrc ? this->pointTgtPoints_ : this->pointSrcPoints_;
2836  const labelList& pointOtherPatchEdges =
2837  isSrc ? this->pointTgtEdges_ : this->pointSrcEdges_;
2838  const labelList& pointOtherPatchFaces =
2839  isSrc ? this->pointTgtFaces_ : this->pointSrcFaces_;
2840 
2841  forAll(triPoints_[trii], triPointi)
2842  {
2843  const label pointi = triPoint(trii, triPointi);
2844 
2845  if
2846  (
2847  // If the point does not intersect an other patch point ...
2848  findIndex
2849  (
2850  otherPatchFacePatchPoints,
2851  pointOtherPatchPoints[pointi]
2852  ) == -1
2853  // ... and does not intersect an other patch edge ...
2854  && findIndex
2855  (
2856  otherPatchFacePatchEdges,
2857  pointOtherPatchEdges[pointi]
2858  ) == -1
2859  // ... and not intersect the other patch face ...
2860  && pointOtherPatchFaces[pointi] != otherPatchFacei
2861  )
2862  {
2863  // ... then this triangle does not entirely intersect the
2864  // other patch face.
2865  return false;
2866  }
2867  }
2868 
2869  return true;
2870  };
2871 
2872  // Get the polygons which comprise all of the intersected triangles
2873  auto getIntersectionPolygon = [&]
2874  (
2875  const label patchFacei,
2876  const label otherPatchFacei,
2877  const bool isSrc,
2878  DynamicList<label>& triis,
2879  DynamicList<label>& edgeis
2880  )
2881  {
2882  const labelList& patchFaceTris =
2883  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
2884  const labelList& triPatchFaces =
2885  isSrc ? triSrcFace_ : triTgtFace_;
2886 
2887  const triFace otherPatchFacePatchPoints
2888  (
2889  patchFacePatchPoints(otherPatchFacei, !isSrc)
2890  );
2891  const triFace otherPatchFacePatchEdges
2892  (
2893  patchFacePatchEdges(otherPatchFacei, !isSrc)
2894  );
2895 
2896  // Get an initial intersected tri
2897  label trii0 = -1;
2898  forAll(patchFaceTris, patchFaceTrii)
2899  {
2900  const label trii = patchFaceTris[patchFaceTrii];
2901  if
2902  (
2903  triIsIntersected
2904  (
2905  trii,
2906  otherPatchFacei,
2907  otherPatchFacePatchPoints,
2908  otherPatchFacePatchEdges,
2909  isSrc
2910  )
2911  )
2912  {
2913  trii0 = trii;
2914  break;
2915  }
2916  }
2917  if (trii0 == -1) return;
2918 
2919  // Populate the star
2920  star::context starContext = star_.populate
2921  (
2922  trii0,
2923  true,
2924  [&](const label edgei, const label trii)
2925  {
2926  return
2927  triPatchFaces[trii] == patchFacei
2928  && triIsIntersected
2929  (
2930  trii,
2931  otherPatchFacei,
2932  otherPatchFacePatchPoints,
2933  otherPatchFacePatchEdges,
2934  isSrc
2935  );
2936  },
2937  triEdges_,
2938  edgeTris_
2939  );
2940 
2941  // Add triangles to the list
2942  forAllStarFaces(star_, starTrii, trii)
2943  {
2944  triis.append(trii);
2945  }
2946 
2947  // Walk around the star polygon, adding edges to the list
2948  forAllStarEdges(star_, i, starEdgei, edgei)
2949  {
2950  edgeis.append(edgei);
2951  }
2952  };
2953  DynamicList<label> srcPolyTris, srcPolyEdges;
2954  DynamicList<label> tgtPolyTris, tgtPolyEdges;
2955  getIntersectionPolygon
2956  (
2957  srcFacei,
2958  tgtFacei,
2959  true,
2960  srcPolyTris,
2961  srcPolyEdges
2962  );
2963  getIntersectionPolygon
2964  (
2965  tgtFacei,
2966  srcFacei,
2967  false,
2968  tgtPolyTris,
2969  tgtPolyEdges
2970  );
2971 
2972  // Return if there is nothing to do
2973  if (srcPolyEdges.size() == 0 && tgtPolyEdges.size() == 0)
2974  {
2975  return true;
2976  }
2977 
2978  // Align the source and target polygon edges
2979  if (srcPolyEdges.size() && tgtPolyEdges.size())
2980  {
2981  inplaceReverseList(tgtPolyEdges);
2982 
2983  label tgtPolyEdgei0 = 0;
2984  forAll(tgtPolyEdges, tgtPolyEdgei)
2985  {
2986  const edge srcE = edgePoints(srcPolyEdges[0]);
2987  const edge tgtE = edgePoints(tgtPolyEdges[tgtPolyEdgei]);
2988 
2989  if (edge::compare(srcE, tgtE) != 0)
2990  {
2991  tgtPolyEdgei0 = tgtPolyEdgei;
2992  break;
2993  }
2994  }
2995 
2996  inplaceRotateList
2997  (
2998  static_cast<labelList&>(tgtPolyEdges),
2999  - tgtPolyEdgei0
3000  );
3001  }
3002  bool aligned;
3003  if (srcPolyEdges.size() == tgtPolyEdges.size())
3004  {
3005  aligned = true;
3006  forAll(srcPolyEdges, polyEdgei)
3007  {
3008  const edge srcE = edgePoints(srcPolyEdges[polyEdgei]);
3009  const edge tgtE = edgePoints(tgtPolyEdges[polyEdgei]);
3010 
3011  if (edge::compare(srcE, tgtE) == 0)
3012  {
3013  aligned = false;
3014  break;
3015  }
3016  }
3017  }
3018  else
3019  {
3020  aligned = false;
3021  }
3022  if (!aligned)
3023  {
3024  if (this->debug > 1)
3025  {
3026  WarningInFunction
3027  << indent << "Failed to combine intersected parts of source "
3028  << "face #" << srcFacei << " and target face #" << tgtFacei
3029  << endl;
3030  writePatchFace(srcFacei, true);
3031  writePatchFace(tgtFacei, false);
3032  }
3033  return false;
3034  }
3035 
3036  // Mark triangles that belong to either polygon as candidates
3037  forAll(srcPolyTris, srcPolyTrii)
3038  {
3039  const label srcTrii = srcPolyTris[srcPolyTrii];
3040  markedTriTris_.append(srcTrii);
3041  triMarkedTris_[srcTrii] = markedTriTris_.size() - 1;
3042  }
3043  forAll(tgtPolyTris, tgtPolyTrii)
3044  {
3045  const label tgtTrii = tgtPolyTris[tgtPolyTrii];
3046  markedTriTris_.append(tgtTrii);
3047  triMarkedTris_[tgtTrii] = markedTriTris_.size() - 1;
3048  }
3049 
3050  // Change the aligned edges so that they connect opposite sides, thereby
3051  // disconnecting the intersected triangles from the rest of their surfaces.
3052  // Tris in the polygon are connected to the source edge and non-star tris
3053  // are connected to the target edge. If there are no non-star tris
3054  // connected to a pair of aligned edges then the target edge will not be
3055  // connected to any triangles as a result of this operation. However, we do
3056  // not remove the target edge in this case because we need it later to
3057  // attach the new intersected face to.
3058  forAll(srcPolyEdges, polyEdgei)
3059  {
3060  const label srcEdgei = srcPolyEdges[polyEdgei];
3061  const label tgtEdgei = tgtPolyEdges[polyEdgei];
3062 
3063  if (srcEdgei == tgtEdgei) continue;
3064 
3065  const label srcEdgeTrii =
3066  edgeTris_[srcEdgei][0] == -1
3067  || triMarkedTris_[edgeTris_[srcEdgei][0]] == -1
3068  || triSrcFace_[edgeTris_[srcEdgei][0]] == -1;
3069  const label tgtEdgeTrii =
3070  edgeTris_[tgtEdgei][0] == -1
3071  || triMarkedTris_[edgeTris_[tgtEdgei][0]] == -1
3072  || triTgtFace_[edgeTris_[tgtEdgei][0]] == -1;
3073 
3074  const label tgtTrii = edgeTris_[tgtEdgei][tgtEdgeTrii];
3075  const label srcTrij = edgeTris_[srcEdgei][!srcEdgeTrii];
3076  const label tgtTrij = edgeTris_[tgtEdgei][!tgtEdgeTrii];
3077 
3078  edgeTris_[srcEdgei][!srcEdgeTrii] = tgtTrii;
3079  edgeTris_[tgtEdgei][tgtEdgeTrii] = srcTrij;
3080 
3081  const label tgtTriEdgei = findIndex(triEdges_[tgtTrii], tgtEdgei);
3082  triEdges_[tgtTrii][tgtTriEdgei] = srcEdgei;
3083 
3084  if (srcTrij != -1)
3085  {
3086  const label srcTriEdgej = findIndex(triEdges_[srcTrij], srcEdgei);
3087  triEdges_[srcTrij][srcTriEdgej] = tgtEdgei;
3088  }
3089 
3090  if (srcTrij != -1 && tgtTrij != -1)
3091  {
3092  edgeFrontEdges_[tgtEdgei] = frontEdgeEdges_.size();
3093  frontEdgeEdges_.append(tgtEdgei);
3094  }
3095  }
3096 
3097  // Circle the edges to create the intersection face
3098  const label facei = this->faces().size();
3099  this->faces_.append(face(srcPolyEdges.size()));
3100  faceEdges_.append(labelList(srcPolyEdges.size()));
3101  this->srcFaceFaces_[srcFacei].append(facei);
3102  this->tgtFaceFaces_[tgtFacei].append(facei);
3103  this->faceSrcFaces_.append(srcFacei);
3104  this->faceTgtFaces_.append(tgtFacei);
3105  forAll(srcPolyEdges, polyEdgei)
3106  {
3107  const label srcEdgei = srcPolyEdges[polyEdgei];
3108  const label srcEdgeTrii =
3109  edgeTris_[srcEdgei][0] == -1
3110  || triMarkedTris_[edgeTris_[srcEdgei][0]] == -1
3111  || triSrcFace_[edgeTris_[srcEdgei][0]] == -1;
3112 
3113  const label srcTrii = edgeTris_[srcEdgei][srcEdgeTrii];
3114  const label srcTriEdgei = findIndex(triEdges_[srcTrii], srcEdgei);
3115 
3116  this->faces_.last()[polyEdgei] =
3117  triEdgePoints(srcTrii, srcTriEdgei).start();
3118 
3119  const label tgtEdgei = tgtPolyEdges[polyEdgei];
3120 
3121  faceEdges_.last()[polyEdgei] = tgtEdgei;
3122  intersectEdgeFaces_[tgtEdgei][intersectEdgeFaces_[tgtEdgei][0] != -1] =
3123  facei;
3124  }
3125 
3126  // Clear the marked tris
3127  forAll(markedTriTris_, candidateTrii)
3128  {
3129  const label trii = markedTriTris_[candidateTrii];
3130  if (trii != -1)
3131  {
3132  triMarkedTris_[trii] = -1;
3133  }
3134  }
3135  markedTriTris_.clear();
3136 
3137  // Check
3138  checkPatchFace(srcFacei, true);
3139  checkPatchFace(tgtFacei, false);
3140 
3141  // Remove the intersection triangles
3142  forAll(srcPolyTris, srcPolyTrii)
3143  {
3144  removeTri(srcPolyTris[srcPolyTrii]);
3145  }
3146  forAll(tgtPolyTris, tgtPolyTrii)
3147  {
3148  removeTri(tgtPolyTris[tgtPolyTrii]);
3149  }
3150 
3151  // Check
3152  checkPatchFace(srcFacei, true);
3153  checkPatchFace(tgtFacei, false);
3154 
3155  // Update the propagation front
3156  forAll(tgtPolyEdges, polyEdgei)
3157  {
3158  const label edgei = tgtPolyEdges[polyEdgei];
3159 
3160  const label trii0 = edgeTris_[edgei][0];
3161  const label trii1 = edgeTris_[edgei][1];
3162 
3163  const bool isFront =
3164  trii0 != -1
3165  && trii1 != -1
3166  && (triSrcFace_[trii0] == -1) != (triSrcFace_[trii1] == -1);
3167 
3168  if (isFront && edgeFrontEdges_[edgei] == -1)
3169  {
3170  edgeFrontEdges_[edgei] = frontEdgeEdges_.size();
3171  frontEdgeEdges_.append(edgei);
3172  }
3173 
3174  if (!isFront && edgeFrontEdges_[edgei] != -1)
3175  {
3176  frontEdgeEdges_[edgeFrontEdges_[edgei]] = -1;
3177  edgeFrontEdges_[edgei] = -1;
3178  }
3179  }
3180 
3181  // Check
3182  checkPatchFace(srcFacei, true);
3183  checkPatchFace(tgtFacei, false);
3184 
3185  return true;
3186 }
3187 
3188 
3189 template<class SrcPatchType, class TgtPatchType>
3190 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::initialise
3191 (
3192  const vectorField& srcPointNormals
3193 )
3194 {
3195  // Clear the base class data ...
3196 
3197  this->points_.clear();
3198 
3199  this->srcPointPoints_ = -1;
3200  this->tgtPointPoints_ = -1;
3201  this->pointSrcPoints_.clear();
3202  this->pointTgtPoints_.clear();
3203 
3204  forAll(this->srcEdgePoints_, srcEdgei)
3205  {
3206  this->srcEdgePoints_[srcEdgei].clear();
3207  }
3208  forAll(this->tgtEdgePoints_, tgtEdgei)
3209  {
3210  this->tgtEdgePoints_[tgtEdgei].clear();
3211  }
3212  this->pointSrcEdges_.clear();
3213  this->pointTgtEdges_.clear();
3214 
3215  this->pointSrcFaces_.clear();
3216  this->pointTgtFaces_.clear();
3217 
3218  this->faces_.clear();
3219 
3220  forAll(this->srcFaceFaces_, srcFacei)
3221  {
3222  this->srcFaceFaces_[srcFacei].clear();
3223  }
3224  forAll(this->tgtFaceFaces_, tgtFacei)
3225  {
3226  this->tgtFaceFaces_[tgtFacei].clear();
3227  }
3228  this->faceSrcFaces_.clear();
3229  this->faceTgtFaces_.clear();
3230 
3231  // Initialise with the source and target patch data ...
3232 
3233  // Points
3234  const label nPoints =
3235  this->srcPatch_.nPoints() + this->tgtPatch_.nPoints();
3236  srcPoints_.resize(nPoints, point::uniform(NaN));
3237  srcPointNormals_.resize(nPoints, vector::uniform(NaN));
3238  tgtPoints_.resize(nPoints, point::uniform(NaN));
3239  pointPoints_.resize(nPoints, -1);
3240  this->pointSrcPoints_.resize(nPoints, -1);
3241  this->pointTgtPoints_.resize(nPoints, -1);
3242  this->pointSrcEdges_.resize(nPoints, -1);
3243  this->pointTgtEdges_.resize(nPoints, -1);
3244  this->pointSrcFaces_.resize(nPoints, -1);
3245  this->pointTgtFaces_.resize(nPoints, -1);
3246  forAll(this->srcPatch_.localPoints(), srcPointi)
3247  {
3248  const label pointi = srcPointi;
3249 
3250  srcPoints_[pointi] = this->srcPatch_.localPoints()[srcPointi];
3251  tgtPoints_[pointi] = this->srcPatch_.localPoints()[srcPointi];
3252  srcPointNormals_[pointi] = srcPointNormals[srcPointi];
3253  pointPoints_[pointi] = pointi;
3254  this->srcPointPoints_[srcPointi] = pointi;
3255  this->pointSrcPoints_[pointi] = srcPointi;
3256  }
3257  forAll(this->tgtPatch_.localPoints(), tgtPointi)
3258  {
3259  const label pointi = this->srcPatch_.nPoints() + tgtPointi;
3260 
3261  srcPoints_[pointi] = this->tgtPatch_.localPoints()[tgtPointi];
3262  tgtPoints_[pointi] = this->tgtPatch_.localPoints()[tgtPointi];
3263  pointPoints_[pointi] = pointi;
3264  this->tgtPointPoints_[tgtPointi] = pointi;
3265  this->pointTgtPoints_[pointi] = tgtPointi;
3266  }
3267 
3268  // Edges
3269  const label nEdges = this->srcPatch_.nEdges() + this->tgtPatch_.nEdges();
3270  edgeTris_.resize(nEdges, labelPair(-1, -1));
3271  intersectEdgeFaces_.resize(nEdges, labelPair(-1, -1));
3272  forAll(this->srcPatch_.faceEdges(), srcFacei)
3273  {
3274  const label trii = srcFacei;
3275 
3276  forAll(this->srcPatch_.faceEdges()[srcFacei], srcFaceEdgei)
3277  {
3278  const label srcEdgei =
3279  this->srcPatch_.faceEdges()[srcFacei][srcFaceEdgei];
3280  const label edgei = srcEdgei;
3281 
3282  const edge& e = this->srcPatch_.edges()[srcEdgei];
3283  const edge fe =
3284  this->srcPatch_.localFaces()[srcFacei].faceEdge(srcFaceEdgei);
3285 
3286  edgeTris_[edgei][edge::compare(e, fe) < 0] = trii;
3287  }
3288  }
3289  forAll(this->srcPatch_.edges(), srcEdgei)
3290  {
3291  const edge& e = this->srcPatch_.edges()[srcEdgei];
3292 
3293  this->srcEdgePoints_[srcEdgei].append(this->srcPointPoints_[e[0]]);
3294  this->srcEdgePoints_[srcEdgei].append(this->srcPointPoints_[e[1]]);
3295  }
3296  forAll(this->tgtPatch_.faceEdges(), tgtFacei)
3297  {
3298  const label trii = this->srcPatch_.size() + tgtFacei;
3299 
3300  forAll(this->tgtPatch_.faceEdges()[tgtFacei], tgtFaceEdgei)
3301  {
3302  const label tgtEdgei =
3303  this->tgtPatch_.faceEdges()[tgtFacei][tgtFaceEdgei];
3304  const label edgei = this->srcPatch_.nEdges() + tgtEdgei;
3305 
3306  const edge& e = this->tgtPatch_.edges()[tgtEdgei];
3307  const edge fe =
3308  this->tgtPatch_.localFaces()[tgtFacei].faceEdge(tgtFaceEdgei);
3309 
3310  edgeTris_[edgei][edge::compare(e, fe) < 0] = trii;
3311  }
3312  }
3313  forAll(this->tgtPatch_.edges(), tgtEdgei)
3314  {
3315  const edge& e = this->tgtPatch_.edges()[tgtEdgei];
3316 
3317  this->tgtEdgePoints_[tgtEdgei].append(this->tgtPointPoints_[e[0]]);
3318  this->tgtEdgePoints_[tgtEdgei].append(this->tgtPointPoints_[e[1]]);
3319  }
3320 
3321  // Tris
3322  const label nTris = this->srcPatch_.size() + this->tgtPatch_.size();
3323  triPoints_.resize(nTris);
3324  triEdges_.resize(nTris);
3325  triSrcFace_.resize(nTris, -1);
3326  triTgtFace_.resize(nTris, -1);
3327  forAll(this->srcPatch_.localFaces(), srcFacei)
3328  {
3329  const label trii = srcFacei;
3330 
3331  forAll(this->srcPatch_.localFaces()[srcFacei], i)
3332  {
3333  triPoints_[trii][i] = this->srcPatch_.localFaces()[srcFacei][i];
3334  triEdges_[trii][i] = this->srcPatch_.faceEdges()[srcFacei][i];
3335  }
3336  triSrcFace_[trii] = srcFacei;
3337  srcFaceTris_[srcFacei].resize(1, trii);
3338  }
3339  forAll(this->tgtPatch_.localFaces(), tgtFacei)
3340  {
3341  const label trii = this->srcPatch_.size() + tgtFacei;
3342 
3343  forAll(this->tgtPatch_.localFaces()[tgtFacei], i)
3344  {
3345  triPoints_[trii][i] =
3346  this->srcPatch_.nPoints()
3347  + this->tgtPatch_.localFaces()[tgtFacei][i];
3348  triEdges_[trii][i] =
3349  this->srcPatch_.nEdges()
3350  + this->tgtPatch_.faceEdges()[tgtFacei][i];
3351  }
3352  triTgtFace_[trii] = tgtFacei;
3353  tgtFaceTris_[tgtFacei].resize(1, trii);
3354  }
3355 
3356  // Removal
3357  removedEdges_.clear();
3358  removedTris_.clear();
3359 
3360  // Front propagation
3361  frontEdgeEdges_.clear();
3362  edgeFrontEdges_ = DynamicList<label>(edgeTris_.size(), -1);
3363 
3364  // Insertion
3365  candidateTriTris_.clear();
3366  triCandidateTris_ = DynamicList<label>(triPoints_.size(), -1);
3367 
3368  // Marking
3369  markedTriTris_.clear();
3370  triMarkedTris_ = DynamicList<label>(triPoints_.size(), -1);
3371 
3372  checkPatchFaces(true);
3373  checkPatchFaces(false);
3374 }
3375 
3376 
3377 template<class SrcPatchType, class TgtPatchType>
3378 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::clean()
3379 {
3380  // Resolve point-points
3381  inplaceRenumber(pointPoints_, triPoints_);
3382  inplaceRenumber(pointPoints_, this->srcPointPoints_);
3383  inplaceRenumber(pointPoints_, this->tgtPointPoints_);
3384  inplaceRenumber(pointPoints_, this->srcEdgePoints_);
3385  inplaceRenumber(pointPoints_, this->tgtEdgePoints_);
3386 
3387  // Remove points by shuffling up
3388  labelList oldPointNewPoints(pointPoints_.size(), -1);
3389  label pointi = 0;
3390  forAll(pointPoints_, pointj)
3391  {
3392  if (pointPoints_[pointj] == pointj)
3393  {
3394  oldPointNewPoints[pointj] = pointi;
3395  srcPoints_[pointi] = srcPoints_[pointj];
3396  srcPointNormals_[pointi] = srcPointNormals_[pointj];
3397  tgtPoints_[pointi] = tgtPoints_[pointj];
3398  pointPoints_[pointi] = pointi;
3399  this->pointSrcPoints_[pointi] = this->pointSrcPoints_[pointj];
3400  this->pointTgtPoints_[pointi] = this->pointTgtPoints_[pointj];
3401  this->pointSrcEdges_[pointi] = this->pointSrcEdges_[pointj];
3402  this->pointTgtEdges_[pointi] = this->pointTgtEdges_[pointj];
3403  this->pointSrcFaces_[pointi] = this->pointSrcFaces_[pointj];
3404  this->pointTgtFaces_[pointi] = this->pointTgtFaces_[pointj];
3405  ++ pointi;
3406  }
3407  }
3408  srcPoints_.resize(pointi);
3409  srcPointNormals_.resize(pointi);
3410  tgtPoints_.resize(pointi);
3411  pointPoints_.resize(pointi);
3412  this->pointSrcPoints_.resize(pointi);
3413  this->pointTgtPoints_.resize(pointi);
3414  this->pointSrcEdges_.resize(pointi);
3415  this->pointTgtEdges_.resize(pointi);
3416  this->pointSrcFaces_.resize(pointi);
3417  this->pointTgtFaces_.resize(pointi);
3418 
3419  // Remove edges by shuffling up
3420  labelList oldEdgeNewEdges(edgeTris_.size(), -1);
3421  label edgei = 0;
3422  forAll(edgeTris_, edgej)
3423  {
3424  if
3425  (
3426  edgeTris_[edgej] != labelPair(-1, -1)
3427  || intersectEdgeFaces_[edgej] != labelPair(-1, -1)
3428  )
3429  {
3430  oldEdgeNewEdges[edgej] = edgei;
3431  edgeTris_[edgei] = edgeTris_[edgej];
3432  intersectEdgeFaces_[edgei] = intersectEdgeFaces_[edgej];
3433  edgeFrontEdges_[edgei] = edgeFrontEdges_[edgej];
3434  ++ edgei;
3435  }
3436  }
3437  edgeTris_.resize(edgei);
3438  intersectEdgeFaces_.resize(edgei);
3439  edgeFrontEdges_.resize(edgei);
3440 
3441  // Remove tris by shuffling up
3442  labelList oldTriNewTris(triPoints_.size(), -1);
3443  label trii = 0;
3444  forAll(triPoints_, trij)
3445  {
3446  if (triPoints_[trij] != FixedList<label, 3>({-1, -1, -1}))
3447  {
3448  oldTriNewTris[trij] = trii;
3449  triPoints_[trii] = triPoints_[trij];
3450  triEdges_[trii] = triEdges_[trij];
3451  triSrcFace_[trii] = triSrcFace_[trij];
3452  triTgtFace_[trii] = triTgtFace_[trij];
3453  ++ trii;
3454  }
3455  }
3456  triPoints_.resize(trii);
3457  triEdges_.resize(trii);
3458  triSrcFace_.resize(trii);
3459  triTgtFace_.resize(trii);
3460 
3461  // Map
3462  inplaceRenumber(oldTriNewTris, edgeTris_);
3463  inplaceRenumber(oldPointNewPoints, triPoints_);
3464  inplaceRenumber(oldEdgeNewEdges, triEdges_);
3465  inplaceRenumber(oldTriNewTris, srcFaceTris_);
3466  inplaceRenumber(oldTriNewTris, tgtFaceTris_);
3467  inplaceRenumber(oldEdgeNewEdges, frontEdgeEdges_);
3468  inplaceRenumber(oldPointNewPoints, this->faces_);
3469  inplaceRenumber(oldEdgeNewEdges, faceEdges_);
3470  inplaceRenumber(oldPointNewPoints, this->srcPointPoints_);
3471  inplaceRenumber(oldPointNewPoints, this->tgtPointPoints_);
3472  inplaceRenumber(oldPointNewPoints, this->srcEdgePoints_);
3473  inplaceRenumber(oldPointNewPoints, this->tgtEdgePoints_);
3474 
3475  // Removal
3476  removedEdges_.clear();
3477  removedTris_.clear();
3478 
3479  // Insertion
3480  candidateTriTris_.clear();
3481  triCandidateTris_ = DynamicList<label>(triPoints_.size(), -1);
3482 
3483  // Marking
3484  markedTriTris_.clear();
3485  triMarkedTris_ = DynamicList<label>(triPoints_.size(), -1);
3486 
3487  checkPatchFaces(true);
3488  checkPatchFaces(false);
3489 }
3490 
3491 
3492 template<class SrcPatchType, class TgtPatchType>
3493 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::finalise()
3494 {
3495  clean();
3496 
3497  // Add the remaining triangles as uncoupled faces ...
3498 
3499  const label nFaces = this->faces_.size();
3500 
3501  this->faces_.resize(nFaces + triPoints_.size());
3502  faceEdges_.resize(nFaces + triPoints_.size());
3503  this->faceSrcFaces_.resize(nFaces + triPoints_.size());
3504  this->faceTgtFaces_.resize(nFaces + triPoints_.size());
3505 
3506  forAll(triPoints_, trii)
3507  {
3508  const label facei = nFaces + trii;
3509 
3510  if (triSrcFace_[trii] != -1)
3511  {
3512  this->faces_[facei] = triPoints(trii);
3513  faceEdges_[facei] = labelList(triEdges_[trii]);
3514  this->srcFaceFaces_[triSrcFace_[trii]].append(facei);
3515  this->faceSrcFaces_[facei] = triSrcFace_[trii];
3516  this->faceTgtFaces_[facei] = -1;
3517  }
3518  else
3519  {
3520  this->faces_[facei] = triPoints(trii).reverseFace();
3521  faceEdges_[facei] = labelList(reverseList(triEdges_[trii]));
3522  this->tgtFaceFaces_[triTgtFace_[trii]].append(facei);
3523  this->faceSrcFaces_[facei] = -1;
3524  this->faceTgtFaces_[facei] = triTgtFace_[trii];
3525  }
3526  }
3527 
3528  nonIntersectEdgeFaces_.resize(edgeTris_.size());
3529  forAll(edgeTris_, edgei)
3530  {
3531  forAll(edgeTris_[edgei], edgeTrii)
3532  {
3533  nonIntersectEdgeFaces_[edgei][edgeTrii] =
3534  edgeTris_[edgei][edgeTrii] + nFaces;
3535  }
3536  }
3537 
3538  checkPatchFaces(true);
3539  checkPatchFaces(false);
3540 }
3541 
3542 
3543 template<class SrcPatchType, class TgtPatchType>
3544 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::unFinalise()
3545 {
3546  clean();
3547 
3548  // Remove all uncoupled faces and store them as triangles ...
3549 
3550  label nFaces = 0;
3551  while
3552  (
3553  nFaces < this->faces_.size()
3554  && this->faceSrcFaces_[nFaces] != -1
3555  && this->faceTgtFaces_[nFaces] != -1
3556  )
3557  {
3558  ++ nFaces;
3559  }
3560 
3561  this->faces_.resize(nFaces);
3562  faceEdges_.resize(nFaces);
3563  this->faceSrcFaces_.resize(nFaces);
3564  this->faceTgtFaces_.resize(nFaces);
3565 
3566  forAll(triPoints_, trii)
3567  {
3568  DynamicList<label>& faceis =
3569  triSrcFace_[trii] != -1
3570  ? this->srcFaceFaces_[triSrcFace_[trii]]
3571  : this->tgtFaceFaces_[triTgtFace_[trii]];
3572 
3573  label n = 0;
3574  while (n < faceis.size() && faceis[n] < nFaces)
3575  {
3576  ++ n;
3577  }
3578 
3579  faceis.resize(n);
3580  }
3581 
3582  nonIntersectEdgeFaces_.clear();
3583 
3584  checkPatchFaces(true);
3585  checkPatchFaces(false);
3586 }
3587 
3588 
3589 template<class SrcPatchType, class TgtPatchType>
3590 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::write()
3591 {
3592  if (this->debug > 2)
3593  {
3594  finalise();
3595 
3596  // Use base class to write the patch
3597  this->report(name(writei_));
3598 
3599  unFinalise();
3600 
3601  // Write the edge front
3602  const fileName frontFileName =
3603  type() + "_front_" + name(writei_) + ".vtk";
3604  Info<< indent << "Writing front to " << frontFileName << endl;
3605 
3606  DynamicList<point> writePoints(frontEdgeEdges_.size()*2);
3607  DynamicList<labelPair> writeLines(frontEdgeEdges_.size()*2);
3608  forAll(frontEdgeEdges_, frontEdgei)
3609  {
3610  const label edgei = frontEdgeEdges_[frontEdgei];
3611 
3612  if (edgei == -1) continue;
3613 
3614  const edge e = edgePoints(edgei);
3615  writePoints.append(tgtPoints_[e.start()]);
3616  writePoints.append(tgtPoints_[e.end()]);
3617 
3618  const label i = writePoints.size() - 2;
3619  writeLines.append(labelPair(i, i + 1));
3620  }
3621 
3622  vtkWritePolyData::write
3623  (
3624  frontFileName,
3625  "front",
3626  false,
3627  writePoints,
3628  labelList(),
3629  writeLines,
3630  faceList()
3631  );
3632 
3633  writei_ ++;
3634  }
3635 }
3636 
3637 
3638 template<class SrcPatchType, class TgtPatchType>
3639 void Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::writePatchFace
3640 (
3641  const label patchFacei,
3642  const bool isSrc
3643 ) const
3644 {
3645  OFstream os
3646  (
3647  word(isSrc ? "src" : "tgt") + "Face_" + name(patchFacei) + ".obj"
3648  );
3649 
3650  OFstream tos
3651  (
3652  word(isSrc ? "src" : "tgt") + "FaceTris_" + name(patchFacei) + ".obj"
3653  );
3654 
3655  Info<< indent << "Writing patch face to " << os.name()
3656  << " and patch face triangulation to " << tos.name()
3657  << incrIndent << endl;
3658 
3659  forAll(patchFacePatchPoints(patchFacei, isSrc), patchFacePatchPointi)
3660  {
3661  const label patchPointi =
3662  patchFacePatchPoints(patchFacei, isSrc)[patchFacePatchPointi];
3663  const point& p =
3664  isSrc
3665  ? this->srcPatch_.localPoints()[patchPointi]
3666  : this->tgtPatch_.localPoints()[patchPointi];
3667  os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
3668  }
3669  os << "f 1 2 3" << nl;
3670 
3671  const labelList& patchFaceTris =
3672  isSrc ? srcFaceTris_[patchFacei] : tgtFaceTris_[patchFacei];
3673 
3674  forAll(patchFaceTris, patchFaceTrii)
3675  {
3676  const label trii = patchFaceTris[patchFaceTrii];
3677 
3678  Info<< indent << "tri #" << trii << " points=" << triPoints(trii)
3679  << " edges=" << triEdges_[trii] << endl;
3680 
3681  forAll(triPoints_[trii], triPointi)
3682  {
3683  const label pointi = triPoint(trii, triPointi);
3684  const point& p = isSrc ? srcPoints_[pointi] : tgtPoints_[pointi];
3685  tos << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
3686  }
3687  tos << "f";
3688  forAll(triPoints_[trii], triPointi)
3689  {
3690  tos << " " << 1 + 3*patchFaceTrii + triPointi;
3691  }
3692  tos << nl;
3693  }
3694 
3695  Info<< decrIndent;
3696 }
3697 
3698 
3699 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
3700 
3701 template<class SrcPatchType, class TgtPatchType>
3702 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::TriPatchIntersection
3703 (
3704  const SrcPatchType& srcPatch,
3705  const TgtPatchType& tgtPatch,
3706  const scalar snapTol
3707 )
3708 :
3709  TriPatchIntersection(srcPatch, srcPatch.pointNormals(), tgtPatch, snapTol)
3710 {}
3711 
3712 
3713 template<class SrcPatchType, class TgtPatchType>
3714 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::TriPatchIntersection
3715 (
3716  const SrcPatchType& srcPatch,
3717  const vectorField& srcPointNormals,
3718  const TgtPatchType& tgtPatch,
3719  const scalar snapTol
3720 )
3721 :
3722  PatchIntersection<SrcPatchType, TgtPatchType>(srcPatch, tgtPatch),
3723 
3724  srcPoints_(),
3725  srcPointNormals_(),
3726  tgtPoints_(this->points_),
3727  pointPoints_(),
3728 
3729  edgeTris_(),
3730  intersectEdgeFaces_(),
3731  nonIntersectEdgeFaces_(),
3732 
3733  triPoints_(),
3734  triEdges_(),
3735  triSrcFace_(),
3736  triTgtFace_(),
3737  srcFaceTris_(srcPatch.size()),
3738  tgtFaceTris_(tgtPatch.size()),
3739 
3740  faceEdges_(),
3741 
3742  removedTris_(),
3743  removedEdges_(),
3744 
3745  frontEdgeEdges_(),
3746  edgeFrontEdges_(),
3747 
3748  candidateTriTris_(),
3749  triCandidateTris_(),
3750 
3751  markedTriTris_(),
3752  triMarkedTris_(),
3753 
3754  polygonTriangulate_(),
3755 
3756  star_(),
3757 
3758  writei_(0)
3759 {
3760  cpuTime time;
3761 
3762  Info<< indent << type() << ": Intersecting "
3763  << this->srcPatch_.size() << " source tri faces and "
3764  << this->tgtPatch_.size() << " target tri faces" << incrIndent << endl;
3765 
3766  if (this->debug)
3767  {
3768  Info<< indent << "Writing tri patches" << incrIndent << endl;
3769  const fileName srcFileName = type() + "_srcPatch.vtk";
3770  Info<< indent << "Writing patch to " << srcFileName << endl;
3771  vtkWritePolyData::write
3772  (
3773  srcFileName,
3774  "source",
3775  false,
3776  this->srcPatch_.localPoints(),
3777  labelList(),
3778  labelListList(),
3779  this->srcPatch_.localFaces(),
3780  "normals",
3781  true,
3782  srcPointNormals
3783  );
3784  const fileName tgtFileName = type() + "_tgtPatch.vtk";
3785  Info<< indent << "Writing patch to " << tgtFileName << endl;
3786  vtkWritePolyData::write
3787  (
3788  tgtFileName,
3789  "target",
3790  false,
3791  this->tgtPatch_.localPoints(),
3792  labelList(),
3793  labelListList(),
3794  this->tgtPatch_.localFaces()
3795  );
3796  Info<< decrIndent;
3797  }
3798 
3799  // Create bound spheres for patch faces for proximity testing
3800  List<Tuple2<point, scalar>> srcFaceSpheres(this->srcPatch_.size());
3801  List<Tuple2<point, scalar>> tgtFaceSpheres(this->tgtPatch_.size());
3802  forAll(this->srcPatch_, srcFacei)
3803  {
3804  const triFace& srcFace = this->srcPatch_.localFaces()[srcFacei];
3805  srcFaceSpheres[srcFacei] =
3806  trivialBoundSphere
3807  (
3808  this->srcPatch_.localPoints(),
3809  {srcFace[0], srcFace[1], srcFace[2], -1},
3810  3
3811  );
3812  srcFaceSpheres[srcFacei].second() *= 1 + snapTol;
3813  }
3814  forAll(this->tgtPatch_, tgtFacei)
3815  {
3816  const triFace& tgtFace = this->tgtPatch_.localFaces()[tgtFacei];
3817  tgtFaceSpheres[tgtFacei] =
3818  trivialBoundSphere
3819  (
3820  this->tgtPatch_.localPoints(),
3821  {tgtFace[0], tgtFace[1], tgtFace[2], -1},
3822  3
3823  );
3824  tgtFaceSpheres[tgtFacei].second() *= 1 + snapTol;
3825  }
3826 
3827  // Construct table to store what faces have been snapped
3828  HashSet<labelPair, labelPair::Hash<>> srcFaceTgtFaceSnaps
3829  (
3830  12*(this->srcPatch_.size() + this->tgtPatch_.size())
3831  );
3832 
3833  // Count the number of successes and failures
3834  label nIntersections = 0, nIntersectionFailures = 0;
3835 
3836  // Function for intersecting two patch faces
3837  auto intersect = [&](const label srcFacei, const label tgtFacei)
3838  {
3839  // Single snapping stage
3840  /*
3841  snapPatchFaceTris(srcFacei, tgtFacei, snapTol);
3842  */
3843 
3844  // Propagate out and snap everything in advance of the intersections
3845  const List<triFace>& srcLocalFaces = this->srcPatch_.localFaces();
3846  const List<triFace>& tgtLocalFaces = this->tgtPatch_.localFaces();
3847  const labelListList& srcPointFaces = this->srcPatch_.pointFaces();
3848  const labelListList& tgtPointFaces = this->tgtPatch_.pointFaces();
3849  forAll(srcLocalFaces[srcFacei], srcFacePointi)
3850  {
3851  const label srcPointi =
3852  srcLocalFaces[srcFacei][srcFacePointi];
3853 
3854  forAll(srcPointFaces[srcPointi], srcPointFacei)
3855  {
3856  const label srcFacej =
3857  srcPointFaces[srcPointi][srcPointFacei];
3858 
3859  forAll(tgtLocalFaces[tgtFacei], tgtFacePointi)
3860  {
3861  const label tgtPointi =
3862  tgtLocalFaces[tgtFacei][tgtFacePointi];
3863 
3864  forAll(tgtPointFaces[tgtPointi], tgtPointFacei)
3865  {
3866  const label tgtFacej =
3867  tgtPointFaces[tgtPointi][tgtPointFacei];
3868 
3869  const point& srcC = srcFaceSpheres[srcFacej].first();
3870  const scalar srcR = srcFaceSpheres[srcFacej].second();
3871  const point& tgtC = tgtFaceSpheres[tgtFacej].first();
3872  const scalar tgtR = tgtFaceSpheres[tgtFacej].second();
3873 
3874  if
3875  (
3876  magSqr(srcC - tgtC) < sqr(srcR + tgtR)
3877  && !srcFaceTgtFaceSnaps.found({srcFacej, tgtFacej})
3878  )
3879  {
3880  srcFaceTgtFaceSnaps.insert({srcFacej, tgtFacej});
3881  snapPatchFaceTris(srcFacej, tgtFacej, snapTol);
3882  }
3883  }
3884  }
3885  }
3886  }
3887 
3888  intersectPatchFaceTris(srcFacei, tgtFacei);
3889 
3890  write();
3891 
3892  const bool conformFailure = !conformPatchFaceTris(srcFacei, tgtFacei);
3893 
3894  const bool combineFailure = !combinePatchFaceTris(srcFacei, tgtFacei);
3895 
3896  nIntersections ++;
3897  nIntersectionFailures += conformFailure || combineFailure;
3898 
3899  write();
3900  };
3901 
3902  // Target search tree. Used to find an initial source and target face to
3903  // intersect. Once the first intersection has been done the rest follow via
3904  // a front-propagation algorithm, without reference to this tree.
3905  indexedOctree<treeDataPrimitivePatch<TgtPatchType>> tgtTree
3906  (
3907  treeDataPrimitivePatch<TgtPatchType>
3908  (
3909  false,
3910  this->tgtPatch_,
3911  indexedOctree<treeDataPrimitivePatch<TgtPatchType>>::perturbTol()
3912  ),
3913  treeBoundBox(this->tgtPatch_.localPoints()).extend(1e-4),
3914  8,
3915  10,
3916  3
3917  );
3918 
3919  // Find intersection operation, excluding a hash set of previously obtained
3920  // intersections. Used to get multiple target triangles intersected by a
3921  // single source point and point normal (i.e., two target triangles if the
3922  // source point intersects a target edge).
3923  class findIntersectExcludingOp
3924  {
3925  public:
3926 
3927  const indexedOctree<treeDataPrimitivePatch<TgtPatchType>>& tree_;
3928 
3929  const labelHashSet& excludingIndices_;
3930 
3931  public:
3932 
3933  //- Construct from components
3934  findIntersectExcludingOp
3935  (
3936  const indexedOctree<treeDataPrimitivePatch<TgtPatchType>>& tree,
3937  const labelHashSet& hitIndices
3938  )
3939  :
3940  tree_(tree),
3941  excludingIndices_(hitIndices)
3942  {}
3943 
3944  //- Calculate intersection of triangle with ray. Sets result
3945  // accordingly
3946  bool operator()
3947  (
3948  const label index,
3949  const point& start,
3950  const point& end,
3951  point& intersectionPoint
3952  ) const
3953  {
3954  if (excludingIndices_.found(index))
3955  {
3956  return false;
3957  }
3958  else
3959  {
3960  typename
3961  treeDataPrimitivePatch<TgtPatchType>::findIntersectOp
3962  iop(tree_);
3963 
3964  return
3965  iop
3966  (
3967  index,
3968  start,
3969  end,
3970  intersectionPoint
3971  );
3972  }
3973  }
3974  };
3975 
3976  // Populate local data from the source and target patches
3977  initialise(srcPointNormals);
3978  write();
3979 
3980  // Loop the source points, looking for ones that have not been intersected
3981  forAll(this->points_, pointi)
3982  {
3983  // Get the next source patch point
3984  const label srcPointi = this->pointSrcPoints_[pointi];
3985 
3986  // Continue if this point is already intersected
3987  if
3988  (
3989  srcPointi == -1
3990  || this->pointTgtPoints_[pointi] != -1
3991  || this->pointTgtEdges_[pointi] != -1
3992  || this->pointTgtFaces_[pointi] != -1
3993  ) continue;
3994 
3995  // Get the projection geometry for this source point
3996  const point& srcP = srcPoints_[pointi];
3997  const vector& srcN = srcPointNormals_[pointi];
3998  scalar srcL = 0;
3999  forAll(this->srcPatch_.pointEdges()[srcPointi], srcPointEdgei)
4000  {
4001  const label srcEdgei =
4002  this->srcPatch_.pointEdges()[srcPointi][srcPointEdgei];
4003  srcL +=
4004  this->srcPatch_.edges()[srcEdgei].mag
4005  (
4006  this->srcPatch_.localPoints()
4007  );
4008  }
4009  srcL /= this->srcPatch_.pointEdges()[srcPointi].size();
4010 
4011  // Find all the target faces that this source point projects to
4012  labelHashSet tgtFaceis;
4013  while (true)
4014  {
4015  pointIndexHit hit =
4016  tgtTree.findLine
4017  (
4018  srcP - srcL*srcN,
4019  srcP + srcL*srcN,
4020  findIntersectExcludingOp(tgtTree, tgtFaceis)
4021  );
4022 
4023  if (!hit.hit()) break;
4024 
4025  tgtFaceis.insert(hit.index());
4026  }
4027 
4028  // Continue if this point does not project to the opposite patch
4029  if (tgtFaceis.empty()) continue;
4030 
4031  // Loop all the potential source/target face intersections until an
4032  // edge front is generated
4033  forAll(this->srcPatch_.pointFaces()[srcPointi], srcPointFacei)
4034  {
4035  const label srcFacei =
4036  this->srcPatch_.pointFaces()[srcPointi][srcPointFacei];
4037 
4038  forAllConstIter(labelHashSet, tgtFaceis, tgtFaceiIter)
4039  {
4040  intersect(srcFacei, tgtFaceiIter.key());
4041 
4042  if (frontEdgeEdges_.size()) break;
4043  }
4044 
4045  if (frontEdgeEdges_.size()) break;
4046  }
4047 
4048  // Propagate until the edge front is empty
4049  while (frontEdgeEdges_.size())
4050  {
4051  const label edgei = frontEdgeEdges_.remove();
4052 
4053  if (edgei == -1) continue;
4054 
4055  edgeFrontEdges_[edgei] = -1;
4056 
4057  label srcTrii = edgeTris_[edgei][0];
4058  label tgtTrii = edgeTris_[edgei][1];
4059 
4060  if (triSrcFace_[srcTrii] == -1)
4061  {
4062  Swap(srcTrii, tgtTrii);
4063  }
4064 
4065  intersect(triSrcFace_[srcTrii], triTgtFace_[tgtTrii]);
4066  }
4067  }
4068 
4069  // Populate data in the base class which was not generated as part of the
4070  // intersection process
4071  finalise();
4072 
4073  this->report();
4074 
4075  // Warn about any failures
4076  if (nIntersectionFailures)
4077  {
4078  Info<< indent << "*** Topology could not be generated in "
4079  << nIntersectionFailures << "/" << nIntersections << " cases"
4080  << endl << indent << " The intersection may be incomplete"
4081  << endl;
4082  }
4083 
4084  Info<< indent << this->faces_.size() << " faces generated in "
4085  << time.cpuTimeIncrement() << 's' << endl;
4086 
4087  Info<< decrIndent;
4088 }
4089 
4090 
4091 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
4092 
4093 template<class SrcPatchType, class TgtPatchType>
4094 Foam::TriPatchIntersection<SrcPatchType, TgtPatchType>::
4095 ~TriPatchIntersection()
4096 {}
4097 
4098 
4099 // ************************************************************************* //
A
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
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:434
forAllConstIter
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
boundSphere.H
Functions for constructing bounding spheres of lists of points.
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::empty
bool empty() const
Return true if the hash table is empty.
Definition: HashTableI.H:72
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:294
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::PointIndexHit
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:54
Foam::PointIndexHit::index
label index() const
Return index.
Definition: PointIndexHit.H:111
Foam::PointIndexHit::hit
bool hit() const
Is there a hit.
Definition: PointIndexHit.H:105
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:4095
Foam::TriPatchIntersection::TriPatchIntersection
TriPatchIntersection(const SrcPatchType &srcPatch, const TgtPatchType &tgtPatch, const scalar snapTol)
Construct from a source and a target patch.
Definition: TriPatchIntersection.C:3703
Foam::TriPatchIntersection::faceEdges
const DynamicList< labelList > & faceEdges() const
...
Definition: TriPatchIntersection.H:450
Foam::UList::first
T & first()
Return the first element of the list.
Definition: UListI.H:114
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
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:251
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
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:1803
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:1555
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:1866
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:241
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::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:49
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
Foam::name
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
Foam::reverseList
ListType reverseList(const ListType &list)
Reverse a list. First element becomes last element etc.
Definition: ListOpsTemplates.C:760
Foam::inplaceReverseList
void inplaceReverseList(ListType &list)
Inplace reversal of a list using Swap.
Definition: ListOpsTemplates.C:777
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:817
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:234
Foam::point
vector point
Point is a vector.
Definition: point.H:41
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::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
Definition: ListOpsTemplates.C:73
Foam::max
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:104
Foam::trivialBoundSphere
Tuple2< point, scalar > trivialBoundSphere(const PointField &ps, const FixedList< label, 4 > &pis, const label nPs)
Compute a bounding sphere of four points or less.
Foam::triPointRef
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
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::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:227
Foam::nl
static const char nl
Definition: Ostream.H:266
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:55
Foam::magSqr
dimensioned< scalar > magSqr(const dimensioned< Type > &)
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....