Click here to Skip to main content
15,940,921 members
Please Sign up or sign in to vote.
0.00/5 (No votes)
See more:
I am trying to construct the domain for lattice boltzmann simulation. I have a binary stl file and I'm trying to generate the block mesh inside the domain. here is my code. I just keep getting the block and vertices of size 0 or 1. here is the code. I appreciate your help to get this bug.

What I have tried:

#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;
};


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 });
    }
}

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;
}

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.001);
    return 0;
}
Posted
Comments
Rick York 17-Mar-23 21:46pm    
This does not appear to be the standard STL file format. Is it a proprietary format? The standard is documented at : https://en.wikipedia.org/wiki/STL_(file_format).

Personally, I would make the triangle structure out of vertexes like this :
struct Triangle
{
    Vertex v1;
    Vertex v2;
    Vertex v3;
};
and then I would implement basic math operators for a Vertex, starting with subtraction. The first few lines of the distanceToTriangle function would look like this :
double distanceToTriangle( const Vertex & pt, const Triangle & triangle )
{
   Vertex u = triangle.v2 - triangle.v1;
   Vertex v = triangle.v3 - triangle.v1;
   Vertex w = point - Triangle.v1;
   // and so on

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)



CodeProject, 20 Bay Street, 11th Floor Toronto, Ontario, Canada M5J 2N8 +1 (416) 849-8900