Rotating cuboid not working correctly in 1.8

[stonepreston]

Rotating a cuboid does not produce the expected result in 1.8 but worked fine in 1.7. Its possible I am doing something wrong with the new version of the lib, so any guidance is appreciated

Output in 1.8 (notice how its bowing out at each edge)
https://ibb.co/YT4JYw8k

Correct output in 1.7:
https://ibb.co/NdxDjQ1v

Minimal reproducer in 1.8:

#include <olb.h>

using namespace olb;

using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = descriptors::D2Q9<>;

const T physDeltaX = 0.0078125;  // grid spacing [m]
const T physDeltaT = 0.00078125; // temporal spacing [s]
const T physLength = 10.0;       // length of the squared cuboid [m]
const T physLidVelocity = 10.0;  // velocity imposed on lid [m/s]
const T physViscosity = 0.001;   // kinetic viscosity of fluid [m*m/s]
const T physDensity = 1.0;       // fluid density [kg/(m*m*m)]
const T physMaxT = 30.0;         // maximal simulation time [s]

void prepareGeometry(const UnitConverter<T, DESCRIPTOR> &converter,
                     SuperGeometry<T, 2> &superGeometry) {
    // Set material numbers to assign physics to lattice nodes
    superGeometry.rename(0, 2);
    superGeometry.rename(2, 1, {1, 1});
    Vector<T, 2> extend{3, 3};
    Vector<T, 2> origin{4, 3};
    IndicatorCuboid2D<T> cube(extend, origin, 0.174533);
    superGeometry.rename(1, 2, cube);

    superGeometry.getStatistics().print();
}

int main(int argc, char *argv[]) {
    // === 1st Step: Initialization ===
    initialize(&argc, &argv);
    OstreamManager clout(std::cout, "main");

    // Provide the unit converter the characteristic entities
    const UnitConverter<T, DESCRIPTOR> converter(
        physDeltaX,      // physDeltaX: spacing between two lattice cells in [m]
        physDeltaT,      // physDeltaT: time step in [s]
        physLength,      // charPhysLength: reference length of simulation geometry in [m]
        physLidVelocity, // charPhysVelocity: highest expected velocity during simulation in [m/s]
        physViscosity,   // physViscosity: physical kinematic viscosity in [m^2/s]
        physDensity      // physDensity: physical density [kg/m^3]
    );
    converter.print();

    // === 2nd Step: Prepare Geometry ===
    Vector extend{physLength, physLength};
    Vector origin{0, 0};
    IndicatorCuboid2D<T> cuboid(extend, origin);
    CuboidDecomposition2D<T> cuboidDecomposition(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());
    HeuristicLoadBalancer<T> loadBalancer(cuboidDecomposition);

    SuperGeometry<T, 2> superGeometry(cuboidDecomposition, loadBalancer);
    prepareGeometry(converter, superGeometry);

    SuperVTMwriter2D<T> vtmWriter("1-8-rotate-cuboid");
    SuperGeometryF<T, DESCRIPTOR::d> geometryF(superGeometry);
    geometryF.getName() = "geometry";
    vtmWriter.write(geometryF);
}

Minimal Reproducer in 1.7:

#include "olb2D.h"
#include "olb2D.hh"

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;

using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D2Q9<>;

const T physDeltaX = 0.0078125;  // grid spacing [m]
const T physDeltaT = 0.00078125; // temporal spacing [s]
const T physLength = 10.0;       // length of the squared cuboid [m]
const T physLidVelocity = 10.0;  // velocity imposed on lid [m/s]
const T physViscosity = 0.001;   // kinetic viscosity of fluid [m*m/s]
const T physDensity = 1.0;       // fluid density [kg/(m*m*m)]
const T physMaxT = 30.0;         // maximal simulation time [s]

void prepareGeometry(UnitConverter<T, DESCRIPTOR> const &converter,
                     SuperGeometry<T, 2> &superGeometry) {
    superGeometry.rename(0, 2);
    superGeometry.rename(2, 1, {1, 1});
    Vector<T, 2> extend{3, 3};
    Vector<T, 2> origin{4, 3};
    IndicatorCuboid2D<T> cube(extend, origin, 0.174533);
    superGeometry.rename(1, 2, cube);

    superGeometry.getStatistics().print();
}

int main(int argc, char *argv[]) {

    // === 1st Step: Initialization ===
    olbInit(&argc, &argv);
    OstreamManager clout(std::cout, "main");

    // Provide the unit converter the characteristic entities
    const UnitConverter<T, DESCRIPTOR> converter(
        physDeltaX,      // physDeltaX: spacing between two lattice cells in [m]
        physDeltaT,      // physDeltaT: time step in [s]
        physLength,      // charPhysLength: reference length of simulation geometry in [m]
        physLidVelocity, // charPhysVelocity: highest expected velocity during simulation in [m/s]
        physViscosity,   // physViscosity: physical kinematic viscosity in [m^2/s]
        physDensity      // physDensity: physical density [kg/m^3]
    );
    converter.print();

    // === 2nd Step: Prepare Geometry ===
    Vector<T, 2> extend(physLength, physLength);
    Vector<T, 2> origin(0, 0);
    IndicatorCuboid2D<T> cuboid(extend, origin);

    CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), 1);
    HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
    SuperGeometry<T, 2> superGeometry(cuboidGeometry, loadBalancer);
    prepareGeometry(converter, superGeometry);

    // === 3rd Step: Prepare Lattice ===

    SuperLattice<T, DESCRIPTOR> sLattice(superGeometry);
    SuperVTMwriter2D<T> vtmWriter("1-7-rotate-cuboid");
    SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
    vtmWriter.write(geometry);
}

