Complete Cylinder implemetation and rename scaling function in geometries

apps/gpu/RotatingGrid/RotatingGrid.cpp

@@ -77,7 +77,7 @@ int main()
         // Simulation parameters
         //////////////////////////////////////////////////////////////////////////
         enum RotationOrInterpolation {Rot, Int};
-        const RotationOrInterpolation rotOrInt = Rot;
+        const RotationOrInterpolation rotOrInt = Int;
 
         const std::string path("./output/RotatingGrid");
         const std::string simulationName = rotOrInt == Int ? "RotatingGridInterpolationTest" : "RotatingGrid";
@@ -197,7 +197,7 @@ int main()
 
         Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory, &scalingFactory);
 
-        const std::string gridName = rotOrInt == Rot ? "rot_grid" : "rot";
+        const std::string gridName = rotOrInt == Rot ? "rot_grid" : "grid";
         gridBuilder->writeGridsToVtk(para->getOutputPath() + gridName);
         // NeighborDebugWriter::writeNeighborLinkLinesDebug(para.get());
 

src/gpu/GridGenerator/geometries/Conglomerate/Conglomerate.cpp

@@ -162,13 +162,13 @@ bool Conglomerate::isPointInObject(const double& x1, const double& x2, const dou
     return false;
 }
 
-void Conglomerate::scale(double delta)
+void Conglomerate::changeSizeByDelta(double delta)
 {
     for (uint i = 0; i < numberOfAddObjects; i++)
-        addObjects[i]->scale(delta);
+        addObjects[i]->changeSizeByDelta(delta);
 
     for (uint i = 0; i < numberOfSubtractObjects; i++)
-        subtractObjects[i]->scale(delta);
+        subtractObjects[i]->changeSizeByDelta(delta);
 }
 
 void Conglomerate::findInnerNodes(SPtr<GridImp> grid)

src/gpu/GridGenerator/geometries/Conglomerate/Conglomerate.h

@@ -63,7 +63,7 @@ public:
     double getX3Minimum() override;
     double getX3Maximum() override;
 
-    void scale(double delta) override;
+    void changeSizeByDelta(double delta) override;
 
     bool isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset, const double& maxOffset) override;
 

src/gpu/GridGenerator/geometries/Cuboid/Cuboid.cpp

@@ -131,7 +131,7 @@ bool Cuboid::isBetween(const real& coord, const real& start, const real& end)
 }
 
 
-void Cuboid::scale(double delta)
+void Cuboid::changeSizeByDelta(double delta)
 {
     this->minX1 -= delta;
     this->minX2 -= delta;

src/gpu/GridGenerator/geometries/Cuboid/Cuboid.h

@@ -54,7 +54,7 @@ public:
     double getX3Minimum() override;
     double getX3Maximum() override;
 
-    void scale(double delta) override;
+    void changeSizeByDelta(double delta) override;
 
     bool isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset, const double& maxOffset) override;
 

src/gpu/GridGenerator/geometries/Cylinder/Cylinder.cpp

@@ -80,14 +80,49 @@ double Cylinder::getX3Maximum()
     return getMaximumCoordinate(z);
 }
 
-////////////////// #TODO
+double Cylinder::getRadius() const
+{
+    return radius;
+}
+
+double Cylinder::getHeight() const
+{
+    return height;
+}
 
-bool Cylinder::isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset, const double& maxOffset)
+bool Cylinder::isInCircle(double delta1, double delta2, double offset) const
 {
- return false;
+    return (delta1 * delta1 + delta2 * delta2) < ((this->radius - offset) * (this->radius - offset));
 }
 
+bool Cylinder::isPointInObject(const double &x1, const double &x2, const double &x3, const double &minOffset,
+                               const double &maxOffset)
+{
+    double offset = maxOffset;
+    if (x1 < center.at(x) || x2 < center.at(y) || x3 < center.at(z)) offset = minOffset;
+
+    const double deltaX1 = x1 - center.at(x);
+    const double deltaX2 = x2 - center.at(y);
+    const double deltaX3 = x3 - center.at(z);
+
+    switch (principalAxis) {
+        case x:
+            if (deltaX1 > 0.5 * height || deltaX1 < -0.5 * height) return false;
+            return isInCircle(deltaX2, deltaX3, offset);
+        case y:
+            if (deltaX2 > 0.5 * height || deltaX2 < -0.5 * height) return false;
+            return isInCircle(deltaX1, deltaX3, offset);
+        case z:
+            if (deltaX3 > 0.5 * height || deltaX3 < -0.5 * height) return false;
+            return isInCircle(deltaX1, deltaX2, offset);
+    }
+
+    VF_LOG_CRITICAL("Unknown principal axis in Cylinder.");
+    return false;
+}
 
