blockMesh generator

[Cplusc]

Following my question regarding hexahedral mesh generation in a domain using the STL file representing the surface triangulation here https://cplusplus.com/forum/general/285328/ I ended up writing my own mesh generator using block mesh. I made a simple geometry like cylinder in comsol and export it in binary format. The problem is the generated blocks are all of size zero. Here is the code, I appreciate your help.

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <queue>
using namespace std;

typedef std::vector<double>          VecDbl_t;
typedef std::vector<int>             VecInt_t;
typedef std::vector<VecDbl_t>        VecVecDbl_t;
typedef std::vector<VecInt_t>        VecVecInt_t;

struct Vertex {
    double x, y, z;
};


struct Triangle {
    double x1, y1, z1;
    double x2, y2, z2;
    double x3, y3, z3;
};

struct BlockMesh {
    std::vector<Vertex> vertices;
    VecVecInt_t blocks;
};

std::vector<Triangle> readSTL(const std::string& stlFile) {
    std::vector<Triangle> triangles;

    std::ifstream file(stlFile, std::ios::in | std::ios::binary);

    if (!file) {
        std::cerr << "Error opening file " << stlFile << std::endl;
        exit(EXIT_FAILURE);
    }

    file.seekg(80, std::ios::beg);
    uint32_t numTriangles;
    file.read((char*)&numTriangles, sizeof(numTriangles));
    for (uint32_t i = 0; i < numTriangles; i++) {

        file.seekg(3 * sizeof(float), std::ios::cur);
        float x1, y1, z1, x2, y2, z2, x3, y3, z3;
        file.read((char*)&x1, sizeof(float));
        file.read((char*)&y1, sizeof(float));
        file.read((char*)&z1, sizeof(float));
        file.read((char*)&x2, sizeof(float));
        file.read((char*)&y2, sizeof(float));
        file.read((char*)&z2, sizeof(float));
        file.read((char*)&x3, sizeof(float));
        file.read((char*)&y3, sizeof(float));
        file.read((char*)&z3, sizeof(float));
        Triangle tri = { x1, y1, z1, x2, y2, z2, x3, y3, z3 };
        triangles.push_back(tri);
        file.seekg(sizeof(uint16_t), std::ios::cur);
    }
    file.close();

    std::cout << "Number of triangles: " << triangles.size() << std::endl;
    for (uint32_t i = 0; i < triangles.size(); i++) {
        std::cout << "Triangle " << i + 1 << ":" << std::endl;
        std::cout << "Vertex 1: " << triangles[i].x1 << ", " << triangles[i].y1 << ", " << triangles[i].z1 << std::endl;
        std::cout << "Vertex 2: " << triangles[i].x2 << ", " << triangles[i].y2 << ", " << triangles[i].z2 << std::endl;
        std::cout << "Vertex 3: " << triangles[i].x3 << ", " << triangles[i].y3 << ", " << triangles[i].z3 << std::endl;
    }

    return triangles;
}


double distanceToTriangle(double x, double y, double z, const Triangle& triangle) {
    double ux = triangle.x2 - triangle.x1;
    double uy = triangle.y2 - triangle.y1;
    double uz = triangle.z2 - triangle.z1;
    double vx = triangle.x3 - triangle.x1;
    double vy = triangle.y3 - triangle.y1;
    double vz = triangle.z3 - triangle.z1;
    double wx = x - triangle.x1;
    double wy = y - triangle.y1;
    double wz = z - triangle.z1;

    double uu = ux * ux + uy * uy + uz * uz;
    double uv = ux * vx + uy * vy + uz * vz;
    double vv = vx * vx + vy * vy + vz * vz;
    double uw = ux * wx + uy * wy + uz * wz;
    double vw = vx * wx + vy * wy + vz * wz;

    double denom = uu * vv - uv * uv;
    double s = (uv * vw - vv * uw) / denom;
    double t = (uv * uw - uu * vw) / denom;

    double dist = std::sqrt(uu * vv - uv * uv) / std::sqrt(denom);

    if (s >= 0 && t >= 0 && s + t <= 1) {
        return dist;
    }
    else {
        double dist1 = std::sqrt((x - triangle.x1) * (x - triangle.x1) + (y - triangle.y1) * (y - triangle.y1) + (z - triangle.z1) * (z - triangle.z1));
        double dist2 = std::sqrt((x - triangle.x2) * (x - triangle.x2) + (y - triangle.y2) * (y - triangle.y2) + (z - triangle.z2) * (z - triangle.z2));
        double dist3 = std::sqrt((x - triangle.x3) * (x - triangle.x3) + (y - triangle.y3) * (y - triangle.y3) + (z - triangle.z3) * (z - triangle.z3));
        return std::min({ dist1, dist2, dist3 });
    }
}

