15,940,921 members
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> 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;
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;
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) {

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
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```