38 template<
class SourcePatch,
class TargetPatch>
46 { interpolationMethod::imDirect,
"directAMI" },
47 { interpolationMethod::imMapNearest,
"mapNearestAMI" },
48 { interpolationMethod::imFaceAreaWeight,
"faceAreaWeightAMI" },
49 { interpolationMethod::imPartialFaceAreaWeight,
"partialFaceAreaWeightAMI" }
52 template<
class SourcePatch,
class TargetPatch>
56 template<
class SourcePatch,
class TargetPatch>
74 faceAreaIntersect::triangulate
76 patch.localFaces()[patchFacei],
87 patchPoints[patchFaceTris[i][0]],
88 patchPoints[patchFaceTris[i][1]],
89 patchPoints[patchFaceTris[i][2]]
100 template<
class SourcePatch,
class TargetPatch>
107 addProfiling(ami,
"AMIInterpolation::projectPointsToSurface");
122 pts[i] =
pi.hitPoint();
134 <<
"Error projecting points to surface: "
135 << nMiss <<
" faces could not be determined"
141 template<
class SourcePatch,
class TargetPatch>
145 const word& patchName,
149 const bool conformal,
151 const scalar lowWeightTol
154 addProfiling(ami,
"AMIInterpolation::normaliseWeights");
157 wghtSum.setSize(wght.size(), 0.0);
158 label nLowWeight = 0;
166 scalar denom = patchAreas[facei];
183 if (t < lowWeightTol)
202 <<
"AMI: Patch " << patchName
204 <<
" min:" <<
gMin(wghtSum)
205 <<
" max:" <<
gMax(wghtSum)
213 <<
"AMI: Patch " << patchName
214 <<
" identified " << nLow
215 <<
" faces with weights less than " << lowWeightTol
223 template<
class SourcePatch,
class TargetPatch>
226 const autoPtr<mapDistribute>& targetMapPtr,
231 const labelList& sourceRestrictAddressing,
232 const labelList& targetRestrictAddressing,
238 autoPtr<mapDistribute>& tgtMap
243 label sourceCoarseSize =
245 sourceRestrictAddressing.size()
246 ?
max(sourceRestrictAddressing)+1
250 label targetCoarseSize =
252 targetRestrictAddressing.size()
253 ?
max(targetRestrictAddressing)+1
259 srcMagSf.setSize(sourceRestrictAddressing.size(), 0.0);
260 forAll(sourceRestrictAddressing, facei)
262 label coarseFacei = sourceRestrictAddressing[facei];
263 srcMagSf[coarseFacei] += fineSrcMagSf[facei];
269 if (targetMapPtr.valid())
271 const mapDistribute& map = targetMapPtr();
274 labelList allRestrict(targetRestrictAddressing);
275 map.distribute(allRestrict);
307 tgtSubMap[Pstream::myProcNo()] =
identity(targetCoarseSize);
310 forAll(map.subMap(), proci)
312 if (proci != Pstream::myProcNo())
322 const labelList& elems = map.subMap()[proci];
324 map.constructMap()[Pstream::myProcNo()];
326 newSubMap.
setSize(elems.size());
328 labelList oldToNew(targetCoarseSize, -1);
331 for (
const label elemi : elems)
333 label fineElem = elemsMap[elemi];
334 label coarseElem = allRestrict[fineElem];
335 if (oldToNew[coarseElem] == -1)
337 oldToNew[coarseElem] = newi;
338 newSubMap[newi] = coarseElem;
342 newSubMap.setSize(newi);
354 tgtConstructMap[Pstream::myProcNo()] =
identity(targetCoarseSize);
357 labelList tgtCompactMap(map.constructSize());
366 const labelList& elemsMap = map.constructMap()[Pstream::myProcNo()];
367 for (
const label fineElem : elemsMap)
369 label coarseElem = allRestrict[fineElem];
370 tgtCompactMap[fineElem] = coarseElem;
374 label compacti = targetCoarseSize;
377 forAll(map.constructMap(), proci)
379 if (proci != Pstream::myProcNo())
384 const labelList& elems = map.constructMap()[proci];
386 labelList& newConstructMap = tgtConstructMap[proci];
387 newConstructMap.
setSize(elems.size());
394 for (
const label elemi : elems)
396 remoteTargetCoarseSize =
max
398 remoteTargetCoarseSize,
402 remoteTargetCoarseSize += 1;
405 labelList oldToNew(remoteTargetCoarseSize, -1);
408 for (
const label fineElem : elems)
411 label coarseElem = allRestrict[fineElem];
412 if (oldToNew[coarseElem] == -1)
414 oldToNew[coarseElem] = newi;
415 tgtCompactMap[fineElem] = compacti;
416 newConstructMap[newi] = compacti++;
422 label compacti = oldToNew[coarseElem];
423 tgtCompactMap[fineElem] = newConstructMap[compacti];
426 newConstructMap.setSize(newi);
431 srcAddress.setSize(sourceCoarseSize);
432 srcWeights.setSize(sourceCoarseSize);
434 forAll(fineSrcAddress, facei)
438 const labelList& elems = fineSrcAddress[facei];
439 const scalarList& weights = fineSrcWeights[facei];
440 const scalar fineArea = fineSrcMagSf[facei];
442 label coarseFacei = sourceRestrictAddressing[facei];
444 labelList& newElems = srcAddress[coarseFacei];
445 scalarList& newWeights = srcWeights[coarseFacei];
449 label elemi = elems[i];
450 label coarseElemi = tgtCompactMap[elemi];
452 label index = newElems.find(coarseElemi);
455 newElems.append(coarseElemi);
456 newWeights.append(fineArea*weights[i]);
460 newWeights[index] += fineArea*weights[i];
470 std::move(tgtSubMap),
471 std::move(tgtConstructMap)
477 srcAddress.setSize(sourceCoarseSize);
478 srcWeights.setSize(sourceCoarseSize);
480 forAll(fineSrcAddress, facei)
484 const labelList& elems = fineSrcAddress[facei];
485 const scalarList& weights = fineSrcWeights[facei];
486 const scalar fineArea = fineSrcMagSf[facei];
488 label coarseFacei = sourceRestrictAddressing[facei];
490 labelList& newElems = srcAddress[coarseFacei];
491 scalarList& newWeights = srcWeights[coarseFacei];
495 label elemi = elems[i];
496 label coarseElemi = targetRestrictAddressing[elemi];
498 label index = newElems.find(coarseElemi);
501 newElems.append(coarseElemi);
502 newWeights.append(fineArea*weights[i]);
506 newWeights[index] += fineArea*weights[i];
527 template<
class SourcePatch,
class TargetPatch>
530 const SourcePatch& srcPatch,
531 const TargetPatch& tgtPatch,
532 const autoPtr<searchableSurface>& surfPtr
539 SourcePatch srcPatch0
552 OFstream os(
"amiSrcPoints.obj");
553 for (
const point& pt : srcPoints)
560 TargetPatch tgtPatch0
573 OFstream os(
"amiTgtPoints.obj");
574 for (
const point& pt : tgtPoints)
582 projectPointsToSurface(surfPtr(), srcPoints);
583 projectPointsToSurface(surfPtr(), tgtPoints);
587 update(srcPatch0, tgtPatch0);
591 update(srcPatch, tgtPatch);
598 template<
class SourcePatch,
class TargetPatch>
601 const SourcePatch& srcPatch,
602 const TargetPatch& tgtPatch,
604 const bool requireMatch,
606 const scalar lowWeightCorrection,
607 const bool reverseTarget
610 methodName_(interpolationMethodNames_[method]),
611 reverseTarget_(reverseTarget),
612 requireMatch_(requireMatch),
613 singlePatchProc_(-999),
614 lowWeightCorrection_(lowWeightCorrection),
627 update(srcPatch, tgtPatch);
631 template<
class SourcePatch,
class TargetPatch>
634 const SourcePatch& srcPatch,
635 const TargetPatch& tgtPatch,
637 const bool requireMatch,
638 const word& methodName,
639 const scalar lowWeightCorrection,
640 const bool reverseTarget
643 methodName_(methodName),
644 reverseTarget_(reverseTarget),
645 requireMatch_(requireMatch),
646 singlePatchProc_(-999),
647 lowWeightCorrection_(lowWeightCorrection),
660 update(srcPatch, tgtPatch);
664 template<
class SourcePatch,
class TargetPatch>
667 const SourcePatch& srcPatch,
668 const TargetPatch& tgtPatch,
671 const bool requireMatch,
673 const scalar lowWeightCorrection,
674 const bool reverseTarget
677 methodName_(interpolationMethodNames_[method]),
678 reverseTarget_(reverseTarget),
679 requireMatch_(requireMatch),
680 singlePatchProc_(-999),
681 lowWeightCorrection_(lowWeightCorrection),
694 constructFromSurface(srcPatch, tgtPatch, surfPtr);
698 template<
class SourcePatch,
class TargetPatch>
701 const SourcePatch& srcPatch,
702 const TargetPatch& tgtPatch,
705 const bool requireMatch,
706 const word& methodName,
707 const scalar lowWeightCorrection,
708 const bool reverseTarget
711 methodName_(methodName),
712 reverseTarget_(reverseTarget),
713 requireMatch_(requireMatch),
714 singlePatchProc_(-999),
715 lowWeightCorrection_(lowWeightCorrection),
728 constructFromSurface(srcPatch, tgtPatch, surfPtr);
732 template<
class SourcePatch,
class TargetPatch>
736 const labelList& sourceRestrictAddressing,
737 const labelList& targetRestrictAddressing
740 methodName_(fineAMI.methodName_),
741 reverseTarget_(fineAMI.reverseTarget_),
742 requireMatch_(fineAMI.requireMatch_),
743 singlePatchProc_(fineAMI.singlePatchProc_),
744 lowWeightCorrection_(-1.0),
753 triMode_(fineAMI.triMode_),
757 label sourceCoarseSize =
759 sourceRestrictAddressing.size()
760 ?
max(sourceRestrictAddressing)+1
764 label neighbourCoarseSize =
766 targetRestrictAddressing.size()
767 ?
max(targetRestrictAddressing)+1
773 Pout<<
"AMI: Creating addressing and weights as agglomeration of AMI :"
776 <<
" coarse source size:" << sourceCoarseSize
777 <<
" neighbour source size:" << neighbourCoarseSize
783 fineAMI.
srcAddress().size() != sourceRestrictAddressing.size()
784 || fineAMI.
tgtAddress().size() != targetRestrictAddressing.size()
788 <<
"Size mismatch." <<
nl
789 <<
"Source patch size:" << fineAMI.
srcAddress().size() <<
nl
790 <<
"Source agglomeration size:"
791 << sourceRestrictAddressing.size() <<
nl
792 <<
"Target patch size:" << fineAMI.
tgtAddress().size() <<
nl
793 <<
"Target agglomeration size:"
794 << targetRestrictAddressing.size()
808 sourceRestrictAddressing,
809 targetRestrictAddressing,
825 targetRestrictAddressing,
826 sourceRestrictAddressing,
839 template<
class SourcePatch,
class TargetPatch>
846 template<
class SourcePatch,
class TargetPatch>
849 const SourcePatch& srcPatch,
850 const TargetPatch& tgtPatch
858 if (srcTotalSize == 0)
860 DebugInfo<<
"AMI: no source faces present - no addressing constructed"
867 <<
"AMI: Creating addressing and weights between "
868 << srcTotalSize <<
" source faces and "
869 << tgtTotalSize <<
" target faces"
873 singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
875 if (singlePatchProc_ == -1)
896 distributeAndMergePatches
927 requireMatch_ && (lowWeightCorrection_ < 0)
942 AMIPtr->setMagSf(tgtPatch, map, srcMagSf_, tgtMagSf_);
954 writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
962 for (
labelList& addressing : srcAddress_)
964 for (
label& addr : addressing)
966 addr = tgtFaceIDs[addr];
970 for (
labelList& addressing : tgtAddress_)
972 globalSrcFaces.inplaceToGlobal(addressing);
978 mapDistributeBase::distribute
980 Pstream::commsTypes::nonBlocking,
993 mapDistributeBase::distribute
995 Pstream::commsTypes::nonBlocking,
1009 AMIPtr->normaliseWeights(
true, *
this);
1013 srcMapPtr_.reset(
new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
1014 tgtMapPtr_.reset(
new mapDistribute(globalTgtFaces, srcAddress_, cMap));
1028 requireMatch_ && (lowWeightCorrection_ < 0)
1040 srcMagSf_.transfer(AMIPtr->srcMagSf());
1041 tgtMagSf_.transfer(AMIPtr->tgtMagSf());
1043 AMIPtr->normaliseWeights(
true, *
this);
1048 Info<<
"AMIInterpolation : Constructed addressing and weights" <<
nl
1050 << faceAreaIntersect::triangulationModeNames_[triMode_] <<
nl
1051 <<
" singlePatchProc:" << singlePatchProc_ <<
nl
1052 <<
" srcMagSf :" <<
gSum(srcMagSf_) <<
nl
1053 <<
" tgtMagSf :" <<
gSum(tgtMagSf_) <<
nl
1059 template<
class SourcePatch,
class TargetPatch>
1062 const SourcePatch& srcPatch,
1063 const TargetPatch& tgtPatch
1078 lowWeightCorrection_,
1084 if (singlePatchProc_ == -1)
1087 labelListList& srcConstructMap = srcMapPtr_->constructMap();
1090 labelListList& tgtConstructMap = tgtMapPtr_->constructMap();
1093 labelListList& newSrcConstructMap = newPtr->srcMapPtr_->constructMap();
1096 labelListList& newTgtConstructMap = newPtr->tgtMapPtr_->constructMap();
1107 srcConstructMap[proci].size(),
1108 (mapMap.size() + newMapMap.size())
1115 newSrcConstructMap[proci].size(),
1116 (mapMap.size() + newMapMap.size())
1123 forAll(srcConstructMap[proci], srci)
1125 srcConstructMap[proci][srci] =
1126 mapMap[srcConstructMap[proci][srci]];
1132 forAll(newSrcConstructMap[proci], srci)
1134 newSrcConstructMap[proci][srci] =
1135 newMapMap[newSrcConstructMap[proci][srci]];
1139 forAll(tgtAddress_, tgti)
1141 forAll(tgtAddress_[tgti], tgtj)
1143 tgtAddress_[tgti][tgtj] =
1144 mapMap[tgtAddress_[tgti][tgtj]];
1148 forAll(newPtr->tgtAddress_, tgti)
1150 forAll(newPtr->tgtAddress_[tgti], tgtj)
1152 newPtr->tgtAddress_[tgti][tgtj] =
1153 newMapMap[newPtr->tgtAddress_[tgti][tgtj]];
1167 tgtConstructMap[proci].size(),
1168 (mapMap.size() + newMapMap.size())
1175 newTgtConstructMap[proci].size(),
1176 (mapMap.size() + newMapMap.size())
1183 forAll(tgtConstructMap[proci], tgti)
1185 tgtConstructMap[proci][tgti] =
1186 mapMap[tgtConstructMap[proci][tgti]];
1192 forAll(newTgtConstructMap[proci], tgti)
1194 newTgtConstructMap[proci][tgti] =
1195 newMapMap[newTgtConstructMap[proci][tgti]];
1199 forAll(srcAddress_, srci)
1201 forAll(srcAddress_[srci], srcj)
1203 srcAddress_[srci][srcj] =
1204 mapMap[srcAddress_[srci][srcj]];
1208 forAll(newPtr->srcAddress_, srci)
1210 forAll(newPtr->srcAddress_[srci], srcj)
1212 newPtr->srcAddress_[srci][srcj] =
1213 newMapMap[newPtr->srcAddress_[srci][srcj]];
1219 srcMapPtr_->constructSize() += newPtr->srcMapPtr_->constructSize();
1220 tgtMapPtr_->constructSize() += newPtr->tgtMapPtr_->constructSize();
1225 srcSubMap[proci].
append(newSrcSubMap[proci]);
1226 srcConstructMap[proci].
append(newSrcConstructMap[proci]);
1228 tgtSubMap[proci].
append(newTgtSubMap[proci]);
1229 tgtConstructMap[proci].
append(newTgtConstructMap[proci]);
1234 forAll(srcMagSf_, srcFacei)
1236 srcAddress_[srcFacei].append(newPtr->srcAddress()[srcFacei]);
1237 srcWeights_[srcFacei].append(newPtr->srcWeights()[srcFacei]);
1238 srcWeightsSum_[srcFacei] += newPtr->srcWeightsSum()[srcFacei];
1242 forAll(tgtMagSf_, tgtFacei)
1244 tgtAddress_[tgtFacei].append(newPtr->tgtAddress()[tgtFacei]);
1245 tgtWeights_[tgtFacei].append(newPtr->tgtWeights()[tgtFacei]);
1246 tgtWeightsSum_[tgtFacei] += newPtr->tgtWeightsSum()[tgtFacei];
1251 template<
class SourcePatch,
class TargetPatch>
1254 const bool conformal,
1267 lowWeightCorrection_
1279 lowWeightCorrection_
1284 template<
class SourcePatch,
class TargetPatch>
1285 template<
class Type,
class CombineOp>
1289 const CombineOp& cop,
1294 addProfiling(ami,
"AMIInterpolation::interpolateToTarget");
1296 if (
fld.size() != srcAddress_.size())
1299 <<
"Supplied field size is not equal to source patch size" <<
nl
1300 <<
" source patch = " << srcAddress_.size() <<
nl
1301 <<
" target patch = " << tgtAddress_.size() <<
nl
1302 <<
" supplied field = " <<
fld.size()
1306 if (lowWeightCorrection_ > 0)
1308 if (defaultValues.
size() != tgtAddress_.size())
1311 <<
"Employing default values when sum of weights falls below "
1312 << lowWeightCorrection_
1313 <<
" but supplied default field size is not equal to target "
1314 <<
"patch size" <<
nl
1315 <<
" default values = " << defaultValues.
size() <<
nl
1316 <<
" target patch = " << tgtAddress_.size() <<
nl
1321 result.
setSize(tgtAddress_.size());
1323 if (singlePatchProc_ == -1)
1332 if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1334 result[facei] = defaultValues[facei];
1338 const labelList& faces = tgtAddress_[facei];
1339 const scalarList& weights = tgtWeights_[facei];
1343 cop(result[facei], facei, work[faces[i]], weights[i]);
1352 if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1354 result[facei] = defaultValues[facei];
1358 const labelList& faces = tgtAddress_[facei];
1359 const scalarList& weights = tgtWeights_[facei];
1363 cop(result[facei], facei,
fld[faces[i]], weights[i]);
1371 template<
class SourcePatch,
class TargetPatch>
1372 template<
class Type,
class CombineOp>
1376 const CombineOp& cop,
1381 addProfiling(ami,
"AMIInterpolation::interpolateToSource");
1383 if (
fld.size() != tgtAddress_.size())
1386 <<
"Supplied field size is not equal to target patch size" <<
nl
1387 <<
" source patch = " << srcAddress_.size() <<
nl
1388 <<
" target patch = " << tgtAddress_.size() <<
nl
1389 <<
" supplied field = " <<
fld.size()
1393 if (lowWeightCorrection_ > 0)
1395 if (defaultValues.
size() != srcAddress_.size())
1398 <<
"Employing default values when sum of weights falls below "
1399 << lowWeightCorrection_
1400 <<
" but supplied default field size is not equal to target "
1401 <<
"patch size" <<
nl
1402 <<
" default values = " << defaultValues.
size() <<
nl
1403 <<
" source patch = " << srcAddress_.size() <<
nl
1408 result.
setSize(srcAddress_.size());
1410 if (singlePatchProc_ == -1)
1419 if (srcWeightsSum_[facei] < lowWeightCorrection_)
1421 result[facei] = defaultValues[facei];
1425 const labelList& faces = srcAddress_[facei];
1426 const scalarList& weights = srcWeights_[facei];
1430 cop(result[facei], facei, work[faces[i]], weights[i]);
1439 if (srcWeightsSum_[facei] < lowWeightCorrection_)
1441 result[facei] = defaultValues[facei];
1445 const labelList& faces = srcAddress_[facei];
1446 const scalarList& weights = srcWeights_[facei];
1450 cop(result[facei], facei,
fld[faces[i]], weights[i]);
1458 template<
class SourcePatch,
class TargetPatch>
1459 template<
class Type,
class CombineOp>
1464 const CombineOp& cop,
1489 template<
class SourcePatch,
class TargetPatch>
1490 template<
class Type,
class CombineOp>
1495 const CombineOp& cop,
1499 return interpolateToSource(tFld(), cop, defaultValues);
1503 template<
class SourcePatch,
class TargetPatch>
1504 template<
class Type,
class CombineOp>
1509 const CombineOp& cop,
1534 template<
class SourcePatch,
class TargetPatch>
1535 template<
class Type,
class CombineOp>
1540 const CombineOp& cop,
1544 return interpolateToTarget(tFld(), cop, defaultValues);
1548 template<
class SourcePatch,
class TargetPatch>
1549 template<
class Type>
1561 template<
class SourcePatch,
class TargetPatch>
1562 template<
class Type>
1570 return interpolateToSource(tFld(),
plusEqOp<Type>(), defaultValues);
1574 template<
class SourcePatch,
class TargetPatch>
1575 template<
class Type>
1587 template<
class SourcePatch,
class TargetPatch>
1588 template<
class Type>
1596 return interpolateToTarget(tFld(),
plusEqOp<Type>(), defaultValues);
1600 template<
class SourcePatch,
class TargetPatch>
1603 const SourcePatch& srcPatch,
1604 const TargetPatch& tgtPatch,
1606 const label tgtFacei,
1611 const pointField& srcPoints = srcPatch.points();
1614 const labelList& addr = tgtAddress_[tgtFacei];
1618 label nearestFacei = -1;
1620 for (
const label srcFacei : addr)
1622 const face&
f = srcPatch[srcFacei];
1624 pointHit ray =
f.ray(tgtPoint,
n, srcPoints);
1634 nearestFacei = srcFacei;
1641 return nearestFacei;
1648 template<
class SourcePatch,
class TargetPatch>
1651 const SourcePatch& srcPatch,
1652 const TargetPatch& tgtPatch,
1654 const label srcFacei,
1659 const pointField& tgtPoints = tgtPatch.points();
1663 label nearestFacei = -1;
1666 const labelList& addr = srcAddress_[srcFacei];
1668 for (
const label tgtFacei : addr)
1670 const face&
f = tgtPatch[tgtFacei];
1672 pointHit ray =
f.ray(srcPoint,
n, tgtPoints);
1682 nearestFacei = tgtFacei;
1689 return nearestFacei;
1696 template<
class SourcePatch,
class TargetPatch>
1699 const SourcePatch& srcPatch,
1700 const TargetPatch& tgtPatch,
1712 const point& srcPt = srcPatch.faceCentres()[i];
1714 for (
const label tgtPti : addr)
1716 const point& tgtPt = tgtPatch.faceCentres()[tgtPti];
1721 os <<
"l " << pti <<
" " << pti + 1 <<
endl;