54 void Foam::hierarchGeomDecomp::setDecompOrder()
58 if (order.size() != 3)
62 "hierarchGeomDecomp::hierarchGeomDecomp"
63 "(const dictionary& decompositionDict)",
65 ) <<
"number of characters in order (" << order <<
") != 3"
69 for (label i = 0; i < 3; i++)
75 else if (order[i] ==
'y')
79 else if (order[i] ==
'z')
87 "hierarchGeomDecomp::hierarchGeomDecomp"
88 "(const dictionary& decompositionDict)",
90 ) <<
"Illegal decomposition order " << order <<
endl
99 const List<scalar>& l,
105 if (initHigh <= initLow)
111 label high = initHigh;
113 while ((high - low) > 1)
115 label mid = (low + high)/2;
147 void Foam::hierarchGeomDecomp::calculateSortedWeightedSizes
152 const label globalCurrentSize,
158 sortedWeightedSizes[0] = 0;
161 label pointI = current[indices[i]];
162 sortedWeightedSizes[i + 1] = sortedWeightedSizes[i] + weights[pointI];
167 sortedWeightedSizes[current.size()],
172 sortedWeightedSizes *= (globalCurrentSize/globalCurrentLength);
178 void Foam::hierarchGeomDecomp::findBinary
181 const List<scalar>& values,
182 const label minIndex,
183 const scalar minValue,
184 const scalar maxValue,
185 const scalar wantedSize,
192 label low = minIndex;
193 scalar lowValue = minValue;
195 scalar highValue = maxValue;
197 label high = values.size();
200 scalar midValuePrev = VGREAT;
204 label size =
returnReduce(mid-minIndex, sumOp<label>());
208 Pout<<
" low:" << low <<
" lowValue:" << lowValue
209 <<
" high:" << high <<
" highValue:" << highValue
210 <<
" mid:" << mid <<
" midValue:" << midValue <<
endl
211 <<
" globalSize:" << size <<
" wantedSize:" << wantedSize
212 <<
" sizeTol:" << sizeTol <<
endl;
215 if (wantedSize < size - sizeTol)
218 highValue = midValue;
220 else if (wantedSize > size + sizeTol)
231 midValue = 0.5*(lowValue+highValue);
232 mid =
findLower(values, midValue, low, high);
235 bool hasNotChanged = (
mag(midValue-midValuePrev) < SMALL);
239 WarningIn(
"hierarchGeomDecomp::findBinary(..)")
240 <<
"unable to find desired deomposition split, making do!"
245 midValuePrev = midValue;
252 void Foam::hierarchGeomDecomp::findBinary
255 const List<scalar>& sortedWeightedSizes,
256 const List<scalar>& values,
257 const label minIndex,
258 const scalar minValue,
259 const scalar maxValue,
260 const scalar wantedSize,
267 label low = minIndex;
268 scalar lowValue = minValue;
270 scalar highValue = maxValue;
272 label high = values.size();
275 scalar midValuePrev = VGREAT;
281 sortedWeightedSizes[mid] - sortedWeightedSizes[minIndex],
287 Pout<<
" low:" << low <<
" lowValue:" << lowValue
288 <<
" high:" << high <<
" highValue:" << highValue
289 <<
" mid:" << mid <<
" midValue:" << midValue <<
endl
290 <<
" globalSize:" << weightedSize
291 <<
" wantedSize:" << wantedSize
292 <<
" sizeTol:" << sizeTol <<
endl;
295 if (wantedSize < weightedSize - sizeTol)
298 highValue = midValue;
300 else if (wantedSize > weightedSize + sizeTol)
311 midValue = 0.5*(lowValue+highValue);
312 mid =
findLower(values, midValue, low, high);
315 bool hasNotChanged = (
mag(midValue-midValuePrev) < SMALL);
319 WarningIn(
"hierarchGeomDecomp::findBinary(..)")
320 <<
"unable to find desired deomposition split, making do!"
325 midValuePrev = midValue;
331 void Foam::hierarchGeomDecomp::sortComponent
342 label compI = decompOrder_[componentIndex];
346 Pout<<
"sortComponent : Sorting slice of size " << current.size()
347 <<
" in component " << compI <<
endl;
351 SortableList<scalar> sortedCoord(current.size());
355 label pointI = current[i];
357 sortedCoord[i] = points[pointI][compI];
361 label globalCurrentSize =
returnReduce(current.size(), sumOp<label>());
377 ? sortedCoord[sortedCoord.size()-1]
385 Pout<<
"sortComponent : minCoord:" << minCoord
386 <<
" maxCoord:" << maxCoord <<
endl;
392 scalar leftCoord = minCoord;
395 for (label bin = 0; bin < n_[compI]; bin++)
403 label localSize = -1;
406 scalar rightCoord = -GREAT;
408 if (bin == n_[compI]-1)
411 localSize = current.size()-leftIndex;
412 rightCoord = maxCoord;
417 localSize = label(current.size()/n_[compI]);
418 rightCoord = sortedCoord[leftIndex+localSize];
426 label rightIndex = current.size();
427 rightCoord = maxCoord;
437 globalCurrentSize/n_[compI],
442 localSize = rightIndex - leftIndex;
447 Pout<<
"For component " << compI <<
", bin " << bin
448 <<
" copying" <<
endl
449 <<
"from " << leftCoord <<
" at local index "
451 <<
"to " << rightCoord <<
" localSize:"
462 label pointI = current[sortedCoord.indices()[leftIndex+i]];
465 finalDecomp[pointI] += bin*mult;
472 if (componentIndex < 2)
498 leftIndex += localSize;
499 leftCoord = rightCoord;
505 void Foam::hierarchGeomDecomp::sortComponent
517 label compI = decompOrder_[componentIndex];
521 Pout<<
"sortComponent : Sorting slice of size " << current.size()
522 <<
" in component " << compI <<
endl;
526 SortableList<scalar> sortedCoord(current.size());
530 label pointI = current[i];
532 sortedCoord[i] = points[pointI][compI];
536 label globalCurrentSize =
returnReduce(current.size(), sumOp<label>());
540 scalarField sortedWeightedSizes(current.size()+1, 0);
541 calculateSortedWeightedSizes
544 sortedCoord.indices(),
564 ? sortedCoord[sortedCoord.size()-1]
572 Pout<<
"sortComponent : minCoord:" << minCoord
573 <<
" maxCoord:" << maxCoord <<
endl;
579 scalar leftCoord = minCoord;
582 for (label bin = 0; bin < n_[compI]; bin++)
590 label localSize = -1;
593 scalar rightCoord = -GREAT;
595 if (bin == n_[compI]-1)
598 localSize = current.size()-leftIndex;
599 rightCoord = maxCoord;
608 label rightIndex = current.size();
609 rightCoord = maxCoord;
620 globalCurrentSize/n_[compI],
625 localSize = rightIndex - leftIndex;
630 Pout<<
"For component " << compI <<
", bin " << bin
631 <<
" copying" <<
endl
632 <<
"from " << leftCoord <<
" at local index "
634 <<
"to " << rightCoord <<
" localSize:"
645 label pointI = current[sortedCoord.indices()[leftIndex+i]];
648 finalDecomp[pointI] += bin*mult;
655 if (componentIndex < 2)
682 leftIndex += localSize;
683 leftCoord = rightCoord;
690 Foam::hierarchGeomDecomp::hierarchGeomDecomp
702 Foam::hierarchGeomDecomp::hierarchGeomDecomp
708 geomDecomp(decompositionDict, hierarchGeomDecomp::typeName),
738 label allSize = points.
size();
741 const label sizeTol =
max(1, label(1
E-3*allSize/nProcessors_));
779 label allSize = points.
size();
782 const label sizeTol =
max(1, label(1
E-3*allSize/nProcessors_));