#include <map> #include <string> #include "Message.h" #include "GmshDefines.h" #include "gmshRegion.h" #include "gmshFace.h" #include "gmshEdge.h" #include "MElement.h" #include "SBoundingBox3d.h" static int getNumVerticesForElementTypeMSH(int type) { switch (type) { case PNT : return 1; case LGN1: return 2; case LGN2: return 2 + 1; case TRI1: return 3; case TRI2: return 3 + 3; case QUA1: return 4; case QUA2: return 4 + 4 + 1; case TET1: return 4; case TET2: return 4 + 6; case HEX1: return 8; case HEX2: return 8 + 12 + 6 + 1; case PRI1: return 6; case PRI2: return 6 + 9 + 3; case PYR1: return 5; case PYR2: return 5 + 8 + 1; default: return 0; } } template<class T> void copyElements(std::vector<T*> &dst, const std::vector<MElement*> &src) { dst.resize(src.size()); for(unsigned int i = 0; i < src.size(); i++) dst[i] = (T*)src[i]; } static void storeElementsInEntities(GModel *m, int type, std::map<int, std::vector<MElement*> > &map) { std::map<int, std::vector<MElement*> >::const_iterator it = map.begin(); std::map<int, std::vector<MElement*> >::const_iterator ite = map.end(); for(; it != ite; ++it){ switch(type){ case LGN1: { GEdge *e = m->edgeByTag(it->first); if(!e){ e = new gmshEdge(m, it->first); m->add(e); } copyElements(e->lines, it->second); } break; case TRI1: case QUA1: { GFace *f = m->faceByTag(it->first); if(!f){ f = new gmshFace(m, it->first); m->add(f); } if(type == TRI1) copyElements(f->triangles, it->second); else if(type == QUA1) copyElements(f->quadrangles, it->second); } break; case TET1: case HEX1: case PRI1: case PYR1: { GRegion *r = m->regionByTag(it->first); if(!r){ r = new gmshRegion(m, it->first); m->add(r); } if(type == TET1) copyElements(r->tetrahedra, it->second); else if(type == HEX1) copyElements(r->hexahedra, it->second); else if(type == PRI1) copyElements(r->prisms, it->second); else if(type == PYR1) copyElements(r->pyramids, it->second); } break; } } } template<class T> static void associateEntityWithVertices(GEntity *ge, std::vector<T*> &elements) { for(unsigned int i = 0; i < elements.size(); i++) for(int j = 0; j < elements[i]->getNumVertices(); j++) elements[i]->getVertex(j)->setEntity(ge); } static void associateEntityWithVertices(GModel *m) { // loop on regions, then on faces, edges and vertices and store the // entity pointer in the the elements' vertices (this way we // associate the entity of lowest geometrical degree with each // vertex) for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){ associateEntityWithVertices(*it, (*it)->tetrahedra); associateEntityWithVertices(*it, (*it)->hexahedra); associateEntityWithVertices(*it, (*it)->prisms); associateEntityWithVertices(*it, (*it)->pyramids); } for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){ associateEntityWithVertices(*it, (*it)->triangles); associateEntityWithVertices(*it, (*it)->quadrangles); } for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){ associateEntityWithVertices(*it, (*it)->lines); } for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it){ (*it)->mesh_vertices[0]->setEntity(*it); } } static void storeVerticesInEntities(std::map<int, MVertex*> &vertices) { std::map<int, MVertex*>::const_iterator it = vertices.begin(); std::map<int, MVertex*>::const_iterator ite = vertices.end(); for(; it != ite; ++it){ MVertex *v = it->second; GEntity *ge = v->onWhat(); if(ge) ge->mesh_vertices.push_back(v); else delete v; // we delete all unused vertices } } static void storePhysicalTagsInEntities(GModel *m, int dim, std::map<int, std::map<int, std::string> > &map) { std::map<int, std::map<int, std::string> >::const_iterator it = map.begin(); std::map<int, std::map<int, std::string> >::const_iterator ite = map.end(); for(; it != ite; ++it){ GEntity *ge = 0; switch(dim){ case 0: ge = m->vertexByTag(it->first); break; case 1: ge = m->edgeByTag(it->first); break; case 2: ge = m->faceByTag(it->first); break; case 3: ge = m->regionByTag(it->first); break; } if(ge){ std::map<int, std::string>::const_iterator it2 = it->second.begin(); std::map<int, std::string>::const_iterator ite2 = it->second.end(); for(; it2 != ite2; ++it2) ge->physicals.push_back(it2->first); } } } int GModel::readMSH(const std::string &name) { FILE *fp = fopen(name.c_str(), "r"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } int elementTypes[7] = {LGN1, TRI1, QUA1, TET1, HEX1, PRI1, PYR1}; double version = 1.0; char str[256]; SBoundingBox3d bbox; std::map<int, MVertex*> vertices; std::map<int, std::vector<MVertex*> > points; std::map<int, std::vector<MElement*> > elements[7]; std::map<int, std::map<int, std::string> > physicals[4]; while(1) { do { if(!fgets(str, sizeof(str), fp) || feof(fp)) break; } while(str[0] != '$'); if(feof(fp)) break; if(!strncmp(&str[1], "MeshFormat", 10)) { int format, size; fscanf(fp, "%lf %d %d\n", &version, &format, &size); } else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5)) { int numVertices; fscanf(fp, "%d", &numVertices); Msg(INFO, "%d vertices", numVertices); int progress = (numVertices > 100000) ? numVertices / 50 : numVertices / 10; for(int i = 0; i < numVertices; i++) { int num; double x, y, z; fscanf(fp, "%d %lf %lf %lf", &num, &x, &y, &z); bbox += SPoint3(x, y, z); if(vertices.count(num)) Msg(WARNING, "Skipping duplicate vertex %d", num); else vertices[num] = new MVertex(x, y, z); if(progress && (i % progress == progress - 1)) Msg(PROGRESS, "Read %d vertices", i + 1); } Msg(PROGRESS, ""); } else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) { int numElements; fscanf(fp, "%d", &numElements); Msg(INFO, "%d elements", numElements); int progress = (numElements > 100000) ? numElements / 50 : numElements / 10; for(int i = 0; i < numElements; i++) { int num, type, physical = 1, elementary = 1, partition = 1, numVertices; if(version <= 1.0){ fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices); int check = getNumVerticesForElementTypeMSH(type); if(!check){ Msg(GERROR, "Unknown type for element %d", num); continue; } if(numVertices != check){ Msg(GERROR, "Wrong number of vertices (%d) for element %d", numVertices, num); continue; } } else{ int numTags; fscanf(fp, "%d %d %d", &num, &type, &numTags); for(int j = 0; j < numTags; j++){ int tag; fscanf(fp, "%d", &tag); if(j == 0) physical = tag; else if(j == 1) elementary = tag; else if(j == 2) partition = tag; // ignore any other tags for now } numVertices = getNumVerticesForElementTypeMSH(type); if(!numVertices){ Msg(GERROR, "Unknown type (%d) for element %d", type, num); continue; } } int n[30]; for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &n[j]); int dim = 0; switch (type) { case PNT: points[elementary].push_back(vertices[n[0]]); dim = 0; break; case LGN1: elements[0][elementary].push_back (new MLine(vertices[n[0]], vertices[n[1]], num, partition)); dim = 1; break; case TRI1: elements[1][elementary].push_back (new MTriangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], num, partition)); dim = 2; break; case QUA1: elements[2][elementary].push_back (new MQuadrangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], vertices[n[3]], num, partition)); dim = 2; break; case TET1: elements[3][elementary].push_back (new MTetrahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]], vertices[n[3]], num, partition)); dim = 3; break; case HEX1: elements[4][elementary].push_back (new MHexahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]], vertices[n[3]], vertices[n[4]], vertices[n[5]], vertices[n[6]], vertices[n[7]], num, partition)); dim = 3; break; case PRI1: elements[5][elementary].push_back (new MPrism(vertices[n[0]], vertices[n[1]], vertices[n[2]], vertices[n[3]], vertices[n[4]], vertices[n[5]], num, partition)); dim = 3; break; case PYR1: elements[6][elementary].push_back (new MPyramid(vertices[n[0]], vertices[n[1]], vertices[n[2]], vertices[n[3]], vertices[n[4]], num, partition)); dim = 3; break; default: Msg(GERROR, "Unknown type (%d) for element %d", type, num); break; } if(!physicals[dim].count(elementary) || !physicals[dim][elementary].count(physical)) physicals[dim][elementary][physical] = "unnamed"; if(progress && (i % progress == progress - 1)) Msg(PROGRESS, "Read %d elements", i + 1); } Msg(PROGRESS, ""); } do { if(!fgets(str, sizeof(str), fp) || feof(fp)) Msg(GERROR, "Prematured end of mesh file"); } while(str[0] != '$'); } // store the elements in their associated elementary entity. If the // entity does not exist, create a new one. for(int i = 0; i < 7; i++) storeElementsInEntities(this, elementTypes[i], elements[i]); // treat points separately { std::map<int, std::vector<MVertex*> >::const_iterator it = points.begin(); std::map<int, std::vector<MVertex*> >::const_iterator ite = points.end(); for(; it != ite; ++it){ GVertex *v = vertexByTag(it->first); if(!v){ v = new gmshVertex(this, it->first); add(v); } v->mesh_vertices.push_back(it->second[0]); } } // associate the correct geometrical entity with each mesh vertex associateEntityWithVertices(this); // special case for geometry vertices: now that the correct // geometrical entity has been associated with the vertices, we // reset mesh_vertices so that it can be filled again below for(viter it = firstVertex(); it != lastVertex(); ++it) (*it)->mesh_vertices.clear(); // store the vertices in their associated geometrical entity storeVerticesInEntities(vertices); // store the physical tags for(int i = 0; i < 4; i++) storePhysicalTagsInEntities(this, i, physicals[i]); boundingBox += bbox; fclose(fp); return 1; } template<class T> static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, int saveAll, double version, int &num, int elementary, std::vector<int> &physicals) { for(unsigned int i = 0; i < ele.size(); i++) if(saveAll) ele[i]->writeMSH(fp, version, ++num, elementary, elementary); else for(unsigned int j = 0; j < physicals.size(); j++) ele[i]->writeMSH(fp, version, ++num, elementary, physicals[j]); } int GModel::writeMSH(const std::string &name, double version, bool saveAll, double scalingFactor) { FILE *fp = fopen(name.c_str(), "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } // if there are no physicals we save all the elements if(noPhysicals()) saveAll = true; // get the number of vertices and renumber the vertices in a // continuous sequence int numVertices = renumberMeshVertices(); // get the number of elements int numElements = 0; for(viter it = firstVertex(); it != lastVertex(); ++it) numElements += (saveAll ? 1 : (*it)->physicals.size()) * (*it)->mesh_vertices.size(); for(eiter it = firstEdge(); it != lastEdge(); ++it) numElements += (saveAll ? 1 : (*it)->physicals.size()) * (*it)->lines.size(); for(fiter it = firstFace(); it != lastFace(); ++it) numElements += (saveAll ? 1 : (*it)->physicals.size()) * ((*it)->triangles.size() + (*it)->quadrangles.size()); for(riter it = firstRegion(); it != lastRegion(); ++it) numElements += (saveAll ? 1 : (*it)->physicals.size()) * ((*it)->tetrahedra.size() + (*it)->hexahedra.size() + (*it)->prisms.size() + (*it)->pyramids.size()); if(version > 2.0){ fprintf(fp, "$MeshFormat\n"); fprintf(fp, "%g %d %d\n", version, 0, (int)sizeof(double)); fprintf(fp, "$EndMeshFormat\n"); fprintf(fp, "$Nodes\n"); } else fprintf(fp, "$NOD\n"); fprintf(fp, "%d\n", numVertices); for(viter it = firstVertex(); it != lastVertex(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor); for(eiter it = firstEdge(); it != lastEdge(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor); for(riter it = firstRegion(); it != lastRegion(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor); if(version > 2.0){ fprintf(fp, "$EndNodes\n"); fprintf(fp, "$Elements\n"); } else{ fprintf(fp, "$ENDNOD\n"); fprintf(fp, "$ELM\n"); } fprintf(fp, "%d\n", numElements); int num = 0; for(viter it = firstVertex(); it != lastVertex(); ++it){ writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, num, (*it)->tag(), (*it)->physicals); } for(eiter it = firstEdge(); it != lastEdge(); ++it){ writeElementsMSH(fp, (*it)->lines, saveAll, version, num, (*it)->tag(), (*it)->physicals); } for(fiter it = firstFace(); it != lastFace(); ++it){ writeElementsMSH(fp, (*it)->triangles, saveAll, version, num, (*it)->tag(), (*it)->physicals); writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, num, (*it)->tag(), (*it)->physicals); } for(riter it = firstRegion(); it != lastRegion(); ++it){ writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, num, (*it)->tag(), (*it)->physicals); writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, num, (*it)->tag(), (*it)->physicals); writeElementsMSH(fp, (*it)->prisms, saveAll, version, num, (*it)->tag(), (*it)->physicals); writeElementsMSH(fp, (*it)->pyramids, saveAll, version, num, (*it)->tag(), (*it)->physicals); } if(version > 2.0){ fprintf(fp, "$EndElements\n"); } else{ fprintf(fp, "$ENDELM\n"); } return 1; } int GModel::writePOS(const std::string &name, double scalingFactor) { FILE *fp = fopen(name.c_str(), "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } if(numRegion()){ fprintf(fp, "View \"Volumes\" {\n"); fprintf(fp, "T2(1.e5,30,%d){\"Elementary Entity\", \"Element Number\", " "\"Gamma\", \"Eta\", \"Rho\"};\n", (1<<16)|(4<<8)); for(riter it = firstRegion(); it != lastRegion(); ++it) { for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++) (*it)->tetrahedra[i]->writePOS(fp, scalingFactor); for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++) (*it)->hexahedra[i]->writePOS(fp, scalingFactor); for(unsigned int i = 0; i < (*it)->prisms.size(); i++) (*it)->prisms[i]->writePOS(fp, scalingFactor); for(unsigned int i = 0; i < (*it)->pyramids.size(); i++) (*it)->pyramids[i]->writePOS(fp, scalingFactor); } fprintf(fp, "};\n"); } if(numFace()){ fprintf(fp, "View \"Surfaces\" {\n"); fprintf(fp, "T2(1.e5,30,%d){\"Elementary Entity\", \"Element Number\", " "\"Gamma\", \"Eta\", \"Rho\"};\n", (1<<16)|(4<<8)); for(fiter it = firstFace(); it != lastFace(); ++it) { for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writePOS(fp, scalingFactor); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writePOS(fp, scalingFactor); } fprintf(fp, "};\n"); } if(numEdge()){ fprintf(fp, "View \"Lines\" {\n"); fprintf(fp, "T2(1.e5,30,%d){\"Elementary Entity\", \"Element Number\", " "\"Gamma\", \"Eta\", \"Rho\"};\n", (1<<16)|(4<<8)); for(eiter it = firstEdge(); it != lastEdge(); ++it) { for(unsigned int i = 0; i < (*it)->lines.size(); i++) (*it)->lines[i]->writePOS(fp, scalingFactor); } fprintf(fp, "};\n"); } fclose(fp); return 1; } static void swapBytes(char *array, int size, int n) { char *x = new char[size]; for(int i = 0; i < n; i++) { char *a = &array[i * size]; memcpy(x, a, size); for(int c = 0; c < size; c++) a[size - 1 - c] = x[c]; } delete [] x; } int GModel::readSTL(const std::string &name, double tolerance) { FILE *fp = fopen(name.c_str(), "rb"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } // read all facets and store triplets of points std::vector<SPoint3> points; SBoundingBox3d bbox; // "solid" or binary data header char buffer[256]; fgets(buffer, sizeof(buffer), fp); if(!strncmp(buffer, "solid", 5)) { // ASCII STL while(!feof(fp)) { // "facet normal x y z" or "endsolid" if(!fgets(buffer, sizeof(buffer), fp)) break; if(!strncmp(buffer, "endsolid", 8)) break; char s1[256], s2[256]; float x, y, z; sscanf(buffer, "%s %s %f %f %f", s1, s2, &x, &y, &z); // "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; sscanf(buffer, "%s %f %f %f", s1, &x, &y, &z); SPoint3 p(x, y, z); points.push_back(p); bbox += p; } // "endloop" if(!fgets(buffer, sizeof(buffer), fp)) break; // "endfacet" if(!fgets(buffer, sizeof(buffer), fp)) break; } } else{ // Binary STL rewind(fp); char header[80]; if(fread(header, sizeof(char), 80, fp)){ unsigned int nfacets = 0; size_t ret = fread(&nfacets, sizeof(unsigned int), 1, fp); bool swap = false; if(nfacets > 10000000){ Msg(INFO, "Swapping bytes from binary file"); swap = true; swapBytes((char*)&nfacets, sizeof(unsigned int), 1); } if(ret && nfacets){ char *data = new char[nfacets * 50 * sizeof(char)]; ret = fread(data, sizeof(char), nfacets * 50, fp); if(ret == nfacets * 50){ for(unsigned int 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]); points.push_back(p); bbox += p; } } } delete data; } } } if(!points.size()){ Msg(GERROR, "No facets found in STL file"); return 0; } if(points.size() % 3){ Msg(GERROR, "Wrong number of points in STL file"); return 0; } Msg(INFO, "%d facets", points.size() / 3); // create face GFace *face = new gmshFace(this, numFace() + 1); add(face); // create (unique) vertices and triangles SPoint3 min = bbox.min(); SPoint3 max = bbox.max(); double dx = max.x() - min.x(), dy = max.y() - min.y(), dz = max.z() - min.z(); double lc = sqrt(dx * dx + dy * dy + dz * dz); MVertexLessThanLexicographic::tolerance = lc * tolerance; std::set<MVertex*, MVertexLessThanLexicographic> vertices; for(unsigned int i = 0; i < points.size(); i += 3){ MVertex *v[3]; for(int j = 0; j < 3; j++){ double x = points[i + j].x(), y = points[i + j].y(), z = points[i + j].z(); MVertex w(x, y, z); std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = vertices.find(&w); if(it != vertices.end()) { v[j] = *it; } else { v[j] = new MVertex(x, y, z); vertices.insert(v[j]); face->mesh_vertices.push_back(v[j]); } } face->triangles.push_back(new MTriangle(v[0], v[1], v[2])); } boundingBox += bbox; fclose(fp); return 1; } int GModel::writeSTL(const std::string &name, double scalingFactor) { FILE *fp = fopen(name.c_str(), "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } fprintf(fp, "solid Created by Gmsh\n"); for(fiter it = firstFace(); it != lastFace(); ++it) { for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeSTL(fp, scalingFactor); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writeSTL(fp, scalingFactor); } fprintf(fp, "endsolid Created by Gmsh\n"); fclose(fp); return 1; } int GModel::writeVRML(const std::string &name, double scalingFactor) { FILE *fp = fopen(name.c_str(), "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } renumberMeshVertices(); fprintf(fp, "#VRML V1.0 ascii\n"); fprintf(fp, "#created by Gmsh\n"); fprintf(fp, "Coordinate3 {\n"); fprintf(fp, " point [\n"); for(viter it = firstVertex(); it != lastVertex(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor); for(eiter it = firstEdge(); it != lastEdge(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor); fprintf(fp, " ]\n"); fprintf(fp, "}\n"); for(eiter it = firstEdge(); it != lastEdge(); ++it){ fprintf(fp, "DEF Curve%d IndexedLineSet {\n", (*it)->tag()); fprintf(fp, " coordIndex [\n"); for(unsigned int i = 0; i < (*it)->lines.size(); i++) (*it)->lines[i]->writeVRML(fp); fprintf(fp, " ]\n"); fprintf(fp, "}\n"); } for(fiter it = firstFace(); it != lastFace(); ++it){ fprintf(fp, "DEF Surface%d IndexedFaceSet {\n", (*it)->tag()); fprintf(fp, " coordIndex [\n"); for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeVRML(fp); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writeVRML(fp); fprintf(fp, " ]\n"); fprintf(fp, "}\n"); } fclose(fp); return 1; } static void addInGroupOfNodesUNV(FILE *fp, std::set<int> &nodes, int num) { std::set<int>::iterator it = nodes.find(num); if(it == nodes.end()){ nodes.insert(num); fprintf(fp, "%10d%10d%2d%2d%2d%2d%2d%2d\n", num, 1, 0, 1, 0, 0, 0, 0); fprintf(fp, " 0.0000000000000000D+00 1.0000000000000000D+00" " 0.0000000000000000D+00\n"); fprintf(fp, " 0.0000000000000000D+00 0.0000000000000000D+00" " 0.0000000000000000D+00\n"); fprintf(fp, "%10d%10d%10d%10d%10d%10d\n", 0, 0, 0, 0, 0, 0); } } template<class T> static void addInGroupOfNodesUNV(FILE *fp, std::set<int> &nodes, std::vector<T*> &elements) { for(unsigned int i = 0; i < elements.size(); i++) for(int j = 0; j < elements[i]->getNumVertices(); j++) addInGroupOfNodesUNV(fp, nodes, elements[i]->getVertex(j)->getNum()); } int GModel::writeUNV(const std::string &name, double scalingFactor) { FILE *fp = fopen(name.c_str(), "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } // IDEAS records const int NODES=2411, ELEMENTS=2412, GROUPOFNODES=790; // IDEAS elements const int BEAM=21, THINSHLL=91, QUAD=94, SOLIDFEM=111, WEDGE=112, BRICK=115; //const int BEAM2=24, THINSHLL2=92, QUAD2=95/*?*/, SOLIDFEM2=118; renumberMeshVertices(); fprintf(fp, "%6d\n", -1); fprintf(fp, "%6d\n", NODES); for(viter it = firstVertex(); it != lastVertex(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor); for(eiter it = firstEdge(); it != lastEdge(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor); for(riter it = firstRegion(); it != lastRegion(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor); fprintf(fp, "%6d\n", -1); fprintf(fp, "%6d\n", -1); fprintf(fp, "%6d\n", ELEMENTS); for(eiter it = firstEdge(); it != lastEdge(); ++it){ for(unsigned int i = 0; i < (*it)->lines.size(); i++) (*it)->lines[i]->writeUNV(fp, BEAM, (*it)->tag()); } for(fiter it = firstFace(); it != lastFace(); ++it){ for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeUNV(fp, THINSHLL, (*it)->tag()); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writeUNV(fp, QUAD, (*it)->tag()); } for(riter it = firstRegion(); it != lastRegion(); ++it){ for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++) (*it)->tetrahedra[i]->writeUNV(fp, SOLIDFEM, (*it)->tag()); for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++) (*it)->hexahedra[i]->writeUNV(fp, BRICK, (*it)->tag()); for(unsigned int i = 0; i < (*it)->prisms.size(); i++) (*it)->prisms[i]->writeUNV(fp, WEDGE, (*it)->tag()); } fprintf(fp, "%6d\n", -1); std::map<int, std::vector<GEntity*> > physicals[4]; getPhysicalGroups(physicals); for(int dim = 0; dim < 4; dim++){ std::map<int, std::vector<GEntity*> >::const_iterator it = physicals[dim].begin(); std::map<int, std::vector<GEntity*> >::const_iterator ite = physicals[dim].end(); for(; it != ite; ++it){ fprintf(fp, "%6d\n", -1); fprintf(fp, "%6d\n", GROUPOFNODES); fprintf(fp, "%10d%10d\n", it->first, 1); fprintf(fp, "LOAD SET %2d\n", 1); std::set<int> nodes; for(unsigned int i = 0; i < it->second.size(); i++){ // we could also do this using the mesh_vertices of the entity // and all the entities in the closure GVertex *v; switch(dim){ case 0: v = (GVertex*)it->second[i]; for(unsigned int j = 0; j < v->mesh_vertices.size(); j++) addInGroupOfNodesUNV(fp, nodes, v->mesh_vertices[j]->getNum()); break; case 1: addInGroupOfNodesUNV(fp, nodes, ((GEdge*)it->second[i])->lines); break; case 2: addInGroupOfNodesUNV(fp, nodes, ((GFace*)it->second[i])->triangles); addInGroupOfNodesUNV(fp, nodes, ((GFace*)it->second[i])->quadrangles); break; case 3: addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->tetrahedra); addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->hexahedra); addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->prisms); addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->pyramids); break; } } fprintf(fp, "%6d\n", -1); } } fclose(fp); return 1; } int GModel::readMESH(const std::string &name) { FILE *fp = fopen(name.c_str(), "r"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } char buffer[256]; fgets(buffer, sizeof(buffer), fp); char str[256]; int format; sscanf(buffer, "%s %d", str, &format); if(format != 1){ Msg(GERROR, "Non-ASCII INRIA mesh import is not yet available"); return 0; } std::map<int, MVertex*> vertices; int elementTypes[2] = {TRI1, QUA1}; std::map<int, std::vector<MElement*> > elements[2]; SBoundingBox3d bbox; while(!feof(fp)) { fgets(buffer, sizeof(buffer), fp); if(buffer[0] != '#'){ // skip comments sscanf(buffer, "%s", str); if(!strcmp(str, "Dimension")){ fgets(buffer, sizeof(buffer), fp); } else if(!strcmp(str, "Vertices")){ fgets(buffer, sizeof(buffer), fp); int nbv; sscanf(buffer, "%d", &nbv); Msg(INFO, "%d vertices", nbv); for(int i = 0; i < nbv; i++) { fgets(buffer, sizeof(buffer), fp); int cl; double x, y, z; sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &cl); bbox += SPoint3(x, y, z); vertices[i + 1] = new MVertex(x, y, z); } } else if(!strcmp(str, "Triangles")){ fgets(buffer, sizeof(buffer), fp); int nbt; sscanf(buffer, "%d", &nbt); Msg(INFO, "%d triangles", nbt); for(int i = 0; i < nbt; i++) { fgets(buffer, sizeof(buffer), fp); int n1, n2, n3, cl; sscanf(buffer, "%d %d %d %d", &n1, &n2, &n3, &cl); elements[0][cl].push_back (new MTriangle(vertices[n1], vertices[n2], vertices[n3])); } } else if(!strcmp(str, "Quadrilaterals")) { fgets(buffer, sizeof(buffer), fp); int nbq; sscanf(buffer, "%d", &nbq); Msg(INFO, "%d quadrangles", nbq); for(int i = 0; i < nbq; i++) { fgets(buffer, sizeof(buffer), fp); int n1, n2, n3, n4, cl; sscanf(buffer, "%d %d %d %d %d", &n1, &n2, &n3, &n4, &cl); elements[1][cl].push_back (new MQuadrangle(vertices[n1], vertices[n2], vertices[n3], vertices[n4])); } } } } // store the elements in their associated elementary entity. If the // entity does not exist, create a new one. for(int i = 0; i < 2; i++) storeElementsInEntities(this, elementTypes[i], elements[i]); // associate the correct geometrical entity with each mesh vertex associateEntityWithVertices(this); // store the vertices in their associated geometrical entity storeVerticesInEntities(vertices); boundingBox += bbox; fclose(fp); return 1; } int GModel::writeMESH(const std::string &name, double scalingFactor) { FILE *fp = fopen(name.c_str(), "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } fprintf(fp, " MeshVersionFormatted 1\n"); fprintf(fp, " Dimension\n"); fprintf(fp, " 3\n"); int numVertices = renumberMeshVertices(); fprintf(fp, " Vertices\n"); fprintf(fp, " %d\n", numVertices); for(viter it = firstVertex(); it != lastVertex(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor); for(eiter it = firstEdge(); it != lastEdge(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor); for(riter it = firstRegion(); it != lastRegion(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor); int numTriangles = 0, numQuadrangles = 0; for(fiter it = firstFace(); it != lastFace(); ++it){ numTriangles += (*it)->triangles.size(); numQuadrangles += (*it)->quadrangles.size(); } if(numTriangles){ fprintf(fp, " Triangles\n"); fprintf(fp, " %d\n", numTriangles); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeMESH(fp, (*it)->tag()); } if(numQuadrangles){ fprintf(fp, " Quadrilaterals\n"); fprintf(fp, " %d\n", numQuadrangles); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->quadrangles[i]->writeMESH(fp, (*it)->tag()); } fprintf(fp, " End\n"); fclose(fp); return 1; }