35 void Foam::starMesh::mergeCoupleFacePoints()
67 const faceList& curFaces = cellFaces_[celli];
72 label nPointsInCell = 0;
74 scalar pointMergeTol = great;
78 nPointsInCell += curFaces[facei].size();
80 edgeList curEdges = curFaces[facei].edges();
84 scalar length = curEdges[edgeI].mag(points_);
86 if (length < pointMergeTol)
88 pointMergeTol = length;
94 pointMergeTol /= 100.0;
98 label nAddedPoints = 0;
102 const face&
f = curFaces[facei];
106 cellPoints[nAddedPoints] =
f[fI];
112 for (
label firstPI = 0; firstPI < cellPoints.size() - 1; firstPI++)
116 label otherPI = firstPI;
117 otherPI < cellPoints.size();
121 if (cellPoints[firstPI] != cellPoints[otherPI])
123 label a = cellPoints[firstPI];
124 label b = cellPoints[otherPI];
126 if (edge (a,
b).
mag(points_) < pointMergeTol)
130 Info<<
"Merging points " << a <<
" and " <<
b <<
endl;
134 label mergeSetA = -1;
135 label mergeSetB = -1;
137 if (renumberPoints[a] > -1)
139 mergeSetA = renumberPoints[a];
142 if (renumberPoints[
b] > -1)
144 mergeSetB = renumberPoints[
b];
147 if (mergeSetA == -1 && mergeSetB == -1)
151 Info<<
"adding new merge group " << nMergeSets
156 renumberPoints[a] = nMergeSets;
157 renumberPoints[
b] = nMergeSets;
159 mergeSetID[nMergeSets] =
min(a,
b);
162 if (nMergeSets >= mergeSetID.size())
164 Info<<
"Resizing mergeSetID" <<
endl;
167 (mergeSetID.size() + mergeIncrement);
170 else if (mergeSetA == -1 && mergeSetB != -1)
173 Info<<
"adding point a into the merge set of b. "
174 <<
"a: " << a <<
endl;
178 renumberPoints[a] = mergeSetB;
181 mergeSetID[mergeSetB] =
182 min(a, mergeSetID[mergeSetB]);
184 else if (mergeSetA != -1 && mergeSetB == -1)
187 Info<<
"adding point b into the merge set of a. "
188 <<
"b: " <<
b <<
endl;
192 renumberPoints[
b] = mergeSetA;
195 mergeSetID[mergeSetA] =
196 min(
b, mergeSetID[mergeSetA]);
198 else if (mergeSetA != mergeSetB)
202 label minMerge =
min(mergeSetA, mergeSetB);
203 label maxMerge =
max(mergeSetA, mergeSetB);
206 Info<<
"Points already belong to two "
207 <<
"different merge sets. "
208 <<
"Eliminate the higher merge set. Sets: "
209 << minMerge <<
" and " << maxMerge <<
endl;
212 forAll(renumberPoints, elimI)
214 if (renumberPoints[elimI] == maxMerge)
216 renumberPoints[elimI] = minMerge;
221 mergeSetID[minMerge] =
222 min(mergeSetID[minMerge], mergeSetID[maxMerge]);
230 mergeSetID.setSize(nMergeSets);
232 Info<<
"Finished creating merge sets. Number of merge sets: "
233 << nMergeSets <<
"." <<
endl;
237 forAll(renumberPoints, pointi)
239 if (renumberPoints[pointi] < 0)
242 renumberPoints[pointi] = pointi;
246 renumberPoints[pointi] = mergeSetID[renumberPoints[pointi]];
257 faceList& prelimFaces = cellFaces_[celli];
259 forAll(prelimFaces, facei)
261 face oldFacePoints = prelimFaces[facei];
263 face& prelimFacePoints = prelimFaces[facei];
265 forAll(prelimFacePoints, pointi)
267 if (renumberPoints[oldFacePoints[pointi]] < 0)
270 <<
"Error in point renumbering. Old face: "
271 << oldFacePoints <<
endl
272 <<
"prelim face: " << prelimFacePoints
276 prelimFacePoints[pointi] =
277 renumberPoints[oldFacePoints[pointi]];
288 const faceList& curFaces = cellFaces_[celli];
292 const face& curFacePoints = curFaces[facei];
294 forAll(curFacePoints, pointi)
296 renumberPoints[curFacePoints[pointi]]++;
301 forAll(cellShapes_, celli)
303 const labelList& curLabels = cellShapes_[celli];
307 if (renumberPoints[curLabels[pointi]] == 0)
310 <<
"Error in point merging for cell "
311 << celli <<
". STAR index: " << starCellIndex_[celli]
313 <<
"Point index: " << curLabels[pointi] <<
" STAR index "
314 << starPointIndex_[curLabels[pointi]] <<
endl
315 <<
"Please check the geometry for the cell." <<
endl;
320 label nUsedPoints = 0;
322 forAll(renumberPoints, pointi)
324 if (renumberPoints[pointi] > 0)
330 renumberPoints[pointi] = nUsedPoints;
331 points_[nUsedPoints] = points_[pointi];
337 renumberPoints[pointi] = -1;
341 Info<<
"Total number of points: " << points_.
size() <<
endl
342 <<
"Number of used points: " << nUsedPoints <<
endl;
347 Info<<
"Renumbering all faces" <<
endl;
351 faceList& newFaces = cellFaces_[celli];
355 face oldFacePoints = newFaces[facei];
357 face& newFacePoints = newFaces[facei];
359 forAll(newFacePoints, pointi)
361 if (renumberPoints[oldFacePoints[pointi]] < 0)
364 <<
"Error in point renumbering for point "
365 << oldFacePoints[pointi]
366 <<
". Renumbering index is -1." <<
endl
367 <<
"Old face: " << oldFacePoints <<
endl
371 newFacePoints[pointi] = renumberPoints[oldFacePoints[pointi]];
376 Info<<
"Renumbering all cell shapes" <<
endl;
378 forAll(cellShapes_, celli)
380 labelList oldLabels = cellShapes_[celli];
382 labelList& curLabels = cellShapes_[celli];
386 if (renumberPoints[curLabels[pointi]] < 0)
389 <<
"Error in point renumbering for cell "
390 << celli <<
". STAR index: " << starCellIndex_[celli]
392 <<
"Point index: " << curLabels[pointi] <<
" STAR index "
393 << starPointIndex_[curLabels[pointi]] <<
" returns invalid "
394 <<
"renumbering index: "
395 << renumberPoints[curLabels[pointi]] <<
"." <<
endl
396 <<
"Old cellShape: " << oldLabels <<
endl
400 curLabels[pointi] = renumberPoints[oldLabels[pointi]];
404 Info<<
"Renumbering STAR point lookup" <<
endl;
406 labelList oldStarPointID = starPointIndex_;
408 starPointIndex_ = -1;
410 forAll(starPointIndex_, pointi)
412 if (renumberPoints[pointi] > -1)
414 starPointIndex_[renumberPoints[pointi]] = oldStarPointID[pointi];
#define forAll(list, i)
Loop across all elements in list.
void size(const label)
Override size to be inconsistent with allocated storage.
void setSize(const label)
Reset size of List.
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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.
Ostream & endl(Ostream &os)
Add newline and flush stream.
void deleteDemandDrivenData(DataType *&dataPtr)
errorManip< error > abort(error &err)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)