AMIInterpolation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2017 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 "AMIInterpolation.H"
27 #include "AMIMethod.H"
28 #include "meshTools.H"
29 #include "mapDistribute.H"
30 #include "flipOp.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 template<class SourcePatch, class TargetPatch>
37 (
38  const interpolationMethod& im
39 )
40 {
41  word method = "unknown-interpolationMethod";
42 
43  switch (im)
44  {
45  case imDirect:
46  {
47  method = "directAMI";
48  break;
49  }
50  case imMapNearest:
51  {
52  method = "mapNearestAMI";
53  break;
54  }
55  case imFaceAreaWeight:
56  {
57  method = "faceAreaWeightAMI";
58  break;
59  }
60  case imPartialFaceAreaWeight:
61  {
62  method = "partialFaceAreaWeightAMI";
63  break;
64  }
65  default:
66  {
68  << "Unhandled interpolationMethod enumeration " << method
69  << abort(FatalError);
70  }
71  }
72 
73  return method;
74 }
75 
76 
77 template<class SourcePatch, class TargetPatch>
80 (
81  const word& im
82 )
83 {
84  interpolationMethod method = imDirect;
85 
86  wordList methods
87  (
89  (
90  "("
91  "directAMI "
92  "mapNearestAMI "
93  "faceAreaWeightAMI "
94  "partialFaceAreaWeightAMI"
95  ")"
96  )()
97  );
98 
99  if (im == "directAMI")
100  {
101  method = imDirect;
102  }
103  else if (im == "mapNearestAMI")
104  {
105  method = imMapNearest;
106  }
107  else if (im == "faceAreaWeightAMI")
108  {
109  method = imFaceAreaWeight;
110  }
111  else if (im == "partialFaceAreaWeightAMI")
112  {
113  method = imPartialFaceAreaWeight;
114  }
115  else
116  {
118  << "Invalid interpolationMethod " << im
119  << ". Valid methods are:" << methods
120  << exit(FatalError);
121  }
122 
123  return method;
124 }
125 
126 
127 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
128 
129 template<class SourcePatch, class TargetPatch>
131 (
132  const searchableSurface& surf,
133  pointField& pts
134 ) const
135 {
136  if (debug)
137  {
138  Info<< "AMI: projecting points to surface" << endl;
139  }
140 
142 
143  surf.findNearest(pts, scalarField(pts.size(), GREAT), nearInfo);
144 
145  label nMiss = 0;
146  forAll(nearInfo, i)
147  {
148  const pointIndexHit& pi = nearInfo[i];
149 
150  if (pi.hit())
151  {
152  pts[i] = pi.hitPoint();
153  }
154  else
155  {
156  pts[i] = pts[i];
157  nMiss++;
158  }
159  }
160 
161  if (nMiss > 0)
162  {
164  << "Error projecting points to surface: "
165  << nMiss << " faces could not be determined"
166  << abort(FatalError);
167  }
168 }
169 
170 
171 template<class SourcePatch, class TargetPatch>
173 (
174  const scalarField& patchAreas,
175  const word& patchName,
176  const labelListList& addr,
177  scalarListList& wght,
178  scalarField& wghtSum,
179  const bool conformal,
180  const bool output,
181  const scalar lowWeightTol
182 )
183 {
184  // Normalise the weights
185  wghtSum.setSize(wght.size(), 0.0);
186  label nLowWeight = 0;
187 
188  forAll(wght, facei)
189  {
190  scalarList& w = wght[facei];
191 
192  if (w.size())
193  {
194  scalar denom = patchAreas[facei];
195 
196  scalar s = sum(w);
197  scalar t = s/denom;
198 
199  if (conformal)
200  {
201  denom = s;
202  }
203 
204  forAll(w, i)
205  {
206  w[i] /= denom;
207  }
208 
209  wghtSum[facei] = t;
210 
211  if (t < lowWeightTol)
212  {
213  nLowWeight++;
214  }
215  }
216  else
217  {
218  wghtSum[facei] = 0;
219  }
220  }
221 
222 
223  if (output)
224  {
225  const label nFace = returnReduce(wght.size(), sumOp<label>());
226 
227  if (nFace)
228  {
229  Info<< indent
230  << "AMI: Patch " << patchName
231  << " sum(weights) min/max/average = "
232  << gMin(wghtSum) << ", "
233  << gMax(wghtSum) << ", "
234  << gAverage(wghtSum) << endl;
235 
236  const label nLow = returnReduce(nLowWeight, sumOp<label>());
237 
238  if (nLow)
239  {
240  Info<< indent
241  << "AMI: Patch " << patchName
242  << " identified " << nLow
243  << " faces with weights less than " << lowWeightTol
244  << endl;
245  }
246  }
247  }
248 }
249 
250 
251 template<class SourcePatch, class TargetPatch>
253 (
254  const autoPtr<mapDistribute>& targetMapPtr,
255  const scalarField& fineSrcMagSf,
256  const labelListList& fineSrcAddress,
257  const scalarListList& fineSrcWeights,
258 
259  const labelList& sourceRestrictAddressing,
260  const labelList& targetRestrictAddressing,
261 
262  scalarField& srcMagSf,
263  labelListList& srcAddress,
264  scalarListList& srcWeights,
265  scalarField& srcWeightsSum,
266  autoPtr<mapDistribute>& tgtMap
267 )
268 {
269  label sourceCoarseSize =
270  (
271  sourceRestrictAddressing.size()
272  ? max(sourceRestrictAddressing)+1
273  : 0
274  );
275 
276  label targetCoarseSize =
277  (
278  targetRestrictAddressing.size()
279  ? max(targetRestrictAddressing)+1
280  : 0
281  );
282 
283  // Agglomerate face areas
284  {
285  srcMagSf.setSize(sourceRestrictAddressing.size(), 0.0);
286  forAll(sourceRestrictAddressing, facei)
287  {
288  label coarseFacei = sourceRestrictAddressing[facei];
289  srcMagSf[coarseFacei] += fineSrcMagSf[facei];
290  }
291  }
292 
293 
294  // Agglomerate weights and indices
295  if (targetMapPtr.valid())
296  {
297  const mapDistribute& map = targetMapPtr();
298 
299  // Get all restriction addressing.
300  labelList allRestrict(targetRestrictAddressing);
301  map.distribute(allRestrict);
302 
303  // So now we have agglomeration of the target side in
304  // allRestrict:
305  // 0..size-1 : local agglomeration (= targetRestrictAddressing)
306  // size.. : agglomeration data from other processors
307 
308  labelListList tgtSubMap(Pstream::nProcs());
309 
310  // Local subMap is just identity
311  {
312  tgtSubMap[Pstream::myProcNo()] = identity(targetCoarseSize);
313  }
314 
315  forAll(map.subMap(), proci)
316  {
317  if (proci != Pstream::myProcNo())
318  {
319  // Combine entries that point to the same coarse element. All
320  // the elements refer to local data so index into
321  // targetRestrictAddressing or allRestrict (since the same
322  // for local data).
323  const labelList& elems = map.subMap()[proci];
324  labelList& newSubMap = tgtSubMap[proci];
325  newSubMap.setSize(elems.size());
326 
327  labelList oldToNew(targetCoarseSize, -1);
328  label newI = 0;
329 
330  forAll(elems, i)
331  {
332  label fineElem = elems[i];
333  label coarseElem = allRestrict[fineElem];
334  if (oldToNew[coarseElem] == -1)
335  {
336  oldToNew[coarseElem] = newI;
337  newSubMap[newI] = coarseElem;
338  newI++;
339  }
340  }
341  newSubMap.setSize(newI);
342  }
343  }
344 
345  // Reconstruct constructMap by combining entries. Note that order
346  // of handing out indices should be the same as loop above to compact
347  // the sending map
348 
349  labelListList tgtConstructMap(Pstream::nProcs());
350  labelList tgtCompactMap;
351 
352  // Local constructMap is just identity
353  {
354  tgtConstructMap[Pstream::myProcNo()] =
355  identity(targetCoarseSize);
356 
357  tgtCompactMap = targetRestrictAddressing;
358  }
359  tgtCompactMap.setSize(map.constructSize());
360  label compactI = targetCoarseSize;
361 
362  // Compact data from other processors
363  forAll(map.constructMap(), proci)
364  {
365  if (proci != Pstream::myProcNo())
366  {
367  // Combine entries that point to the same coarse element. All
368  // elements now are remote data so we cannot use any local
369  // data here - use allRestrict instead.
370  const labelList& elems = map.constructMap()[proci];
371 
372  labelList& newConstructMap = tgtConstructMap[proci];
373  newConstructMap.setSize(elems.size());
374 
375  if (elems.size())
376  {
377  // Get the maximum target coarse size for this set of
378  // received data.
379  label remoteTargetCoarseSize = labelMin;
380  forAll(elems, i)
381  {
382  remoteTargetCoarseSize = max
383  (
384  remoteTargetCoarseSize,
385  allRestrict[elems[i]]
386  );
387  }
388  remoteTargetCoarseSize += 1;
389 
390  // Combine locally data coming from proci
391  labelList oldToNew(remoteTargetCoarseSize, -1);
392  label newI = 0;
393 
394  forAll(elems, i)
395  {
396  label fineElem = elems[i];
397  // fineElem now points to section from proci
398  label coarseElem = allRestrict[fineElem];
399  if (oldToNew[coarseElem] == -1)
400  {
401  oldToNew[coarseElem] = newI;
402  tgtCompactMap[fineElem] = compactI;
403  newConstructMap[newI] = compactI++;
404  newI++;
405  }
406  else
407  {
408  // Get compact index
409  label compactI = oldToNew[coarseElem];
410  tgtCompactMap[fineElem] = newConstructMap[compactI];
411  }
412  }
413  newConstructMap.setSize(newI);
414  }
415  }
416  }
417 
418 
419  srcAddress.setSize(sourceCoarseSize);
420  srcWeights.setSize(sourceCoarseSize);
421 
422  forAll(fineSrcAddress, facei)
423  {
424  // All the elements contributing to facei. Are slots in
425  // mapDistribute'd data.
426  const labelList& elems = fineSrcAddress[facei];
427  const scalarList& weights = fineSrcWeights[facei];
428  const scalar fineArea = fineSrcMagSf[facei];
429 
430  label coarseFacei = sourceRestrictAddressing[facei];
431 
432  labelList& newElems = srcAddress[coarseFacei];
433  scalarList& newWeights = srcWeights[coarseFacei];
434 
435  forAll(elems, i)
436  {
437  label elemI = elems[i];
438  label coarseElemI = tgtCompactMap[elemI];
439 
440  label index = findIndex(newElems, coarseElemI);
441  if (index == -1)
442  {
443  newElems.append(coarseElemI);
444  newWeights.append(fineArea*weights[i]);
445  }
446  else
447  {
448  newWeights[index] += fineArea*weights[i];
449  }
450  }
451  }
452 
453  tgtMap.reset
454  (
455  new mapDistribute
456  (
457  compactI,
458  tgtSubMap.xfer(),
459  tgtConstructMap.xfer()
460  )
461  );
462  }
463  else
464  {
465  srcAddress.setSize(sourceCoarseSize);
466  srcWeights.setSize(sourceCoarseSize);
467 
468  forAll(fineSrcAddress, facei)
469  {
470  // All the elements contributing to facei. Are slots in
471  // mapDistribute'd data.
472  const labelList& elems = fineSrcAddress[facei];
473  const scalarList& weights = fineSrcWeights[facei];
474  const scalar fineArea = fineSrcMagSf[facei];
475 
476  label coarseFacei = sourceRestrictAddressing[facei];
477 
478  labelList& newElems = srcAddress[coarseFacei];
479  scalarList& newWeights = srcWeights[coarseFacei];
480 
481  forAll(elems, i)
482  {
483  label elemI = elems[i];
484  label coarseElemI = targetRestrictAddressing[elemI];
485 
486  label index = findIndex(newElems, coarseElemI);
487  if (index == -1)
488  {
489  newElems.append(coarseElemI);
490  newWeights.append(fineArea*weights[i]);
491  }
492  else
493  {
494  newWeights[index] += fineArea*weights[i];
495  }
496  }
497  }
498  }
499 
500  // Weights normalisation
501  normaliseWeights
502  (
503  srcMagSf,
504  "source",
505  srcAddress,
506  srcWeights,
507  srcWeightsSum,
508  true,
509  false,
510  -1
511  );
512 }
513 
514 
515 template<class SourcePatch, class TargetPatch>
517 (
518  const SourcePatch& srcPatch,
519  const TargetPatch& tgtPatch,
520  const autoPtr<searchableSurface>& surfPtr
521 )
522 {
523  if (surfPtr.valid())
524  {
525  // Create new patches for source and target
526  pointField srcPoints = srcPatch.points();
527  SourcePatch srcPatch0
528  (
530  (
531  srcPatch,
532  srcPatch.size(),
533  0
534  ),
535  srcPoints
536  );
537 
538  if (debug)
539  {
540  OFstream os("amiSrcPoints.obj");
541  forAll(srcPoints, i)
542  {
543  meshTools::writeOBJ(os, srcPoints[i]);
544  }
545  }
546 
547  pointField tgtPoints = tgtPatch.points();
548  TargetPatch tgtPatch0
549  (
551  (
552  tgtPatch,
553  tgtPatch.size(),
554  0
555  ),
556  tgtPoints
557  );
558 
559  if (debug)
560  {
561  OFstream os("amiTgtPoints.obj");
562  forAll(tgtPoints, i)
563  {
564  meshTools::writeOBJ(os, tgtPoints[i]);
565  }
566  }
567 
568 
569  // Map source and target patches onto projection surface
570  projectPointsToSurface(surfPtr(), srcPoints);
571  projectPointsToSurface(surfPtr(), tgtPoints);
572 
573 
574  // Calculate AMI interpolation
575  update(srcPatch0, tgtPatch0);
576  }
577  else
578  {
579  update(srcPatch, tgtPatch);
580  }
581 }
582 
583 
584 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
585 
586 template<class SourcePatch, class TargetPatch>
588 (
589  const SourcePatch& srcPatch,
590  const TargetPatch& tgtPatch,
592  const bool requireMatch,
593  const interpolationMethod& method,
594  const scalar lowWeightCorrection,
595  const bool reverseTarget
596 )
597 :
598  methodName_(interpolationMethodToWord(method)),
599  reverseTarget_(reverseTarget),
600  requireMatch_(requireMatch),
601  singlePatchProc_(-999),
602  lowWeightCorrection_(lowWeightCorrection),
603  srcAddress_(),
604  srcWeights_(),
605  srcWeightsSum_(),
606  tgtAddress_(),
607  tgtWeights_(),
608  tgtWeightsSum_(),
609  triMode_(triMode),
610  srcMapPtr_(nullptr),
611  tgtMapPtr_(nullptr)
612 {
613  update(srcPatch, tgtPatch);
614 }
615 
616 
617 template<class SourcePatch, class TargetPatch>
619 (
620  const SourcePatch& srcPatch,
621  const TargetPatch& tgtPatch,
623  const bool requireMatch,
624  const word& methodName,
625  const scalar lowWeightCorrection,
626  const bool reverseTarget
627 )
628 :
629  methodName_(methodName),
630  reverseTarget_(reverseTarget),
631  requireMatch_(requireMatch),
632  singlePatchProc_(-999),
633  lowWeightCorrection_(lowWeightCorrection),
634  srcAddress_(),
635  srcWeights_(),
636  srcWeightsSum_(),
637  tgtAddress_(),
638  tgtWeights_(),
639  tgtWeightsSum_(),
640  triMode_(triMode),
641  srcMapPtr_(nullptr),
642  tgtMapPtr_(nullptr)
643 {
644  update(srcPatch, tgtPatch);
645 }
646 
647 
648 template<class SourcePatch, class TargetPatch>
650 (
651  const SourcePatch& srcPatch,
652  const TargetPatch& tgtPatch,
653  const autoPtr<searchableSurface>& surfPtr,
655  const bool requireMatch,
656  const interpolationMethod& method,
657  const scalar lowWeightCorrection,
658  const bool reverseTarget
659 )
660 :
661  methodName_(interpolationMethodToWord(method)),
662  reverseTarget_(reverseTarget),
663  requireMatch_(requireMatch),
664  singlePatchProc_(-999),
665  lowWeightCorrection_(lowWeightCorrection),
666  srcAddress_(),
667  srcWeights_(),
668  srcWeightsSum_(),
669  tgtAddress_(),
670  tgtWeights_(),
671  tgtWeightsSum_(),
672  triMode_(triMode),
673  srcMapPtr_(nullptr),
674  tgtMapPtr_(nullptr)
675 {
676  constructFromSurface(srcPatch, tgtPatch, surfPtr);
677 }
678 
679 
680 template<class SourcePatch, class TargetPatch>
682 (
683  const SourcePatch& srcPatch,
684  const TargetPatch& tgtPatch,
685  const autoPtr<searchableSurface>& surfPtr,
687  const bool requireMatch,
688  const word& methodName,
689  const scalar lowWeightCorrection,
690  const bool reverseTarget
691 )
692 :
693  methodName_(methodName),
694  reverseTarget_(reverseTarget),
695  requireMatch_(requireMatch),
696  singlePatchProc_(-999),
697  lowWeightCorrection_(lowWeightCorrection),
698  srcAddress_(),
699  srcWeights_(),
700  srcWeightsSum_(),
701  tgtAddress_(),
702  tgtWeights_(),
703  tgtWeightsSum_(),
704  triMode_(triMode),
705  srcMapPtr_(nullptr),
706  tgtMapPtr_(nullptr)
707 {
708  constructFromSurface(srcPatch, tgtPatch, surfPtr);
709 }
710 
711 
712 template<class SourcePatch, class TargetPatch>
714 (
716  const labelList& sourceRestrictAddressing,
717  const labelList& targetRestrictAddressing
718 )
719 :
720  methodName_(fineAMI.methodName_),
721  reverseTarget_(fineAMI.reverseTarget_),
722  requireMatch_(fineAMI.requireMatch_),
723  singlePatchProc_(fineAMI.singlePatchProc_),
724  lowWeightCorrection_(-1.0),
725  srcAddress_(),
726  srcWeights_(),
727  srcWeightsSum_(),
728  tgtAddress_(),
729  tgtWeights_(),
730  tgtWeightsSum_(),
731  triMode_(fineAMI.triMode_),
732  srcMapPtr_(nullptr),
733  tgtMapPtr_(nullptr)
734 {
735  label sourceCoarseSize =
736  (
737  sourceRestrictAddressing.size()
738  ? max(sourceRestrictAddressing)+1
739  : 0
740  );
741 
742  label neighbourCoarseSize =
743  (
744  targetRestrictAddressing.size()
745  ? max(targetRestrictAddressing)+1
746  : 0
747  );
748 
749  if (debug & 2)
750  {
751  Pout<< "AMI: Creating addressing and weights as agglomeration of AMI :"
752  << " source:" << fineAMI.srcAddress().size()
753  << " target:" << fineAMI.tgtAddress().size()
754  << " coarse source size:" << sourceCoarseSize
755  << " neighbour source size:" << neighbourCoarseSize
756  << endl;
757  }
758 
759  if
760  (
761  fineAMI.srcAddress().size() != sourceRestrictAddressing.size()
762  || fineAMI.tgtAddress().size() != targetRestrictAddressing.size()
763  )
764  {
766  << "Size mismatch." << nl
767  << "Source patch size:" << fineAMI.srcAddress().size() << nl
768  << "Source agglomeration size:"
769  << sourceRestrictAddressing.size() << nl
770  << "Target patch size:" << fineAMI.tgtAddress().size() << nl
771  << "Target agglomeration size:"
772  << targetRestrictAddressing.size()
773  << exit(FatalError);
774  }
775 
776 
777  // Agglomerate addresses and weights
778 
779  agglomerate
780  (
781  fineAMI.tgtMapPtr_,
782  fineAMI.srcMagSf(),
783  fineAMI.srcAddress(),
784  fineAMI.srcWeights(),
785 
786  sourceRestrictAddressing,
787  targetRestrictAddressing,
788 
789  srcMagSf_,
790  srcAddress_,
791  srcWeights_,
792  srcWeightsSum_,
793  tgtMapPtr_
794  );
795 
796  agglomerate
797  (
798  fineAMI.srcMapPtr_,
799  fineAMI.tgtMagSf(),
800  fineAMI.tgtAddress(),
801  fineAMI.tgtWeights(),
802 
803  targetRestrictAddressing,
804  sourceRestrictAddressing,
805 
806  tgtMagSf_,
807  tgtAddress_,
808  tgtWeights_,
809  tgtWeightsSum_,
810  srcMapPtr_
811  );
812 }
813 
814 
815 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
816 
817 template<class SourcePatch, class TargetPatch>
819 {}
820 
821 
822 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
823 
824 template<class SourcePatch, class TargetPatch>
826 (
827  const SourcePatch& srcPatch,
828  const TargetPatch& tgtPatch
829 )
830 {
831  label srcTotalSize = returnReduce(srcPatch.size(), sumOp<label>());
832  label tgtTotalSize = returnReduce(tgtPatch.size(), sumOp<label>());
833 
834  if (srcTotalSize == 0)
835  {
836  if (debug)
837  {
838  Info<< "AMI: no source faces present - no addressing constructed"
839  << endl;
840  }
841 
842  return;
843  }
844 
845  Info<< indent
846  << "AMI: Creating addressing and weights between "
847  << srcTotalSize << " source faces and "
848  << tgtTotalSize << " target faces"
849  << endl;
850 
851  // Calculate face areas
852  srcMagSf_.setSize(srcPatch.size());
853  forAll(srcMagSf_, facei)
854  {
855  srcMagSf_[facei] = srcPatch[facei].mag(srcPatch.points());
856  }
857  tgtMagSf_.setSize(tgtPatch.size());
858  forAll(tgtMagSf_, facei)
859  {
860  tgtMagSf_[facei] = tgtPatch[facei].mag(tgtPatch.points());
861  }
862 
863  // Calculate if patches present on multiple processors
864  singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
865 
866  if (singlePatchProc_ == -1)
867  {
868  // Convert local addressing to global addressing
869  globalIndex globalSrcFaces(srcPatch.size());
870  globalIndex globalTgtFaces(tgtPatch.size());
871 
872  // Create processor map of overlapping faces. This map gets
873  // (possibly remote) faces from the tgtPatch such that they (together)
874  // cover all of the srcPatch
875  autoPtr<mapDistribute> mapPtr = calcProcMap(srcPatch, tgtPatch);
876  const mapDistribute& map = mapPtr();
877 
878  // Create new target patch that fully encompasses source patch
879 
880  // Faces and points
881  faceList newTgtFaces;
882  pointField newTgtPoints;
883 
884  // Original faces from tgtPatch (in globalIndexing since might be
885  // remote)
886  labelList tgtFaceIDs;
887  distributeAndMergePatches
888  (
889  map,
890  tgtPatch,
891  globalTgtFaces,
892  newTgtFaces,
893  newTgtPoints,
894  tgtFaceIDs
895  );
896 
897  TargetPatch
898  newTgtPatch
899  (
901  (
902  newTgtFaces,
903  newTgtFaces.size()
904  ),
905  newTgtPoints
906  );
907 
908  // Calculate AMI interpolation
910  (
912  (
913  methodName_,
914  srcPatch,
915  newTgtPatch,
916  srcMagSf_,
917  tgtMagSf_,
918  triMode_,
919  reverseTarget_,
920  requireMatch_ && (lowWeightCorrection_ < 0)
921  )
922  );
923 
924  AMIPtr->calculate
925  (
926  srcAddress_,
927  srcWeights_,
928  tgtAddress_,
929  tgtWeights_
930  );
931 
932  // Now
933  // ~~~
934  // srcAddress_ : per srcPatch face a list of the newTgtPatch (not
935  // tgtPatch) faces it overlaps
936  // tgtAddress_ : per newTgtPatch (not tgtPatch) face a list of the
937  // srcPatch faces it overlaps
938 
939 
940  // Rework newTgtPatch indices into globalIndices of tgtPatch
941  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
942 
943 
944  forAll(srcAddress_, i)
945  {
946  labelList& addressing = srcAddress_[i];
947  forAll(addressing, addrI)
948  {
949  addressing[addrI] = tgtFaceIDs[addressing[addrI]];
950  }
951  }
952 
953  forAll(tgtAddress_, i)
954  {
955  labelList& addressing = tgtAddress_[i];
956  forAll(addressing, addrI)
957  {
958  addressing[addrI] = globalSrcFaces.toGlobal(addressing[addrI]);
959  }
960  }
961 
962  // Send data back to originating procs. Note that contributions
963  // from different processors get added (ListAppendEqOp)
964 
965  mapDistributeBase::distribute
966  (
967  Pstream::commsTypes::nonBlocking,
968  List<labelPair>(),
969  tgtPatch.size(),
970  map.constructMap(),
971  false, // has flip
972  map.subMap(),
973  false, // has flip
974  tgtAddress_,
976  flipOp(), // flip operation
977  labelList()
978  );
979 
980  mapDistributeBase::distribute
981  (
982  Pstream::commsTypes::nonBlocking,
983  List<labelPair>(),
984  tgtPatch.size(),
985  map.constructMap(),
986  false,
987  map.subMap(),
988  false,
989  tgtWeights_,
991  flipOp(),
992  scalarList()
993  );
994 
995  // Weights normalisation
996  normaliseWeights
997  (
998  srcMagSf_,
999  "source",
1000  srcAddress_,
1001  srcWeights_,
1002  srcWeightsSum_,
1003  AMIPtr->conformal(),
1004  true,
1005  lowWeightCorrection_
1006  );
1007  normaliseWeights
1008  (
1009  tgtMagSf_,
1010  "target",
1011  tgtAddress_,
1012  tgtWeights_,
1013  tgtWeightsSum_,
1014  AMIPtr->conformal(),
1015  true,
1016  lowWeightCorrection_
1017  );
1018 
1019  // Cache maps and reset addresses
1020  List<Map<label>> cMap;
1021  srcMapPtr_.reset(new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
1022  tgtMapPtr_.reset(new mapDistribute(globalTgtFaces, srcAddress_, cMap));
1023 
1024  if (debug)
1025  {
1026  writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
1027  }
1028  }
1029  else
1030  {
1031  // Calculate AMI interpolation
1033  (
1035  (
1036  methodName_,
1037  srcPatch,
1038  tgtPatch,
1039  srcMagSf_,
1040  tgtMagSf_,
1041  triMode_,
1042  reverseTarget_,
1043  requireMatch_ && (lowWeightCorrection_ < 0)
1044  )
1045  );
1046 
1047  AMIPtr->calculate
1048  (
1049  srcAddress_,
1050  srcWeights_,
1051  tgtAddress_,
1052  tgtWeights_
1053  );
1054 
1055  normaliseWeights
1056  (
1057  srcMagSf_,
1058  "source",
1059  srcAddress_,
1060  srcWeights_,
1061  srcWeightsSum_,
1062  AMIPtr->conformal(),
1063  true,
1064  lowWeightCorrection_
1065  );
1066  normaliseWeights
1067  (
1068  tgtMagSf_,
1069  "target",
1070  tgtAddress_,
1071  tgtWeights_,
1072  tgtWeightsSum_,
1073  AMIPtr->conformal(),
1074  true,
1075  lowWeightCorrection_
1076  );
1077  }
1078 
1079  if (debug)
1080  {
1081  Info<< "AMIInterpolation : Constructed addressing and weights" << nl
1082  << " triMode :"
1083  << faceAreaIntersect::triangulationModeNames_[triMode_] << nl
1084  << " singlePatchProc:" << singlePatchProc_ << nl
1085  << " srcMagSf :" << gSum(srcMagSf_) << nl
1086  << " tgtMagSf :" << gSum(tgtMagSf_) << nl
1087  << endl;
1088  }
1089 }
1090 
1091 
1092 template<class SourcePatch, class TargetPatch>
1093 template<class Type, class CombineOp>
1096  const UList<Type>& fld,
1097  const CombineOp& cop,
1098  List<Type>& result,
1099  const UList<Type>& defaultValues
1100 ) const
1101 {
1102  if (fld.size() != srcAddress_.size())
1103  {
1105  << "Supplied field size is not equal to source patch size" << nl
1106  << " source patch = " << srcAddress_.size() << nl
1107  << " target patch = " << tgtAddress_.size() << nl
1108  << " supplied field = " << fld.size()
1109  << abort(FatalError);
1110  }
1111 
1112  if (lowWeightCorrection_ > 0)
1113  {
1114  if (defaultValues.size() != tgtAddress_.size())
1115  {
1117  << "Employing default values when sum of weights falls below "
1118  << lowWeightCorrection_
1119  << " but supplied default field size is not equal to target "
1120  << "patch size" << nl
1121  << " default values = " << defaultValues.size() << nl
1122  << " target patch = " << tgtAddress_.size() << nl
1123  << abort(FatalError);
1124  }
1125  }
1126 
1127  result.setSize(tgtAddress_.size());
1128 
1129  if (singlePatchProc_ == -1)
1130  {
1131  const mapDistribute& map = srcMapPtr_();
1132 
1133  List<Type> work(fld);
1134  map.distribute(work);
1135 
1136  forAll(result, facei)
1137  {
1138  if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1139  {
1140  result[facei] = defaultValues[facei];
1141  }
1142  else
1143  {
1144  const labelList& faces = tgtAddress_[facei];
1145  const scalarList& weights = tgtWeights_[facei];
1146 
1147  forAll(faces, i)
1148  {
1149  cop(result[facei], facei, work[faces[i]], weights[i]);
1150  }
1151  }
1152  }
1153  }
1154  else
1155  {
1156  forAll(result, facei)
1157  {
1158  if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1159  {
1160  result[facei] = defaultValues[facei];
1161  }
1162  else
1163  {
1164  const labelList& faces = tgtAddress_[facei];
1165  const scalarList& weights = tgtWeights_[facei];
1166 
1167  forAll(faces, i)
1168  {
1169  cop(result[facei], facei, fld[faces[i]], weights[i]);
1170  }
1171  }
1172  }
1173  }
1174 }
1175 
1176 
1177 template<class SourcePatch, class TargetPatch>
1178 template<class Type, class CombineOp>
1181  const UList<Type>& fld,
1182  const CombineOp& cop,
1183  List<Type>& result,
1184  const UList<Type>& defaultValues
1185 ) const
1186 {
1187  if (fld.size() != tgtAddress_.size())
1188  {
1190  << "Supplied field size is not equal to target patch size" << nl
1191  << " source patch = " << srcAddress_.size() << nl
1192  << " target patch = " << tgtAddress_.size() << nl
1193  << " supplied field = " << fld.size()
1194  << abort(FatalError);
1195  }
1196 
1197  if (lowWeightCorrection_ > 0)
1198  {
1199  if (defaultValues.size() != srcAddress_.size())
1200  {
1202  << "Employing default values when sum of weights falls below "
1203  << lowWeightCorrection_
1204  << " but supplied default field size is not equal to target "
1205  << "patch size" << nl
1206  << " default values = " << defaultValues.size() << nl
1207  << " source patch = " << srcAddress_.size() << nl
1208  << abort(FatalError);
1209  }
1210  }
1211 
1212  result.setSize(srcAddress_.size());
1213 
1214  if (singlePatchProc_ == -1)
1215  {
1216  const mapDistribute& map = tgtMapPtr_();
1217 
1218  List<Type> work(fld);
1219  map.distribute(work);
1220 
1221  forAll(result, facei)
1222  {
1223  if (srcWeightsSum_[facei] < lowWeightCorrection_)
1224  {
1225  result[facei] = defaultValues[facei];
1226  }
1227  else
1228  {
1229  const labelList& faces = srcAddress_[facei];
1230  const scalarList& weights = srcWeights_[facei];
1231 
1232  forAll(faces, i)
1233  {
1234  cop(result[facei], facei, work[faces[i]], weights[i]);
1235  }
1236  }
1237  }
1238  }
1239  else
1240  {
1241  forAll(result, facei)
1242  {
1243  if (srcWeightsSum_[facei] < lowWeightCorrection_)
1244  {
1245  result[facei] = defaultValues[facei];
1246  }
1247  else
1248  {
1249  const labelList& faces = srcAddress_[facei];
1250  const scalarList& weights = srcWeights_[facei];
1251 
1252  forAll(faces, i)
1253  {
1254  cop(result[facei], facei, fld[faces[i]], weights[i]);
1255  }
1256  }
1257  }
1258  }
1259 }
1260 
1261 
1262 template<class SourcePatch, class TargetPatch>
1263 template<class Type, class CombineOp>
1267  const Field<Type>& fld,
1268  const CombineOp& cop,
1269  const UList<Type>& defaultValues
1270 ) const
1271 {
1272  tmp<Field<Type>> tresult
1273  (
1274  new Field<Type>
1275  (
1276  srcAddress_.size(),
1277  Zero
1278  )
1279  );
1280 
1281  interpolateToSource
1282  (
1283  fld,
1285  tresult.ref(),
1286  defaultValues
1287  );
1288 
1289  return tresult;
1290 }
1291 
1292 
1293 template<class SourcePatch, class TargetPatch>
1294 template<class Type, class CombineOp>
1298  const tmp<Field<Type>>& tFld,
1299  const CombineOp& cop,
1300  const UList<Type>& defaultValues
1301 ) const
1302 {
1303  return interpolateToSource(tFld(), cop, defaultValues);
1304 }
1305 
1306 
1307 template<class SourcePatch, class TargetPatch>
1308 template<class Type, class CombineOp>
1312  const Field<Type>& fld,
1313  const CombineOp& cop,
1314  const UList<Type>& defaultValues
1315 ) const
1316 {
1317  tmp<Field<Type>> tresult
1318  (
1319  new Field<Type>
1320  (
1321  tgtAddress_.size(),
1322  Zero
1323  )
1324  );
1325 
1326  interpolateToTarget
1327  (
1328  fld,
1330  tresult.ref(),
1331  defaultValues
1332  );
1333 
1334  return tresult;
1335 }
1336 
1337 
1338 template<class SourcePatch, class TargetPatch>
1339 template<class Type, class CombineOp>
1343  const tmp<Field<Type>>& tFld,
1344  const CombineOp& cop,
1345  const UList<Type>& defaultValues
1346 ) const
1347 {
1348  return interpolateToTarget(tFld(), cop, defaultValues);
1349 }
1350 
1351 
1352 template<class SourcePatch, class TargetPatch>
1353 template<class Type>
1357  const Field<Type>& fld,
1358  const UList<Type>& defaultValues
1359 ) const
1360 {
1361  return interpolateToSource(fld, plusEqOp<Type>(), defaultValues);
1362 }
1363 
1364 
1365 template<class SourcePatch, class TargetPatch>
1366 template<class Type>
1370  const tmp<Field<Type>>& tFld,
1371  const UList<Type>& defaultValues
1372 ) const
1373 {
1374  return interpolateToSource(tFld(), plusEqOp<Type>(), defaultValues);
1375 }
1376 
1377 
1378 template<class SourcePatch, class TargetPatch>
1379 template<class Type>
1383  const Field<Type>& fld,
1384  const UList<Type>& defaultValues
1385 ) const
1386 {
1387  return interpolateToTarget(fld, plusEqOp<Type>(), defaultValues);
1388 }
1389 
1390 
1391 template<class SourcePatch, class TargetPatch>
1392 template<class Type>
1396  const tmp<Field<Type>>& tFld,
1397  const UList<Type>& defaultValues
1398 ) const
1399 {
1400  return interpolateToTarget(tFld(), plusEqOp<Type>(), defaultValues);
1401 }
1402 
1403 
1404 template<class SourcePatch, class TargetPatch>
1407  const SourcePatch& srcPatch,
1408  const TargetPatch& tgtPatch,
1409  const vector& n,
1410  const label tgtFacei,
1411  point& tgtPoint
1412 )
1413 const
1414 {
1415  const pointField& srcPoints = srcPatch.points();
1416 
1417  // Source face addresses that intersect target face tgtFacei
1418  const labelList& addr = tgtAddress_[tgtFacei];
1419 
1420  pointHit nearest;
1421  nearest.setDistance(GREAT);
1422  label nearestFacei = -1;
1423 
1424  forAll(addr, i)
1425  {
1426  const label srcFacei = addr[i];
1427  const face& f = srcPatch[srcFacei];
1428 
1429  pointHit ray = f.ray(tgtPoint, n, srcPoints);
1430 
1431  if (ray.hit())
1432  {
1433  // tgtPoint = ray.rawPoint();
1434  return srcFacei;
1435  }
1436  else if (ray.distance() < nearest.distance())
1437  {
1438  nearest = ray;
1439  nearestFacei = srcFacei;
1440  }
1441  }
1442 
1443  if (nearest.hit() || nearest.eligibleMiss())
1444  {
1445  // tgtPoint = nearest.rawPoint();
1446  return nearestFacei;
1447  }
1448 
1449  return -1;
1450 }
1451 
1452 
1453 template<class SourcePatch, class TargetPatch>
1456  const SourcePatch& srcPatch,
1457  const TargetPatch& tgtPatch,
1458  const vector& n,
1459  const label srcFacei,
1460  point& srcPoint
1461 )
1462 const
1463 {
1464  const pointField& tgtPoints = tgtPatch.points();
1465 
1466  pointHit nearest;
1467  nearest.setDistance(GREAT);
1468  label nearestFacei = -1;
1469 
1470  // Target face addresses that intersect source face srcFacei
1471  const labelList& addr = srcAddress_[srcFacei];
1472 
1473  forAll(addr, i)
1474  {
1475  const label tgtFacei = addr[i];
1476  const face& f = tgtPatch[tgtFacei];
1477 
1478  pointHit ray = f.ray(srcPoint, n, tgtPoints);
1479 
1480  if (ray.hit() || ray.eligibleMiss())
1481  {
1482  // srcPoint = ray.rawPoint();
1483  return tgtFacei;
1484  }
1485  else if (ray.distance() < nearest.distance())
1486  {
1487  nearest = ray;
1488  nearestFacei = tgtFacei;
1489  }
1490  }
1491 
1492  if (nearest.hit() || nearest.eligibleMiss())
1493  {
1494  // srcPoint = nearest.rawPoint();
1495  return nearestFacei;
1496  }
1497 
1498  return -1;
1499 }
1500 
1501 
1502 template<class SourcePatch, class TargetPatch>
1505  const SourcePatch& srcPatch,
1506  const TargetPatch& tgtPatch,
1507  const labelListList& srcAddress
1508 )
1509 const
1510 {
1511  OFstream os("faceConnectivity" + Foam::name(Pstream::myProcNo()) + ".obj");
1512 
1513  label ptI = 1;
1514 
1515  forAll(srcAddress, i)
1516  {
1517  const labelList& addr = srcAddress[i];
1518  const point& srcPt = srcPatch.faceCentres()[i];
1519 
1520  forAll(addr, j)
1521  {
1522  label tgtPtI = addr[j];
1523  const point& tgtPt = tgtPatch.faceCentres()[tgtPtI];
1524 
1525  meshTools::writeOBJ(os, srcPt);
1526  meshTools::writeOBJ(os, tgtPt);
1527 
1528  os << "l " << ptI << " " << ptI + 1 << endl;
1529 
1530  ptI += 2;
1531  }
1532  }
1533 }
1534 
1535 
1536 // ************************************************************************* //
static word interpolationMethodToWord(const interpolationMethod &method)
Convert interpolationMethod to word representation.
Class containing functor to negate primitives. Dummy for all other types.
Definition: flipOp.H:50
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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 & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
const scalarListList & tgtWeights() const
Return const access to target patch weights.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Type gMin(const FieldField< Field, Type > &f)
Xfer< List< T > > xfer()
Transfer contents to the Xfer container.
Definition: ListI.H:177
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute data using default commsType.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
Output to file stream.
Definition: OFstream.H:82
const labelListList & subMap() const
From subsetted data back to original data.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
static interpolationMethod wordTointerpolationMethod(const word &method)
Convert word to interpolationMethod.
const labelListList & tgtAddress() const
Return const access to target patch addressing.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
label tgtPointFace(const SourcePatch &srcPatch, const TargetPatch &tgtPatch, const vector &n, const label srcFacei, point &srcPoint) const
Return target patch face index of point on source patch face.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
label srcPointFace(const SourcePatch &srcPatch, const TargetPatch &tgtPatch, const vector &n, const label tgtFacei, point &tgtPoint) const
Return source patch face index of point on target patch face.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
interpolationMethod
Enumeration specifying interpolation method.
void writeFaceConnectivity(const SourcePatch &srcPatch, const TargetPatch &tgtPatch, const labelListList &srcAddress) const
Write face connectivity as OBJ file.
const scalarField & srcMagSf() const
Return const access to source patch face areas.
const Point & hitPoint() const
Return hit point.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
A List obtained as a section of another List.
Definition: SubList.H:53
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
pointHit ray(const point &p, const vector &n, const pointField &, const intersection::algorithm alg=intersection::FULL_RAY, const intersection::direction dir=intersection::VECTOR) const
Return potential intersection with face with a ray starting.
void reset(T *=0)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
void setDistance(const scalar d)
Definition: PointHit.H:186
Type gSum(const FieldField< Field, Type > &f)
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const =0
A class for handling words, derived from string.
Definition: word.H:59
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:184
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
bool hit() const
Is there a hit.
bool hit() const
Is there a hit.
Definition: PointHit.H:120
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
~AMIInterpolation()
Destructor.
void interpolateToTarget(const UList< Type > &fld, const CombineOp &cop, List< Type > &result, const UList< Type > &defaultValues=UList< Type >::null()) const
Interpolate from source to target with supplied op.
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Base class for Arbitrary Mesh Interface (AMI) methods.
Definition: AMIMethod.H:55
volScalarField scalarField(fieldObject, mesh)
const scalarListList & srcWeights() const
Return const access to source patch weights.
static const char nl
Definition: Ostream.H:262
Type gMax(const FieldField< Field, Type > &f)
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void interpolateToSource(const UList< Type > &fld, const CombineOp &cop, List< Type > &result, const UList< Type > &defaultValues=UList< Type >::null()) const
Interpolate from target to source with supplied op.
label constructSize() const
Constructed data size.
void setSize(const label)
Reset size of List.
Definition: List.C:281
Class containing processor-to-processor mapping information.
void update(const SourcePatch &srcPatch, const TargetPatch &tgtPatch)
Update addressing and weights.
Type gAverage(const FieldField< Field, Type > &f)
Input from memory buffer stream.
Definition: IStringStream.H:49
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
messageStream Info
const labelListList & constructMap() const
From subsetted data to new reconstructed data.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Helper class for list to append y onto the end of x.
Definition: ListOps.H:259
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
A class for managing temporary objects.
Definition: PtrList.H:53
bool eligibleMiss() const
Is this an eligible miss.
Definition: PointHit.H:164
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
const labelListList & srcAddress() const
Return const access to source patch addressing.
static const label labelMin
Definition: label.H:61
const scalarField & tgtMagSf() const
Return const access to target patch face areas.
Tuple2< scalar, label > nearInfo
Private class for finding nearest.