1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright held by original author
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "PointEdgeWave.H"
29 #include "processorPolyPatch.H"
30 #include "cyclicPolyPatch.H"
33 #include "PstreamCombineReduceOps.H"
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 Foam::scalar Foam::PointEdgeWave<Type>::propagationTol_ = 0.01;
43 // Offset labelList. Used for transferring from one cyclic half to the other.
45 void Foam::PointEdgeWave<Type>::offset(const label val, labelList& elems)
54 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
56 // Gets point-point correspondence. Is
57 // - list of halfA points (in cyclic patch points)
58 // - list of halfB points (can overlap with A!)
59 // - for every patchPoint its corresponding point
61 void Foam::PointEdgeWave<Type>::calcCyclicAddressing()
65 forAll(mesh_.boundaryMesh(), patchI)
67 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
69 if (isA<cyclicPolyPatch>(patch))
71 label halfSize = patch.size()/2;
73 SubList<face> halfAFaces
83 new primitivePatch(halfAFaces, mesh_.points())
86 SubList<face> halfBFaces
90 patch.start() + halfSize
96 new primitivePatch(halfBFaces, mesh_.points())
103 // Handle leaving domain. Implementation referred to Type
104 template <class Type>
105 void Foam::PointEdgeWave<Type>::leaveDomain
107 const polyPatch& meshPatch,
108 const primitivePatch& patch,
109 const labelList& patchPointLabels,
110 List<Type>& pointInfo
113 const labelList& meshPoints = patch.meshPoints();
115 forAll(patchPointLabels, i)
117 label patchPointI = patchPointLabels[i];
119 const point& pt = patch.points()[meshPoints[patchPointI]];
121 pointInfo[i].leaveDomain(meshPatch, patchPointI, pt);
126 // Handle entering domain. Implementation referred to Type
127 template <class Type>
128 void Foam::PointEdgeWave<Type>::enterDomain
130 const polyPatch& meshPatch,
131 const primitivePatch& patch,
132 const labelList& patchPointLabels,
133 List<Type>& pointInfo
136 const labelList& meshPoints = patch.meshPoints();
138 forAll(patchPointLabels, i)
140 label patchPointI = patchPointLabels[i];
142 const point& pt = patch.points()[meshPoints[patchPointI]];
144 pointInfo[i].enterDomain(meshPatch, patchPointI, pt);
149 // Transform. Implementation referred to Type
150 template <class Type>
151 void Foam::PointEdgeWave<Type>::transform
153 const tensorField& rotTensor,
154 List<Type>& pointInfo
157 if (rotTensor.size() == 1)
159 const tensor& T = rotTensor[0];
163 pointInfo[i].transform(T);
170 "PointEdgeWave<Type>::transform(const tensorField&, List<Type>&)"
171 ) << "Parallel cyclics not supported" << abort(FatalError);
175 pointInfo[i].transform(rotTensor[i]);
181 // Update info for pointI, at position pt, with information from
182 // neighbouring edge.
184 // - changedPoint_, changedPoints_, nChangedPoints_,
185 // - statistics: nEvals_, nUnvisitedPoints_
186 template <class Type>
187 bool Foam::PointEdgeWave<Type>::updatePoint
190 const label neighbourEdgeI,
191 const Type& neighbourInfo,
198 bool wasValid = pointInfo.valid();
201 pointInfo.updatePoint
212 if (!changedPoint_[pointI])
214 changedPoint_[pointI] = true;
215 changedPoints_[nChangedPoints_++] = pointI;
219 if (!wasValid && pointInfo.valid())
228 // Update info for pointI, at position pt, with information from
231 // - changedPoint_, changedPoints_, nChangedPoints_,
232 // - statistics: nEvals_, nUnvisitedPoints_
233 template <class Type>
234 bool Foam::PointEdgeWave<Type>::updatePoint
237 const Type& neighbourInfo,
244 bool wasValid = pointInfo.valid();
247 pointInfo.updatePoint
257 if (!changedPoint_[pointI])
259 changedPoint_[pointI] = true;
260 changedPoints_[nChangedPoints_++] = pointI;
264 if (!wasValid && pointInfo.valid())
273 // Update info for edgeI, at position pt, with information from
274 // neighbouring point.
276 // - changedEdge_, changedEdges_, nChangedEdges_,
277 // - statistics: nEvals_, nUnvisitedEdge_
278 template <class Type>
279 bool Foam::PointEdgeWave<Type>::updateEdge
282 const label neighbourPointI,
283 const Type& neighbourInfo,
290 bool wasValid = edgeInfo.valid();
304 if (!changedEdge_[edgeI])
306 changedEdge_[edgeI] = true;
307 changedEdges_[nChangedEdges_++] = edgeI;
311 if (!wasValid && edgeInfo.valid())
320 // Check if patches of given type name are present
321 template <class Type>
322 template <class PatchType>
323 Foam::label Foam::PointEdgeWave<Type>::countPatchType() const
327 forAll(mesh_.boundaryMesh(), patchI)
329 if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
338 // Collect changed patch points
339 template <class Type>
340 void Foam::PointEdgeWave<Type>::getChangedPatchPoints
342 const primitivePatch& patch,
344 DynamicList<Type>& patchInfo,
345 DynamicList<label>& patchPoints,
346 DynamicList<label>& owner,
347 DynamicList<label>& ownerIndex
350 const labelList& meshPoints = patch.meshPoints();
351 const faceList& localFaces = patch.localFaces();
352 const labelListList& pointFaces = patch.pointFaces();
354 forAll(meshPoints, patchPointI)
356 label meshPointI = meshPoints[patchPointI];
358 if (changedPoint_[meshPointI])
360 patchInfo.append(allPointInfo_[meshPointI]);
361 patchPoints.append(patchPointI);
363 label patchFaceI = pointFaces[patchPointI][0];
365 const face& f = localFaces[patchFaceI];
367 label index = findIndex(f, patchPointI);
369 owner.append(patchFaceI);
370 ownerIndex.append(index);
375 patchPoints.shrink();
381 // Update overall for changed patch points
382 template <class Type>
383 void Foam::PointEdgeWave<Type>::updateFromPatchInfo
385 const polyPatch& meshPatch,
386 const primitivePatch& patch,
387 const labelList& owner,
388 const labelList& ownerIndex,
389 List<Type>& patchInfo
392 const faceList& localFaces = patch.localFaces();
393 const labelList& meshPoints = patch.meshPoints();
395 // Get patch and mesh points.
396 labelList changedPatchPoints(patchInfo.size());
397 labelList changedMeshPoints(patchInfo.size());
401 label faceI = owner[i];
403 const face& f = localFaces[faceI];
405 label index = (f.size() - ownerIndex[i]) % f.size();
407 changedPatchPoints[i] = f[index];
408 changedMeshPoints[i] = meshPoints[f[index]];
411 // Adapt for entering domain
412 enterDomain(meshPatch, patch, changedPatchPoints, patchInfo);
414 // Merge received info
419 changedMeshPoints[i],
422 allPointInfo_[changedMeshPoints[i]]
429 // Transfer all the information to/from neighbouring processors
430 template <class Type>
431 void Foam::PointEdgeWave<Type>::handleProcPatches()
433 // 1. Send all point info on processor patches. Send as
434 // face label + offset in face.
436 forAll(mesh_.boundaryMesh(), patchI)
438 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
440 if (isA<processorPolyPatch>(patch))
442 // Get all changed points in relative addressing
444 DynamicList<Type> patchInfo(patch.nPoints());
445 DynamicList<label> patchPoints(patch.nPoints());
446 DynamicList<label> owner(patch.nPoints());
447 DynamicList<label> ownerIndex(patch.nPoints());
449 getChangedPatchPoints
458 // Adapt for leaving domain
459 leaveDomain(patch, patch, patchPoints, patchInfo);
461 const processorPolyPatch& procPatch =
462 refCast<const processorPolyPatch>(patch);
466 Pout<< "Processor patch " << patchI << ' ' << patch.name()
467 << " communicating with " << procPatch.neighbProcNo()
468 << " Sending:" << patchInfo.size() << endl;
475 procPatch.neighbProcNo()
478 toNeighbour << owner << ownerIndex << patchInfo;
485 // 2. Receive all point info on processor patches.
488 forAll(mesh_.boundaryMesh(), patchI)
490 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
492 if (isA<processorPolyPatch>(patch))
494 const processorPolyPatch& procPatch =
495 refCast<const processorPolyPatch>(patch);
497 List<Type> patchInfo;
499 labelList ownerIndex;
501 IPstream fromNeighbour
504 procPatch.neighbProcNo()
507 fromNeighbour >> owner >> ownerIndex >> patchInfo;
512 Pout<< "Processor patch " << patchI << ' ' << patch.name()
513 << " communicating with " << procPatch.neighbProcNo()
514 << " Received:" << patchInfo.size() << endl;
517 // Apply transform to received data for non-parallel planes
518 if (!procPatch.parallel())
520 transform(procPatch.reverseT(), patchInfo);
537 // 3. Handle all shared points
538 // (Note:irrespective if changed or not for now)
541 const globalMeshData& pd = mesh_.globalData();
543 List<Type> sharedData(pd.nGlobalPoints());
545 forAll(pd.sharedPointLabels(), i)
547 label meshPointI = pd.sharedPointLabels()[i];
549 // Fill my entries in the shared points
550 sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
553 // Combine on master. Reduce operator has to handle a list and call
554 // Type.updatePoint for all elements
555 combineReduce(sharedData, listUpdateOp<Type>());
557 forAll(pd.sharedPointLabels(), i)
559 label meshPointI = pd.sharedPointLabels()[i];
561 // Retrieve my entries from the shared points
565 sharedData[pd.sharedPointAddr()[i]],
567 allPointInfo_[meshPointI]
573 template <class Type>
574 void Foam::PointEdgeWave<Type>::handleCyclicPatches()
576 // 1. Send all point info on cyclic patches. Send as
577 // face label + offset in face.
581 forAll(mesh_.boundaryMesh(), patchI)
583 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
585 if (isA<cyclicPolyPatch>(patch))
587 const primitivePatch& halfA = cycHalves_[cycHalf++];
588 const primitivePatch& halfB = cycHalves_[cycHalf++];
590 // HalfA : get all changed points in relative addressing
592 DynamicList<Type> halfAInfo(halfA.nPoints());
593 DynamicList<label> halfAPoints(halfA.nPoints());
594 DynamicList<label> halfAOwner(halfA.nPoints());
595 DynamicList<label> halfAIndex(halfA.nPoints());
597 getChangedPatchPoints
606 // HalfB : get all changed points in relative addressing
608 DynamicList<Type> halfBInfo(halfB.nPoints());
609 DynamicList<label> halfBPoints(halfB.nPoints());
610 DynamicList<label> halfBOwner(halfB.nPoints());
611 DynamicList<label> halfBIndex(halfB.nPoints());
613 getChangedPatchPoints
623 // HalfA : adapt for leaving domain
624 leaveDomain(patch, halfA, halfAPoints, halfAInfo);
626 // HalfB : adapt for leaving domain
627 leaveDomain(patch, halfB, halfBPoints, halfBInfo);
630 // Apply rotation for non-parallel planes
631 const cyclicPolyPatch& cycPatch =
632 refCast<const cyclicPolyPatch>(patch);
634 if (!cycPatch.parallel())
636 // received data from half1
637 transform(cycPatch.forwardT(), halfAInfo);
639 // received data from half2
640 transform(cycPatch.reverseT(), halfBInfo);
645 Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
646 << " Changed on first half : " << halfAInfo.size()
647 << " Changed on second half : " << halfBInfo.size()
651 // Half1: update with data from halfB
661 // Half2: update with data from halfA
673 //checkCyclic(patch);
680 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
682 // Iterate, propagating changedPointsInfo across mesh, until no change (or
683 // maxIter reached). Initial point values specified.
684 template <class Type>
685 Foam::PointEdgeWave<Type>::PointEdgeWave
687 const polyMesh& mesh,
688 const labelList& changedPoints,
689 const List<Type>& changedPointsInfo,
691 List<Type>& allPointInfo,
692 List<Type>& allEdgeInfo,
698 allPointInfo_(allPointInfo),
699 allEdgeInfo_(allEdgeInfo),
700 changedPoint_(mesh_.nPoints(), false),
701 changedPoints_(mesh_.nPoints()),
703 changedEdge_(mesh_.nEdges(), false),
704 changedEdges_(mesh_.nEdges()),
706 nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
707 cycHalves_(2*nCyclicPatches_),
709 nUnvisitedPoints_(mesh_.nPoints()),
710 nUnvisitedEdges_(mesh_.nEdges())
712 if (allPointInfo_.size() != mesh_.nPoints())
716 "PointEdgeWave<Type>::PointEdgeWave"
717 "(const polyMesh&, const labelList&, const List<Type>,"
718 " List<Type>&, List<Type>&, const label maxIter)"
719 ) << "size of pointInfo work array is not equal to the number"
720 << " of points in the mesh" << endl
721 << " pointInfo :" << allPointInfo_.size() << endl
722 << " mesh.nPoints:" << mesh_.nPoints()
725 if (allEdgeInfo_.size() != mesh_.nEdges())
729 "PointEdgeWave<Type>::PointEdgeWave"
730 "(const polyMesh&, const labelList&, const List<Type>,"
731 " List<Type>&, List<Type>&, const label maxIter)"
732 ) << "size of edgeInfo work array is not equal to the number"
733 << " of edges in the mesh" << endl
734 << " edgeInfo :" << allEdgeInfo_.size() << endl
735 << " mesh.nEdges:" << mesh_.nEdges()
740 // Calculate cyclic halves addressing.
741 if (nCyclicPatches_ > 0)
743 calcCyclicAddressing();
747 // Set from initial changed points data
748 setPointInfo(changedPoints, changedPointsInfo);
752 Pout<< "Seed points : " << nChangedPoints_ << endl;
755 // Iterate until nothing changes
756 label iter = iterate(maxIter);
758 if ((maxIter > 0) && (iter >= maxIter))
762 "PointEdgeWave<Type>::PointEdgeWave"
763 "(const polyMesh&, const labelList&, const List<Type>,"
764 " List<Type>&, List<Type>&, const label maxIter)"
765 ) << "Maximum number of iterations reached. Increase maxIter." << endl
766 << " maxIter:" << maxIter << endl
767 << " nChangedPoints:" << nChangedPoints_ << endl
768 << " nChangedEdges:" << nChangedEdges_ << endl
774 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
776 template <class Type>
777 Foam::PointEdgeWave<Type>::~PointEdgeWave()
781 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
784 template <class Type>
785 Foam::label Foam::PointEdgeWave<Type>::getUnsetPoints() const
787 return nUnvisitedPoints_;
791 template <class Type>
792 Foam::label Foam::PointEdgeWave<Type>::getUnsetEdges() const
794 return nUnvisitedEdges_;
798 // Copy point information into member data
799 template <class Type>
800 void Foam::PointEdgeWave<Type>::setPointInfo
802 const labelList& changedPoints,
803 const List<Type>& changedPointsInfo
806 forAll(changedPoints, changedPointI)
808 label pointI = changedPoints[changedPointI];
810 bool wasValid = allPointInfo_[pointI].valid();
812 // Copy info for pointI
813 allPointInfo_[pointI] = changedPointsInfo[changedPointI];
815 // Maintain count of unset points
816 if (!wasValid && allPointInfo_[pointI].valid())
821 // Mark pointI as changed, both on list and on point itself.
823 if (!changedPoint_[pointI])
825 changedPoint_[pointI] = true;
826 changedPoints_[nChangedPoints_++] = pointI;
832 // Propagate information from edge to point. Return number of points changed.
833 template <class Type>
834 Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
838 label changedEdgeI = 0;
839 changedEdgeI < nChangedEdges_;
843 label edgeI = changedEdges_[changedEdgeI];
845 if (!changedEdge_[edgeI])
847 FatalErrorIn("PointEdgeWave<Type>::edgeToPoint()")
849 << " not marked as having been changed" << nl
850 << "This might be caused by multiple occurences of the same"
851 << " seed point." << abort(FatalError);
855 const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
857 // Evaluate all connected points (= edge endpoints)
858 const edge& e = mesh_.edges()[edgeI];
862 Type& currentWallInfo = allPointInfo_[e[eI]];
864 if (currentWallInfo != neighbourWallInfo)
877 // Reset status of edge
878 changedEdge_[edgeI] = false;
881 // Handled all changed edges by now
884 if (nCyclicPatches_ > 0)
886 // Transfer changed points across cyclic halves
887 handleCyclicPatches();
889 if (Pstream::parRun())
891 // Transfer changed points from neighbouring processors.
897 Pout<< "Changed points : " << nChangedPoints_ << endl;
900 // Sum nChangedPoints over all procs
901 label totNChanged = nChangedPoints_;
903 reduce(totNChanged, sumOp<label>());
909 // Propagate information from point to edge. Return number of edges changed.
910 template <class Type>
911 Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
913 const labelListList& pointEdges = mesh_.pointEdges();
917 label changedPointI = 0;
918 changedPointI < nChangedPoints_;
922 label pointI = changedPoints_[changedPointI];
924 if (!changedPoint_[pointI])
926 FatalErrorIn("PointEdgeWave<Type>::pointToEdge()")
927 << "Point " << pointI
928 << " not marked as having been changed" << nl
929 << "This might be caused by multiple occurences of the same"
930 << " seed point." << abort(FatalError);
933 const Type& neighbourWallInfo = allPointInfo_[pointI];
935 // Evaluate all connected edges
937 const labelList& edgeLabels = pointEdges[pointI];
938 forAll(edgeLabels, edgeLabelI)
940 label edgeI = edgeLabels[edgeLabelI];
942 Type& currentWallInfo = allEdgeInfo_[edgeI];
944 if (currentWallInfo != neighbourWallInfo)
957 // Reset status of point
958 changedPoint_[pointI] = false;
961 // Handled all changed points by now
966 Pout<< "Changed edges : " << nChangedEdges_ << endl;
969 // Sum nChangedPoints over all procs
970 label totNChanged = nChangedEdges_;
972 reduce(totNChanged, sumOp<label>());
979 template <class Type>
980 Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
982 if (nCyclicPatches_ > 0)
984 // Transfer changed points across cyclic halves
985 handleCyclicPatches();
987 if (Pstream::parRun())
989 // Transfer changed points from neighbouring processors.
997 while(iter < maxIter)
1001 Pout<< "Iteration " << iter << endl;
1004 label nEdges = pointToEdge();
1008 Pout<< "Total changed edges : " << nEdges << endl;
1016 label nPoints = edgeToPoint();
1020 Pout<< "Total changed points : " << nPoints << endl;
1022 Pout<< "Total evaluations : " << nEvals_ << endl;
1024 Pout<< "Remaining unvisited points: " << nUnvisitedPoints_ << endl;
1026 Pout<< "Remaining unvisited edges : " << nUnvisitedEdges_ << endl;
1042 // ************************************************************************* //