34 template<
class SourcePatch,
class TargetPatch>
41 word method =
"unknown-interpolationMethod";
52 method =
"mapNearestAMI";
55 case imFaceAreaWeight:
57 method =
"faceAreaWeightAMI";
60 case imPartialFaceAreaWeight:
62 method =
"partialFaceAreaWeightAMI";
68 <<
"Unhandled interpolationMethod enumeration " << method
77 template<
class SourcePatch,
class TargetPatch>
94 "partialFaceAreaWeightAMI" 99 if (im ==
"directAMI")
103 else if (im ==
"mapNearestAMI")
105 method = imMapNearest;
107 else if (im ==
"faceAreaWeightAMI")
109 method = imFaceAreaWeight;
111 else if (im ==
"partialFaceAreaWeightAMI")
113 method = imPartialFaceAreaWeight;
118 <<
"Invalid interpolationMethod " << im
119 <<
". Valid methods are:" << methods
129 template<
class SourcePatch,
class TargetPatch>
138 Info<<
"AMI: projecting points to surface" <<
endl;
164 <<
"Error projecting points to surface: " 165 << nMiss <<
" faces could not be determined" 171 template<
class SourcePatch,
class TargetPatch>
175 const word& patchName,
179 const bool conformal,
181 const scalar lowWeightTol
186 label nLowWeight = 0;
194 scalar denom = patchAreas[facei];
211 if (t < lowWeightTol)
230 <<
"AMI: Patch " << patchName
231 <<
" sum(weights) min/max/average = " 232 <<
gMin(wghtSum) <<
", " 233 <<
gMax(wghtSum) <<
", " 241 <<
"AMI: Patch " << patchName
242 <<
" identified " << nLow
243 <<
" faces with weights less than " << lowWeightTol
251 template<
class SourcePatch,
class TargetPatch>
259 const labelList& sourceRestrictAddressing,
260 const labelList& targetRestrictAddressing,
269 label sourceCoarseSize =
271 sourceRestrictAddressing.
size()
272 ?
max(sourceRestrictAddressing)+1
276 label targetCoarseSize =
278 targetRestrictAddressing.
size()
279 ?
max(targetRestrictAddressing)+1
285 srcMagSf.
setSize(sourceRestrictAddressing.
size(), 0.0);
286 forAll(sourceRestrictAddressing, facei)
288 label coarseFacei = sourceRestrictAddressing[facei];
289 srcMagSf[coarseFacei] += fineSrcMagSf[facei];
295 if (targetMapPtr.
valid())
300 labelList allRestrict(targetRestrictAddressing);
312 tgtSubMap[Pstream::myProcNo()] =
identity(targetCoarseSize);
317 if (proci != Pstream::myProcNo())
327 labelList oldToNew(targetCoarseSize, -1);
332 label fineElem = elems[i];
333 label coarseElem = allRestrict[fineElem];
334 if (oldToNew[coarseElem] == -1)
336 oldToNew[coarseElem] = newI;
337 newSubMap[newI] = coarseElem;
354 tgtConstructMap[Pstream::myProcNo()] =
357 tgtCompactMap = targetRestrictAddressing;
360 label compactI = targetCoarseSize;
365 if (proci != Pstream::myProcNo())
372 labelList& newConstructMap = tgtConstructMap[proci];
382 remoteTargetCoarseSize =
max 384 remoteTargetCoarseSize,
385 allRestrict[elems[i]]
388 remoteTargetCoarseSize += 1;
391 labelList oldToNew(remoteTargetCoarseSize, -1);
396 label fineElem = elems[i];
398 label coarseElem = allRestrict[fineElem];
399 if (oldToNew[coarseElem] == -1)
401 oldToNew[coarseElem] = newI;
402 tgtCompactMap[fineElem] = compactI;
403 newConstructMap[newI] = compactI++;
409 label compactI = oldToNew[coarseElem];
410 tgtCompactMap[fineElem] = newConstructMap[compactI];
419 srcAddress.
setSize(sourceCoarseSize);
420 srcWeights.
setSize(sourceCoarseSize);
422 forAll(fineSrcAddress, facei)
426 const labelList& elems = fineSrcAddress[facei];
427 const scalarList& weights = fineSrcWeights[facei];
428 const scalar fineArea = fineSrcMagSf[facei];
430 label coarseFacei = sourceRestrictAddressing[facei];
432 labelList& newElems = srcAddress[coarseFacei];
433 scalarList& newWeights = srcWeights[coarseFacei];
437 label elemI = elems[i];
438 label coarseElemI = tgtCompactMap[elemI];
443 newElems.
append(coarseElemI);
444 newWeights.
append(fineArea*weights[i]);
448 newWeights[index] += fineArea*weights[i];
459 tgtConstructMap.
xfer()
465 srcAddress.
setSize(sourceCoarseSize);
466 srcWeights.
setSize(sourceCoarseSize);
468 forAll(fineSrcAddress, facei)
472 const labelList& elems = fineSrcAddress[facei];
473 const scalarList& weights = fineSrcWeights[facei];
474 const scalar fineArea = fineSrcMagSf[facei];
476 label coarseFacei = sourceRestrictAddressing[facei];
478 labelList& newElems = srcAddress[coarseFacei];
479 scalarList& newWeights = srcWeights[coarseFacei];
483 label elemI = elems[i];
484 label coarseElemI = targetRestrictAddressing[elemI];
489 newElems.
append(coarseElemI);
490 newWeights.
append(fineArea*weights[i]);
494 newWeights[index] += fineArea*weights[i];
515 template<
class SourcePatch,
class TargetPatch>
518 const SourcePatch& srcPatch,
519 const TargetPatch& tgtPatch,
527 SourcePatch srcPatch0
548 TargetPatch tgtPatch0
570 projectPointsToSurface(surfPtr(), srcPoints);
571 projectPointsToSurface(surfPtr(), tgtPoints);
575 update(srcPatch0, tgtPatch0);
579 update(srcPatch, tgtPatch);
586 template<
class SourcePatch,
class TargetPatch>
589 const SourcePatch& srcPatch,
590 const TargetPatch& tgtPatch,
592 const bool requireMatch,
594 const scalar lowWeightCorrection,
595 const bool reverseTarget
598 methodName_(interpolationMethodToWord(method)),
599 reverseTarget_(reverseTarget),
600 requireMatch_(requireMatch),
601 singlePatchProc_(-999),
602 lowWeightCorrection_(lowWeightCorrection),
613 update(srcPatch, tgtPatch);
617 template<
class SourcePatch,
class TargetPatch>
620 const SourcePatch& srcPatch,
621 const TargetPatch& tgtPatch,
623 const bool requireMatch,
624 const word& methodName,
625 const scalar lowWeightCorrection,
626 const bool reverseTarget
629 methodName_(methodName),
630 reverseTarget_(reverseTarget),
631 requireMatch_(requireMatch),
632 singlePatchProc_(-999),
633 lowWeightCorrection_(lowWeightCorrection),
644 update(srcPatch, tgtPatch);
648 template<
class SourcePatch,
class TargetPatch>
651 const SourcePatch& srcPatch,
652 const TargetPatch& tgtPatch,
655 const bool requireMatch,
657 const scalar lowWeightCorrection,
658 const bool reverseTarget
661 methodName_(interpolationMethodToWord(method)),
662 reverseTarget_(reverseTarget),
663 requireMatch_(requireMatch),
664 singlePatchProc_(-999),
665 lowWeightCorrection_(lowWeightCorrection),
676 constructFromSurface(srcPatch, tgtPatch, surfPtr);
680 template<
class SourcePatch,
class TargetPatch>
683 const SourcePatch& srcPatch,
684 const TargetPatch& tgtPatch,
687 const bool requireMatch,
688 const word& methodName,
689 const scalar lowWeightCorrection,
690 const bool reverseTarget
693 methodName_(methodName),
694 reverseTarget_(reverseTarget),
695 requireMatch_(requireMatch),
696 singlePatchProc_(-999),
697 lowWeightCorrection_(lowWeightCorrection),
708 constructFromSurface(srcPatch, tgtPatch, surfPtr);
712 template<
class SourcePatch,
class TargetPatch>
716 const labelList& sourceRestrictAddressing,
717 const labelList& targetRestrictAddressing
720 methodName_(fineAMI.methodName_),
721 reverseTarget_(fineAMI.reverseTarget_),
722 requireMatch_(fineAMI.requireMatch_),
723 singlePatchProc_(fineAMI.singlePatchProc_),
724 lowWeightCorrection_(-1.0),
731 triMode_(fineAMI.triMode_),
735 label sourceCoarseSize =
737 sourceRestrictAddressing.
size()
738 ?
max(sourceRestrictAddressing)+1
742 label neighbourCoarseSize =
744 targetRestrictAddressing.
size()
745 ?
max(targetRestrictAddressing)+1
751 Pout<<
"AMI: Creating addressing and weights as agglomeration of AMI :" 754 <<
" coarse source size:" << sourceCoarseSize
755 <<
" neighbour source size:" << neighbourCoarseSize
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()
786 sourceRestrictAddressing,
787 targetRestrictAddressing,
803 targetRestrictAddressing,
804 sourceRestrictAddressing,
817 template<
class SourcePatch,
class TargetPatch>
824 template<
class SourcePatch,
class TargetPatch>
827 const SourcePatch& srcPatch,
828 const TargetPatch& tgtPatch
834 if (srcTotalSize == 0)
838 Info<<
"AMI: no source faces present - no addressing constructed" 846 <<
"AMI: Creating addressing and weights between " 847 << srcTotalSize <<
" source faces and " 848 << tgtTotalSize <<
" target faces" 852 srcMagSf_.setSize(srcPatch.size());
855 srcMagSf_[facei] = srcPatch[facei].mag(srcPatch.points());
857 tgtMagSf_.setSize(tgtPatch.size());
860 tgtMagSf_[facei] = tgtPatch[facei].mag(tgtPatch.points());
864 singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
866 if (singlePatchProc_ == -1)
887 distributeAndMergePatches
920 requireMatch_ && (lowWeightCorrection_ < 0)
949 addressing[addrI] = tgtFaceIDs[addressing[addrI]];
958 addressing[addrI] = globalSrcFaces.toGlobal(addressing[addrI]);
965 mapDistributeBase::distribute
967 Pstream::commsTypes::nonBlocking,
980 mapDistributeBase::distribute
982 Pstream::commsTypes::nonBlocking,
1003 AMIPtr->conformal(),
1005 lowWeightCorrection_
1014 AMIPtr->conformal(),
1016 lowWeightCorrection_
1021 srcMapPtr_.reset(
new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
1022 tgtMapPtr_.reset(
new mapDistribute(globalTgtFaces, srcAddress_, cMap));
1026 writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
1043 requireMatch_ && (lowWeightCorrection_ < 0)
1062 AMIPtr->conformal(),
1064 lowWeightCorrection_
1073 AMIPtr->conformal(),
1075 lowWeightCorrection_
1081 Info<<
"AMIInterpolation : Constructed addressing and weights" <<
nl 1083 << faceAreaIntersect::triangulationModeNames_[triMode_] <<
nl 1084 <<
" singlePatchProc:" << singlePatchProc_ <<
nl 1085 <<
" srcMagSf :" <<
gSum(srcMagSf_) <<
nl 1086 <<
" tgtMagSf :" <<
gSum(tgtMagSf_) <<
nl 1092 template<
class SourcePatch,
class TargetPatch>
1093 template<
class Type,
class CombineOp>
1097 const CombineOp& cop,
1102 if (fld.
size() != srcAddress_.size())
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()
1112 if (lowWeightCorrection_ > 0)
1114 if (defaultValues.
size() != tgtAddress_.size())
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 1127 result.
setSize(tgtAddress_.size());
1129 if (singlePatchProc_ == -1)
1138 if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1140 result[facei] = defaultValues[facei];
1144 const labelList& faces = tgtAddress_[facei];
1145 const scalarList& weights = tgtWeights_[facei];
1149 cop(result[facei], facei, work[faces[i]], weights[i]);
1158 if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1160 result[facei] = defaultValues[facei];
1164 const labelList& faces = tgtAddress_[facei];
1165 const scalarList& weights = tgtWeights_[facei];
1169 cop(result[facei], facei, fld[faces[i]], weights[i]);
1177 template<
class SourcePatch,
class TargetPatch>
1178 template<
class Type,
class CombineOp>
1182 const CombineOp& cop,
1187 if (fld.
size() != tgtAddress_.size())
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()
1197 if (lowWeightCorrection_ > 0)
1199 if (defaultValues.
size() != srcAddress_.size())
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 1212 result.
setSize(srcAddress_.size());
1214 if (singlePatchProc_ == -1)
1223 if (srcWeightsSum_[facei] < lowWeightCorrection_)
1225 result[facei] = defaultValues[facei];
1229 const labelList& faces = srcAddress_[facei];
1230 const scalarList& weights = srcWeights_[facei];
1234 cop(result[facei], facei, work[faces[i]], weights[i]);
1243 if (srcWeightsSum_[facei] < lowWeightCorrection_)
1245 result[facei] = defaultValues[facei];
1249 const labelList& faces = srcAddress_[facei];
1250 const scalarList& weights = srcWeights_[facei];
1254 cop(result[facei], facei, fld[faces[i]], weights[i]);
1262 template<
class SourcePatch,
class TargetPatch>
1263 template<
class Type,
class CombineOp>
1268 const CombineOp& cop,
1293 template<
class SourcePatch,
class TargetPatch>
1294 template<
class Type,
class CombineOp>
1299 const CombineOp& cop,
1303 return interpolateToSource(tFld(), cop, defaultValues);
1307 template<
class SourcePatch,
class TargetPatch>
1308 template<
class Type,
class CombineOp>
1313 const CombineOp& cop,
1338 template<
class SourcePatch,
class TargetPatch>
1339 template<
class Type,
class CombineOp>
1344 const CombineOp& cop,
1348 return interpolateToTarget(tFld(), cop, defaultValues);
1352 template<
class SourcePatch,
class TargetPatch>
1353 template<
class Type>
1361 return interpolateToSource(fld,
plusEqOp<Type>(), defaultValues);
1365 template<
class SourcePatch,
class TargetPatch>
1366 template<
class Type>
1374 return interpolateToSource(tFld(),
plusEqOp<Type>(), defaultValues);
1378 template<
class SourcePatch,
class TargetPatch>
1379 template<
class Type>
1387 return interpolateToTarget(fld,
plusEqOp<Type>(), defaultValues);
1391 template<
class SourcePatch,
class TargetPatch>
1392 template<
class Type>
1400 return interpolateToTarget(tFld(),
plusEqOp<Type>(), defaultValues);
1404 template<
class SourcePatch,
class TargetPatch>
1407 const SourcePatch& srcPatch,
1408 const TargetPatch& tgtPatch,
1410 const label tgtFacei,
1415 const pointField& srcPoints = srcPatch.points();
1418 const labelList& addr = tgtAddress_[tgtFacei];
1422 label nearestFacei = -1;
1426 const label srcFacei = addr[i];
1427 const face&
f = srcPatch[srcFacei];
1439 nearestFacei = srcFacei;
1446 return nearestFacei;
1453 template<
class SourcePatch,
class TargetPatch>
1456 const SourcePatch& srcPatch,
1457 const TargetPatch& tgtPatch,
1459 const label srcFacei,
1464 const pointField& tgtPoints = tgtPatch.points();
1468 label nearestFacei = -1;
1471 const labelList& addr = srcAddress_[srcFacei];
1475 const label tgtFacei = addr[i];
1476 const face&
f = tgtPatch[tgtFacei];
1488 nearestFacei = tgtFacei;
1495 return nearestFacei;
1502 template<
class SourcePatch,
class TargetPatch>
1505 const SourcePatch& srcPatch,
1506 const TargetPatch& tgtPatch,
1518 const point& srcPt = srcPatch.faceCentres()[i];
1522 label tgtPtI = addr[j];
1523 const point& tgtPt = tgtPatch.faceCentres()[tgtPtI];
1528 os <<
"l " << ptI <<
" " << ptI + 1 <<
endl;
static word interpolationMethodToWord(const interpolationMethod &method)
Convert interpolationMethod to word representation.
Class containing functor to negate primitives. Dummy for all other types.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & indent(Ostream &os)
Indent stream.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
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.
Type gMin(const FieldField< Field, Type > &f)
Xfer< List< T > > xfer()
Transfer contents to the Xfer container.
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.
const labelListList & subMap() const
From subsetted data back to original data.
void size(const label)
Override size to be inconsistent with allocated storage.
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.
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.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
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...
A List obtained as a section of another List.
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.
void setDistance(const scalar d)
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.
void append(const T &)
Append an element at the end of the list.
List< scalar > scalarList
A List of scalars.
bool hit() const
Is there a hit.
bool hit() const
Is there a hit.
List< label > labelList
A List of labels.
~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)
errorManip< error > abort(error &err)
Base class for Arbitrary Mesh Interface (AMI) methods.
volScalarField scalarField(fieldObject, mesh)
const scalarListList & srcWeights() const
Return const access to source patch weights.
Type gMax(const FieldField< Field, Type > &f)
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.
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.
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.
prefixOSstream Pout(cout, "Pout")
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
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.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
scalar distance() const
Return distance to hit.
A class for managing temporary objects.
bool eligibleMiss() const
Is this an eligible miss.
label size() const
Return the number of elements in the UList.
const labelListList & srcAddress() const
Return const access to source patch addressing.
static const label labelMin
const scalarField & tgtMagSf() const
Return const access to target patch face areas.
Tuple2< scalar, label > nearInfo
Private class for finding nearest.