-void Cylinder::scale(double delta)
+void Cylinder::changeSizeByDelta(double delta)
 {
+    this->radius += delta;
+    this->height += 2 * delta;
 }
\ No newline at end of file

src/gpu/GridGenerator/geometries/Cylinder/Cylinder.h

@@ -63,15 +63,20 @@ public:
     double getX3Minimum() override;
     double getX3Maximum() override;
 
+    double getRadius() const;
+    double getHeight() const;
+
     bool isPointInObject(const double &x1, const double &x2, const double &x3, const double &minOffset,
                          const double &maxOffset) override;
-    void scale(double delta) override;
+    void changeSizeByDelta(double delta) override;
 
 private:
     double getCentroidCoordinate(PrincipalAxis coordinateDirection) const;
     double getMinimunCoordinate(PrincipalAxis coordinateDirection) const;
     double getMaximumCoordinate(PrincipalAxis coordinateDirection) const;
 
+    bool isInCircle(double delta1, double delta2, double offset) const;
+
     const std::map<PrincipalAxis, std::array<double, 3>> unitVectors{ { x, { 1, 0, 0 } },
                                                                       { y, { 0, 1, 0 } },
                                                                       { z, { 0, 0, 1 } } };

src/gpu/GridGenerator/geometries/Cylinder/CylinderTest.cpp

@@ -2,11 +2,37 @@
 
 #include "Cylinder.h"
 
+
+void testChangeSizeByDelta(Cylinder &cylinder)
+{
+    const double delta = 0.2;
+    const double tolerance = 0.0000000001;
+    const double oldRadius = cylinder.getRadius();
+    const double oldHeight = cylinder.getHeight();
+    const double oldMinimumX1 = cylinder.getX1Minimum();
+    const double oldMinimumX2 = cylinder.getX2Minimum();
+    const double oldMinimumX3 = cylinder.getX3Minimum();
+    const double oldMaximumX1 = cylinder.getX1Maximum();
+    const double oldMaximumX2 = cylinder.getX2Maximum();
+    const double oldMaximumX3 = cylinder.getX3Maximum();
+
+    cylinder.changeSizeByDelta(delta);
+
+    EXPECT_THAT(cylinder.getRadius(), testing::Eq(oldRadius + delta));
+    EXPECT_THAT(cylinder.getHeight(), testing::Eq(oldHeight + 2 * delta));
+    EXPECT_THAT(cylinder.getX1Minimum(), testing::DoubleNear(oldMinimumX1 - delta, tolerance));
+    EXPECT_THAT(cylinder.getX2Minimum(), testing::DoubleNear(oldMinimumX2 - delta, tolerance));
+    EXPECT_THAT(cylinder.getX3Minimum(), testing::DoubleNear(oldMinimumX3 - delta, tolerance));
+    EXPECT_THAT(cylinder.getX1Maximum(), testing::DoubleNear(oldMaximumX1 + delta, tolerance));
+    EXPECT_THAT(cylinder.getX2Maximum(), testing::DoubleNear(oldMaximumX2 + delta, tolerance));
+    EXPECT_THAT(cylinder.getX3Maximum(), testing::DoubleNear(oldMaximumX3 + delta, tolerance));
+}
+
 // CylinderTestAxisNormalToX ////////////////////////////////////////////////
 
 class CylinderTestAxisNormalToX : public testing::Test
 {
-protected:
+public:
     std::array<double, 3> center = { 0.1, 0.2, 0.3 };
     double radius = 2.0;
     double height = 8.0;
@@ -34,6 +60,35 @@ TEST_F(CylinderTestAxisNormalToX, getMaximum)
     EXPECT_THAT(cylinder.getX3Maximum(), testing::Eq(radius + center[2]));
 }
 
+TEST_F(CylinderTestAxisNormalToX, isPointInObject)
+{
+    double epsilon = 0.0001;
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2], 0.0, 0.0));
+
+    // x
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] - 0.5 * height + epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] - 0.5 * height - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] + 0.5 * height - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] + 0.5 * height + epsilon, center[1], center[2], 0.0, 0.0));
+
+    // y
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] - radius + epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] - radius - epsilon, center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] + radius - epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] + radius + epsilon, center[2], 0.0, 0.0));
+
+    // z
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] - radius + epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] - radius - epsilon, 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] + radius - epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] + radius + epsilon, 0.0, 0.0));
+}
+
+TEST_F(CylinderTestAxisNormalToX, changeSizeByDelta)
+{
+    testChangeSizeByDelta(cylinder);
+}
+
 // CylinderTestAxisNormalToY ////////////////////////////////////////////////
 
 class CylinderTestAxisNormalToY : public testing::Test
