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