// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // issues on https://gitlab.onelab.info/gmsh/gmsh/issues. #include <stdio.h> #include <string> #include <algorithm> #include <sstream> #include "GModel.h" #include "OS.h" #include "MLine.h" #include "MTriangle.h" #include "MQuadrangle.h" #include "MVertexRTree.h" #include "discreteFace.h" #include "StringUtils.h" #include "Context.h" static bool invalidChar(char c) { return !(c >= 32 && c <= 126); } int GModel::readSTL(const std::string &name, double tolerance) { FILE *fp = Fopen(name.c_str(), "rb"); if(!fp) { Msg::Error("Unable to open file '%s'", name.c_str()); return 0; } // store triplets of points for each solid found in the file std::vector<std::vector<SPoint3> > points; SBoundingBox3d bbox; std::vector<std::string> names; // "solid", or binary data header char buffer[256]; if(!fgets(buffer, sizeof(buffer), fp)) { fclose(fp); return 0; } //SPoint3 p0(1.9e6, 4e6, 0); bool binary = strncmp(buffer, "solid", 5) && strncmp(buffer, "SOLID", 5); // ASCII STL if(!binary) { if(strlen(buffer) > 6) names.push_back(&buffer[6]); else names.push_back(""); points.resize(1); while(!feof(fp)) { // "facet normal x y z" or "endsolid" if(!fgets(buffer, sizeof(buffer), fp)) break; if(!strncmp(buffer, "endsolid", 8) || !strncmp(buffer, "ENDSOLID", 8)) { // "solid" if(!fgets(buffer, sizeof(buffer), fp)) break; if(!strncmp(buffer, "solid", 5) || !strncmp(buffer, "SOLID", 5)) { if(strlen(buffer) > 6) names.push_back(&buffer[6]); else names.push_back(""); points.resize(points.size() + 1); // "facet normal x y z" if(!fgets(buffer, sizeof(buffer), fp)) break; } } // "outer loop" if(!fgets(buffer, sizeof(buffer), fp)) break; // "vertex x y z" for(int i = 0; i < 3; i++) { if(!fgets(buffer, sizeof(buffer), fp)) break; char s1[256]; double x, y, z; if(sscanf(buffer, "%s %lf %lf %lf", s1, &x, &y, &z) != 4) break; SPoint3 p(x, y, z); //p -= p0; points.back().push_back(p); bbox += p; } // "endloop" if(!fgets(buffer, sizeof(buffer), fp)) break; // "endfacet" if(!fgets(buffer, sizeof(buffer), fp)) break; } } // check if we could parse something bool empty = true; for(std::size_t i = 0; i < points.size(); i++) { if(points[i].size()) { empty = false; break; } } if(empty) points.clear(); // binary STL (we also try to read in binary mode if the header told // us the format was ASCII but we could not read any vertices) if(binary || empty) { if(binary) Msg::Info("Mesh is in binary format"); else Msg::Info("Wrong ASCII header or empty file: trying binary read"); rewind(fp); while(!feof(fp)) { char header[80]; if(!fread(header, sizeof(char), 80, fp)) break; unsigned int nfacets = 0; size_t ret = fread(&nfacets, sizeof(unsigned int), 1, fp); bool swap = false; if(nfacets > 100000000) { Msg::Info("Swapping bytes from binary file"); swap = true; SwapBytes((char *)&nfacets, sizeof(unsigned int), 1); } if(ret && nfacets) { names.push_back(header); points.resize(points.size() + 1); char *data = new char[nfacets * 50 * sizeof(char)]; ret = fread(data, sizeof(char), nfacets * 50, fp); if(ret == nfacets * 50) { for(std::size_t i = 0; i < nfacets; i++) { float *xyz = (float *)&data[i * 50 * sizeof(char)]; if(swap) SwapBytes((char *)xyz, sizeof(float), 12); for(int j = 0; j < 3; j++) { SPoint3 p(xyz[3 + 3 * j], xyz[3 + 3 * j + 1], xyz[3 + 3 * j + 2]); //p -= p0; points.back().push_back(p); bbox += p; } } } delete[] data; } } } // cleanup names if(names.size() != points.size()){ Msg::Debug("Invalid number of names in STL file - should never happen"); names.resize(points.size()); } for(std::size_t i = 0; i < names.size(); i++){ names[i].erase(remove_if(names[i].begin(), names[i].end(), invalidChar), names[i].end()); } std::vector<GFace *> faces; for(std::size_t i = 0; i < points.size(); i++) { if(points[i].empty()) { Msg::Error("No facets found in STL file for solid %d %s", i, names[i].c_str()); fclose(fp); return 0; } if(points[i].size() % 3) { Msg::Error("Wrong number of points (%d) in STL file for solid %d %s", points[i].size(), i, names[i].c_str()); fclose(fp); return 0; } Msg::Info("%d facets in solid %d %s", points[i].size() / 3, i, names[i].c_str()); // create face GFace *face = new discreteFace(this, getMaxElementaryNumber(2) + 1); faces.push_back(face); add(face); if(!names[i].empty()) setElementaryName(2, face->tag(), names[i]); } // create triangles using unique vertices double eps = norm(SVector3(bbox.max(), bbox.min())) * tolerance; std::vector<MVertex *> vertices; for(std::size_t i = 0; i < points.size(); i++) for(std::size_t j = 0; j < points[i].size(); j++) vertices.push_back( new MVertex(points[i][j].x(), points[i][j].y(), points[i][j].z())); MVertexRTree pos(eps); pos.insert(vertices); std::set<MFace, Less_Face> unique; int nbDuplic = 0, nbDegen = 0; for(std::size_t i = 0; i < points.size(); i++) { for(std::size_t j = 0; j < points[i].size(); j += 3) { MVertex *v[3]; for(int k = 0; k < 3; k++) { double x = points[i][j + k].x(); double y = points[i][j + k].y(); double z = points[i][j + k].z(); v[k] = pos.find(x, y, z); } if(CTX::instance()->mesh.stlRemoveDuplicateTriangles){ MFace mf(v[0], v[1], v[2]); if(unique.find(mf) == unique.end()) { faces[i]->triangles.push_back(new MTriangle(v[0], v[1], v[2])); unique.insert(mf); } else { nbDuplic++; } } else{ if(v[0] == v[1] || v[0] == v[2] || v[1] == v[2]){ Msg::Debug("Skipping degenerated triangle %lu %lu %lu", v[0]->getNum(), v[1]->getNum(), v[2]->getNum()); nbDegen++; } else{ faces[i]->triangles.push_back(new MTriangle(v[0], v[1], v[2])); } } } } if(nbDuplic || nbDegen) Msg::Warning("%d duplicate/%d degenerate triangles in STL file", nbDuplic, nbDegen); _associateEntityWithMeshVertices(); _storeVerticesInEntities(vertices); // will delete unused vertices fclose(fp); return 1; } static void writeSTLfaces(FILE *fp, std::vector<GFace*> &faces, bool binary, double scalingFactor, const std::string &name) { bool useGeoSTL = false; unsigned int nfacets = 0; for(std::vector<GFace*>::iterator it = faces.begin(); it != faces.end(); ++it) { nfacets += (*it)->triangles.size() + 2 * (*it)->quadrangles.size(); } if(!nfacets) { // use CAD STL if there is no mesh useGeoSTL = true; for(std::vector<GFace*>::iterator it = faces.begin(); it != faces.end(); ++it) { (*it)->buildSTLTriangulation(); nfacets += (*it)->stl_triangles.size() / 3; } } if(!binary) { fprintf(fp, "solid %s\n", name.c_str()); } else { char header[80]; strncpy(header, name.c_str(), 80); fwrite(header, sizeof(char), 80, fp); fwrite(&nfacets, sizeof(unsigned int), 1, fp); } for(std::vector<GFace*>::iterator it = faces.begin(); it != faces.end(); ++it) { if(useGeoSTL && (*it)->stl_vertices_uv.size()) { for(std::size_t i = 0; i < (*it)->stl_triangles.size(); i += 3) { SPoint2 &p1((*it)->stl_vertices_uv[(*it)->stl_triangles[i]]); SPoint2 &p2((*it)->stl_vertices_uv[(*it)->stl_triangles[i + 1]]); SPoint2 &p3((*it)->stl_vertices_uv[(*it)->stl_triangles[i + 2]]); GPoint gp1 = (*it)->point(p1); GPoint gp2 = (*it)->point(p2); GPoint gp3 = (*it)->point(p3); double x[3] = {gp1.x(), gp2.x(), gp3.x()}; double y[3] = {gp1.y(), gp2.y(), gp3.y()}; double z[3] = {gp1.z(), gp2.z(), gp3.z()}; double n[3]; normal3points(x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], n); if(!binary) { fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]); fprintf(fp, " outer loop\n"); for(int j = 0; j < 3; j++) fprintf(fp, " vertex %g %g %g\n", x[j] * scalingFactor, y[j] * scalingFactor, z[j] * scalingFactor); fprintf(fp, " endloop\n"); fprintf(fp, "endfacet\n"); } else { char data[50]; float *coords = (float *)data; coords[0] = (float)n[0]; coords[1] = (float)n[1]; coords[2] = (float)n[2]; for(int j = 0; j < 3; j++) { coords[3 + 3 * j] = (float)(x[j] * scalingFactor); coords[3 + 3 * j + 1] = (float)(y[j] * scalingFactor); coords[3 + 3 * j + 2] = (float)(z[j] * scalingFactor); } data[48] = data[49] = 0; fwrite(data, sizeof(char), 50, fp); } } } else { for(std::size_t i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeSTL(fp, binary, scalingFactor); for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writeSTL(fp, binary, scalingFactor); } } if(!binary) fprintf(fp, "endsolid %s\n", name.c_str()); } int GModel::writeSTL(const std::string &name, bool binary, bool saveAll, double scalingFactor, int oneSolidPerSurface) { FILE *fp = Fopen(name.c_str(), binary ? "wb" : "w"); if(!fp) { Msg::Error("Unable to open file '%s'", name.c_str()); return 0; } if(noPhysicalGroups()) saveAll = true; if(oneSolidPerSurface == 1){ // one solid per surface for(fiter it = firstFace(); it != lastFace(); ++it) { if(saveAll || (*it)->physicals.size()) { std::vector<GFace*> faces(1, *it); std::string name = getElementaryName(2, (*it)->tag()); if(name.empty()){ std::ostringstream s; s << "Gmsh Surface " << (*it)->tag(); name = s.str(); } writeSTLfaces(fp, faces, binary, scalingFactor, name); } } } else if(oneSolidPerSurface == 2){ // one solid per physical surface std::map<int, std::vector<GEntity *> > phys; getPhysicalGroups(2, phys); for(std::map<int, std::vector<GEntity *> >::iterator it = phys.begin(); it != phys.end(); it++){ std::vector<GFace*> faces; for(std::size_t i = 0; i < it->second.size(); i++){ faces.push_back(static_cast<GFace*>(it->second[i])); } std::string name = getPhysicalName(2, it->first); if(name.empty()){ std::ostringstream s; s << "Gmsh Physical Surface " << it->first; name = s.str(); } writeSTLfaces(fp, faces, binary, scalingFactor, name); } } else{ // one solid std::vector<GFace*> faces; for(fiter it = firstFace(); it != lastFace(); ++it) { if(saveAll || (*it)->physicals.size()) { faces.push_back(*it); } } std::string name = "Created by Gmsh"; writeSTLfaces(fp, faces, binary, scalingFactor, name); } fclose(fp); return 1; }