@@ -66,6 +121,35 @@ TEST_F(CylinderTestAxisNormalToY, getMaximum)
     EXPECT_THAT(cylinder.getX3Maximum(), testing::Eq(radius + center[2]));
 }
 
+TEST_F(CylinderTestAxisNormalToY, isPointInObject)
+{
+    double epsilon = 0.0001;
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2], 0.0, 0.0));
+
+    // x
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] - radius + epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] - radius - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] + radius - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] + radius + epsilon, center[1], center[2], 0.0, 0.0));
+
+    // y
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] - 0.5 * height + epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] - 0.5 * height - epsilon, center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] + 0.5 * height - epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] + 0.5 * height + epsilon, center[2], 0.0, 0.0));
+
+    // z
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] - radius + epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] - radius - epsilon, 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] + radius - epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] + radius + epsilon, 0.0, 0.0));
+}
+
+TEST_F(CylinderTestAxisNormalToY, changeSizeByDelta)
+{
+    testChangeSizeByDelta(cylinder);
+}
+
 // CylinderTestAxisNormalToZ ////////////////////////////////////////////////
 
 class CylinderTestAxisNormalToZ : public testing::Test
@@ -96,4 +180,33 @@ TEST_F(CylinderTestAxisNormalToZ, getMaximum)
     EXPECT_THAT(cylinder.getX1Maximum(), testing::Eq(radius + center[0]));
     EXPECT_THAT(cylinder.getX2Maximum(), testing::Eq(radius + center[1]));
     EXPECT_THAT(cylinder.getX3Maximum(), testing::Eq(0.5 * height + center[2]));
-}
\ No newline at end of file
+}
+
+TEST_F(CylinderTestAxisNormalToZ, isPointInObject)
+{
+    double epsilon = 0.0001;
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2], 0.0, 0.0));
+
+    // x
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] - radius + epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] - radius - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] + radius - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] + radius + epsilon, center[1], center[2], 0.0, 0.0));
+
+    // y
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] - radius + epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] - radius - epsilon, center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] + radius - epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] + radius + epsilon, center[2], 0.0, 0.0));
+
+    // z
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] - 0.5 * height + epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] - 0.5 * height - epsilon, 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] + 0.5 * height - epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] + 0.5 * height + epsilon, 0.0, 0.0));
+}
+
+TEST_F(CylinderTestAxisNormalToZ, changeSizeByDelta)
+{
+    testChangeSizeByDelta(cylinder);
+}

src/gpu/GridGenerator/geometries/Object.h

@@ -59,7 +59,7 @@ public:
     virtual double getX3Maximum()  = 0;
 
 
-    virtual void scale(double delta) = 0;
+    virtual void changeSizeByDelta(double delta) = 0;
 
 
     virtual bool isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset, const double& maxOffset) = 0;

src/gpu/GridGenerator/geometries/Sphere/Sphere.cpp

@@ -115,7 +115,7 @@ bool Sphere::isPointInObject(const double& x1, const double& x2, const double& x
 }
 
 
-void Sphere::scale(double delta)
+void Sphere::changeSizeByDelta(double delta)
 {
     this->radius += delta;
 }

src/gpu/GridGenerator/geometries/Sphere/Sphere.h

@@ -59,7 +59,7 @@ public:
     bool isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset, const double& maxOffset) override;
 
 
-    void scale(double delta) override;
+    void changeSizeByDelta(double delta) override;
     
     int getIntersection(const Vertex &P, const Vertex &direction, Vertex &pointOnObject, real &qVal) override;
 

src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.cpp

@@ -185,7 +185,7 @@ bool intersectPlane(const Vertex &normal, const Vertex &pointOnPlane, const Vert
     return false;
 }
 
