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,
812 targetRestrictAddressing,
813 sourceRestrictAddressing,
835 template<
class SourcePatch,
class TargetPatch>
842 template<
class SourcePatch,
class TargetPatch>
845 const SourcePatch& srcPatch,
846 const TargetPatch& tgtPatch
852 if (srcTotalSize == 0)
856 Info<<
"AMI: no source faces present - no addressing constructed" 864 <<
"AMI: Creating addressing and weights between " 865 << srcTotalSize <<
" source faces and " 866 << tgtTotalSize <<
" target faces" 870 srcMagSf_.setSize(srcPatch.size());
873 srcMagSf_[facei] = srcPatch[facei].mag(srcPatch.points());
875 tgtMagSf_.setSize(tgtPatch.size());
878 tgtMagSf_[facei] = tgtPatch[facei].mag(tgtPatch.points());
882 singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
884 if (singlePatchProc_ == -1)
904 distributeAndMergePatches
937 requireMatch_ && (lowWeightCorrection_ < 0)
966 addressing[addrI] = tgtFaceIDs[addressing[addrI]];
975 addressing[addrI] = globalSrcFaces.toGlobal(addressing[addrI]);
982 mapDistributeBase::distribute
984 Pstream::nonBlocking,
997 mapDistributeBase::distribute
999 Pstream::nonBlocking,
1020 AMIPtr->conformal(),
1022 lowWeightCorrection_
1031 AMIPtr->conformal(),
1033 lowWeightCorrection_
1038 srcMapPtr_.reset(
new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
1039 tgtMapPtr_.reset(
new mapDistribute(globalTgtFaces, srcAddress_, cMap));
1043 writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
1060 requireMatch_ && (lowWeightCorrection_ < 0)
1079 AMIPtr->conformal(),
1081 lowWeightCorrection_
1090 AMIPtr->conformal(),
1092 lowWeightCorrection_
1098 Info<<
"AMIInterpolation : Constructed addressing and weights" <<
nl 1100 << faceAreaIntersect::triangulationModeNames_[triMode_] <<
nl 1101 <<
" singlePatchProc:" << singlePatchProc_ <<
nl 1102 <<
" srcMagSf :" <<
gSum(srcMagSf_) <<
nl 1103 <<
" tgtMagSf :" <<
gSum(tgtMagSf_) <<
nl 1109 template<
class SourcePatch,
class TargetPatch>
1110 template<
class Type,
class CombineOp>
1114 const CombineOp& cop,
1119 if (fld.
size() != srcAddress_.size())
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()
1129 if (lowWeightCorrection_ > 0)
1131 if (defaultValues.
size() != tgtAddress_.size())
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 1144 result.
setSize(tgtAddress_.size());
1146 if (singlePatchProc_ == -1)
1155 if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1157 result[facei] = defaultValues[facei];
1161 const labelList& faces = tgtAddress_[facei];
1162 const scalarList& weights = tgtWeights_[facei];
1166 cop(result[facei], facei, work[faces[i]], weights[i]);
1175 if (tgtWeightsSum_[facei] < lowWeightCorrection_)
1177 result[facei] = defaultValues[facei];
1181 const labelList& faces = tgtAddress_[facei];
1182 const scalarList& weights = tgtWeights_[facei];
1186 cop(result[facei], facei, fld[faces[i]], weights[i]);
1194 template<
class SourcePatch,
class TargetPatch>
1195 template<
class Type,
class CombineOp>
1199 const CombineOp& cop,
1204 if (fld.
size() != tgtAddress_.size())
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()
1214 if (lowWeightCorrection_ > 0)
1216 if (defaultValues.
size() != srcAddress_.size())
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 1229 result.
setSize(srcAddress_.size());
1231 if (singlePatchProc_ == -1)
1240 if (srcWeightsSum_[facei] < lowWeightCorrection_)
1242 result[facei] = defaultValues[facei];
1246 const labelList& faces = srcAddress_[facei];
1247 const scalarList& weights = srcWeights_[facei];
1251 cop(result[facei], facei, work[faces[i]], weights[i]);
1260 if (srcWeightsSum_[facei] < lowWeightCorrection_)
1262 result[facei] = defaultValues[facei];
1266 const labelList& faces = srcAddress_[facei];
1267 const scalarList& weights = srcWeights_[facei];
1271 cop(result[facei], facei, fld[faces[i]], weights[i]);
1279 template<
class SourcePatch,
class TargetPatch>
1280 template<
class Type,
class CombineOp>
1285 const CombineOp& cop,
1310 template<
class SourcePatch,
class TargetPatch>
1311 template<
class Type,
class CombineOp>
1316 const CombineOp& cop,
1320 return interpolateToSource(tFld(), cop, defaultValues);
1324 template<
class SourcePatch,
class TargetPatch>
1325 template<
class Type,
class CombineOp>
1330 const CombineOp& cop,
1355 template<
class SourcePatch,
class TargetPatch>
1356 template<
class Type,
class CombineOp>
1361 const CombineOp& cop,
1365 return interpolateToTarget(tFld(), cop, defaultValues);
1369 template<
class SourcePatch,
class TargetPatch>
1370 template<
class Type>
1378 return interpolateToSource(fld,
plusEqOp<Type>(), defaultValues);
1382 template<
class SourcePatch,
class TargetPatch>
1383 template<
class Type>
1391 return interpolateToSource(tFld(),
plusEqOp<Type>(), defaultValues);
1395 template<
class SourcePatch,
class TargetPatch>
1396 template<
class Type>
1404 return interpolateToTarget(fld,
plusEqOp<Type>(), defaultValues);
1408 template<
class SourcePatch,
class TargetPatch>
1409 template<
class Type>
1417 return interpolateToTarget(tFld(),
plusEqOp<Type>(), defaultValues);
1421 template<
class SourcePatch,
class TargetPatch>
1424 const SourcePatch& srcPatch,
1425 const TargetPatch& tgtPatch,
1427 const label tgtFacei,
1432 const pointField& srcPoints = srcPatch.points();
1435 const labelList& addr = tgtAddress_[tgtFacei];
1439 label srcFacei = addr[i];
1440 const face&
f = srcPatch[srcFacei];
1455 label srcFacei = addr[i];
1456 const face&
f = srcPatch[srcFacei];
1458 vector nFace(-srcPatch.faceNormals()[srcFacei]);
1459 nFace += tgtPatch.faceNormals()[tgtFacei];
1475 template<
class SourcePatch,
class TargetPatch>
1478 const SourcePatch& srcPatch,
1479 const TargetPatch& tgtPatch,
1481 const label srcFacei,
1486 const pointField& tgtPoints = tgtPatch.points();
1489 const labelList& addr = srcAddress_[srcFacei];
1493 label tgtFacei = addr[i];
1494 const face&
f = tgtPatch[tgtFacei];
1509 label tgtFacei = addr[i];
1510 const face&
f = tgtPatch[tgtFacei];
1512 vector nFace(-srcPatch.faceNormals()[srcFacei]);
1513 nFace += tgtPatch.faceNormals()[tgtFacei];
1529 template<
class SourcePatch,
class TargetPatch>
1532 const SourcePatch& srcPatch,
1533 const TargetPatch& tgtPatch,
1545 const point& srcPt = srcPatch.faceCentres()[i];
1548 label tgtPtI = addr[j];
1549 const point& tgtPt = tgtPatch.faceCentres()[tgtPtI];
1554 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)
const labelListList & constructMap() const
From subsetted data to new reconstructed data.
A face is a list of labels corresponding to mesh vertices.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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.
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.
void size(const label)
Override size to be inconsistent with allocated storage.
static interpolationMethod wordTointerpolationMethod(const word &method)
Convert word to interpolationMethod.
const Point & rawPoint() const
Return point with no checking.
Ostream & endl(Ostream &os)
Add newline and flush stream.
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...
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...
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))
void reset(T *=0)
If object pointer already set, delete object and set to given.
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.
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.
~AMIInterpolation()
Destructor.
const labelListList & tgtAddress() const
Return const access to target patch addressing.
errorManip< error > abort(error &err)
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Base class for Arbitrary Mesh Interface (AMI) methods.
volScalarField scalarField(fieldObject, mesh)
prefixOSstream Pout(cout,"Pout")
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,.
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.
label size() const
Return the number of elements in the UList.
void setSize(const label)
Reset size of List.
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.
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.
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...
bool hit() const
Is there a hit.
A class for managing temporary objects.
T & ref() const
Return non-const reference or generate a fatal error.
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
Tuple2< scalar, label > nearInfo
Private class for finding nearest.