nschloe / pygmsh

:spider_web: Gmsh for Python
GNU General Public License v3.0
858 stars 162 forks source link

Generate Hexahedral mesh inside complex geometry using gmsh #575

Open Resacd opened 1 year ago

Resacd commented 1 year ago

I wrote the following code to generate block mesh inside complex geometry but I' getting error in gmsh::model::mesh::generate(3). here is the error

Unhandled exception at 0x00007FFE19C7CD29 in gmsh.exe: Microsoft C++ exception: std::basic_string<char,std::char_traits,std::allocator > at memory location 0x0000002CEF0FEC08.

and here is part of output.

Info : Starting subloop 1525 in curve loop 1 (are you sure about this?) Info : Starting subloop 1526 in curve loop 1 (are you sure about this?) Info : Starting subloop 1527 in curve loop 1 (are you sure about this?) Info : Starting subloop 1528 in curve loop 1 (are you sure about this?) Info : Starting subloop 1529 in curve loop 1 (are you sure about this?) Info : Starting subloop 1530 in curve loop 1 (are you sure about this?) Info : Starting subloop 1531 in curve loop 1 (are you sure about this?) Info : Starting subloop 1532 in curve loop 1 (are you sure about this?) Info : Starting subloop 1533 in curve loop 1 (are you sure about this?) Info : Starting subloop 1534 in curve loop 1 (are you sure about this?) Info : Starting subloop 1535 in curve loop 1 (are you sure about this?) Error : Mesh module not compiled

maybe it helps to mention that I have 1536 triangles for the surface. it seems that it already iterated all over 1536 triangles and trying to exceed the bound.

I reinstalled the library and still the same issue. I installed the library using vcpkg. I strongly appreciate your help.

#include <gmsh.h>
#include <gmshc.h>
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>

struct Point {
    float x;
    float y;
    float z;
};

struct Triangle {
    float x1, y1, z1;
    float x2, y2, z2;
    float x3, y3, z3;
    float nx, ny, nz;
};

std::vector<Triangle> readSTLASCII(const std::string& filename) {
    std::vector<Triangle> triangles;
    std::ifstream file(filename);
    if (!file.is_open()) {
        std::cerr << "Error: cannot open file " << filename << std::endl;
        return triangles;
    }

    std::string line;
    std::stringstream ss;
    std::string facet, normal, outer, loop, vertex, endloop, endfacet;
    float x, y, z;
    Triangle triangle;
    int triangle_count = 0;
    while (getline(file, line)) {
        ss.clear();
        ss.str(line);
        ss >> facet >> normal >> triangle.nx >> triangle.ny >> triangle.nz;
        if (ss.fail() || facet != "facet" || normal != "normal") {
            continue;
        }

        getline(file, line);
        ss.clear();
        ss.str(line);
        ss >> outer >> loop;
        if (ss.fail() || outer != "outer" || loop != "loop") {
            continue;
        }

        //std::cout << "Triangle " << triangle_count << ": ";
        for (int i = 0; i < 3; ++i) {
            getline(file, line);
            ss.clear();
            ss.str(line);
            ss >> vertex >> x >> y >> z;
            if (ss.fail() || vertex != "vertex") {
                break;
            }
            switch (i) {
            case 0:
                triangle.x1 = x;
                triangle.y1 = y;
                triangle.z1 = z;
                // std::cout << x << "," << y << "," << z << ";";
                break;
            case 1:
                triangle.x2 = x;
                triangle.y2 = y;
                triangle.z2 = z;
                //std::cout << x << "," << y << "," << z << ";";
                break;
            case 2:
                triangle.x3 = x;
                triangle.y3 = y;
                triangle.z3 = z;
                //std::cout << x << "," << y << "," << z << ";";
                break;
            }
        }

        // std::cout << triangle.nx << "," << triangle.ny << "," << triangle.nz << std::endl;
        getline(file, line);
        ss.clear();
        ss.str(line);
        ss >> endloop;
        if (ss.fail() || endloop != "endloop") {
            continue;
        }

        getline(file, line);
        ss.clear();
        ss.str(line);
        ss >> endfacet;
        if (ss.fail() || endfacet != "endfacet") {
            continue;
        }

        triangles.push_back(triangle);
        triangle_count++;
    }

    file.close();
    return triangles;
}

void generateMesh(const std::vector<Triangle>& triangles) {
    // Initialize Gmsh
    gmsh::initialize();

    // Create a new model
    gmsh::model::add("hex_mesh");

    // Define the geometry
    std::vector<std::size_t> pointTags;
    std::vector<int> curveTags;
    for (const auto& t : triangles) {
        std::vector<std::size_t> currentPoints;
        currentPoints.push_back(gmsh::model::geo::addPoint(t.x1, t.y1, t.z1, 1e-3));
        currentPoints.push_back(gmsh::model::geo::addPoint(t.x2, t.y2, t.z2, 1e-3));
        currentPoints.push_back(gmsh::model::geo::addPoint(t.x3, t.y3, t.z3, 1e-3));
        std::vector<int> currentCurves;
        for (std::size_t i = 0; i < currentPoints.size(); ++i) {
            std::size_t j = (i + 1) % currentPoints.size();
            currentCurves.push_back(gmsh::model::geo::addLine(currentPoints[i], currentPoints[j]));
        }
        curveTags.insert(curveTags.end(), currentCurves.begin(), currentCurves.end());
        pointTags.insert(pointTags.end(), currentPoints.begin(), currentPoints.end());
    }
    gmsh::model::geo::addCurveLoop(curveTags);
    gmsh::model::geo::addPlaneSurface({ 1 });

    // Define the mesh
    gmsh::model::geo::synchronize();

    // Set the characteristic length of the mesh
    gmsh::vectorpair entities;
    for (auto tag : pointTags) {
        entities.emplace_back(0, tag);
    }
    gmsh::model::mesh::setSize(entities, 0.1);

    // Generate the mesh
    gmsh::model::mesh::generate(3);

    // Save the mesh to a file
    gmsh::write("hex_mesh.msh");

    // Finalize Gmsh
    gmsh::finalize();
}

int main() {

    std::vector<Triangle>Triangles = readSTLASCII("cube.stl");
    generateMesh(Triangles);
    return 0;
}