BlockMesh createBlockMesh(const std::string& stlFile, double dx, double maxDist) {

    std::vector<Triangle> triangles = readSTL(stlFile);

    double minX = std::numeric_limits<double>::max();
    double minY = std::numeric_limits<double>::max();
    double minZ = std::numeric_limits<double>::max();
    double maxX = std::numeric_limits<double>::lowest();
    double maxY = std::numeric_limits<double>::lowest();
    double maxZ = std::numeric_limits<double>::lowest();
    for (const auto& triangle : triangles) {
        minX = std::min(minX, std::min(triangle.x1, std::min(triangle.x2, triangle.x3)));
        minY = std::min(minY, std::min(triangle.y1, std::min(triangle.y2, triangle.y3)));
        minZ = std::min(minZ, std::min(triangle.z1, std::min(triangle.z2, triangle.z3)));
        maxX = std::max(maxX, std::max(triangle.x1, std::max(triangle.x2, triangle.x3)));
        maxY = std::max(maxY, std::max(triangle.y1, std::max(triangle.y2, triangle.y3)));
        maxZ = std::max(maxZ, std::max(triangle.z1, std::max(triangle.z2, triangle.z3)));
    }

    double blockX = std::ceil((maxX - minX) / dx);
    double blockY = std::ceil((maxY - minY) / dx);
    double blockZ = std::ceil((maxZ - minZ) / dx);

    std::vector<Vertex> vertices;
    VecVecInt_t blocks(blockX * blockY * blockZ);
    for (int i = 0; i < blockX; i++) {
        for (int j = 0; j < blockY; j++) {
            for (int k = 0; k < blockZ; k++) {
                double x = minX + i * dx;
                double y = minY + j * dx;
                double z = minZ + k * dx;
                bool inside = false;
                for (const auto& triangle : triangles) {
                    double dist = distanceToTriangle(x, y, z, triangle);
                    if (dist < maxDist) {
                        inside = !inside;
                    }
                }
                if (inside) {
                    int idx = k + j * blockZ + i * blockZ * blockY;
                    blocks[idx].push_back(vertices.size());
                    vertices.push_back({ x, y, z });
                }
            }
        }
    }

    BlockMesh blockMesh = { vertices, blocks };
    return blockMesh;
}


int main() {

    BlockMesh mesh = createBlockMesh("cylinder.stl", 0.1,0.1);
    VecVecInt_t block = mesh.blocks;
    return 0;
} 

[salem c]

> The problem is the generated blocks are all of size zero.
So have you worked out how to use the debugger yet?

Because now would seem to be a good time.

