cuboid3d.hh

open lattice boltzmann project www.openlb.org

                }
            }
        }
    }

    int rest = p - bestIx*bestIy*bestIz;

    // split in one cuboid
    if (rest==0) {
        divide(bestIx, bestIy, bestIz, childrenC);
        return;
    }
    else {

        // add in z than in y direction
        if (nZ>nY && nZ>nX) {
        int restY = rest%bestIy;

            // split in two cuboid
            if (restY==0) {
                int restX = rest/bestIy;
                CuboidGeometry2D<T> helpG(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restX);

                int yN_child = 0;
                T globPosY_child = _globPosY;

                for (int iY=0; iY<bestIy; iY++) {
                    yN_child         = (_nY+bestIy-iY-1)/bestIy;
                    for (int iC=0; iC<helpG.get_nC(); iC++) {
                        int xN_child     = helpG.get_cuboid(iC).get_nX();
                        int zN_child     = helpG.get_cuboid(iC).get_nY();
                        T globPosX_child = helpG.get_cuboid(iC).get_globPosX();
                        T globPosZ_child = helpG.get_cuboid(iC).get_globPosY();

                        Cuboid3D<T> child(globPosX_child, globPosY_child, globPosZ_child,
                                   _delta, xN_child, yN_child, zN_child);
                        childrenC.push_back(child);

                    }
                    globPosY_child += yN_child*_delta;
                }
                return;
            }

            // spilt in four cuboid

            int restX = rest/bestIy+1;
            int yN_child = 0;
            T globPosY_child = _globPosY;
            int splited_nY = (int) (_nY * (T)((bestIx*bestIz+restX)*restY)/(T)p);
            CuboidGeometry2D<T> helpG0(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restX);

            for (int iY=0; iY<restY; iY++) {
                yN_child         = (splited_nY+restY-iY-1)/restY;
                for (int iC=0; iC<helpG0.get_nC(); iC++) {
                    int xN_child     = helpG0.get_cuboid(iC).get_nX();
                    int zN_child     = helpG0.get_cuboid(iC).get_nY();
                    T globPosX_child = helpG0.get_cuboid(iC).get_globPosX();
                    T globPosZ_child = helpG0.get_cuboid(iC).get_globPosY();

                    Cuboid3D<T> child(globPosX_child, globPosY_child, globPosZ_child,
                                   _delta, xN_child, yN_child, zN_child);
                    childrenC.push_back(child);
                }
                globPosY_child += yN_child*_delta;
            }

            splited_nY = _nY - splited_nY;
            restX = rest/bestIy;
            CuboidGeometry2D<T> helpG1(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restX);
            yN_child = 0;

            for (int iY=0; iY<bestIy-restY; iY++) {
                yN_child         = (splited_nY+bestIy-restY-iY-1)/(bestIy-restY);
                for (int iC=0; iC<helpG1.get_nC(); iC++) {
                    int xN_child     = helpG1.get_cuboid(iC).get_nX();
                    int zN_child     = helpG1.get_cuboid(iC).get_nY();
                    T globPosX_child = helpG1.get_cuboid(iC).get_globPosX();
                    T globPosZ_child = helpG1.get_cuboid(iC).get_globPosY();

                    Cuboid3D<T> child(globPosX_child, globPosY_child, globPosZ_child,
                               _delta, xN_child, yN_child, zN_child);
                    childrenC.push_back(child);
                }
                globPosY_child += yN_child*_delta;
            }
            return;
        }

        // add in x than in y direction
        else if (nX>nY && nX>nZ) {
        int restY = rest%bestIy;

            // split in two cuboid
            if (restY==0) {
                int restZ = rest/bestIy;
                CuboidGeometry2D<T> helpG(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restZ);

                int yN_child = 0;
                T globPosY_child = _globPosY;

                for (int iY=0; iY<bestIy; iY++) {
                    yN_child         = (_nY+bestIy-iY-1)/bestIy;
                    for (int iC=0; iC<helpG.get_nC(); iC++) {
                        int xN_child     = helpG.get_cuboid(iC).get_nX();
                        int zN_child     = helpG.get_cuboid(iC).get_nY();
                        T globPosX_child = helpG.get_cuboid(iC).get_globPosX();
                        T globPosZ_child = helpG.get_cuboid(iC).get_globPosY();

                        Cuboid3D<T> child(globPosX_child, globPosY_child, globPosZ_child,
                                   _delta, xN_child, yN_child, zN_child);
                        childrenC.push_back(child);

                    }
                    globPosY_child += yN_child*_delta;
                }
                return;
            }

            // spilt in four cuboid

            int restZ = rest/bestIy+1;

            int yN_child = 0;
            T globPosY_child = _globPosY;
            int splited_nY = (int) (_nY * (T)((bestIx*bestIz+restZ)*restY)/(T)p);
            CuboidGeometry2D<T> helpG0(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restZ);

            for (int iY=0; iY<restY; iY++) {
                yN_child         = (splited_nY+restY-iY-1)/restY;
                for (int iC=0; iC<helpG0.get_nC(); iC++) {
                    int xN_child     = helpG0.get_cuboid(iC).get_nX();
                    int zN_child     = helpG0.get_cuboid(iC).get_nY();
                    T globPosX_child = helpG0.get_cuboid(iC).get_globPosX();
                    T globPosZ_child = helpG0.get_cuboid(iC).get_globPosY();

                    Cuboid3D<T> child(globPosX_child, globPosY_child, globPosZ_child,
                                   _delta, xN_child, yN_child, zN_child);
                    childrenC.push_back(child);
                }
                globPosY_child += yN_child*_delta;
            }

            splited_nY = _nY - splited_nY;
            restZ = rest/bestIy;

            CuboidGeometry2D<T> helpG1(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restZ);
            yN_child = 0;

            for (int iY=0; iY<bestIy-restY; iY++) {
                yN_child         = (splited_nY+bestIy-restY-iY-1)/(bestIy-restY);
                for (int iC=0; iC<helpG1.get_nC(); iC++) {
                    int xN_child     = helpG1.get_cuboid(iC).get_nX();
                    int zN_child     = helpG1.get_cuboid(iC).get_nY();
                    T globPosX_child = helpG1.get_cuboid(iC).get_globPosX();
                    T globPosZ_child = helpG1.get_cuboid(iC).get_globPosY();

                    Cuboid3D<T> child(globPosX_child, globPosY_child, globPosZ_child,
                               _delta, xN_child, yN_child, zN_child);
                    childrenC.push_back(child);
                }
                globPosY_child += yN_child*_delta;
            }
            return;
        }

        // add in y than in x direction
        else {
            int restX = rest%bestIx;

            // split in two cuboid
            if (restX==0) {
                int restZ = rest/bestIx;
                CuboidGeometry2D<T> helpG(_globPosZ, _globPosY, _delta, _nZ, _nY, bestIz*bestIy+restZ);


                int xN_child = 0;
                T globPosX_child = _globPosX;

                for (int iX=0; iX<bestIx; iX++) {
                    xN_child         = (_nX+bestIx-iX-1)/bestIx;
                    for (int iC=0; iC<helpG.get_nC(); iC++) {
                        int zN_child     = helpG.get_cuboid(iC).get_nX();
                        int yN_child     = helpG.get_cuboid(iC).get_nY();
                        T globPosZ_child = helpG.get_cuboid(iC).get_globPosX();
                        T globPosY_child = helpG.get_cuboid(iC).get_globPosY();

                        Cuboid3D<T> child(globPosZ_child, globPosX_child, globPosY_child,
                                   _delta, zN_child, xN_child, yN_child);
                        childrenC.push_back(child);
                    }
                    globPosX_child += xN_child*_delta;
                }
                return;
            }

            // spilt in four cuboid

            int restZ = rest/bestIx+1;
            int xN_child = 0;
            T globPosX_child = _globPosX;
            int splited_nX = (int) (_nX * (T)((bestIz*bestIy+restZ)*restX)/(T)p);
            CuboidGeometry2D<T> helpG0(_globPosZ, _globPosY, _delta, _nZ, _nY, bestIz*bestIy+restZ);

            for (int iX=0; iX<restX; iX++) {
                xN_child         = (splited_nX+restX-iX-1)/restX;
                for (int iC=0; iC<helpG0.get_nC(); iC++) {
                    int zN_child     = helpG0.get_cuboid(iC).get_nX();
                    int yN_child     = helpG0.get_cuboid(iC).get_nY();
                    T globPosZ_child = helpG0.get_cuboid(iC).get_globPosX();
                    T globPosY_child = helpG0.get_cuboid(iC).get_globPosY();

                    Cuboid3D<T> child(globPosZ_child, globPosX_child, globPosY_child,
                                   _delta, zN_child, xN_child, yN_child);
                    childrenC.push_back(child);
                }
                globPosX_child += xN_child*_delta;
            }

            splited_nX = _nX - splited_nX;
            restZ = rest/bestIx;
            CuboidGeometry2D<T> helpG1(_globPosZ, _globPosY, _delta, _nZ, _nY, bestIz*bestIy+restZ);
            xN_child = 0;

            for (int iX=0; iX<bestIx-restX; iX++) {
                xN_child         = (splited_nX+bestIx-restX-iX-1)/(bestIx-restX);
                for (int iC=0; iC<helpG1.get_nC(); iC++) {
                    int zN_child     = helpG1.get_cuboid(iC).get_nX();
                    int yN_child     = helpG1.get_cuboid(iC).get_nY();
                    T globPosZ_child = helpG1.get_cuboid(iC).get_globPosX();
                    T globPosY_child = helpG1.get_cuboid(iC).get_globPosY();

                    Cuboid3D<T> child(globPosZ_child, globPosX_child, globPosY_child,
                               _delta, zN_child, xN_child, yN_child);
                    childrenC.push_back(child);
                }
                globPosX_child += xN_child*_delta;
            }
            return;
        }
    }
}

}  // namespace olb

#endif