-void TriangularMesh::scale(double offset)
+void TriangularMesh::changeSizeByDelta(double offset)
 {
     std::vector<Triangle> triangles = this->triangleVec;
 

src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.h

@@ -97,7 +97,7 @@ public:
     double getX3Centroid() override { throw "Not implemented in TriangularMesh"; }
     double getX3Minimum() override { return minmax.minZ; }
     double getX3Maximum() override { return minmax.maxZ; }
-    void scale(double delta) override;
+    void changeSizeByDelta(double delta) override;
     bool isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset,
         const double& maxOffset) override {
         return false;

src/gpu/GridGenerator/geometries/VerticalCylinder/VerticalCylinder.cpp

@@ -111,10 +111,10 @@ bool VerticalCylinder::isPointInObject(const double& x1, const double& x2, const
 }
 
 
-void VerticalCylinder::scale(double delta)
+void VerticalCylinder::changeSizeByDelta(double delta)
 {
     this->radius += delta;
-    this->height += delta;
+    this->height += 2 * delta;
 }
 
 double VerticalCylinder::getRadius() const

src/gpu/GridGenerator/geometries/VerticalCylinder/VerticalCylinder.h

@@ -62,7 +62,7 @@ public:
     bool isPointInObject(const double& x1, const double& x2, const double& x3, const double& minOffset, const double& maxOffset) override;
 
 
-    void scale(double delta) override;
+    void changeSizeByDelta(double delta) override;
 
 protected:
     double centerX;

src/gpu/GridGenerator/geometries/VerticalCylinder/VerticalCylinderTest.cpp

@@ -32,4 +32,26 @@ TEST_F(VerticalCylinderTest, getMaximum)
     EXPECT_THAT(cylinder.getX3Maximum(), testing::Eq(0.5 * height + center[2]));
 }
 
-// TEST_F()
\ No newline at end of file
+TEST_F(VerticalCylinderTest, isPointInObject)
+{
+    double epsilon = 0.0001;
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2], 0.0, 0.0));
+
+    // x
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] - radius + epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] - radius - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0] + radius - epsilon, center[1], center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0] + radius + epsilon, center[1], center[2], 0.0, 0.0));
+
+    // y
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] - radius + epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] - radius - epsilon, center[2], 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1] + radius - epsilon, center[2], 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1] + radius + epsilon, center[2], 0.0, 0.0));
+
+    // z
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] - 0.5 * height + epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] - 0.5 * height - epsilon, 0.0, 0.0));
+    EXPECT_TRUE(cylinder.isPointInObject(center[0], center[1], center[2] + 0.5 * height - epsilon, 0.0, 0.0));
+    EXPECT_FALSE(cylinder.isPointInObject(center[0], center[1], center[2] + 0.5 * height + epsilon, 0.0, 0.0));
+}
\ No newline at end of file

src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp

@@ -86,7 +86,7 @@ void MultipleGridBuilder::addGeometry(SPtr<Object> solidObject, uint level)
 {
     this->solidObject = solidObject;
     auto gridShape = solidObject->clone();
-    gridShape->scale(4.0);
+    gridShape->changeSizeByDelta(4.0);
 
     this->addGrid(gridShape, level);
 }
@@ -169,7 +169,7 @@ void MultipleGridBuilder::addIntermediateGridsToList(uint levelDifference, uint
         {
             const real scalingFactor = nodesBetweenGrids * spacings[i] * calculateDelta(levelFine);
             SPtr<Object> gridShapeClone = gridShape->clone();
-            gridShapeClone->scale(scalingFactor);
+            gridShapeClone->changeSizeByDelta(scalingFactor);
 
             const auto grid = makeGrid(gridShapeClone, level++, 0);
             grids.push_back(grid);
@@ -267,8 +267,12 @@ SPtr<Grid> MultipleGridBuilder::makeRotatingGrid(SPtr<VerticalCylinder> cylinder
                                                               staggeredCoordinates[4],
                                                               staggeredCoordinates[5], delta, level);
 
-    cylinder = std::make_shared<VerticalCylinder>(cylinder->getX1Centroid(), cylinder->getX2Centroid(), cylinder->getX3Centroid(), cylinder->getRadius()+5*delta, cylinder->getHeight()+6*delta);
-    VF_LOG_WARNING("RotGr: LargerCylinder Size x : {}, {} \n z: {}, {}", cylinder->getX1Minimum(), cylinder->getX1Maximum(), cylinder->getX3Minimum(), cylinder->getX3Maximum());
+    // #TODO: change next line to use staggering
+    cylinder = std::make_shared<VerticalCylinder>(cylinder->getX1Centroid(), cylinder->getX2Centroid(),
+                                                  cylinder->getX3Centroid(), cylinder->getRadius() + 6.0 * delta,
+                                                  cylinder->getHeight() + 6.0 * delta);
+    VF_LOG_WARNING("RotGr: LargerCylinder Size x : {}, {} \n z: {}, {}", cylinder->getX1Minimum(), cylinder->getX1Maximum(),
+                   cylinder->getX3Minimum(), cylinder->getX3Maximum());
 
     SPtr<Grid> newGrid = this->makeGrid(cylinder, staggeredCoordinates[0],
                                                   staggeredCoordinates[1],