• This topic was modified 3 weeks, 6 days ago by stonepreston.

[stonepreston]

I did find a new rotation class, but using that does not result in the expected output either:

    superGeometry.rename(0, 2);
    superGeometry.rename(2, 1, {1, 1});
    Vector<T, 2> extend{2, 2};
    Vector<T, 2> origin{4, 3};
    IndicatorCuboid2D<T> cube(extend, origin);
    IndicatorRotate<T, 2> rotated({5, 4}, 0.174533, cube);
    superGeometry.rename(1, 2, rotated);

[stonepreston]

it seems as if the IndicatorRotate method is actually rotating the indicator into the page (i.e rotating from x to z axis). I would expect the rotation to be from the x to y axis. Is there a way to achieve rotation on the x,y plane instead?

[upqdb]

Hello,

the problem does not actually come from wrong rotation. In the current version, both the IndicatorCuboid2D class and the IndicatorRotate class have their bounding boxes not respecting the rotation correctly, this results in the geometry being cut off. This problem was already noticed and fixed, however it is not yet in the release branch.

For now I will provide you the updated getMin and getMax functions for the IndicatorRotate Class.

template <typename S, unsigned D>
Vector<S,D>& IndicatorRotate<S,D>::getMin()
{
  if constexpr (D == 2){
    const S inMin[2] = {_indicator.getMin()[0], _indicator.getMin()[1]};
    const S inMax[2] = {_indicator.getMax()[0], _indicator.getMax()[1]};
    const S radius = util::sqrt(util::pow(inMax[0] - inMin[0], 2) + util::pow(inMax[1] - inMin[1], 2));
    _min[0] = _rotationPoint[0] - radius;
    _min[1] = _rotationPoint[1] - radius;
    return _min;
  }
  else if constexpr (D == 3){
    const S inMin[3] = {_indicator.getMin()[0], _indicator.getMin()[1], _indicator.getMin()[2]};
    const S inMax[3] = {_indicator.getMax()[0], _indicator.getMax()[1], _indicator.getMax()[2]};
    const S radius = util::sqrt(util::pow(inMax[0] - inMin[0], 2) + util::pow(inMax[1] - inMin[1], 2) + util::pow(inMax[2] - inMin[2], 2));
    _min[0] = _rotationPoint[0] - radius;
    _min[1] = _rotationPoint[1] - radius;
    _min[2] = _rotationPoint[2] - radius;
    return _min;
  }
}

template <typename S, unsigned D>
Vector<S,D>& IndicatorRotate<S,D>::getMax()
{
  if constexpr (D == 2){
    const S inMin[2] = {_indicator.getMin()[0], _indicator.getMin()[1]};
    const S inMax[2] = {_indicator.getMax()[0], _indicator.getMax()[1]};

    const S radius = util::sqrt(util::pow(inMax[0] - inMin[0], 2) + util::pow(inMax[1] - inMin[1], 2));
    _max[0] = _rotationPoint[0] + radius;
    _max[1] = _rotationPoint[1] + radius;
    return _max;
  }
  else if constexpr (D == 3){
    const S inMin[3] = {_indicator.getMin()[0], _indicator.getMin()[1], _indicator.getMin()[2]};
    const S inMax[3] = {_indicator.getMax()[0], _indicator.getMax()[1], _indicator.getMax()[2]};
    const S radius = util::sqrt(util::pow(inMax[0] - inMin[0], 2) + util::pow(inMax[1] - inMin[1], 2) + util::pow(inMax[2] - inMin[2], 2));
    _max[0] = _rotationPoint[0] + radius;
    _max[1] = _rotationPoint[1] + radius;
    _max[2] = _rotationPoint[2] + radius;

    return _max;
  }
}

This should fix your problem, please report if you run into any more issues!

[stonepreston]

Can you provide the functions for the IndicatorCuboid2d class as well?

