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-2016 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_(NULL),
611  tgtMapPtr_(NULL)
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_(NULL),
642  tgtMapPtr_(NULL)
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_(NULL),
674  tgtMapPtr_(NULL)
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_(NULL),
706  tgtMapPtr_(NULL)
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_(NULL),
733  tgtMapPtr_(NULL)
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  //if (tgtMapPtr_.valid())
797  //{
798  // Pout<< "tgtMap:" << endl;
799  // string oldPrefix = Pout.prefix();
800  // Pout.prefix() = oldPrefix + " ";
801  // tgtMapPtr_().printLayout(Pout);
802  // Pout.prefix() = oldPrefix;
803  //}
804 
805  agglomerate
806  (
807  fineAMI.srcMapPtr_,
808  fineAMI.tgtMagSf(),
809  fineAMI.tgtAddress(),
810  fineAMI.tgtWeights(),
811 
812  targetRestrictAddressing,
813  sourceRestrictAddressing,
814 
815  tgtMagSf_,
816  tgtAddress_,
817  tgtWeights_,
818  tgtWeightsSum_,
819  srcMapPtr_
820  );
821 
822  //if (srcMapPtr_.valid())
823  //{
824  // Pout<< "srcMap:" << endl;
825  // string oldPrefix = Pout.prefix();
826  // Pout.prefix() = oldPrefix + " ";
827  // srcMapPtr_().printLayout(Pout);
828  // Pout.prefix() = oldPrefix;
829  //}
830 }
831 
832 
833 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
834 
835 template<class SourcePatch, class TargetPatch>
837 {}
838 
839 
840 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
841 
842 template<class SourcePatch, class TargetPatch>
844 (
845  const SourcePatch& srcPatch,
846  const TargetPatch& tgtPatch
847 )
848 {
849  label srcTotalSize = returnReduce(srcPatch.size(), sumOp<label>());
850  label tgtTotalSize = returnReduce(tgtPatch.size(), sumOp<label>());
851 
852  if (srcTotalSize == 0)
853  {
854  if (debug)
855  {
856  Info<< "AMI: no source faces present - no addressing constructed"
857  << endl;
858  }
859 
860  return;
861  }
862 
863  Info<< indent
864  << "AMI: Creating addressing and weights between "
865  << srcTotalSize << " source faces and "
866  << tgtTotalSize << " target faces"
867  << endl;
868 
869  // Calculate face areas
870  srcMagSf_.setSize(srcPatch.size());
871  forAll(srcMagSf_, facei)
872  {
873  srcMagSf_[facei] = srcPatch[facei].mag(srcPatch.points());
874  }
875  tgtMagSf_.setSize(tgtPatch.size());
876  forAll(tgtMagSf_, facei)
877  {
878  tgtMagSf_[facei] = tgtPatch[facei].mag(tgtPatch.points());
879  }
880 
881  // Calculate if patches present on multiple processors
882  singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
883 
884  if (singlePatchProc_ == -1)
885  {
886  // convert local addressing to global addressing
887  globalIndex globalSrcFaces(srcPatch.size());
888  globalIndex globalTgtFaces(tgtPatch.size());
889 
890  // Create processor map of overlapping faces. This map gets
891  // (possibly remote) faces from the tgtPatch such that they (together)
892  // cover all of the srcPatch
893  autoPtr<mapDistribute> mapPtr = calcProcMap(srcPatch, tgtPatch);
894  const mapDistribute& map = mapPtr();
895 
896  // create new target patch that fully encompasses source patch
897 
898  // faces and points
899  faceList newTgtFaces;
900  pointField newTgtPoints;
901  // original faces from tgtPatch (in globalIndexing since might be
902  // remote)
903  labelList tgtFaceIDs;
904  distributeAndMergePatches
905  (
906  map,
907  tgtPatch,
908  globalTgtFaces,
909  newTgtFaces,
910  newTgtPoints,
911  tgtFaceIDs
912  );
913 
914  TargetPatch
915  newTgtPatch
916  (
918  (
919  newTgtFaces,
920  newTgtFaces.size()
921  ),
922  newTgtPoints
923  );
924 
925  // calculate AMI interpolation
927  (
929  (
930  methodName_,
931  srcPatch,
932  newTgtPatch,
933  srcMagSf_,
934  tgtMagSf_,
935  triMode_,
936  reverseTarget_,
937  requireMatch_ && (lowWeightCorrection_ < 0)
938  )
939  );
940 
941  AMIPtr->calculate
942  (
943  srcAddress_,
944  srcWeights_,
945  tgtAddress_,
946  tgtWeights_
947  );
948 
949  // Now
950  // ~~~
951  // srcAddress_ : per srcPatch face a list of the newTgtPatch (not
952  // tgtPatch) faces it overlaps
953  // tgtAddress_ : per newTgtPatch (not tgtPatch) face a list of the
954  // srcPatch faces it overlaps
955 
956 
957  // Rework newTgtPatch indices into globalIndices of tgtPatch
958  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
959 
960 
961  forAll(srcAddress_, i)
962  {
963  labelList& addressing = srcAddress_[i];
964  forAll(addressing, addrI)
965  {
966  addressing[addrI] = tgtFaceIDs[addressing[addrI]];
967  }
968  }
969 
970  forAll(tgtAddress_, i)
971  {
972  labelList& addressing = tgtAddress_[i];
973  forAll(addressing, addrI)
974  {
975  addressing[addrI] = globalSrcFaces.toGlobal(addressing[addrI]);
976  }
977  }
978 
979  // send data back to originating procs. Note that contributions
980  // from different processors get added (ListAppendEqOp)
981 
982  mapDistributeBase::distribute
983  (
984  Pstream::nonBlocking,
985  List<labelPair>(),
986  tgtPatch.size(),
987  map.constructMap(),
988  false, // has flip
989  map.subMap(),
990  false, // has flip
991  tgtAddress_,
993  flipOp(), // flip operation
994  labelList()
995  );
996 
997  mapDistributeBase::distribute
998  (
999  Pstream::nonBlocking,
1000  List<labelPair>(),
1001  tgtPatch.size(),
1002  map.constructMap(),
1003  false,
1004  map.subMap(),
1005  false,
1006  tgtWeights_,
1008  flipOp(),
1009  scalarList()
1010  );
1011 
1012  // weights normalisation
1013  normaliseWeights
1014  (
1015  srcMagSf_,
1016  "source",
1017  srcAddress_,
1018  srcWeights_,
1019  srcWeightsSum_,
1020  AMIPtr->conformal(),
1021  true,
1022  lowWeightCorrection_
1023  );
1024  normaliseWeights
1025  (
1026  tgtMagSf_,
1027  "target",
1028  tgtAddress_,
1029  tgtWeights_,
1030  tgtWeightsSum_,
1031  AMIPtr->conformal(),
1032  true,
1033  lowWeightCorrection_
1034  );
1035 
1036  // cache maps and reset addresses
1037  List<Map<label>> cMap;
1038  srcMapPtr_.reset(new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
1039  tgtMapPtr_.reset(new mapDistribute(globalTgtFaces, srcAddress_, cMap));
1040 
1041  if (debug)
1042  {
1043  writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
1044  }
1045  }
1046  else
1047  {
1048  // calculate AMI interpolation
1050  (
1052  (
1053  methodName_,
1054  srcPatch,
1055  tgtPatch,
1056  srcMagSf_,
1057  tgtMagSf_,
1058  triMode_,
1059  reverseTarget_,
1060  requireMatch_ && (lowWeightCorrection_ < 0)
1061  )
1062  );
1063 
1064  AMIPtr->calculate
1065  (
1066  srcAddress_,
1067  srcWeights_,
1068  tgtAddress_,
1069  tgtWeights_
1070  );
1071 
1072  normaliseWeights
1073  (
1074  srcMagSf_,
1075  "source",
1076  srcAddress_,
1077  srcWeights_,
1078  srcWeightsSum_,
1079  AMIPtr->conformal(),
1080  true,
1081  lowWeightCorrection_
1082  );
1083  normaliseWeights
1084  (
1085  tgtMagSf_,
1086  "target",
1087  tgtAddress_,
1088  tgtWeights_,
1089  tgtWeightsSum_,
1090  AMIPtr->conformal(),
1091  true,
1092  lowWeightCorrection_
1093  );
1094  }
1095 
1096  if (debug)
1097  {
1098  Info<< "AMIInterpolation : Constructed addressing and weights" << nl
1099  << " triMode :"
1100  << faceAreaIntersect::triangulationModeNames_[triMode_] << nl
1101  << " singlePatchProc:" << singlePatchProc_ << nl
1102  << " srcMagSf :" << gSum(srcMagSf_) << nl
1103  << " tgtMagSf :" << gSum(tgtMagSf_) << nl
1104  << endl;
1105  }
1106 }
1107 
1108 
1109 template<class SourcePatch, class TargetPatch>
1110 template<class Type, class CombineOp>
1113  const UList<Type>& fld,
1114  const CombineOp& cop,
1115  List<Type>& result,
1116  const UList<Type>& defaultValues
1117 ) const
1118 {
1119  if (fld.size() != srcAddress_.size())
1120  {
1122  << "Supplied field size is not equal to source patch size" << nl
1123  << " source patch = " << srcAddress_.size() << nl
1124  << " target patch = " << tgtAddress_.size() << nl
1125  << " supplied field = " << fld.size()
1126  << abort(FatalError);
1127  }
1128 
1129  if (lowWeightCorrection_ > 0)
1130  {
1131  if (defaultValues.size() != tgtAddress_.size())
1132  {
1134  << "Employing default values when sum of weights falls below "
1135  << lowWeightCorrection_
1136  << " but supplied default field size is not equal to target "
1137  << "patch size" << nl
1138  << " default values = " << defaultValues.size() << nl
1139  << " target patch = " << tgtAddress_.size() << nl
1140  << abort(FatalError);
1141  }
1142  }
1143 
1144  result.setSize(tgtAddress_.size());
1145 
1146  if (singlePatchProc_ == -1)
1147  {
1148  const mapDistribute& map = srcMapPtr_();
1149 
1150  List<Type> work(fld);
1151  map.distribute(work);
1152 
1153  forAll(result, facei)
1154  {
1155  if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1156  {
1157  result[facei] = defaultValues[facei];
1158  }
1159  else
1160  {
1161  const labelList& faces = tgtAddress_[facei];
1162  const scalarList& weights = tgtWeights_[facei];
1163 
1164  forAll(faces, i)
1165  {
1166  cop(result[facei], facei, work[faces[i]], weights[i]);
1167  }
1168  }
1169  }
1170  }
1171  else
1172  {
1173  forAll(result, facei)
1174  {
1175  if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1176  {
1177  result[facei] = defaultValues[facei];
1178  }
1179  else
1180  {
1181  const labelList& faces = tgtAddress_[facei];
1182  const scalarList& weights = tgtWeights_[facei];
1183 
1184  forAll(faces, i)
1185  {
1186  cop(result[facei], facei, fld[faces[i]], weights[i]);
1187  }
1188  }
1189  }
1190  }
1191 }
1192 
1193 
1194 template<class SourcePatch, class TargetPatch>
1195 template<class Type, class CombineOp>
1198  const UList<Type>& fld,
1199  const CombineOp& cop,
1200  List<Type>& result,
1201  const UList<Type>& defaultValues
1202 ) const
1203 {
1204  if (fld.size() != tgtAddress_.size())
1205  {
1207  << "Supplied field size is not equal to target patch size" << nl
1208  << " source patch = " << srcAddress_.size() << nl
1209  << " target patch = " << tgtAddress_.size() << nl
1210  << " supplied field = " << fld.size()
1211  << abort(FatalError);
1212  }
1213 
1214  if (lowWeightCorrection_ > 0)
1215  {
1216  if (defaultValues.size() != srcAddress_.size())
1217  {
1219  << "Employing default values when sum of weights falls below "
1220  << lowWeightCorrection_
1221  << " but supplied default field size is not equal to target "
1222  << "patch size" << nl
1223  << " default values = " << defaultValues.size() << nl
1224  << " source patch = " << srcAddress_.size() << nl
1225  << abort(FatalError);
1226  }
1227  }
1228 
1229  result.setSize(srcAddress_.size());
1230 
1231  if (singlePatchProc_ == -1)
1232  {
1233  const mapDistribute& map = tgtMapPtr_();
1234 
1235  List<Type> work(fld);
1236  map.distribute(work);
1237 
1238  forAll(result, facei)
1239  {
1240  if (srcWeightsSum_[facei] < lowWeightCorrection_)
1241  {
1242  result[facei] = defaultValues[facei];
1243  }
1244  else
1245  {
1246  const labelList& faces = srcAddress_[facei];
1247  const scalarList& weights = srcWeights_[facei];
1248 
1249  forAll(faces, i)
1250  {
1251  cop(result[facei], facei, work[faces[i]], weights[i]);
1252  }
1253  }
1254  }
1255  }
1256  else
1257  {
1258  forAll(result, facei)
1259  {
1260  if (srcWeightsSum_[facei] < lowWeightCorrection_)
1261  {
1262  result[facei] = defaultValues[facei];
1263  }
1264  else
1265  {
1266  const labelList& faces = srcAddress_[facei];
1267  const scalarList& weights = srcWeights_[facei];
1268 
1269  forAll(faces, i)
1270  {
1271  cop(result[facei], facei, fld[faces[i]], weights[i]);
1272  }
1273  }
1274  }
1275  }
1276 }
1277 
1278 
1279 template<class SourcePatch, class TargetPatch>
1280 template<class Type, class CombineOp>
1284  const Field<Type>& fld,
1285  const CombineOp& cop,
1286  const UList<Type>& defaultValues
1287 ) const
1288 {
1289  tmp<Field<Type>> tresult
1290  (
1291  new Field<Type>
1292  (
1293  srcAddress_.size(),
1294  Zero
1295  )
1296  );
1297 
1298  interpolateToSource
1299  (
1300  fld,
1302  tresult.ref(),
1303  defaultValues
1304  );
1305 
1306  return tresult;
1307 }
1308 
1309 
1310 template<class SourcePatch, class TargetPatch>
1311 template<class Type, class CombineOp>
1315  const tmp<Field<Type>>& tFld,
1316  const CombineOp& cop,
1317  const UList<Type>& defaultValues
1318 ) const
1319 {
1320  return interpolateToSource(tFld(), cop, defaultValues);
1321 }
1322 
1323 
1324 template<class SourcePatch, class TargetPatch>
1325 template<class Type, class CombineOp>
1329  const Field<Type>& fld,
1330  const CombineOp& cop,
1331  const UList<Type>& defaultValues
1332 ) const
1333 {
1334  tmp<Field<Type>> tresult
1335  (
1336  new Field<Type>
1337  (
1338  tgtAddress_.size(),
1339  Zero
1340  )
1341  );
1342 
1343  interpolateToTarget
1344  (
1345  fld,
1347  tresult.ref(),
1348  defaultValues
1349  );
1350 
1351  return tresult;
1352 }
1353 
1354 
1355 template<class SourcePatch, class TargetPatch>
1356 template<class Type, class CombineOp>
1360  const tmp<Field<Type>>& tFld,
1361  const CombineOp& cop,
1362  const UList<Type>& defaultValues
1363 ) const
1364 {
1365  return interpolateToTarget(tFld(), cop, defaultValues);
1366 }
1367 
1368 
1369 template<class SourcePatch, class TargetPatch>
1370 template<class Type>
1374  const Field<Type>& fld,
1375  const UList<Type>& defaultValues
1376 ) const
1377 {
1378  return interpolateToSource(fld, plusEqOp<Type>(), defaultValues);
1379 }
1380 
1381 
1382 template<class SourcePatch, class TargetPatch>
1383 template<class Type>
1387  const tmp<Field<Type>>& tFld,
1388  const UList<Type>& defaultValues
1389 ) const
1390 {
1391  return interpolateToSource(tFld(), plusEqOp<Type>(), defaultValues);
1392 }
1393 
1394 
1395 template<class SourcePatch, class TargetPatch>
1396 template<class Type>
1400  const Field<Type>& fld,
1401  const UList<Type>& defaultValues
1402 ) const
1403 {
1404  return interpolateToTarget(fld, plusEqOp<Type>(), defaultValues);
1405 }
1406 
1407 
1408 template<class SourcePatch, class TargetPatch>
1409 template<class Type>
1413  const tmp<Field<Type>>& tFld,
1414  const UList<Type>& defaultValues
1415 ) const
1416 {
1417  return interpolateToTarget(tFld(), plusEqOp<Type>(), defaultValues);
1418 }
1419 
1420 
1421 template<class SourcePatch, class TargetPatch>
1424  const SourcePatch& srcPatch,
1425  const TargetPatch& tgtPatch,
1426  const vector& n,
1427  const label tgtFacei,
1428  point& tgtPoint
1429 )
1430 const
1431 {
1432  const pointField& srcPoints = srcPatch.points();
1433 
1434  // source face addresses that intersect target face tgtFacei
1435  const labelList& addr = tgtAddress_[tgtFacei];
1436 
1437  forAll(addr, i)
1438  {
1439  label srcFacei = addr[i];
1440  const face& f = srcPatch[srcFacei];
1441 
1442  pointHit ray = f.ray(tgtPoint, n, srcPoints);
1443 
1444  if (ray.hit())
1445  {
1446  tgtPoint = ray.rawPoint();
1447 
1448  return srcFacei;
1449  }
1450  }
1451 
1452  // no hit registered - try with face normal instead of input normal
1453  forAll(addr, i)
1454  {
1455  label srcFacei = addr[i];
1456  const face& f = srcPatch[srcFacei];
1457 
1458  vector nFace(-srcPatch.faceNormals()[srcFacei]);
1459  nFace += tgtPatch.faceNormals()[tgtFacei];
1460 
1461  pointHit ray = f.ray(tgtPoint, nFace, srcPoints);
1462 
1463  if (ray.hit())
1464  {
1465  tgtPoint = ray.rawPoint();
1466 
1467  return srcFacei;
1468  }
1469  }
1470 
1471  return -1;
1472 }
1473 
1474 
1475 template<class SourcePatch, class TargetPatch>
1478  const SourcePatch& srcPatch,
1479  const TargetPatch& tgtPatch,
1480  const vector& n,
1481  const label srcFacei,
1482  point& srcPoint
1483 )
1484 const
1485 {
1486  const pointField& tgtPoints = tgtPatch.points();
1487 
1488  // target face addresses that intersect source face srcFacei
1489  const labelList& addr = srcAddress_[srcFacei];
1490 
1491  forAll(addr, i)
1492  {
1493  label tgtFacei = addr[i];
1494  const face& f = tgtPatch[tgtFacei];
1495 
1496  pointHit ray = f.ray(srcPoint, n, tgtPoints);
1497 
1498  if (ray.hit())
1499  {
1500  srcPoint = ray.rawPoint();
1501 
1502  return tgtFacei;
1503  }
1504  }
1505 
1506  // no hit registered - try with face normal instead of input normal
1507  forAll(addr, i)
1508  {
1509  label tgtFacei = addr[i];
1510  const face& f = tgtPatch[tgtFacei];
1511 
1512  vector nFace(-srcPatch.faceNormals()[srcFacei]);
1513  nFace += tgtPatch.faceNormals()[tgtFacei];
1514 
1515  pointHit ray = f.ray(srcPoint, n, tgtPoints);
1516 
1517  if (ray.hit())
1518  {
1519  srcPoint = ray.rawPoint();
1520 
1521  return tgtFacei;
1522  }
1523  }
1524 
1525  return -1;
1526 }
1527 
1528 
1529 template<class SourcePatch, class TargetPatch>
1532  const SourcePatch& srcPatch,
1533  const TargetPatch& tgtPatch,
1534  const labelListList& srcAddress
1535 )
1536 const
1537 {
1538  OFstream os("faceConnectivity" + Foam::name(Pstream::myProcNo()) + ".obj");
1539 
1540  label ptI = 1;
1541 
1542  forAll(srcAddress, i)
1543  {
1544  const labelList& addr = srcAddress[i];
1545  const point& srcPt = srcPatch.faceCentres()[i];
1546  forAll(addr, j)
1547  {
1548  label tgtPtI = addr[j];
1549  const point& tgtPt = tgtPatch.faceCentres()[tgtPtI];
1550 
1551  meshTools::writeOBJ(os, srcPt);
1552  meshTools::writeOBJ(os, tgtPt);
1553 
1554  os << "l " << ptI << " " << ptI + 1 << endl;
1555 
1556  ptI += 2;
1557  }
1558  }
1559 }
1560 
1561 
1562 // ************************************************************************* //
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
const labelListList & constructMap() const
From subsetted data to new reconstructed data.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const scalarListList & srcWeights() const
Return const access to source patch weights.
Type gMin(const FieldField< Field, Type > &f)
Xfer< List< T > > xfer()
Transfer contents to the Xfer container.
Definition: ListI.H:90
const scalarField & tgtMagSf() const
Return const access to target patch face areas.
bool hit() const
Is there a hit.
const labelListList & subMap() const
From subsetted data back to original data.
Output to file stream.
Definition: OFstream.H:81
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
static interpolationMethod wordTointerpolationMethod(const word &method)
Convert word to interpolationMethod.
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
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
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.
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
void writeFaceConnectivity(const SourcePatch &srcPatch, const TargetPatch &tgtPatch, const labelListList &srcAddress) const
Write face connectivity as OBJ file.
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.
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))
void reset(T *=0)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
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:97
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
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.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
~AMIInterpolation()
Destructor.
const labelListList & tgtAddress() const
Return const access to target patch addressing.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Definition: autoPtrI.H:83
Base class for Arbitrary Mesh Interface (AMI) methods.
Definition: AMIMethod.H:55
volScalarField scalarField(fieldObject, mesh)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
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,.
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute data using default commsType.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
void setSize(const label)
Reset size of List.
Definition: List.C:295
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.
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
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.
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
const Point & hitPoint() const
Return hit point.
const labelListList & srcAddress() const
Return const access to source patch addressing.
messageStream Info
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:53
bool hit() const
Is there a hit.
Definition: PointHit.H:120
A class for managing temporary objects.
Definition: PtrList.H:54
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
label constructSize() const
Constructed data size.
const scalarField & srcMagSf() const
Return const access to source patch face areas.
const scalarListList & tgtWeights() const
Return const access to target patch weights.
static const label labelMin
Definition: label.H:61
Tuple2< scalar, label > nearInfo
Private class for finding nearest.