From ba95676b4a7eae51317c027ac0828fa4b55fdda0 Mon Sep 17 00:00:00 2001
From: Soeren Peters <peters@irmb.tu-bs.de>
Date: Tue, 29 May 2018 14:47:21 +0200
Subject: [PATCH] - find all Qs

---
 .../grid/BoundaryConditions/Side.cpp          |  14 +-
 src/GridGenerator/grid/Field.cu               |   2 +-
 .../grid/GridBuilder/GridBuilder.h            |   2 +-
 .../grid/GridBuilder/LevelGridBuilder.cpp     | 144 ++++++++++++++----
 .../grid/GridBuilder/LevelGridBuilder.h       |   6 +-
 src/GridGenerator/grid/GridImp.cu             |  51 ++++++-
 src/GridGenerator/grid/GridImp.h              |   1 +
 src/GridGenerator/grid/GridTest.cpp           |  15 ++
 targets/apps/HULC/main.cpp                    |   4 +-
 9 files changed, 190 insertions(+), 49 deletions(-)

diff --git a/src/GridGenerator/grid/BoundaryConditions/Side.cpp b/src/GridGenerator/grid/BoundaryConditions/Side.cpp
index 413f5cb58..df51adefa 100644
--- a/src/GridGenerator/grid/BoundaryConditions/Side.cpp
+++ b/src/GridGenerator/grid/BoundaryConditions/Side.cpp
@@ -15,7 +15,7 @@ void Side::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition
             const uint index = getIndex(grid, coord, constant, v1, v2);
             if (grid->getFieldEntry(index) == FLUID)
             {
-                boundaryCondition->indices.push_back(grid->getSparseIndex(index) + 1);
+                boundaryCondition->indices.push_back(index);
                 setPressureNeighborIndices(boundaryCondition, grid, index);
             }
         }
@@ -42,8 +42,7 @@ void Side::setPressureNeighborIndices(SPtr<BoundaryCondition> boundaryCondition,
             nz = boundaryCondition->side->getDirection() * grid->getDelta() + z;
 
         int neighborIndex = grid->transCoordToIndex(nx, ny, nz);
-        int sparseIndex = grid->getSparseIndex(neighborIndex);
-        pressureBoundaryCondition->neighborIndices.push_back(sparseIndex);
+        pressureBoundaryCondition->neighborIndices.push_back(neighborIndex);
     }
 }
 