The ones provided for the rotate class do work for fixing the rotation issues, but I am needing to be able to combine indicators using union, intersection, difference. I cant seem to work out how to get the rotate indicator to be able to be combined with the cuboid2ds. This is what I did in 1.7:

  Vector<T,2> straightCubeExtend(6.43, 1.85);
  Vector<T,2> straightCubeOrigin;
  straightCubeOrigin[0] = 2;
  Vector<T,2> tiltedCubeExtend(8.43, 1.27);
  Vector<T,2> tiltedCubeOrigin;
  tiltedCubeOrigin[0] = 1;
  tiltedCubeOrigin[1] = 0;
  Vector<T,2> lowerCubeExtend(6.43, .71);
  Vector<T,2> lowerCubeOrigin;
  lowerCubeOrigin[0] = 2;
  std::shared_ptr<IndicatorCuboid2D<T>> straightCube( new IndicatorCuboid2D<T>(straightCubeExtend, straightCubeOrigin ));
  std::shared_ptr<IndicatorCuboid2D<T>> tiltedCube( new IndicatorCuboid2D<T>(tiltedCubeExtend, tiltedCubeOrigin, 0.174533));
  std::shared_ptr<IndicatorCuboid2D<T>> lowerCube( new IndicatorCuboid2D<T>(lowerCubeExtend, lowerCubeOrigin ));
  IndicatorIdentity2D<T> ramp((straightCube * tiltedCube) + lowerCube);
  // Note: have to rename both the 2 material and the 1 material in this indicator region
  superGeometry.rename( 2,2,ramp );
  superGeometry.rename( 1,2,ramp );

but it seems this does not work when I mix IndicatorCuboids with IndicatorRotates

[stonepreston]

this is what I am trying in 1.8 using indicator rotate instead of rotating the cube on construction:

// ramp
Vector<T, 2> straightCubeExtend(6.43, 1.85);
Vector<T, 2> straightCubeOrigin;
straightCubeOrigin[0] = spillwayGeometry.upperReservoirLength;
Vector<T, 2> tiltedCubeExtend(8.43, 1.27);
Vector<T, 2> tiltedCubeOrigin;
tiltedCubeOrigin[0] = 1;
tiltedCubeOrigin[1] = 0;
Vector<T, 2> lowerCubeExtend(6.43, .71);
Vector<T, 2> lowerCubeOrigin;
lowerCubeOrigin[0] = spillwayGeometry.upperReservoirLength;
std::shared_ptr<IndicatorCuboid2D<T>> straightCube(new IndicatorCuboid2D<T>(straightCubeExtend, straightCubeOrigin));
IndicatorCuboid2D<T> tiltedCube(tiltedCubeExtend, tiltedCubeOrigin);
auto rotX = tiltedCubeOrigin[0] + tiltedCubeExtend[0] / 2.0;
auto rotY = tiltedCubeOrigin[1] + tiltedCubeExtend[1] / 2.0;
Vector<T, 2> rotationPoint(rotX, rotY);
std::shared_ptr<IndicatorRotate<T, 2>> rotatedCube(new IndicatorRotate<T, 2>(rotationPoint, 0.174533, tiltedCube));
std::shared_ptr<IndicatorCuboid2D<T>> lowerCube(new IndicatorCuboid2D<T>(lowerCubeExtend, lowerCubeOrigin));
// * is intersection + is union
IndicatorIdentity2D<T> straightId(straightCube);
IndicatorIdentity2D<T> rotatedId(rotatedCube);
IndicatorIdentity2D<T> lowerId(lowerCube);

IndicatorIdentity2D<T> ramp((straightCube * rotatedCube) + lowerCube);

[upqdb]

Hello,

When using the IndicatorRotate class, at the moment it is not possible to use the union and intersection operators. This, however can be solved by upcasting the IndicatorRotate to an IndicatorF2D when creating the pointer for the rotated geometry. I have added a comment in the following code.

  // ramp
  Vector<T, 2> straightCubeExtend(6.43, 1.85);
  Vector<T, 2> straightCubeOrigin;
  straightCubeOrigin[0] = 2;
  Vector<T, 2> tiltedCubeExtend(8.43, 1.27);
  Vector<T, 2> tiltedCubeOrigin;
  tiltedCubeOrigin[0] = 1;
  tiltedCubeOrigin[1] = 0;
  Vector<T, 2> lowerCubeExtend(6.43, .71);
  Vector<T, 2> lowerCubeOrigin;
  lowerCubeOrigin[0] = 2;
  std::shared_ptr<IndicatorCuboid2D<T>> straightCube(new IndicatorCuboid2D<T>(straightCubeExtend, straightCubeOrigin));
  IndicatorCuboid2D<T> tiltedCube(tiltedCubeExtend, tiltedCubeOrigin);

  auto                 rotX = tiltedCubeOrigin[0] + tiltedCubeExtend[0] / 2.0;
  auto                 rotY = tiltedCubeOrigin[1] + tiltedCubeExtend[1] / 2.0;
  Vector<T, 2>         rotationPoint(rotX, rotY);

  // Use IndicatorF2D as template argument for Pointer, thereby upcasting IndicatorRotate to IndicatorF2D
  std::shared_ptr<IndicatorF2D<T>> rotatedCube(new IndicatorRotate<T, 2>(rotationPoint, 0.174533, tiltedCube));
  std::shared_ptr<IndicatorCuboid2D<T>> lowerCube(new IndicatorCuboid2D<T>(lowerCubeExtend, lowerCubeOrigin));

  // * is intersection + is union
  std::shared_ptr<IndicatorF2D<T>> ramp((straightCube * rotatedCube) + lowerCube);

  // Set Material Numbers
  superGeometry.rename(0, 2, ramp);

  superGeometry.getStatistics().print();

This should result in the same geometry that you created in the 1.7 version.

[stonepreston]

thank you!

[Author]

Posts