29#ifndef CUBOID_GEOMETRY_3D_HH
30#define CUBOID_GEOMETRY_3D_HH
46 : _motherCuboid(0,0,0,0,0,0,0), _periodicityOn(false), clout(std::cout,
"CuboidGeometry3D")
55 int nX,
int nY,
int nZ,
int nC)
56 : _motherCuboid(originX, originY, originZ, deltaR, nX, nY, nZ),
57 _periodicityOn(false), clout(std::cout,
"CuboidGeometry3D")
59 _cuboids.reserve(nC+2);
66 std::vector<int> extent,
int nC)
68 extent[0], extent[1], extent[2], nC)
75 extent[0], extent[1], extent[2], nC)
81 motherCuboid.getOrigin(), motherCuboid.getDeltaR(),
82 motherCuboid.getExtent(), nC)
87 : _motherCuboid( indicatorF.getMin()[0], indicatorF.getMin()[1], indicatorF.getMin()[2], voxelSize,
88 (int)((indicatorF.getMax()[0] - indicatorF.getMin()[0]) / voxelSize + 1.5),
89 (int)((indicatorF.getMax()[1] - indicatorF.getMin()[1]) / voxelSize + 1.5),
90 (int)((indicatorF.getMax()[2] - indicatorF.getMin()[2]) / voxelSize + 1.5)),
91 _periodicityOn(false), clout(std::cout,
"CuboidGeometry3D")
93 _cuboids.reserve(nC+2);
109 : _motherCuboid( indicatorF.getMin()[0], indicatorF.getMin()[1], indicatorF.getMin()[2], voxelSize,
110 (int)((indicatorF.getMax()[0] - indicatorF.getMin()[0]) / voxelSize + 1.5),
111 (int)((indicatorF.getMax()[1] - indicatorF.getMin()[1]) / voxelSize + 1.5),
112 (int)((indicatorF.getMax()[2] - indicatorF.getMin()[2]) / voxelSize + 1.5)),
113 _periodicityOn(false), clout(std::cout,
"CuboidGeometry3D")
115 _cuboids.reserve(nC+2);
118 if ( minimizeBy ==
"volume" ) {
121 else if ( minimizeBy ==
"weight" ) {
151 return _motherCuboid;
157 return _motherCuboid;
163 _periodicityOn[0] = periodicityX;
164 _periodicityOn[1] = periodicityY;
165 _periodicityOn[2] = periodicityZ;
173 for (i = 0; i < _cuboids.size(); i++) {
174 if (_cuboids[i].checkPoint(x, y, z, offset)) {
185 return get_iC(coords[0], coords[1], coords[2], offset);
190 int orientationZ)
const
193 for (i = 0; i < _cuboids.size(); i++) {
194 if (_cuboids[i].checkPoint(x, y, z) &&
195 _cuboids[i].checkPoint(x + orientationX / _cuboids[i].getDeltaR(),
196 y + orientationY / _cuboids[i].getDeltaR(),
197 z + orientationZ / _cuboids[i].getDeltaR())) {
207 const int iCtmp = get_iC(physR[0], physR[1], physR[2]);
208 if (iCtmp < getNc()) {
220 const int iCtmp = get_iC(physR);
221 if (iCtmp < getNc()) {
233 const int iCtmp = get_iC(physR);
234 if (iCtmp < getNc()) {
246 int iCtmp = get_iC(physR[0], physR[1], physR[2]);
247 if (iCtmp < getNc()) {
249 latticeR[1] = (int)
util::floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
250 latticeR[2] = (int)
util::floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
251 latticeR[3] = (int)
util::floor( (physR[2] - _cuboids[latticeR[0]].getOrigin()[2] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
262 int iCtmp = get_iC(physR[0], physR[1], physR[2]);
263 if (iCtmp < getNc()) {
265 latticeR[1] = (int)
util::floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() );
266 latticeR[2] = (int)
util::floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() );
267 latticeR[3] = (int)
util::floor( (physR[2] - _cuboids[latticeR[0]].getOrigin()[2] ) / _cuboids[latticeR[0]].getDeltaR() );
279 int iCtmp = get_iC(physR[0], physR[1], physR[2]);
280 if (iCtmp < getNc()) {
282 latticeR[1] = (int)
util::floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() );
283 latticeR[2] = (int)
util::floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() );
284 latticeR[3] = (int)
util::floor( (physR[2] - _cuboids[latticeR[0]].getOrigin()[2] ) / _cuboids[latticeR[0]].getDeltaR() );
295 _cuboids[iCglob].getPhysR(physR, iX, iY, iZ);
296 for (
int iDim = 0; iDim < 3; iDim++) {
297 if (_periodicityOn[iDim]) {
299 physR[iDim] = remainder( physR[iDim] - _motherCuboid.getOrigin()[iDim]
300 + _motherCuboid.getDeltaR() * (_motherCuboid.getExtent()[iDim]),
301 _motherCuboid.getDeltaR() * (_motherCuboid.getExtent()[iDim]));
303 if ( physR[iDim]*physR[iDim] < 0.001 * _motherCuboid.getDeltaR()*_motherCuboid.getDeltaR() ) {
307 if ( physR[iDim] < 0 ) {
308 physR[iDim] += _motherCuboid.getDeltaR() *( _motherCuboid.getExtent()[iDim]);
311 physR[iDim] += _motherCuboid.getOrigin()[iDim];
320 getPhysR(physR, latticeR[0], latticeR[1], latticeR[2], latticeR[3]);
326 getPhysR(physR, latticeR[0], latticeR[1], latticeR[2], latticeR[3]);
332 getPhysR(physR, iCglob, latticeR[0], latticeR[1], latticeR[2]);
342 for (
int iC = 0; iC < getNc(); iC++) {
346 T globX = get(iC).getOrigin()[0];
347 T globY = get(iC).getOrigin()[1];
348 T globZ = get(iC).getOrigin()[2];
349 T nX = get(iC).getNx();
350 T nY = get(iC).getNy();
351 T nZ = get(iC).getNz();
352 T deltaR = get(iC).getDeltaR();
353 if (get(cuboid).checkInters(globX,
354 globX + (nX + overlap - 1)*deltaR,
355 globY - overlap * deltaR,
356 globY + (nY + overlap - 1)*deltaR,
357 globZ - overlap * deltaR,
358 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
363 if (_periodicityOn[0]) {
364 if (get(cuboid).getOrigin()[0] + (get(cuboid).getNx() + overlap - 1)*get(cuboid).getDeltaR() > getMaxPhysR()[0]) {
365 Cuboid3D<T> cub(get(cuboid).getOrigin()[0]-getMaxPhysR()[0],
366 get(cuboid).getOrigin()[1],
367 get(cuboid).getOrigin()[2],
368 get(cuboid).getDeltaR(),
371 get(cuboid).getNz());
373 globX + (nX + overlap - 1)*deltaR,
374 globY - overlap * deltaR,
375 globY + (nY + overlap - 1)*deltaR,
376 globZ - overlap * deltaR,
377 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
381 if (get(cuboid).getOrigin()[0] - overlap*get(cuboid).getDeltaR() < getMinPhysR()[0]) {
382 Cuboid3D<T> cub(get(cuboid).getOrigin()[0]+getMaxPhysR()[0],
383 get(cuboid).getOrigin()[1],
384 get(cuboid).getOrigin()[2],
385 get(cuboid).getDeltaR(),
388 get(cuboid).getNz());
390 globX + (nX + overlap - 1)*deltaR,
391 globY - overlap * deltaR,
392 globY + (nY + overlap - 1)*deltaR,
393 globZ - overlap * deltaR,
394 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
400 if (_periodicityOn[1]) {
401 if (get(cuboid).getOrigin()[1] + (get(cuboid).getNy() + overlap - 1)*get(cuboid).getDeltaR() > getMaxPhysR()[1]) {
403 get(cuboid).getOrigin()[1]-getMaxPhysR()[1],
404 get(cuboid).getOrigin()[2],
405 get(cuboid).getDeltaR(),
408 get(cuboid).getNz());
410 globX + (nX + overlap - 1)*deltaR,
411 globY - overlap * deltaR,
412 globY + (nY + overlap - 1)*deltaR,
413 globZ - overlap * deltaR,
414 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
418 if (get(cuboid).getOrigin()[1] - overlap*get(cuboid).getDeltaR() < getMinPhysR()[1]) {
420 get(cuboid).getOrigin()[1]+getMaxPhysR()[1],
421 get(cuboid).getOrigin()[2],
422 get(cuboid).getDeltaR(),
425 get(cuboid).getNz());
427 globX + (nX + overlap - 1)*deltaR,
428 globY - overlap * deltaR,
429 globY + (nY + overlap - 1)*deltaR,
430 globZ - overlap * deltaR,
431 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
437 if (_periodicityOn[2]) {
438 if (get(cuboid).getOrigin()[2] + (get(cuboid).getNz() + overlap - 1)*get(cuboid).getDeltaR() > getMaxPhysR()[2]) {
440 get(cuboid).getOrigin()[1],
441 get(cuboid).getOrigin()[2]-getMaxPhysR()[2],
442 get(cuboid).getDeltaR(),
445 get(cuboid).getNz());
447 globX + (nX + overlap - 1)*deltaR,
448 globY - overlap * deltaR,
449 globY + (nY + overlap - 1)*deltaR,
450 globZ - overlap * deltaR,
451 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
455 if (get(cuboid).getOrigin()[2] - overlap*get(cuboid).getDeltaR() < getMinPhysR()[2]) {
457 get(cuboid).getOrigin()[1],
458 get(cuboid).getOrigin()[2]+getMaxPhysR()[2],
459 get(cuboid).getDeltaR(),
462 get(cuboid).getNz());
464 globX + (nX + overlap - 1)*deltaR,
465 globY - overlap * deltaR,
466 globY + (nY + overlap - 1)*deltaR,
467 globZ - overlap * deltaR,
468 globZ + (nZ + overlap - 1)*deltaR, overlap)) {
475 std::set<int>::iterator it = dummy.begin();
476 for (; it != dummy.end(); ++it) {
477 neighbours.push_back(*it);
484 return _cuboids.size();
491 for (
unsigned i = 0; i < _cuboids.size(); i++) {
492 if ((T)_cuboids[i].getNx() / (T)_cuboids[i].getNy() < minRatio) {
493 minRatio = (T)_cuboids[i].getNx() / (T)_cuboids[i].getNy();
495 if ((T)_cuboids[i].getNy() / (T)_cuboids[i].getNz() < minRatio) {
496 minRatio = (T)_cuboids[i].getNy() / (T)_cuboids[i].getNz();
498 if ((T)_cuboids[i].getNz() / (T)_cuboids[i].getNx() < minRatio) {
499 minRatio = (T)_cuboids[i].getNz() / (T)_cuboids[i].getNx();
509 for (
unsigned i = 0; i < _cuboids.size(); i++) {
510 if ((T)_cuboids[i].getNx() / (T)_cuboids[i].getNy() > maxRatio) {
511 maxRatio = (T)_cuboids[i].getNx() / (T)_cuboids[i].getNy();
513 if ((T)_cuboids[i].getNy() / (T)_cuboids[i].getNz() > maxRatio) {
514 maxRatio = (T)_cuboids[i].getNy() / (T)_cuboids[i].getNz();
516 if ((T)_cuboids[i].getNz() / (T)_cuboids[i].getNx() > maxRatio) {
517 maxRatio = (T)_cuboids[i].getNz() / (T)_cuboids[i].getNx();
527 for (
unsigned i = 0; i < _cuboids.size(); i++) {
528 if (_cuboids[i].getOrigin()[0] < output[0]) {
529 output[0] = _cuboids[i].getOrigin()[0];
531 if (_cuboids[i].getOrigin()[1] < output[1]) {
532 output[1] = _cuboids[i].getOrigin()[1];
534 if (_cuboids[i].getOrigin()[2] < output[2]) {
535 output[2] = _cuboids[i].getOrigin()[2];
545 output[0] += _cuboids[0].getNx()*_cuboids[0].getDeltaR();
546 output[1] += _cuboids[0].getNy()*_cuboids[0].getDeltaR();
547 output[2] += _cuboids[0].getNz()*_cuboids[0].getDeltaR();
548 for (
unsigned i = 0; i < _cuboids.size(); i++) {
549 if (_cuboids[i].getOrigin()[0] + _cuboids[i].getNx()*_cuboids[i].getDeltaR() > output[0]) {
550 output[0] = _cuboids[i].getOrigin()[0] + _cuboids[i].getNx()*_cuboids[i].getDeltaR();
552 if (_cuboids[i].getOrigin()[1] + _cuboids[i].getNy()*_cuboids[i].getDeltaR() > output[1]) {
553 output[1] = _cuboids[i].getOrigin()[1] + _cuboids[i].getNy()*_cuboids[i].getDeltaR();
555 if (_cuboids[i].getOrigin()[2] + _cuboids[i].getNz()*_cuboids[i].getDeltaR() > output[2]) {
556 output[2] = _cuboids[i].getOrigin()[2] + _cuboids[i].getNz()*_cuboids[i].getDeltaR();
565 T minVolume = _cuboids[0].getPhysVolume();
566 for (
unsigned i = 0; i < _cuboids.size(); i++) {
567 if (_cuboids[i].getPhysVolume() < minVolume) {
568 minVolume = _cuboids[i].getPhysVolume();
577 T maxVolume = _cuboids[0].getPhysVolume();
578 for (
unsigned i = 0; i < _cuboids.size(); i++) {
579 if (_cuboids[i].getPhysVolume() > maxVolume) {
580 maxVolume = _cuboids[i].getPhysVolume();
589 size_t minNodes = _cuboids[0].getLatticeVolume();
590 for (
unsigned i = 0; i < _cuboids.size(); i++) {
591 if (_cuboids[i].getLatticeVolume() < minNodes) {
592 minNodes = _cuboids[i].getLatticeVolume();
601 size_t maxNodes = _cuboids[0].getLatticeVolume();
602 for (
unsigned i = 0; i < _cuboids.size(); i++) {
603 if (_cuboids[i].getLatticeVolume() > maxNodes) {
604 maxNodes = _cuboids[i].getLatticeVolume();
613 size_t minNodes = _cuboids[0].getWeight();
614 for (
unsigned i = 0; i < _cuboids.size(); i++) {
615 if (_cuboids[i].getWeight() < minNodes) {
616 minNodes = _cuboids[i].getWeight();
625 size_t maxNodes = _cuboids[0].getWeight();
626 for (
unsigned i = 0; i < _cuboids.size(); i++) {
627 if (_cuboids[i].getWeight() > maxNodes) {
628 maxNodes = _cuboids[i].getWeight();
637 T minDelta = _cuboids[0].getDeltaR();
638 for (
unsigned i = 0; i < _cuboids.size(); i++) {
639 if (_cuboids[i].getDeltaR() < minDelta) {
640 minDelta = _cuboids[i].getDeltaR();
649 T maxDelta = _cuboids[0].getDeltaR();
650 for (
unsigned i = 0; i < _cuboids.size(); i++) {
651 if (_cuboids[i].getDeltaR() > maxDelta) {
652 maxDelta = _cuboids[i].getDeltaR();
662 return _motherCuboid == rhs._motherCuboid
663 && _periodicityOn == rhs._periodicityOn
664 && _cuboids == rhs._cuboids;
673 _cuboids.push_back(cuboid);
680 _cuboids.erase(_cuboids.begin() + iC);
687 std::vector<bool> allZero;
690 for (
unsigned iC = 0; iC < _cuboids.size(); iC++) {
692 allZero.push_back(
true);
693 for (
int iX = 0; iX < _cuboids[iC].getNx(); iX++) {
694 for (
int iY = 0; iY < _cuboids[iC].getNy(); iY++) {
695 for (
int iZ = 0; iZ < _cuboids[iC].getNz(); iZ++) {
699 getPhysR(physR,latticeR);
701 indicatorF(inside,physR);
709 for (
int iC = _cuboids.size() - 1; iC >= 0; iC--) {
719 std::vector<bool> allZero(_cuboids.size(),
false);
720 for (
unsigned iC = 0; iC < _cuboids.size(); iC++) {
721 allZero[iC] = (_cuboids[iC].getWeight() == 0);
723 for (
int iC = _cuboids.size() - 1; iC >= 0; iC--) {
738 size_t fullCells = 0;
739 int xN = get(iC).getNx();
740 int yN = get(iC).getNy();
741 int zN = get(iC).getNz();
748 for (
int iX = 0; iX < xN; iX++) {
749 for (
int iY = 0; iY < yN; iY++) {
750 for (
int iZ = 0; iZ < zN; iZ++) {
754 getPhysR(physR,latticeR);
755 indicatorF(inside,physR);
777 get(iC).setWeight(fullCells);
778 _cuboids[iC].resize(newX, newY, newZ, maxX - newX + 1, maxY - newY + 1, maxZ - newZ + 1);
788 for (
int iC = getNc() - 1; iC >= 0; iC--) {
789 shrink(iC, indicatorF);
794 T minDelataR = getMinDeltaR();
795 _motherCuboid =
Cuboid3D<T>(minPhysR[0], minPhysR[1], minPhysR[2], minDelataR,
796 (
int)((maxPhysR[0]-minPhysR[0])/minDelataR + 0.5),
797 (
int)((maxPhysR[1]-minPhysR[1])/minDelataR + 0.5),
798 (
int)((maxPhysR[2]-minPhysR[2])/minDelataR + 0.5));
804 _motherCuboid.refine(factor);
805 for (
auto& cuboid : _cuboids) {
806 cuboid.refine(factor);
813 const T tolerance = std::numeric_limits<T>::epsilon();
814 const T currDeltaR = _motherCuboid.getDeltaR();
815 const int factor = std::ceil(currDeltaR / goalDeltaR);
816 if (
util::fabs(currDeltaR / factor - goalDeltaR) < tolerance) {
827 Cuboid3D<T> temp(_cuboids[iC].getOrigin()[0], _cuboids[iC].getOrigin()[1],
828 _cuboids[iC].getOrigin()[2], _cuboids[iC].getDeltaR(),
829 _cuboids[iC].getNx(), _cuboids[iC].getNy(), _cuboids[iC].getNz());
837 Cuboid3D<T> temp(_cuboids[iC].getOrigin()[0], _cuboids[iC].getOrigin()[1],
838 _cuboids[iC].getOrigin()[2], _cuboids[iC].getDeltaR(),
839 _cuboids[iC].getNx(), _cuboids[iC].getNy(), _cuboids[iC].getNz());
840 const int p = std::max(1, temp.
getNx() / width);
841 const int q = std::max(1, temp.
getNy() / width);
842 const int r = std::max(1, temp.
getNz() / width);
843 temp.
divide(p, q, r, _cuboids);
850 T averageWeight = get(iC).getWeight() / (T) p;
852 Cuboid3D<T> temp(_cuboids[iC].getOrigin()[0], _cuboids[iC].getOrigin()[1],
853 _cuboids[iC].getOrigin()[2], _cuboids[iC].getDeltaR(),
854 _cuboids[iC].getNx(), _cuboids[iC].getNy(), _cuboids[iC].getNz());
859 int xN = get(iC).getNx();
860 int yN = get(iC).getNy();
861 int zN = get(iC).getNz();
862 T deltaR = get(iC).getDeltaR();
866 std::vector<int> extend_child = {xN, yN, zN};
867 int localPos_child = 0;
870 if ( get(iC).getNx() >= get(iC).getNy() && get(iC).getNx() >= get(iC).getNz()) {
874 for (
int iChild = 0; iChild < p - 1; iChild++) {
876 int fullCells_minusOne = 0;
878 for (
int iX = localPos_child; iX < xN; iX++) {
879 fullCells_minusOne = fullCells;
880 for (
int iY = 0; iY < yN; iY++) {
881 for (
int iZ = 0; iZ < zN; iZ++) {
885 getPhysR(physR,latticeR);
887 indicatorF(inside,physR);
894 if ( fullCells >= averageWeight ) {
895 if ( (fullCells - averageWeight) > (averageWeight - fullCells_minusOne) ) {
899 extend_child[0] = iX - localPos_child + 1;
901 Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]);
902 _cuboids.push_back(child);
904 globPos_child[0] += extend_child[0]*deltaR;
905 localPos_child += extend_child[0] + 1;
914 extend_child[0] = xN - localPos_child + p - 1;
916 Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]);
917 _cuboids.push_back(child);
922 else if ( get(iC).getNy() >= get(iC).getNx() && get(iC).getNy() >= get(iC).getNz()) {
925 for (
int iChild = 0; iChild < p - 1; iChild++) {
927 int fullCells_minusOne = 0;
929 for (
int iY = localPos_child; iY < yN; iY++) {
930 fullCells_minusOne = fullCells;
931 for (
int iX = 0; iX < xN; iX++) {
932 for (
int iZ = 0; iZ < zN; iZ++) {
936 getPhysR(physR,latticeR);
938 indicatorF(inside,physR);
944 if ( fullCells >= averageWeight ) {
945 if ( (fullCells - averageWeight) > (averageWeight - fullCells_minusOne) ) {
949 extend_child[1] = iY - localPos_child + 1;
951 Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]);
952 _cuboids.push_back(child);
954 globPos_child[1] += extend_child[1]*deltaR;
955 localPos_child += extend_child[1] + 1;
964 extend_child[1] = yN - localPos_child + p - 1;
966 Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]);
967 _cuboids.push_back(child);
974 for (
int iChild = 0; iChild < p - 1; iChild++) {
976 int fullCells_minusOne = 0;
978 for (
int iZ = localPos_child; iZ < zN; iZ++) {
979 fullCells_minusOne = fullCells;
980 for (
int iY = 0; iY < yN; iY++) {
981 for (
int iX = 0; iX < xN; iX++) {
985 getPhysR(physR,latticeR);
987 indicatorF(inside,physR);
993 if ( fullCells >= averageWeight ) {
994 if ( (fullCells - averageWeight) > (averageWeight - fullCells_minusOne) ) {
998 extend_child[2] = iZ - localPos_child + 1;
1000 Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]);
1001 _cuboids.push_back(child);
1003 globPos_child[2] += extend_child[2]*deltaR;
1004 localPos_child += extend_child[2] + 1;
1013 extend_child[2] = zN - localPos_child + p - 1;
1015 Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]);
1016 _cuboids.push_back(child);
1033 std::swap(this->_cuboids, rhs._cuboids);
1034 std::swap(this->_motherCuboid, rhs._motherCuboid);
1035 std::swap(this->_periodicityOn[0], rhs._periodicityOn[0]);
1036 std::swap(this->_periodicityOn[1], rhs._periodicityOn[1]);
1037 std::swap(this->_periodicityOn[2], rhs._periodicityOn[2]);
1038 std::swap(this->clout, rhs.clout);
1044 _cuboids.swap(cuboids);
1050 this->_cuboids.clear();
1051 for (
unsigned iC = 0; iC < cuboids.size(); iC++) {
1059 #ifdef PARALLEL_MODE_OMP
1060 #pragma omp parallel for schedule(dynamic,1)
1062 for (
int iC=0; iC < getNc(); ++iC) {
1063 int latticeR[4] { iC, 0, 0, 0 };
1065 int xN = get(iC).getNx();
1066 int yN = get(iC).getNy();
1067 int zN = get(iC).getNz();
1068 size_t fullCells = 0;
1069 for (
int iX = 0; iX < xN; iX++) {
1070 for (
int iY = 0; iY < yN; iY++) {
1071 for (
int iZ = 0; iZ < zN; iZ++) {
1075 getPhysR(physR,latticeR);
1077 indicatorF(inside,physR);
1084 get(iC).setWeight(fullCells);
1093 + _motherCuboid.getNblock()
1094 + _cuboids.size() > 0 ? 1 + _cuboids.size() * _cuboids[0].getNblock() : 0;
1101 return 3 *
sizeof(bool)
1102 + _motherCuboid.getSerializableSize()
1103 + (_cuboids.size() > 0 ?
1104 sizeof(size_t) + _cuboids.size() * _cuboids[0].getSerializableSize() :
1111 std::size_t currentBlock = 0;
1112 size_t sizeBufferIndex = 0;
1113 bool* dataPtr =
nullptr;
1115 registerVar<bool> (iBlock, sizeBlock, currentBlock, dataPtr, _periodicityOn[0], 3);
1116 registerSerializableOfConstSize (iBlock, sizeBlock, currentBlock, dataPtr, _motherCuboid, loadingMode);
1117 registerStdVectorOfSerializablesOfConstSize (iBlock, sizeBlock, currentBlock, sizeBufferIndex, dataPtr,
1118 _cuboids, loadingMode);
1127 clout <<
"---Cuboid Stucture Statistics---" << std::endl;
1128 clout <<
" Number of Cuboids: " <<
"\t" << getNc() << std::endl;
1129 clout <<
" Delta (min): " <<
"\t" <<
"\t" << getMinDeltaR() << std::endl;
1130 clout <<
" (max): " <<
"\t" <<
"\t" << getMaxDeltaR() << std::endl;
1131 clout <<
" Ratio (min): " <<
"\t" <<
"\t" << getMinRatio() << std::endl;
1132 clout <<
" (max): " <<
"\t" <<
"\t" << getMaxRatio() << std::endl;
1133 clout <<
" Nodes (min): " <<
"\t" <<
"\t" << getMinLatticeVolume() << std::endl;
1134 clout <<
" (max): " <<
"\t" <<
"\t" << getMaxLatticeVolume() << std::endl;
1135 clout <<
" Weight (min): " <<
"\t" <<
"\t" << getMinLatticeWeight() << std::endl;
1136 clout <<
" (max): " <<
"\t" <<
"\t" << getMaxLatticeWeight() << std::endl;
1137 clout <<
"--------------------------------" << std::endl;
1143 clout <<
"Mothercuboid :" << std::endl;
1144 getMotherCuboid().print();
1146 for (
int iC = 0; iC < getNc(); iC++) {
1147 clout <<
"Cuboid #" << iC <<
": " << std::endl;
1159 fout.open(completeFileName.c_str(), std::ios::app);
1161 clout <<
"Error: could not open " << completeFileName << std::endl;
1165 fout <<
"<CuboidGeometry dimension=\"3\" " << _cuboidParameters(getMotherCuboid()) <<
">\n";
1168 for (
int iC = 0; iC < getNc(); ++iC) {
1169 fout <<
"<Cuboid " << _cuboidParameters(get(iC)) <<
" />\n";
1172 fout <<
"</CuboidGeometry>\n";
1186 fout.open(fname.c_str(), std::ios::trunc);
1187 fout <<
"<?xml version=\"1.0\"?>\n";
1188 fout <<
"<XMLContent>\n";
1193 writeToExistingFile(fname, loadBalancer);
1196 fout.open(fname.c_str(), std::ios::app);
1197 fout <<
"</XMLContent>\n";
1208 std::stringstream ss;
1209 ss.flags(std::ios::scientific);
1210 ss.precision (std::numeric_limits<double>::digits10 + 1);
1212 for (
int i = 0; i<3; i++) {
1216 ss <<
"\" origin=\"";
1217 for (
int i = 0; i<3; i++) {
1221 ss <<
"\" deltaR=\"" << cub.
getDeltaR();
A regular single 3D cuboid is the basic component of a 3D cuboid structure which defines the grid.
void divideFractional(int iD, std::vector< T > fractions, std::vector< Cuboid3D< T > > &childrenC) const
Divides the cuboid into fractions along the iDth dimension.
Vector< T, 3 > getOrigin() const
Read only access to left lower corner coordinates.
Vector< int, 3 > const getExtent() const
Read only access to the number of voxels in every dimension.
int getWeightValue() const
Returns the actual value of weight (-1 for getLatticeVolume())
bool checkInters(T globX0, T globX1, T globY0, T globY1, T globZ0, T globZ1, int overlap=0) const
Checks whether there is an intersection with the cuboid extended with an layer of size overlap*delta.
T getDeltaR() const
Read only access to the distance of cuboid nodes.
void divide(int p, int q, int r, std::vector< Cuboid3D< T > > &childrenC) const
Divides the cuboid in p*q*r cuboids of equal volume and add them to the given vector.
int getNz() const
Read access to cuboid depth.
int getNy() const
Read access to cuboid height.
int getNx() const
Read access to cuboid width.
A cuboid geometry represents a voxel mesh.
void splitByWeight(int iC, int p, IndicatorF3D< T > &indicatorF)
Splits cuboid iC, removes it and adds p cuboids of same weight.
T getMinRatio() const
Returns the minimum of the ratio nX/NY in the structure.
Cuboid3D< T > & get(int iC)
Read and write access to a single cuboid.
void swap(CuboidGeometry3D< T > &rhs)
Swaps data from input into object.
size_t getSerializableSize() const override
Binary size for the serializer interface.
bool operator==(CuboidGeometry3D< T > &rhs)
Compares two CuboidGeometries.
T getMaxRatio() const
Returns the maximum of the ratio nX/NY in the structure.
void refine(int factor)
Refines mesh by splitting each cell into factor^3 cells.
bool getFloorLatticeR(const std::vector< T > &physR, std::vector< int > &latticeR) const
Returns true and the util::floor lattice position to the given physical position if the physical posi...
void writeToFile(std::string fileName, LoadBalancer< T > &loadBalancer)
Save CuboidGeometry into XML File.
void shrink(int iC, IndicatorF3D< T > &indicatorF)
Shrink cuboid iC so that no empty planes are left.
size_t getNblock() const override
Number of data blocks for the serializer interface.
void swapCuboids(std::vector< Cuboid3D< T > > &cuboids)
Swaps the vector of cuboids.
void replaceCuboids(std::vector< Cuboid3D< T > > &cuboids)
Replace the vector of cuboids.
int get_iC(T globX, T globY, T globZ, int offset=0) const
Gives for a given point (globX/globY/globZ) the related cuboidID and _p if the point is not in any of...
bool * getBlock(std::size_t iBlock, std::size_t &sizeBlock, bool loadingMode) override
Return a pointer to the memory of the current block and its size for the serializable interface.
T getMaxDeltaR() const
Returns the maximum delta in the structure.
void print() const
Prints cuboid geometry details.
bool getC(T physR[3], int &iC) const
Returns true and the cuboid number of the nearest lattice position to the given physical position if ...
void split(int iC, int p)
Splits cuboid iC, removes it and adds p cuboids of same volume.
bool tryRefineTo(T deltaR)
Tries to refine mesh to given deltaR.
Vector< T, 3 > getMinPhysR() const
Returns the minimum coordinate in the structure.
virtual ~CuboidGeometry3D()
Destructs CuboidGeometry.
size_t getMaxLatticeVolume() const
Returns the maximum number of nodes in the structure.
bool getLatticeR(int latticeR[4], const T physR[3]) const
Returns true and the nearest lattice position to the given physical position if the physical position...
void remove(int iC)
Removes the cuboid iC.
void writeToExistingFile(std::string completeFileName, LoadBalancer< T > &loadBalancer)
Save CuboidGeometry into an existing XML File.
T getMinPhysVolume() const
Returns the minimum volume in the structure.
int getNc() const
Returns the number of cuboids in the structure.
void setWeights(IndicatorF3D< T > &indicatorF)
Sets the number of full cells of each cuboid.
Vector< T, 3 > getMaxPhysR() const
Returns the maximum coordinate in the structure.
void splitRegular(int iC, int width)
Splits cuboid iC, removes it, adds approx. width^3 sized new cuboids.
Cuboid3D< T > getMotherCuboid()
Returns the smallest cuboid that includes all cuboids of the structure.
CuboidGeometry3D()
Constructs empty Geometry.
void add(Cuboid3D< T > cuboid)
Adds a cuboid.
size_t getMinLatticeVolume() const
Returns the minimum number of nodes in the structure.
void removeByWeight()
Removes all cuboids where weight = 0.
void printExtended()
Prints cuboid geometry details plus details of all cuboids.
size_t getMaxLatticeWeight() const
Returns the maximum number of nodes in the structure inside the indicator.
T getMinDeltaR() const
Returns the minimum delta in the structure.
size_t getMinLatticeWeight() const
Returns the minimum number of nodes in the structure inside the indicator.
void getPhysR(T physR[3], const int &iCglob, const int &iX, const int &iY, const int &iZ) const
Returns the physical position to the given lattice position respecting periodicity for the overlap no...
void splitFractional(int iC, int iD, std::vector< T > fractions)
Splits cuboid iC along dimension iD into cuboids of fractions.
T getMaxPhysVolume() const
Returns the maximum volume in the structure.
void getNeighbourhood(int cuboid, std::vector< int > &neighbours, int overlap=0)
Stores the iC of the neighbouring cuboids in a vector;.
void setPeriodicity(bool periodicityX, bool periodicityY, bool periodicityZ)
Set flag to enable/disable periodicity depending of direction. Be aware that not all directions are t...
IndicatorF3D is an application from .
Base class for all LoadBalancer.
std::string getLogOutDir() const
The description of a vector of 3D cuboid – header file.
Directories & directories()
ADf< T, DIM > floor(const ADf< T, DIM > &a)
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Top level namespace for all of OpenLB.
void minimizeByVolume(CuboidGeometry3D< T > &cGeometry, IndicatorF3D< T > &indicatorF, int nC)
Splits into nC cuboids by-volume.
void minimizeByWeight(CuboidGeometry3D< T > &cGeometry, IndicatorF3D< T > &indicatorF, int nC)