patchToPatch.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) 2011-2023 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 "patchToPatch.H"
27 #include "patchToPatchTools.H"
28 #include "cpuTime.H"
29 #include "distributionMap.H"
30 #include "globalIndex.H"
31 #include "indexedOctree.H"
32 #include "treeDataPrimitivePatch.H"
33 #include "vtkWritePolyData.H"
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
40  {
41  labelList f(p.size());
42  forAll(p, i)
43  {
44  f[i] = p[i].first();
45  }
46  return f;
47  }
48 
50  {
51  labelList s(p.size());
52  forAll(p, i)
53  {
54  s[i] = p[i].second();
55  }
56  return s;
57  }
58 
59  template<class ListType>
61  (
62  const ListType& l,
63  typename ListType::const_reference t,
64  const label start=0
65  )
66  {
67  for (label i = start; i < l.size(); i++)
68  {
69  if (l[i] != t)
70  {
71  return i;
72  }
73  }
74 
75  return -1;
76  }
77 
79  {
80  return treeBoundBox(min(a.min(), b.min()), max(a.max(), b.max()));
81  }
82 }
83 
84 
85 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
86 
87 namespace Foam
88 {
91 }
92 
93 
94 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
95 
97 (
98  const primitiveOldTimePatch& srcPatch,
99  const vectorField& srcPointNormals,
100  const vectorField& srcPointNormals0,
101  const label srcFacei
102 ) const
103 {
104  const face& srcFace = srcPatch.localFaces()[srcFacei];
105  const pointField& srcPoints = srcPatch.localPoints();
106  treeBoundBox box = srcBox(srcFace, srcPoints, srcPointNormals);
107 
108  if (srcPatch.has0())
109  {
110  const pointField& srcPoints0 = srcPatch.localPoints();
111  box = combine(box, srcBox(srcFace, srcPoints0, srcPointNormals0));
112  }
113 
114  return box;
115 }
116 
117 
119 (
120  const primitiveOldTimePatch& srcPatch,
121  const vectorField& srcPointNormals,
122  const vectorField& srcPointNormals0
123 ) const
124 {
126 
127  forAll(srcPatch, srcFacei)
128  {
129  box =
130  combine
131  (
132  box,
133  srcBox(srcPatch, srcPointNormals, srcPointNormals0, srcFacei)
134  );
135  }
136 
137  box.inflate(rootSmall);
138 
139  return box;
140 }
141 
142 
144 (
145  const primitiveOldTimePatch& tgtPatch,
146  const label tgtFacei
147 ) const
148 {
149  const face& tgtFace = tgtPatch.localFaces()[tgtFacei];
150  const pointField& tgtPoints = tgtPatch.localPoints();
151  treeBoundBox box(tgtPoints, tgtFace);
152 
153  if (tgtPatch.has0())
154  {
155  const pointField& tgtPoints0 = tgtPatch.localPoints();
156  box = combine(box, treeBoundBox(tgtPoints0, tgtFace));
157  }
158 
159  return box;
160 }
161 
162 
164 (
165  const primitiveOldTimePatch& tgtPatch
166 ) const
167 {
169 
170  forAll(tgtPatch, tgtFacei)
171  {
172  box = combine(box, tgtBox(tgtPatch, tgtFacei));
173  }
174 
175  box.inflate(rootSmall);
176 
177  return box;
178 }
179 
180 
182 (
183  const primitiveOldTimePatch& srcPatch,
184  const vectorField& srcPointNormals,
185  const vectorField& srcPointNormals0,
186  const primitiveOldTimePatch& tgtPatch,
187  const label srcFacei,
188  const label tgtFacei
189 )
190 {
191  // Return if these faces have already been intersected
192  auto intersected = []
193  (
194  const List<DynamicList<label>>& faceLocalFaces,
195  const label facei,
196  const label otherFacei
197  )
198  {
199  forAll(faceLocalFaces[facei], i)
200  {
201  if (faceLocalFaces[facei][i] == otherFacei)
202  {
203  return true;
204  }
205  }
206  return false;
207  };
208  if
209  (
210  intersected(srcLocalTgtFaces_, srcFacei, tgtFacei)
211  || intersected(tgtLocalSrcFaces_, tgtFacei, srcFacei)
212  )
213  {
214  return true;
215  }
216 
217  // Try to intersect these faces
218  return
219  intersectFaces
220  (
221  srcPatch,
222  srcPointNormals,
223  srcPointNormals0,
224  tgtPatch,
225  srcFacei,
226  tgtFacei
227  );
228 }
229 
230 
232 (
233  const primitiveOldTimePatch& srcPatch,
234  const vectorField& srcPointNormals,
235  const vectorField& srcPointNormals0,
236  const primitiveOldTimePatch& tgtPatch,
237  const bool isSrc,
238  const DynamicList<labelPair>& queue,
239  labelList& faceComplete,
240  DynamicList<labelPair>& otherQueue,
241  const labelList& otherFaceComplete,
242  boolList& otherFaceQueued,
243  boolList& otherFaceVisited
244 )
245 {
246  const primitivePatch& otherPatch = isSrc ? tgtPatch : srcPatch;
247 
248  const faceList& otherLocalFaces = otherPatch.localFaces();
249  //const labelListList& otherFaceFaces = otherPatch.faceFaces();
250  const labelListList& otherPointFaces = otherPatch.pointFaces();
251 
252  DynamicList<label> otherQueuedFaces;
253  label nFaceComplete = 0;
254  forAll(queue, queuei)
255  {
256  const label facei = queue[queuei].first();
257  const label otherFacei = queue[queuei].second();
258 
259  // Wave out from the initial target face until all intersections fail
260  DynamicList<label> otherCurrentFaces(1, otherFacei);
261  DynamicList<label> otherVisitedFaces(1, otherFacei);
262  otherFaceVisited[otherFacei] = true;
263  label otherPerimeterReached = false;
264  while (otherCurrentFaces.size())
265  {
266  labelList otherNextFaces;
267  forAll(otherCurrentFaces, otheri)
268  {
269  const label otherFacej = otherCurrentFaces[otheri];
270 
271  if
272  (
273  findOrIntersectFaces
274  (
275  srcPatch,
276  srcPointNormals,
277  srcPointNormals0,
278  tgtPatch,
279  isSrc ? facei : otherFacej,
280  isSrc ? otherFacej : facei
281  )
282  )
283  {
284  /*
285  // Propagate to edge neighbours
286  forAll(otherFaceFaces[otherFacej], otherFaceFacej)
287  {
288  const label otherFacek =
289  otherFaceFaces[otherFacej][otherFaceFacej];
290  if (!otherFaceVisited[otherFacek])
291  {
292  otherFaceVisited[otherFacek] = true;
293  otherVisitedFaces.append(otherFacek);
294  otherNextFaces.append(otherFacek);
295  }
296  }
297  */
298 
299  // Propagate to point neighbours
300  forAll(otherPatch[otherFacej], otherFacePointj)
301  {
302  const label otherPointj =
303  otherLocalFaces[otherFacej][otherFacePointj];
304  forAll(otherPointFaces[otherPointj], otherPointFacej)
305  {
306  const label otherFacek =
307  otherPointFaces[otherPointj][otherPointFacej];
308  if (!otherFaceVisited[otherFacek])
309  {
310  otherFaceVisited[otherFacek] = true;
311  otherVisitedFaces.append(otherFacek);
312  otherNextFaces.append(otherFacek);
313  }
314  }
315  }
316 
317  // Check if this face is connected to the perimeter
318  forAll(otherPatch[otherFacej], otherFaceEdgej)
319  {
320  const label otherEdgej =
321  otherPatch.faceEdges()[otherFacej][otherFaceEdgej];
322 
323  // !!! Two face-edges is not a sufficient condition for
324  // manifoldness. The edges also need to be numbered in
325  // opposite directions in the two faces.
326  if (otherPatch.edgeFaces()[otherEdgej].size() != 2)
327  {
328  otherPerimeterReached = true;
329  }
330  }
331 
332  // If this face is not complete, then add it to the next
333  // iteration's queue
334  if
335  (
336  otherFaceComplete[otherFacej] == 0
337  && !otherFaceQueued[otherFacej]
338  )
339  {
340  otherFaceQueued[otherFacej] = true;
341  otherQueuedFaces.append(otherFacej);
342  otherQueue.append(labelPair(otherFacej, facei));
343  }
344  }
345  }
346 
347  otherCurrentFaces.transfer(otherNextFaces);
348  }
349 
350  // Reset the visited array
351  forAll(otherVisitedFaces, otherVisitedFacei)
352  {
353  otherFaceVisited[otherVisitedFaces[otherVisitedFacei]] = false;
354  }
355 
356  // Set this face as complete
357  nFaceComplete -= faceComplete[facei] == 2;
358  faceComplete[facei] =
359  max(faceComplete[facei], otherPerimeterReached ? 1 : 2);
360  nFaceComplete += faceComplete[facei] == 2;
361  }
362 
363  // Reset the queued array
364  forAll(otherQueuedFaces, otherQueuedFacei)
365  {
366  otherFaceQueued[otherQueuedFaces[otherQueuedFacei]] = false;
367  }
368 
369  return nFaceComplete;
370 }
371 
372 
374 (
375  const primitiveOldTimePatch& srcPatch,
376  const vectorField& srcPointNormals,
377  const vectorField& srcPointNormals0,
378  const primitiveOldTimePatch& tgtPatch
379 )
380 {
381  if (srcPatch.empty() || tgtPatch.empty()) return;
382 
383  if (debug)
384  {
385  Info<< indent << "Writing patches" << incrIndent << endl;
386 
387  const fileName srcFileName = typeName + "_srcPatch.vtk";
388  Info<< indent << "Writing patch to " << srcFileName << endl;
390  (
391  srcFileName,
392  "source",
393  false,
394  srcPatch.localPoints(),
395  labelList(),
396  labelListList(),
397  srcPatch.localFaces()
398  );
399 
400  const fileName tgtFileName = typeName + "_tgtPatch.vtk";
401  Info<< indent << "Writing patch to " << tgtFileName << endl;
403  (
404  tgtFileName,
405  "target",
406  false,
407  tgtPatch.localPoints(),
408  labelList(),
409  labelListList(),
410  tgtPatch.localFaces()
411  );
412 
413  Info<< decrIndent;
414  }
415 
416  auto writeQueue = [&]
417  (
418  const label restarti,
419  const label iterationi,
420  const bool isSrc,
421  const DynamicList<labelPair>& queue
422  )
423  {
424  const primitivePatch& patch = isSrc ? srcPatch : tgtPatch;
425 
426  const fileName queueFileNamePart0 =
427  typeName + "_restarti=" + name(restarti) + "_iterationi=";
428  const fileName queueFileNamePart1 =
429  word(isSrc ? "src" : "tgt") + "Queue";
430 
431  Info<< incrIndent;
432 
433  if (iterationi == 0 || iterationi % 2 == !isSrc)
434  {
435  const fileName queueFileName =
436  queueFileNamePart0 + name(iterationi) + '_'
437  + queueFileNamePart1 + ".vtk";
438  Info<< indent << "Writing queue to " << queueFileName << endl;
440  (
441  queueFileName,
442  queueFileNamePart1,
443  false,
444  patch.localPoints(),
445  labelList(),
446  labelListList(),
447  IndirectList<face>(patch.localFaces(), first(queue))
448  );
449  }
450  else
451  {
452  ln
453  (
454  queueFileNamePart0 + name(iterationi - 1) + '_'
455  + queueFileNamePart1 + ".vtk",
456  queueFileNamePart0 + name(iterationi) + '_'
457  + queueFileNamePart1 + ".vtk"
458  );
459  }
460 
461  Info<< decrIndent;
462  };
463 
464  auto writeQueues = [&]
465  (
466  const label restarti,
467  const label iterationi,
468  const DynamicList<labelPair>& srcQueue,
469  const DynamicList<labelPair>& tgtQueue
470  )
471  {
472  if (debug)
473  {
474  Info<< indent << "Iteration #" << iterationi
475  << ": processing " << srcQueue.size() << '/'
476  << srcPatch.size() << " source and "
477  << tgtQueue.size() << '/' << tgtPatch.size()
478  << " target faces " << endl;
479  }
480 
481  if (debug > 1)
482  {
483  writeQueue(restarti, iterationi, true, srcQueue);
484  writeQueue(restarti, iterationi, false, tgtQueue);
485  }
486  };
487 
488  // Build a search tree for the target patch
490  const treeBoundBox tgtTreeBox =
491  treeBoundBox(tgtPatch.localPoints()).extend(1e-4);
493  (
494  treeType
495  (
496  false,
497  tgtPatch,
499  ),
500  tgtTreeBox,
501  8,
502  10,
503  3
504  );
505 
506  // Set up complete arrays and loop until they are full. Note that the
507  // *FaceComplete lists can take three values; 0 is incomplete, 1 is
508  // complete for this iteration only, 2 is fully complete.
509  if (debug)
510  {
511  Info<< indent << "Calculating coupling" << incrIndent << endl;
512  }
513 
514  label nSrcFaceComplete = 0, nTgtFaceComplete = 0;
515  labelList srcFaceComplete(srcPatch.size(), 0);
516  labelList tgtFaceComplete(tgtPatch.size(), 0);
517  boolList srcFaceQueued(srcPatch.size(), false);
518  boolList tgtFaceQueued(tgtPatch.size(), false);
519  boolList srcFaceVisited(srcPatch.size(), false);
520  boolList tgtFaceVisited(tgtPatch.size(), false);
521  label srcFacei = 0;
522  label restarti = 0;
523  while (srcFacei < srcPatch.size() && srcFacei != -1)
524  {
525  // Consider this face only once
526  srcFaceComplete[srcFacei] = 2;
527  nSrcFaceComplete ++;
528 
529  // Find target faces that overlap this source face's bound box
530  const labelList seedTgtFaces =
531  tgtTree.findBox
532  (
533  srcBox(srcPatch, srcPointNormals, srcPointNormals0, srcFacei)
534  );
535 
536  if (!seedTgtFaces.empty())
537  {
538  if (debug)
539  {
540  Info<< indent << "Restart #" << restarti
541  << " from at source face at "
542  << srcPatch.faceCentres()[srcFacei]
543  << incrIndent << endl;
544  }
545 
546  // Initialise queues with the target faces identified
547  DynamicList<labelPair> srcQueue, tgtQueue;
548  forAll(seedTgtFaces, seedTgtFacei)
549  {
550  const label tgtFacei = seedTgtFaces[seedTgtFacei];
551  srcQueue.append(labelPair(srcFacei, tgtFacei));
552  tgtQueue.append(labelPair(tgtFacei, srcFacei));
553  }
554 
555  if (debug)
556  {
557  writeQueues(restarti, 0, srcQueue, tgtQueue);
558  }
559 
560  // Do intersections until queues are empty
561  label iterationi = 0;
562  while (true)
563  {
564  tgtQueue.clear();
565 
566  nSrcFaceComplete +=
567  intersectPatchQueue
568  (
569  srcPatch,
570  srcPointNormals,
571  srcPointNormals0,
572  tgtPatch,
573  true,
574  srcQueue,
575  srcFaceComplete,
576  tgtQueue,
577  tgtFaceComplete,
578  tgtFaceQueued,
579  tgtFaceVisited
580  );
581 
582  if (debug)
583  {
584  writeQueues(restarti, 2*iterationi + 1, srcQueue, tgtQueue);
585  }
586 
587  if (!tgtQueue.size()) break;
588 
589  srcQueue.clear();
590 
591  nTgtFaceComplete +=
592  intersectPatchQueue
593  (
594  srcPatch,
595  srcPointNormals,
596  srcPointNormals0,
597  tgtPatch,
598  false,
599  tgtQueue,
600  tgtFaceComplete,
601  srcQueue,
602  srcFaceComplete,
603  srcFaceQueued,
604  srcFaceVisited
605  );
606 
607  if (debug)
608  {
609  writeQueues(restarti, 2*iterationi + 2, srcQueue, tgtQueue);
610  }
611 
612  if (!srcQueue.size()) break;
613 
614  ++ iterationi;
615  }
616 
617  if (debug)
618  {
619  Info<< indent << "Completed " << nSrcFaceComplete << '/'
620  << srcPatch.size() << " source faces and "
621  << nTgtFaceComplete << '/' << tgtPatch.size()
622  << " target faces " << decrIndent << endl;
623  }
624  }
625 
626  // Find the next incomplete face
627  srcFacei = findNotIndex(srcFaceComplete, 2, srcFacei);
628 
629  ++ restarti;
630  }
631 
632  if (debug)
633  {
634  Info<< decrIndent;
635  }
636 }
637 
638 
640 (
641  const primitiveOldTimePatch& srcPatch,
642  const vectorField& srcPointNormals,
643  const vectorField& srcPointNormals0,
644  const primitiveOldTimePatch& tgtPatch
645 )
646 {
647  srcLocalTgtFaces_.resize(srcPatch.size());
648  forAll(srcLocalTgtFaces_, i)
649  {
650  srcLocalTgtFaces_[i].clear();
651  }
652 
653  tgtLocalSrcFaces_.resize(tgtPatch.size());
654  forAll(tgtLocalSrcFaces_, i)
655  {
656  tgtLocalSrcFaces_[i].clear();
657  }
658 }
659 
660 
662 (
663  const primitiveOldTimePatch& srcPatch,
664  const vectorField& srcPointNormals,
665  const vectorField& srcPointNormals0,
666  const primitiveOldTimePatch& tgtPatch
667 )
668 {
669  // Determine which local target faces are actually used
670  boolList oldLocalTgtFaceIsUsed(tgtPatch.size(), false);
671  forAll(srcLocalTgtFaces_, srcFacei)
672  {
673  forAll(srcLocalTgtFaces_[srcFacei], i)
674  {
675  oldLocalTgtFaceIsUsed[srcLocalTgtFaces_[srcFacei][i]] = true;
676  }
677  }
678 
679  // Trim the target map
680  labelList oldToNewLocalTgtFace, newToOldLocalTgtFace;
682  (
683  oldLocalTgtFaceIsUsed,
684  tgtMapPtr_(),
685  oldToNewLocalTgtFace,
686  newToOldLocalTgtFace
687  );
688 
689  // Renumber the addressing
690  forAll(srcLocalTgtFaces_, srcFacei)
691  {
692  forAll(srcLocalTgtFaces_[srcFacei], i)
693  {
694  srcLocalTgtFaces_[srcFacei][i] =
695  oldToNewLocalTgtFace[srcLocalTgtFaces_[srcFacei][i]];
696  }
697  }
698 
699  // Trim the local target faces
700  tgtLocalSrcFaces_ =
701  List<DynamicList<label>>(tgtLocalSrcFaces_, newToOldLocalTgtFace);
702  localTgtProcFacesPtr_() =
703  List<remote>(localTgtProcFacesPtr_(), newToOldLocalTgtFace);
704 
705  return newToOldLocalTgtFace;
706 }
707 
708 
710 (
711  const primitiveOldTimePatch& srcPatch
712 )
713 {
714  localSrcProcFacesPtr_.reset
715  (
716  new List<remote>
717  (
719  )
720  );
721 }
722 
723 
725 {
727  (
728  tgtPatch.size(),
729  tgtMapPtr_(),
730  localSrcProcFacesPtr_(),
731  tgtLocalSrcFaces_
732  );
733 }
734 
735 
737 (
738  const primitiveOldTimePatch& srcPatch,
739  const vectorField& srcPointNormals,
740  const vectorField& srcPointNormals0,
741  const primitiveOldTimePatch& tgtPatch,
742  const transformer& tgtToSrc
743 )
744 {
745  label nCouples = 0;
746 
747  forAll(srcLocalTgtFaces_, srcFacei)
748  {
749  nCouples += srcLocalTgtFaces_[srcFacei].size();
750  }
751  forAll(tgtLocalSrcFaces_, tgtFacei)
752  {
753  nCouples += tgtLocalSrcFaces_[tgtFacei].size();
754  }
755 
756  reduce(nCouples, sumOp<label>());
757 
758  return nCouples;
759 }
760 
761 
762 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
763 
765 :
766  reverse_(reverse),
767  singleProcess_(-labelMax),
768  srcLocalTgtFaces_(),
769  tgtLocalSrcFaces_(),
770  srcMapPtr_(nullptr),
771  tgtMapPtr_(nullptr),
772  localSrcProcFacesPtr_(nullptr),
773  localTgtProcFacesPtr_(nullptr)
774 {}
775 
776 
777 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
778 
780 {}
781 
782 
783 // * * * * * * * * * * * * * * * * Selector * * * * * * * * * * * * * * * * //
784 
786 (
787  const word& patchToPatchType,
788  const bool reverse
789 )
790 {
791  boolConstructorTable::iterator cstrIter =
792  boolConstructorTablePtr_->find(patchToPatchType);
793 
794  if (cstrIter == boolConstructorTablePtr_->end())
795  {
797  << "Unknown " << typeName << " type "
798  << patchToPatchType << endl << endl
799  << "Valid " << typeName << " types are : " << endl
800  << boolConstructorTablePtr_->sortedToc()
801  << exit(FatalError);
802  }
803 
804  return cstrIter()(reverse);
805 }
806 
807 
808 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
809 
811 {
812  PackedBoolList result(srcLocalTgtFaces_.size());
813  forAll(srcLocalTgtFaces_, srcFacei)
814  {
815  result[srcFacei] = !srcLocalTgtFaces_[srcFacei].empty();
816  }
817  return result;
818 }
819 
820 
822 {
823  PackedBoolList result(tgtLocalSrcFaces_.size());
824  forAll(tgtLocalSrcFaces_, tgtFacei)
825  {
826  result[tgtFacei] = !tgtLocalSrcFaces_[tgtFacei].empty();
827  }
828  return result;
829 }
830 
831 
834 {
835  return
836  isSingleProcess()
837  ? patchToPatchTools::localToRemote(srcLocalTgtFaces_)
839  (
840  srcLocalTgtFaces_,
841  localTgtProcFacesPtr_()
842  );
843 }
844 
845 
848 {
849  return
850  isSingleProcess()
851  ? patchToPatchTools::localToRemote(tgtLocalSrcFaces_)
853  (
854  tgtLocalSrcFaces_,
855  localSrcProcFacesPtr_()
856  );
857 }
858 
859 
861 (
862  const primitiveOldTimePatch& srcPatch,
863  const vectorField& srcPointNormals,
864  const vectorField& srcPointNormals0,
865  const primitiveOldTimePatch& tgtPatch,
866  const transformer& tgtToSrc
867 )
868 {
869  cpuTime time;
870 
871  // Determine numbers of faces on both sides, report, and quit if either
872  // side is empty
873  const label srcTotalSize = returnReduce(srcPatch.size(), sumOp<label>());
874  const label tgtTotalSize = returnReduce(tgtPatch.size(), sumOp<label>());
875  if (srcTotalSize == 0 || tgtTotalSize == 0)
876  {
877  return;
878  }
879 
880  // If a transformation is given then transform the target to the source
881  tmpNrc<primitiveOldTimePatch> tTgtPatchPtr(tgtPatch);
882  tmpNrc<pointField> tTgtPointsPtr(tgtPatch.localPoints());
883  tmpNrc<pointField> tTgtPoints0Ptr
884  (
885  tgtPatch.has0() ? tgtPatch.localPoints0() : NullObjectRef<pointField>()
886  );
887  if (!isNull(tgtToSrc))
888  {
889  tTgtPointsPtr = new pointField(tgtPatch.localPoints());
890  tgtToSrc.transformPosition
891  (
892  tTgtPointsPtr.ref(),
893  tTgtPointsPtr.ref()
894  );
895 
896  if (tgtPatch.has0())
897  {
898  tTgtPoints0Ptr = new pointField(tgtPatch.localPoints0());
899  tgtToSrc.transformPosition
900  (
901  tTgtPoints0Ptr.ref(),
902  tTgtPoints0Ptr.ref()
903  );
904  }
905 
906  tTgtPatchPtr =
907  tgtPatch.has0()
909  (
910  SubList<face>(tgtPatch.localFaces(), tgtPatch.size()),
911  tTgtPointsPtr(),
912  tTgtPoints0Ptr()
913  )
915  (
916  SubList<face>(tgtPatch.localFaces(), tgtPatch.size()),
917  tTgtPointsPtr()
918  );
919  }
920  const primitiveOldTimePatch& tTgtPatch = tTgtPatchPtr();
921 
922  Info<< indent << typeName << ": Calculating couplings between "
923  << srcTotalSize << " source faces and " << tgtTotalSize
924  << " target faces" << incrIndent << endl;
925 
926  // Determine if patches are present on multiple processors
927  singleProcess_ =
929  (
930  srcPatch.size(),
931  tTgtPatch.size()
932  );
933 
934  // Do intersection in serial or parallel as appropriate
935  if (isSingleProcess())
936  {
937  // Initialise the workspace
938  initialise(srcPatch, srcPointNormals, srcPointNormals0, tTgtPatch);
939 
940  // Intersect the patches
941  const treeBoundBox srcPatchBox =
942  srcBox(srcPatch, srcPointNormals, srcPointNormals0);
943  const treeBoundBox tgtPatchBox = tgtBox(tTgtPatch);
944  if (srcPatchBox.overlaps(tgtPatchBox))
945  {
946  intersectPatches
947  (
948  srcPatch,
949  srcPointNormals,
950  srcPointNormals0,
951  tTgtPatch
952  );
953  }
954  }
955  else
956  {
957  // Distribute the target patch so that everything is locally available
958  // to the source. This is done based on bound boxes, so quite a lot of
959  // faces will get distributed that ultimately are not used. These will
960  // be filtered out after the intersection has been completed.
961  tgtMapPtr_ =
963  (
964  tgtPatchSendFaces
965  (
966  srcPatch,
967  srcPointNormals,
968  srcPointNormals0,
969  tTgtPatch
970  )
971  );
973  localTgtProcFacesPtr_.reset
974  (
975  new List<remote>
976  (
977  distributePatch
978  (
979  tgtMapPtr_(),
980  tTgtPatch,
981  localTTgtPatchPtr
982  )
983  )
984  );
985 
986  // Massage target patch into form that can be used by the serial
987  // intersection interface
988  const primitiveOldTimePatch localTTgtPatch =
989  tgtPatch.has0()
991  (
992  SubList<face>(localTTgtPatchPtr(), localTTgtPatchPtr().size()),
993  localTTgtPatchPtr().points(),
994  localTTgtPatchPtr().points0()
995  )
997  (
998  SubList<face>(localTTgtPatchPtr(), localTTgtPatchPtr().size()),
999  localTTgtPatchPtr().points()
1000  );
1001 
1002  // Initialise the workspace
1003  initialise(srcPatch, srcPointNormals, srcPointNormals, localTTgtPatch);
1004 
1005  // Intersect the patches
1006  if (localTTgtPatch.size())
1007  {
1008  intersectPatches
1009  (
1010  srcPatch,
1011  srcPointNormals,
1012  srcPointNormals,
1013  localTTgtPatch
1014  );
1015  }
1016 
1017  // Trim the local target patch
1018  finaliseLocal
1019  (
1020  srcPatch,
1021  srcPointNormals,
1022  srcPointNormals,
1023  localTTgtPatch
1024  );
1025 
1026  // Distribute the source patch
1027  srcMapPtr_ =
1029  (
1031  (
1032  tgtLocalSrcFaces_,
1033  localTgtProcFacesPtr_()
1034  )
1035  );
1036  distributeSrc(srcPatch);
1037 
1038  // Reverse distribute coupling data back to the target
1039  rDistributeTgt(tgtPatch);
1040  }
1041 
1042  // Finalise the intersection
1043  const label nCouples =
1044  finalise
1045  (
1046  srcPatch,
1047  srcPointNormals,
1048  srcPointNormals,
1049  tgtPatch,
1050  tgtToSrc
1051  );
1052 
1053  if (nCouples != 0)
1054  {
1055  Info<< indent << nCouples << " couplings calculated in "
1056  << time.cpuTimeIncrement() << 's' << endl;
1057  }
1058  else
1059  {
1060  Info<< indent << "No couplings found" << endl;
1061  }
1062 
1063  Info<< decrIndent;
1064 }
1065 
1066 
1068 (
1069  const primitivePatch& srcPatch,
1070  const vectorField& srcPointNormals,
1071  const primitivePatch& tgtPatch,
1072  const transformer& tgtToSrc
1073 )
1074 {
1075  update
1076  (
1077  primitiveOldTimePatch(srcPatch),
1078  srcPointNormals,
1079  NullObjectRef<vectorField>(),
1080  primitiveOldTimePatch(tgtPatch),
1081  tgtToSrc
1082  );
1083 }
1084 
1085 
1086 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:78
void transfer(List< T > &)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:285
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
void max(const dimensioned< Type > &)
void min(const dimensioned< Type > &)
A List with indirect addressing.
Definition: IndirectList.H:105
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
A bit-packed bool list.
bool empty() const
Return true if the list is empty (ie, size() is zero).
Definition: PackedListI.H:718
bool has0() const
Return whether or not old-time geometry is available.
const Field< PointType > & localPoints0() const
Return pointField of old-time points in patch.
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
const labelListList & pointFaces() const
Return point-face addressing.
const labelListList & edgeFaces() const
Return edge-face addressing.
const labelListList & faceEdges() const
Return face-edge addressing.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
const Field< PointType > & faceCentres() const
Return face centres for patch.
A List obtained as a section of another List.
Definition: SubList.H:56
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:74
T & first()
Return the first element of the list.
Definition: UListI.H:114
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:60
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:66
void inflate(const scalar s)
Inflate box by factor*mag(span) in all dimensions.
Definition: boundBox.C:199
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTime.H:55
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:60
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
A class for handling file names.
Definition: fileName.H:82
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:72
Class to generate coupling geometry between two primitive patches.
Definition: patchToPatch.H:56
virtual treeBoundBox srcBox(const face &srcFace, const pointField &srcPoints, const vectorField &srcPointNormals) const =0
Get the bound box for a source face.
virtual ~patchToPatch()
Destructor.
Definition: patchToPatch.C:779
virtual void rDistributeTgt(const primitiveOldTimePatch &tgtPatch)
Send data that resulted from an intersection between the source.
Definition: patchToPatch.C:724
virtual labelList finaliseLocal(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch)
Finalise the intersection locally. Trims the local target patch.
Definition: patchToPatch.C:662
virtual void distributeSrc(const primitiveOldTimePatch &srcPatch)
Distribute the source patch so that any information needed by.
Definition: patchToPatch.C:710
treeBoundBox tgtBox(const primitiveOldTimePatch &tgtPatch, const label tgtFacei) const
Get the bound box for a target face.
Definition: patchToPatch.C:144
static autoPtr< patchToPatch > New(const word &patchToPatchType, const bool reverse)
Select from name.
Definition: patchToPatch.C:786
List< List< remote > > tgtSrcProcFaces() const
For each target face, the coupled source procs and faces.
Definition: patchToPatch.C:847
PackedBoolList srcCoupled() const
Return a list indicating which source faces are coupled.
Definition: patchToPatch.C:810
bool findOrIntersectFaces(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch, const label srcFacei, const label tgtFacei)
Intersect two faces.
Definition: patchToPatch.C:182
List< List< remote > > srcTgtProcFaces() const
For each source face, the coupled target procs and faces.
Definition: patchToPatch.C:833
virtual void initialise(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch)
Initialisation.
Definition: patchToPatch.C:640
void intersectPatches(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch)
Intersect the patches.
Definition: patchToPatch.C:374
label intersectPatchQueue(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch, const bool isSrc, const DynamicList< labelPair > &queue, labelList &faceComplete, DynamicList< labelPair > &otherQueue, const labelList &otherFaceComplete, boolList &otherFaceQueued, boolList &otherFaceVisited)
Intersect a queue of source-target face pairs. Update completion.
Definition: patchToPatch.C:232
void update(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch, const transformer &tgtToSrc=NullObjectRef< transformer >())
Update addressing and weights for the given patches.
Definition: patchToPatch.C:861
PackedBoolList tgtCoupled() const
Return a list indicating which target faces are coupled.
Definition: patchToPatch.C:821
patchToPatch(const bool reverse)
Construct from components.
Definition: patchToPatch.C:764
virtual label finalise(const primitiveOldTimePatch &srcPatch, const vectorField &srcPointNormals, const vectorField &srcPointNormals0, const primitiveOldTimePatch &tgtPatch, const transformer &tgtToSrc)
Finalising.
Definition: patchToPatch.C:737
A class for managing temporary objects without reference counting.
Definition: tmpNrc.H:53
T & ref()
Return non-const reference or generate a fatal error.
Definition: tmpNrcI.H:149
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space.
Definition: transformer.H:84
vector transformPosition(const vector &v) const
Transform the given position.
Definition: transformerI.H:153
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:90
treeBoundBox extend(const scalar s) const
Return asymmetrically extended bounding box, with guaranteed.
bool overlaps(const boundBox &) const
Overlaps other bounding box?
Definition: boundBoxI.H:126
static const treeBoundBox invertedBox
As per boundBox::invertedBox, but with great instead of vGreat.
Definition: treeBoundBox.H:108
Encapsulation of data needed to search on PrimitivePatches.
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const pointField & points
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField & b
Definition: createFields.H:25
label singleProcess(const label srcSize, const label tgtSize)
Determine whether this intersection is confined to a single processor or.
autoPtr< distributionMap > constructDistributionMap(const labelListList &procSendIndices)
Turn a list of send-to-processor indices into a distribution map.
List< remote > distributeAddressing(const distributionMap &map)
Construct local addressing from the given distribution map. The result is a.
void rDistributeTgtAddressing(const label tgtSize, const distributionMap &tgtMap, const List< remote > &localSrcProcFaces, labelListList &tgtLocalSrcFaces)
Reverse distribute a set of target addressing.
labelListList procSendIndices(const labelListList &tgtLocalSrcFaces, const List< remote > &localTgtProcFaces)
Given a local intersection addressing, determine what elements need sending.
void trimDistributionMap(const boolList &oldIsUsed, distributionMap &map, labelList &oldToNew, labelList &newToOld)
Trim unused elements from a distribution map. Return the map and maps.
List< List< remote > > localToRemote(const labelListList &indices, const List< remote > &indexToProcIndex=NullObjectRef< List< remote >>())
Map local to remote indices.
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.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
const doubleScalar e
Definition: doubleScalar.H:106
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:241
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
label findNotIndex(const ListType &l, typename ListType::const_reference t, const label start=0)
Definition: patchToPatch.C:61
messageStream Info
labelList second(const UList< labelPair > &p)
Definition: patchToPatch.C:49
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:234
void reverse(UList< T > &, const label n)
Definition: UListI.H:334
labelList first(const UList< labelPair > &p)
Definition: patchToPatch.C:39
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
PrimitiveOldTimePatch< SubList< face >, const pointField & > primitiveOldTimePatch
Addressing for a faceList slice.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
bool isNull(const T &t)
Return true if t is a reference to the nullObject of type T.
Definition: nullObjectI.H:58
error FatalError
static const label labelMax
Definition: label.H:62
treeBoundBox combine(const treeBoundBox &a, const treeBoundBox &b)
Definition: patchToPatch.C:78
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:227
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
Definition: POSIX.C:908
labelList f(nPoints)
volScalarField & p