40 const Foam::scalar Foam::meshToMesh0::directHitTol = 1
e-5;
47 void Foam::meshToMesh0::calcAddressing()
52 <<
"Calculating mesh-to-mesh cell addressing" <<
endl;
75 List<bool> boundaryCell(fromCells.size(),
false);
87 boundaryCell[bCells[facei]] =
true;
91 treeBoundBox
meshBb(fromPoints);
93 scalar typDim =
meshBb.avgDim()/(2.0*
cbrt(scalar(fromCells.size())));
95 treeBoundBox shiftedBb
104 <<
" bounding box : " <<
meshBb <<
nl
105 <<
" bounding box (shifted) : " << shiftedBb <<
nl
106 <<
" typical dimension :" << shiftedBb.typDim() <<
endl;
109 indexedOctree<treeDataCell> oc
120 oc.print(
Pout,
false, 0);
136 if (cuttingPatches_.
found(toPatch.name()))
142 boundaryAddressing_[
patchi],
143 toPatch.faceCentres(),
151 patchMap_.found(toPatch.name())
152 && fromMeshPatches_.
found(patchMap_.find(toPatch.name())())
157 fromMeshPatches_.
find(patchMap_.find(toPatch.name())())()
160 if (fromPatch.empty())
163 <<
"Source patch " << fromPatch.name()
164 <<
" has no faces. Not performing mapping for it."
167 boundaryAddressing_[
patchi] = -1;
171 treeBoundBox wallBb(fromPatch.localPoints());
173 wallBb.avgDim()/(2.0*
sqrt(scalar(fromPatch.size())));
175 treeBoundBox shiftedBb
178 wallBb.max() +
vector(typDim, typDim, typDim)
183 indexedOctree<treeDataFace> oc
185 treeDataFace(
false, fromPatch),
193 toPatch.faceCentres();
197 const scalar distSqr =
sqr(wallBb.mag());
201 boundaryAddressing_[
patchi][toi] = oc.findNearest
203 centresToBoundary[toi],
214 <<
"Finished calculating mesh-to-mesh cell addressing" <<
endl;
219 void Foam::meshToMesh0::cellAddresses
223 const fvMesh& fromMesh,
224 const List<bool>& boundaryCell,
225 const indexedOctree<treeDataCell>& oc
240 const vectorField& centresFrom = fromMesh.cellCentres();
249 scalar distSqr =
magSqr(
p - centresFrom[curCell]);
258 const labelList& neighbours = cc[curCell];
262 const scalar curDistSqr =
263 magSqr(
p - centresFrom[neighbours[ni]]);
267 if (curDistSqr < (1 - small)*distSqr)
269 curCell = neighbours[ni];
270 distSqr = curDistSqr;
276 cellAddressing_[toi] = -1;
279 if (fromMesh.pointInCell(
p, curCell))
281 cellAddressing_[toi] = curCell;
288 if (boundaryCell[curCell])
290 cellAddressing_[toi] = oc.findInside(
p);
291 if (cellAddressing_[toi] != -1)
293 curCell = cellAddressing_[toi];
302 const labelList& neighbours = cc[curCell];
308 if (fromMesh.pointInCell(
p, neighbours[ni]))
310 cellAddressing_[toi] = neighbours[ni];
321 const labelList& neighbours = cc[curCell];
326 const labelList& nn = cc[neighbours[ni]];
332 if (fromMesh.pointInCell(
p, nn[ni]))
334 cellAddressing_[toi] = nn[ni];
346 cellAddressing_[toi] = oc.findInside(
p);
348 if (cellAddressing_[toi] != -1)
350 curCell = cellAddressing_[toi];
359 void Foam::meshToMesh0::calculateInverseDistanceWeights()
const
364 <<
"Calculating inverse distance weighting factors" <<
endl;
367 if (inverseDistanceWeightsPtr_)
370 <<
"weighting factors already calculated"
377 inverseDistanceWeightsPtr_ =
new scalarListList(toMesh_.nCells());
385 forAll(cellAddressing_, celli)
387 if (cellAddressing_[celli] != -1)
389 const vector& target = centreTo[celli];
390 scalar m =
mag(target - centreFrom[cellAddressing_[celli]]);
392 const labelList& neighbours = cc[cellAddressing_[celli]];
396 label directCelli = -1;
397 if (m < directHitTol || neighbours.empty())
405 scalar nm =
mag(target - centreFrom[neighbours[ni]]);
406 if (nm < directHitTol)
408 directCelli = neighbours[ni];
415 if (directCelli != -1)
418 invDistCoeffs[directCelli].setSize(1);
419 invDistCoeffs[directCelli][0] = 1.0;
420 V_ += fromMesh_.V()[cellAddressing_[directCelli]];
424 invDistCoeffs[celli].setSize(neighbours.size() + 1);
428 scalar invDist = 1.0/m;
429 invDistCoeffs[celli][0] = invDist;
430 scalar sumInvDist = invDist;
435 invDist = 1.0/
mag(target - centreFrom[neighbours[ni]]);
436 invDistCoeffs[celli][ni + 1] = invDist;
437 sumInvDist += invDist;
441 forAll(invDistCoeffs[celli], i)
443 invDistCoeffs[celli][i] /= sumInvDist;
448 invDistCoeffs[celli][0]
449 *fromMesh_.V()[cellAddressing_[celli]];
450 for (
label i = 1; i < invDistCoeffs[celli].size(); i++)
453 invDistCoeffs[celli][i]*fromMesh_.V()[neighbours[i-1]];
461 void Foam::meshToMesh0::calculateInverseVolumeWeights()
const
466 <<
"Calculating inverse volume weighting factors" <<
endl;
469 if (inverseVolumeWeightsPtr_)
472 <<
"weighting factors already calculated"
484 tetOverlapVolume overlapEngine;
488 const labelList& overlapCells = cellToCell[celli];
490 if (overlapCells.size() > 0)
492 invVolCoeffs[celli].
setSize(overlapCells.size());
496 label cellFrom = overlapCells[j];
497 treeBoundBox bbFromMesh
502 fromMesh_.cellPoints()[cellFrom]
506 scalar v = overlapEngine.cellCellOverlapVolumeMinDecomp
515 invVolCoeffs[celli][j] = v/toMesh_.V()[celli];
524 void Foam::meshToMesh0::calculateCellToCellAddressing()
const
529 <<
"Calculating cell to cell addressing" <<
endl;
532 if (cellToCellAddressingPtr_)
535 <<
"addressing already calculated"
542 tetOverlapVolume overlapEngine;
544 cellToCellAddressingPtr_ =
new labelListList(toMesh_.nCells());
551 overlapEngine.overlappingCells(fromMesh_, toMesh_, iTo);
552 if (overLapCells.size() > 0)
554 cellToCell[iTo].
setSize(overLapCells.size());
557 cellToCell[iTo][j] = overLapCells[j];
558 V_ += fromMesh_.V()[overLapCells[j]];
569 const fvMesh& meshFrom,
570 const fvMesh& meshTo,
571 const HashTable<word>& patchMap,
578 cellAddressing_(toMesh_.nCells()),
579 boundaryAddressing_(toMesh_.boundaryMesh().size()),
580 inverseDistanceWeightsPtr_(nullptr),
581 inverseVolumeWeightsPtr_(nullptr),
582 cellToCellAddressingPtr_(nullptr),
587 fromMeshPatches_.insert
589 fromMesh_.boundaryMesh()[
patchi].name(),
596 toMeshPatches_.insert
598 toMesh_.boundaryMesh()[
patchi].name(),
603 forAll(cuttingPatchNames, i)
605 if (toMeshPatches_.found(cuttingPatchNames[i]))
607 cuttingPatches_.insert
609 cuttingPatchNames[i],
610 toMeshPatches_.find(cuttingPatchNames[i])()
616 <<
"Cannot find cutting-patch " << cuttingPatchNames[i]
617 <<
" in destination mesh" <<
endl;
624 if (isA<processorPolyPatch>(toMesh_.boundaryMesh()[
patchi]))
626 cuttingPatches_.insert
628 toMesh_.boundaryMesh()[
patchi].name(),
640 const fvMesh& meshFrom,
646 cellAddressing_(toMesh_.nCells()),
647 boundaryAddressing_(toMesh_.boundaryMesh().size()),
648 inverseDistanceWeightsPtr_(nullptr),
649 inverseVolumeWeightsPtr_(nullptr),
650 cellToCellAddressingPtr_(nullptr),
655 if (fromMesh_.boundary().size() != toMesh_.boundary().size())
658 <<
"Incompatible meshes: different number of patches, "
659 <<
"fromMesh = " << fromMesh_.boundary().size()
660 <<
", toMesh = " << toMesh_.boundary().size()
668 fromMesh_.boundaryMesh()[
patchi].name()
669 != toMesh_.boundaryMesh()[
patchi].name()
673 <<
"Incompatible meshes: different patch names for patch "
675 <<
", fromMesh = " << fromMesh_.boundary()[
patchi].name()
676 <<
", toMesh = " << toMesh_.boundary()[
patchi].name()
682 fromMesh_.boundaryMesh()[
patchi].type()
683 != toMesh_.boundaryMesh()[
patchi].type()
687 <<
"Incompatible meshes: different patch types for patch "
689 <<
", fromMesh = " << fromMesh_.boundary()[
patchi].type()
690 <<
", toMesh = " << toMesh_.boundary()[
patchi].type()
694 fromMeshPatches_.insert
696 fromMesh_.boundaryMesh()[
patchi].name(),
700 toMeshPatches_.insert
702 toMesh_.boundaryMesh()[
patchi].name(),
708 toMesh_.boundaryMesh()[
patchi].name(),
709 fromMesh_.boundaryMesh()[
patchi].name()
731 if (!inverseDistanceWeightsPtr_)
733 calculateInverseDistanceWeights();
736 return *inverseDistanceWeightsPtr_;
742 if (!inverseVolumeWeightsPtr_)
744 calculateInverseVolumeWeights();
747 return *inverseVolumeWeightsPtr_;
753 if (!cellToCellAddressingPtr_)
755 calculateCellToCellAddressing();
758 return *cellToCellAddressingPtr_;
#define forAll(list, i)
Loop across all elements in list.
SubField< vector > subField
Declare type of subField.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
bool found(const Key &) const
Return true if hashedEntry is found in table.
void setSize(const label)
Reset size of List.
meshToMesh0(const fvMesh &fromMesh, const fvMesh &toMesh, const HashTable< word > &patchMap, const wordList &cuttingPatchNames)
Construct from the two meshes, the patch name map for the patches.
~meshToMesh0()
Destructor.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const pointField & points() const
Return raw points.
const vectorField & cellCentres() const
const cellList & cells() const
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
#define InfoInFunction
Report an information message using Foam::Info.
const dimensionedScalar e
Elementary charge.
List< scalarList > scalarListList
errorManipArg< error, int > exit(error &err, const int errNo=1)
List< word > wordList
A List of words.
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
List< cell > cellList
list of cells
Ostream & endl(Ostream &os)
Add newline and flush stream.
void deleteDemandDrivenData(DataType *&dataPtr)
vectorField pointField
pointField is a vectorField.
Vector< scalar > vector
A scalar version of the templated Vector.
dimensionedScalar sqrt(const dimensionedScalar &ds)
List< labelList > labelListList
A List of labelList.
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
PtrList< polyPatch > polyPatchList
container classes for polyPatch
Field< vector > vectorField
Specialisation of Field<T> for vector.
prefixOSstream Pout(cout, "Pout")
dimensionedScalar cbrt(const dimensionedScalar &ds)
UList< label > labelUList
dimensioned< scalar > magSqr(const dimensioned< Type > &)
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(1e-3))