@@ -68,6 +67,9 @@ void Geometry::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondi
 
     for (int i = 0; i < grid->getSize(); i++)
     {
+        if (grid->getFieldEntry(i) != Q)
+            continue;
+
         for (int dir = 0; dir < grid->getEndDirection(); dir++)
         {
             const int qIndex = dir * grid->getSize() + i;
@@ -76,13 +78,15 @@ void Geometry::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondi
             qNode[dir] = q;
             if (q > 0)
                 qFound = true;
+            else
+                qNode[dir] = -1.0;
+
 
         }
 
         if (qFound)
         {
-            const int sparseIndex = grid->getSparseIndex(i);
-            geometryBoundaryCondition->indices.push_back(sparseIndex);
+            geometryBoundaryCondition->indices.push_back(i);
             geometryBoundaryCondition->qs.push_back(qNode);
         }
 
diff --git a/src/GridGenerator/grid/Field.cu b/src/GridGenerator/grid/Field.cu
index 9d95db80b..6c2f533ff 100644
--- a/src/GridGenerator/grid/Field.cu
+++ b/src/GridGenerator/grid/Field.cu
@@ -63,7 +63,7 @@ HOSTDEVICE bool Field::isFineToCoarseNode(uint index) const
 HOSTDEVICE bool Field::isFluid(uint index) const
 {
     const char type = field[index];
-    return type == FLUID || type == FLUID_CFC || type == FLUID_CFF || type == FLUID_FCC || type == FLUID_FCF;
+    return type == FLUID || type == FLUID_CFC || type == FLUID_CFF || type == FLUID_FCC || type == FLUID_FCF || type == Q;
 }
 
 HOSTDEVICE bool Field::isSolid(uint index) const
diff --git a/src/GridGenerator/grid/GridBuilder/GridBuilder.h b/src/GridGenerator/grid/GridBuilder/GridBuilder.h
index b7fff238a..9a41a5db5 100644
--- a/src/GridGenerator/grid/GridBuilder/GridBuilder.h
+++ b/src/GridGenerator/grid/GridBuilder/GridBuilder.h
@@ -42,7 +42,7 @@ public:
     virtual VF_PUBLIC uint getNumberOfGridLevels() = 0;
 
 
-    virtual void writeArrows(std::string fileName, std::shared_ptr<ArrowTransformator> trans) const = 0;
+    virtual void writeArrows(std::string fileName) const = 0;
 
 	virtual SPtr<Grid> getGrid(uint level) = 0;
 
diff --git a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
index 9d9558077..d4c45aeb6 100644
--- a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
+++ b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
@@ -30,6 +30,8 @@
 #include "grid/BoundaryConditions/BoundaryCondition.h"
 #include "grid/BoundaryConditions/Side.h"
 
+#include <basics/utilities/UbTuple.h>
+
 
 #define GEOFLUID 19
 #define GEOSOLID 16
@@ -264,7 +266,7 @@ void LevelGridBuilder::getVelocityValues(real* vx, real* vy, real* vz, int* indi
     {
         for(int i = 0; i < boundaryCondition->indices.size(); i++)
         {
-            indices[allIndicesCounter] = boundaryCondition->indices[i];  
+            indices[allIndicesCounter] = grids[level]->getSparseIndex(boundaryCondition->indices[i]) +1;  
 
             vx[allIndicesCounter] = boundaryCondition->vx;
             vy[allIndicesCounter] = boundaryCondition->vy;
@@ -291,9 +293,9 @@ void LevelGridBuilder::getPressureValues(real* rho, int* indices, int* neighborI
     {
         for (int i = 0; i < boundaryCondition->indices.size(); i++)
         {
-            indices[allIndicesCounter] = boundaryCondition->indices[i];
+            indices[allIndicesCounter] = grids[level]->getSparseIndex(boundaryCondition->indices[i]) + 1;
 
-            neighborIndices[allIndicesCounter] = boundaryCondition->neighborIndices[i];;
+            neighborIndices[allIndicesCounter] = grids[level]->getSparseIndex(boundaryCondition->neighborIndices[i]) + 1;
 
             rho[allIndicesCounter] = boundaryCondition->rho;
             allIndicesCounter++;
@@ -354,7 +356,7 @@ void LevelGridBuilder::getGeometryIndices(int* indices, int level) const
 {
     for (uint i = 0; i < geometryBoundaryCondition->indices.size(); i++)
     {
-        indices[i] = geometryBoundaryCondition->indices[i];
+        indices[i] = grids[level]->getSparseIndex(geometryBoundaryCondition->indices[i]) + 1;
     }
 }
 
@@ -485,36 +487,128 @@ void LevelGridBuilder::fillRBForNode(int index, int direction, int directionSign
     qNode.clear();
 }
 
-void LevelGridBuilder::writeArrows(std::string fileName, std::shared_ptr<ArrowTransformator> trans) const
+
+
+
+
+
+#include <fstream>
+using namespace std;
+void writeLines(std::string filename, std::vector<UbTupleFloat3> nodes, std::vector<UbTupleInt2> lines)
 {
-    Grid* grid = grids[0].get();
+    string vtkfilename = filename + ".bin.vtu";
+
+    ofstream out(vtkfilename.c_str(), ios::out | ios::binary);
+
+    int nofNodes = (int)nodes.size();
+    int nofCells = (int)lines.size();
+
+    int bytesPerByteVal = 4; //==sizeof(int)
+    int bytesPoints = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
+    int bytesCellConnectivty = 2 /*nodes per line */ * nofCells * sizeof(int);
+    int bytesCellOffsets = 1 /*offset per line */ * nofCells * sizeof(int);
+    int bytesCellTypes = 1 /*type of line */ * nofCells * sizeof(unsigned char);
+
+    int offset = 0;
+    //VTK FILE
+    out << "<?xml version=\"1.0\"?>\n";
+    out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" >" << "\n";
+    out << "   <UnstructuredGrid>" << "\n";
+    out << "      <Piece NumberOfPoints=\"" << nofNodes << "\" NumberOfCells=\"" << nofCells << "\">\n";
+
+    //POINTS SECTION
+    out << "         <Points>\n";
+    out << "            <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"appended\" offset=\"" << offset << "\"  />\n";
+    out << "         </Points>\n";
+    offset += (bytesPerByteVal + bytesPoints);
+
+    //CELLS SECTION
+    out << "         <Cells>\n";
+    out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset << "\" />\n";
+    offset += (bytesPerByteVal + bytesCellConnectivty);
+    out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset << "\" />\n";
+    offset += (bytesPerByteVal + bytesCellOffsets);
+    out << "            <DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\"" << offset << "\" />\n ";
+    offset += (bytesPerByteVal + bytesCellTypes);
+    out << "         </Cells>\n";
+
+    out << "      </Piece>\n";
+    out << "   </UnstructuredGrid>\n";
+
+    // AppendedData SECTION
+    out << "   <AppendedData encoding=\"raw\">\n";
+    out << "_";
+
+    //POINTS SECTION
+    out.write((char*)&bytesPoints, bytesPerByteVal);
+    for (int n = 0; n < nofNodes; n++)
+    {
+        out.write((char*)&val<1>(nodes[n]), sizeof(float));
+        out.write((char*)&val<2>(nodes[n]), sizeof(float));
+        out.write((char*)&val<3>(nodes[n]), sizeof(float));
+    }
 
-    //std::shared_ptr<PolyDataWriterWrapper> writer = std::shared_ptr<PolyDataWriterWrapper>(new PolyDataWriterWrapper());
-    for (int index = 0; index < Qs[GEOMQS].size(); index++)
+    //CELLS SECTION
+    //cellConnectivity
+    out.write((char*)&bytesCellConnectivty, bytesPerByteVal);
+    for (int c = 0; c < nofCells; c++)
     {
-        Vertex startNode = getVertex(getMatrixIndex(index));
-        //for (int qi = grid.d.dir_start; qi <= grid.d.dir_end; qi++)
-            //writeArrow(index, qi, startNode, trans, writer);
+        out.write((char*)&val<1>(lines[c]), sizeof(int));
+        out.write((char*)&val<2>(lines[c]), sizeof(int));
+
     }
-    //writer->writePolyDataToFile(fileName);
+
+    //cellOffsets
+    out.write((char*)&bytesCellOffsets, bytesPerByteVal);
+    int itmp;
+    for (int c = 1; c <= nofCells; c++)
+    {
+        itmp = 2 * c;
+        out.write((char*)&itmp, sizeof(int));
+    }
+
+    //cellTypes
+    out.write((char*)&bytesCellTypes, bytesPerByteVal);
+    unsigned char vtkCellType = 3;
+    for (int c = 0; c < nofCells; c++)
+    {
+        out.write((char*)&vtkCellType, sizeof(unsigned char));
+    }
+    out << "\n</AppendedData>\n";
+    out << "</VTKFile>";
+    out << endl;
+    out.close();
 }
 
-void LevelGridBuilder::writeArrow(const int i, const int qi, const Vertex& startNode, std::shared_ptr<const ArrowTransformator> trans/*, std::shared_ptr<PolyDataWriterWrapper> writer*/) const
+void LevelGridBuilder::writeArrows(std::string fileName) const
 {
-    Grid* grid = grids[0].get();
+    std::vector<UbTupleFloat3> nodes;
+    std::vector<UbTupleInt2> cells;
 
-    real qval = Qs[GEOMQS][i][qi + 1];
-    if (qval > 0.0f)
+    int actualNodeNumber = 0;
+    for (int index = 0; index < geometryBoundaryCondition->indices.size(); index++)
     {
-        int qReverse = grid->getEndDirection() - qi;
-        Vertex dir((real)grid->getDirection()[qReverse * DIMENSION + 0], (real)grid->getDirection()[qReverse* DIMENSION + 1], (real)grid->getDirection()[qReverse * DIMENSION + 2]);
-        Vertex nodeOnGeometry(startNode + (dir * qval));
-        std::shared_ptr<Arrow> arrow = ArrowImp::make(startNode, nodeOnGeometry);
-        trans->transformGridToWorld(arrow);
-        //writer->addVectorArrow(arrow);
+        Vertex startNode = getVertex(geometryBoundaryCondition->indices[index]);
+        for (int qi = 0; qi <= 26; qi++)
+        {
+            real qval = geometryBoundaryCondition->qs[index][qi];
+            if (qval > 0.0f)
+            {
+                Vertex dir((real)grids[0]->getDirection()[qi * DIMENSION + 0], (real)grids[0]->getDirection()[qi* DIMENSION + 1], (real)grids[0]->getDirection()[qi * DIMENSION + 2]);
+                Vertex nodeOnGeometry(startNode + (dir * qval));
+
+                nodes.push_back(makeUbTuple(float(startNode.x), float(startNode.y), float(startNode.z)));
+                nodes.push_back(makeUbTuple(float(nodeOnGeometry.x), float(nodeOnGeometry.y), float(nodeOnGeometry.z)));
+                actualNodeNumber += 2;
+                cells.push_back(makeUbTuple(actualNodeNumber - 2, actualNodeNumber - 1));
+            }
+        }
     }
+
+    writeLines(fileName, nodes, cells);
 }
 
+
 Vertex LevelGridBuilder::getVertex(int matrixIndex) const
 {
     real x, y, z;
@@ -522,9 +616,3 @@ Vertex LevelGridBuilder::getVertex(int matrixIndex) const
     return Vertex(x,y,z);
 }
 
-int LevelGridBuilder::getMatrixIndex(int i) const
-{
-    int index = (int)Qs[GEOMQS][i][0];
-    return this->grids[0]->getSparseIndex(index);
-}
-
diff --git a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
index eb7e43f9a..635e06033 100644
--- a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
+++ b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
@@ -83,7 +83,7 @@ public:
 
     VF_PUBLIC virtual void setPressValues(real* RhoBC, int* kN, int channelSide, int level) const;
 
-    VF_PUBLIC void writeArrows(std::string fileName, std::shared_ptr<ArrowTransformator> trans) const;
+    VF_PUBLIC void writeArrows(std::string fileName) const;
 
 protected:
 
@@ -110,11 +110,7 @@ protected:
     void addQsToVector(int index);
     void fillRBForNode(int index, int direction, int directionSign, int rb);
 
-    int getMatrixIndex(const int i) const;
     Vertex getVertex(const int matrixIndex) const;
-    void writeArrow(const int i, const int qi, const Vertex& startNode,
-                    std::shared_ptr<const ArrowTransformator> trans/*, std::shared_ptr<PolyDataWriterWrapper> writer*/)
-    const;
 
 private:
     Device device;
diff --git a/src/GridGenerator/grid/GridImp.cu b/src/GridGenerator/grid/GridImp.cu
index c2fffa281..057457ad1 100644
--- a/src/GridGenerator/grid/GridImp.cu
+++ b/src/GridGenerator/grid/GridImp.cu
@@ -208,7 +208,46 @@ HOSTDEVICE bool GridImp::isValidInnerStopper(uint index) const
         || this->field.is(index, SOLID) && (nodeInPreviousCellIs(index, FLUID) || nodeInPreviousCellIs(index, FLUID_CFF));
 }
 
+HOSTDEVICE bool GridImp::hasNeighbor(uint index, char type) const
+{
+    real x, y, z;
+    this->transIndexToCoords(index, x, y, z);
+
+    const real neighborX = x + this->delta > endX ? endX : x + this->delta;
+    const real neighborY = y + this->delta > endY ? endY : y + this->delta;
+    const real neighborZ = z + this->delta > endZ ? endZ : z + this->delta;
+
+    const real neighborMinusX = x - this->delta < startX ? startX : x - this->delta;
+    const real neighborMinusY = y - this->delta < startY ? startY : y - this->delta;
+    const real neighborMinusZ = z - this->delta < startZ ? startZ : z - this->delta;
+
+
+    const uint indexMXY = transCoordToIndex(-neighborX, neighborY, z);
+    const uint indexMYZ = transCoordToIndex(x, -neighborY, neighborZ);
+    const uint indexMXZ = transCoordToIndex(-neighborX, y, neighborZ);
 
+    const uint indexXMY = transCoordToIndex(neighborX, -neighborY, z);
+    const uint indexYMZ = transCoordToIndex(x, neighborY, -neighborZ);
+    const uint indexXMZ = transCoordToIndex(neighborX, y, -neighborZ);
+
+    const uint indexMXYMZ = transCoordToIndex(-neighborX, neighborY, -neighborZ);
+    const uint indexMXYZ  = transCoordToIndex(-neighborX, neighborY, neighborZ);
+    const uint indexMXMYZ = transCoordToIndex(-neighborX, -neighborY, neighborZ);
+    const uint indexXMYMZ = transCoordToIndex(neighborX, -neighborY, -neighborZ);
+    const uint indexXMYZ  = transCoordToIndex(neighborX, -neighborY, neighborZ);
+    const uint indexXYMZ  = transCoordToIndex(neighborX, neighborY, -neighborZ);
+
+
+
+    return nodeInNextCellIs(index, type) || nodeInPreviousCellIs(index, type) || indexMXY || indexMYZ || indexMXZ || indexXMY || indexYMZ || indexXMZ 
+        || indexMXYMZ
+        || indexMXYZ
+        || indexMXMYZ
+        || indexXMYMZ
+        || indexXMYZ
+        || indexXYMZ;
+
+}
 
 HOSTDEVICE bool GridImp::nodeInNextCellIs(int index, char type) const
 {
@@ -596,13 +635,11 @@ HOSTDEVICE void GridImp::findQs(Triangle &triangle)
                     continue;
 
                 const Vertex point(x, y, z);
-                const int value = triangle.isUnderFace(point);
 
-                if (value == Q)
+                if(hasNeighbor(index, STOPPER_OVERLAP_GRID))
                 {
                     field.setFieldEntry(index, Q);
                     calculateQs(point, triangle);
-
                 }
             }
         }
@@ -635,13 +672,13 @@ HOSTDEVICE void GridImp::calculateQs(const Vertex &point, const Triangle &triang
 
         error = triangle.getTriangleIntersection(point, direction, pointOnTriangle, subdistance);
 
-        if (error != 0 && subdistance <= 1.0f)
+        if (error != 0 && subdistance < 1.0 && subdistance > 0.0)
         {
             //solid_node = VertexInteger(actualPoint.x + direction.x, actualPoint.y + direction.y, actualPoint.z + direction.z);
             distribution.f[i*size + transCoordToIndex(point.x, point.y, point.z)] = subdistance;
             //printf("Q%d %d: %2.8f \n", i, grid.transCoordToIndex(actualPoint), grid.d.f[index]);
-        } else
-            distribution.f[i*size + transCoordToIndex(point.x, point.y, point.z)] = -1.0;
+        } /*else
+            distribution.f[i*size + transCoordToIndex(point.x, point.y, point.z)] = -1.0;*/
     }
 }
 
@@ -717,7 +754,7 @@ HOSTDEVICE real GridImp::getMaximumOnNodes(const real& maxExact, const real& dec
     real maxNode = ceil(maxExact - 1.0);
     maxNode += decimalStart;
 
-    while (maxNode < maxExact)
+    while (maxNode <= maxExact)
         maxNode += delta;
     return maxNode;
 }
diff --git a/src/GridGenerator/grid/GridImp.h b/src/GridGenerator/grid/GridImp.h
index f5aa081de..ffe7c9f23 100644
--- a/src/GridGenerator/grid/GridImp.h
+++ b/src/GridGenerator/grid/GridImp.h
@@ -104,6 +104,7 @@ public:
     HOSTDEVICE void setNodeTo(uint index, char type);
     HOSTDEVICE bool isNode(uint index, char type) const;
     HOSTDEVICE bool nodeInNextCellIs(int index, char type) const;
+    HOSTDEVICE bool hasNeighbor(uint index, char type)const;
 
     HOSTDEVICE Field getField() const;
     HOSTDEVICE char getFieldEntry(uint index) const override;
diff --git a/src/GridGenerator/grid/GridTest.cpp b/src/GridGenerator/grid/GridTest.cpp
index bcf0c4b92..166acacc2 100644
--- a/src/GridGenerator/grid/GridTest.cpp
+++ b/src/GridGenerator/grid/GridTest.cpp
@@ -55,7 +55,22 @@ TEST(GridTest, getMaximumOnNode)
     EXPECT_TRUE(actual == expected);
 }
 
+TEST(GridTest, getBoundingBoxOnNodes)
+{
+    auto gridStrategyDummy = SPtr<GridStrategy>(new GridStrategyDummy);
+
+    const real startX = -1.0;
+    const real startY = -1.0;
+    const real startZ = -1.0;
+    const real delta = 0.25;
+    const auto sut = GridImp::makeShared(NULL, startX, startY, startZ, 10, 10, 10, delta, gridStrategyDummy, Distribution());
 
+    Triangle t = Triangle(Vertex(0,0,0), Vertex(1,0,0), Vertex(1,1,0), Vertex(0.0f, 0.0f, 0.0f));
+
+    const auto actual = sut->getBoundingBoxOnNodes(t);
+
+    EXPECT_THAT(actual.maxX, RealEq(1.25));
+}
 
 
 
diff --git a/targets/apps/HULC/main.cpp b/targets/apps/HULC/main.cpp
index 2047c39e2..5c194bd85 100644
--- a/targets/apps/HULC/main.cpp
+++ b/targets/apps/HULC/main.cpp
@@ -276,7 +276,7 @@ void multipleLevel(const std::string& configPath)
     //gridBuilder->addCoarseGrid(-16, -14, -14, 59, 28, 29, 1.0);
     //TriangularMesh* triangularMesh = TriangularMesh::make("D:/GRIDGENERATION/STL/input/local_input/bruecke.stl");
 
-    gridBuilder->addCoarseGrid(-10, -10, -10, 10, 10, 10, 0.25);
+    gridBuilder->addCoarseGrid(-10, -10, -10, 9, 9, 9, 0.25);
     TriangularMesh* triangularMesh = TriangularMesh::make("D:/GRIDGENERATION/STL/cubeBinaer1x1.stl");
 
     gridBuilder->addGeometry(triangularMesh);
@@ -319,7 +319,7 @@ void multipleLevel(const std::string& configPath)
     //SimulationFileWriter::write("D:/GRIDGENERATION/files/", gridBuilder, FILEFORMAT::ASCII);
 
     gridBuilder->writeGridsToVtk("D:/GRIDGENERATION/");
-
+    gridBuilder->writeArrows("D:/arrows");
 
 
     SPtr<Parameter> para = Parameter::make();
-- 
GitLab