Just a few basic commands will tell you a lot
- breakpoints to stop at a point in the code (I'd suggest line 102 to begin)
- print to examine the values of variables.
- step to execute the next line

[seeplus]

Have you checked that readSTL() is working as expected? Incorrect read routines are a common source of issues.

[Cplusc]

@seeplus yes, the reader is working properly. But I'll do a double check.

[Cplusc]

I just have one doubt on this code. I used Möller–Trumbore ray-triangle intersection algorithm, but I'm not if it is correctly check whether the node falls within the bounded region or not.

[Cplusc]

@salem c There was one bug with stl reader which I fixed it in the post but still I'm getting the size of inernal vector for blocks as zeros. In the stl that I'm getting a 2d vector for blocks of size 8000 but all internal vector are empty.

[seeplus]

Are you sure L148 is what you want? How far have you got with the debugging as suggested by salem_c above?

[salem c]

> In the stl that I'm getting a 2d vector for blocks of size 8000 but all internal vector are empty.
And you're only finding this out now?

Your compile/edit/test cycle needs to be a lot tighter.

You start small - like writing 5 lines of code. You compile and test (eg, step through with the debugger). If it doesn't work, you fix it before you bury the problem under more code. It's easier to fix, because you only just wrote the code. It isn't days old by the time you try to run it, by which time you've forgotten half of what you intended to begin with.

Sure, as your experience grows, you can decide how many lines to write at once, but you always need to be prepared for that big fat "NOPE!" when you try to run it.

If it's some GUI boilerplate code, maybe you can get by with say 50 lines at once. If it's some tricky algorithm, maybe it's down to one line.


Writing 150+ lines without ever trying to run it, finding it doesn't work and dumping the whole lot on a forum isn't a long term strategy.

[Cplusc]

@salem c thanks, but how you can write a code just in 5 lines or so to test the code for what it's been written for. Actually this is not a big code at all, there is 2 structs for points and triangles. One function for reading stl file and other two functions check the distance of a point from triangle in space. All of them have been tested separately. The main part of the code is create block function. The problem should be here.

[Cplusc]

For future readers, the bug in the code is the createBlockMesh() just generate those nodes that fall within a distance less that maxDist, but the internal node will not be generated.

[salem c]

> Actually this is not a big code at all
Big enough for you to dump on a forum with no means for anyone else to replicate the problem except by eyeballing it.

Without cylinder.stl, it's just a bunch of characters on screen.

> but how you can write a code just in 5 lines or so to test the code for what it's been written for.
Your first cut would have been.

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <queue>
using namespace std;

typedef std::vector<double>          VecDbl_t;
typedef std::vector<int>             VecInt_t;
typedef std::vector<VecDbl_t>        VecVecDbl_t;
typedef std::vector<VecInt_t>        VecVecInt_t;

struct Triangle {
    double x1, y1, z1;
    double x2, y2, z2;
    double x3, y3, z3;
};

std::vector<Triangle> readSTL(const std::string& stlFile) {
    std::vector<Triangle> triangles;

    std::ifstream file(stlFile, std::ios::in | std::ios::binary);

    if (!file) {
        std::cerr << "Error opening file " << stlFile << std::endl;
        exit(EXIT_FAILURE);
    }

    file.seekg(80, std::ios::beg);
    uint32_t numTriangles;
    file.read((char*)&numTriangles, sizeof(numTriangles));
    for (uint32_t i = 0; i < numTriangles; i++) {

        file.seekg(3 * sizeof(float), std::ios::cur);
        float x1, y1, z1, x2, y2, z2, x3, y3, z3;
        file.read((char*)&x1, sizeof(float));
        file.read((char*)&y1, sizeof(float));
        file.read((char*)&z1, sizeof(float));
        file.read((char*)&x2, sizeof(float));
        file.read((char*)&y2, sizeof(float));
        file.read((char*)&z2, sizeof(float));
        file.read((char*)&x3, sizeof(float));
        file.read((char*)&y3, sizeof(float));
        file.read((char*)&z3, sizeof(float));
        Triangle tri = { x1, y1, z1, x2, y2, z2, x3, y3, z3 };
        triangles.push_back(tri);
        file.seekg(sizeof(uint16_t), std::ios::cur);
    }
    file.close();

    std::cout << "Number of triangles: " << triangles.size() << std::endl;
    for (uint32_t i = 0; i < triangles.size(); i++) {
        std::cout << "Triangle " << i + 1 << ":" << std::endl;
        std::cout << "Vertex 1: " << triangles[i].x1 << ", " << triangles[i].y1 << ", " << triangles[i].z1 << std::endl;
        std::cout << "Vertex 2: " << triangles[i].x2 << ", " << triangles[i].y2 << ", " << triangles[i].z2 << std::endl;
        std::cout << "Vertex 3: " << triangles[i].x3 << ", " << triangles[i].y3 << ", " << triangles[i].z3 << std::endl;
    }

    return triangles;
}

int main() {
    std::vector<Triangle> triangles = readSTL("cylinder.stl");
    return 0;
} 

The 10 lines of interest are 35 to 45.
The rest is just boilerplate to make everything work and just print some stuff out.

Ideally, cylinder.stl is a really small model. There's no point flooding your debug with 1000's of triangles when it only takes one triangle to prove whether you can read a triangle.

You make sure it works, then you move on to say adding
double distanceToTriangle(double x, double y, double z, const Triangle& triangle)

Which you test by doing this.

void test_distanceToTriangle(const std::vector<Triangle> &triangles ) {
    struct TestDistanceToTriangle {
        double x, y, z;
        Triangle triangle;
        // double expectedAnswer;
    }
    TestDistanceToTriangle tests[] = {
        { 0, 0, 0, { {0,0,0}, {0,1,0}, {1,1,1} } },
        // add more tests here as you think of them.
        // Especially when you find a bug, you add a test to find that bug
        // so you know when you fixed it, and it will stay fixed.
    };
    for ( auto t : tests ) {
        double d = distanceToTriangle(t.x, t.y, t.z, &t.triangle);
        cout << "Distance = " << d << endl;
        // if ( d != t.expectedAnswer )
    }
}

int main() {
    std::vector<Triangle> triangles = readSTL("cylinder.stl");
    test_distanceToTriangle();
    return 0;
} 

For bonus points, you can start to make it self-checking every time you run the code.


After you're happy your distance calculations work, you can just do

int main() {
    std::vector<Triangle> triangles = readSTL("cylinder.stl");
    //test_distanceToTriangle();
    return 0;
} 

knowing that it's just a simple edit to go back to testing distance calculations if you later discover a problem.

https://en.wikipedia.org/wiki/Test-driven_development

The rest is rinse and repeat.