// $Id: GModelIO.cpp,v 1.53 2006-09-11 23:11:08 geuzaine Exp $ // // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // // This program is free software; you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation; either version 2 of the License, or // (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 // USA. // // Please report all bugs and problems to <gmsh@geuz.org>. #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 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; } template<class T> static void addElements(std::vector<T*> &dst, const std::vector<MElement*> &src) { for(unsigned int i = 0; i < src.size(); i++) dst.push_back((T*)src[i]); } static void storeElementsInEntities(GModel *m, std::map<int, std::vector<MElement*> > &map) { std::map<int, std::vector<MElement*> >::const_iterator it = map.begin(); for(; it != map.end(); ++it){ if(!it->second.size()) continue; int numEdges = it->second[0]->getNumEdges(); switch(numEdges){ case 1: { GEdge *e = m->edgeByTag(it->first); if(!e){ e = new gmshEdge(m, it->first); m->add(e); } addElements(e->lines, it->second); } break; case 3: case 4: { GFace *f = m->faceByTag(it->first); if(!f){ f = new gmshFace(m, it->first); m->add(f); } if(numEdges == 3) addElements(f->triangles, it->second); else addElements(f->quadrangles, it->second); } break; case 6: case 12: case 9: case 8: { GRegion *r = m->regionByTag(it->first); if(!r){ r = new gmshRegion(m, it->first); m->add(r); } if(numEdges == 6) addElements(r->tetrahedra, it->second); else if(numEdges == 12) addElements(r->hexahedra, it->second); else if(numEdges == 9) addElements(r->prisms, it->second); else addElements(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(); for(; it != vertices.end(); ++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 storeVerticesInEntities(std::vector<MVertex*> &vertices) { for(unsigned int i = 0; i < vertices.size(); i++){ MVertex *v = vertices[i]; if(v){ // the vector can have null entries (first or last element) 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(); for(; it != map.end(); ++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(); for(; it2 != it->second.end(); ++it2) ge->physicals.push_back(it2->first); } } } static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map, std::vector<MVertex*> &vertices) { for(int i = 0; i < num; i++){ if(!map.count(indices[i])){ Msg(GERROR, "Wrong vertex index %d", indices[i]); return false; } else vertices.push_back(map[indices[i]]); } return true; } static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec, std::vector<MVertex*> &vertices) { for(int i = 0; i < num; i++){ if(indices[i] < 0 || indices[i] > (int)(vec.size() - 1)){ Msg(GERROR, "Wrong vertex index %d", indices[i]); return false; } else vertices.push_back(vec[indices[i]]); } return true; } static int getNumVerticesForElementTypeMSH(int type) { switch (type) { case PNT_1 : return 1; case LIN_2 : return 2; case LIN_3 : return 2 + 1; case TRI_3 : return 3; case TRI_6 : return 3 + 3; case QUA_4 : return 4; case QUA_8 : return 4 + 4; case QUA_9 : return 4 + 4 + 1; case TET_4 : return 4; case TET_10 : return 4 + 6; case HEX_8 : return 8; case HEX_20 : return 8 + 12; case HEX_27 : return 8 + 12 + 6 + 1; case PRI_6 : return 6; case PRI_15 : return 6 + 9; case PRI_18 : return 6 + 9 + 3; case PYR_5 : return 5; case PYR_13 : return 5 + 8; case PYR_14 : return 5 + 8 + 1; default: Msg(GERROR, "Unknown type of element %d", type); return 0; } } static void createElementMSH(GModel *m, int num, int type, int physical, int reg, int part, std::vector<MVertex*> &v, std::map<int, std::vector<MVertex*> > &points, std::map<int, std::vector<MElement*> > elem[7], std::map<int, std::map<int, std::string> > physicals[4]) { int dim = 0; switch (type) { case PNT_1: points[reg].push_back(v[0]); dim = 0; break; case LIN_2: elem[0][reg].push_back(new MLine(v, num, part)); dim = 1; break; case LIN_3: elem[0][reg].push_back(new MLine3(v, num, part)); dim = 1; break; case TRI_3: elem[1][reg].push_back(new MTriangle(v, num, part)); dim = 2; break; case TRI_6: elem[1][reg].push_back(new MTriangle6(v, num, part)); dim = 2; break; case QUA_4: elem[2][reg].push_back(new MQuadrangle(v, num, part)); dim = 2; break; case QUA_8: elem[2][reg].push_back(new MQuadrangle8(v, num, part)); dim = 2; break; case QUA_9: elem[2][reg].push_back(new MQuadrangle9(v, num, part)); dim = 2; break; case TET_4: elem[3][reg].push_back(new MTetrahedron(v, num, part)); dim = 3; break; case TET_10: elem[3][reg].push_back(new MTetrahedron10(v, num, part)); dim = 3; break; case HEX_8: elem[4][reg].push_back(new MHexahedron(v, num, part)); dim = 3; break; case HEX_20: elem[4][reg].push_back(new MHexahedron20(v, num, part)); dim = 3; break; case HEX_27: elem[4][reg].push_back(new MHexahedron27(v, num, part)); dim = 3; break; case PRI_6: elem[5][reg].push_back(new MPrism(v, num, part)); dim = 3; break; case PRI_15: elem[5][reg].push_back(new MPrism15(v, num, part)); dim = 3; break; case PRI_18: elem[5][reg].push_back(new MPrism18(v, num, part)); dim = 3; break; case PYR_5: elem[6][reg].push_back(new MPyramid(v, num, part)); dim = 3; break; case PYR_13: elem[6][reg].push_back(new MPyramid13(v, num, part)); dim = 3; break; case PYR_14: elem[6][reg].push_back(new MPyramid14(v, num, part)); dim = 3; break; default: Msg(GERROR, "Unknown type (%d) for element %d", type, num); break; } if(physical && (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical))) physicals[dim][reg][physical] = "unnamed"; if(part) m->getMeshPartitions().insert(part); } int GModel::readMSH(const std::string &name) { FILE *fp = fopen(name.c_str(), "rb"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } double version = 1.0; bool binary = false, swap = false; char str[256]; std::map<int, MVertex*> vertexMap; std::vector<MVertex*> vertexVector; 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)) { if(!fgets(str, sizeof(str), fp)) return 0; int format, size; if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) return 0; if(format){ binary = true; Msg(INFO, "Mesh is in binary format"); int one; if(fread(&one, sizeof(int), 1, fp) != 1) return 0; if(one != 1){ swap = true; Msg(INFO, "Swapping bytes from binary file"); } } } else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5)) { if(!fgets(str, sizeof(str), fp)) return 0; int numVertices; if(sscanf(str, "%d", &numVertices) != 1) return 0; Msg(INFO, "%d vertices", numVertices); int progress = (numVertices > 100000) ? numVertices / 25 : 0; int minVertex = numVertices + 1, maxVertex = -1; for(int i = 0; i < numVertices; i++) { int num; double xyz[3]; if(!binary){ if(fscanf(fp, "%d %lf %lf %lf", &num, &xyz[0], &xyz[1], &xyz[2]) != 4) return 0; } else{ if(fread(&num, sizeof(int), 1, fp) != 1) return 0; if(swap) swapBytes((char*)&num, sizeof(int), 1); if(fread(xyz, sizeof(double), 3, fp) != 3) return 0; if(swap) swapBytes((char*)xyz, sizeof(double), 3); } minVertex = std::min(minVertex, num); maxVertex = std::max(maxVertex, num); if(vertexMap.count(num)) Msg(WARNING, "Skipping duplicate vertex %d", num); else vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2]); if(progress && (i % progress == progress - 1)) Msg(PROGRESS, "Read %d vertices", i + 1); } if(progress) Msg(PROGRESS, ""); // If the vertex numbering is dense, tranfer the map into a // vector to speed up element creation if((int)vertexMap.size() == numVertices && ((minVertex == 1 && maxVertex == numVertices) || (minVertex == 0 && maxVertex == numVertices - 1))){ Msg(INFO, "Vertex numbering is dense"); vertexVector.resize(vertexMap.size() + 1); if(minVertex == 1) vertexVector[0] = 0; else vertexVector[numVertices] = 0; std::map<int, MVertex*>::const_iterator it = vertexMap.begin(); for(; it != vertexMap.end(); ++it) vertexVector[it->first] = it->second; vertexMap.clear(); } } else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) { if(!fgets(str, sizeof(str), fp)) return 0; int numElements; sscanf(str, "%d", &numElements); Msg(INFO, "%d elements", numElements); int progress = (numElements > 100000) ? numElements / 25 : 0; if(!binary){ for(int i = 0; i < numElements; i++) { int num, type, physical = 0, elementary = 0, partition = 0, numVertices; if(version <= 1.0){ fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices); if(numVertices != getNumVerticesForElementTypeMSH(type)) return 0; } 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 } if(!(numVertices = getNumVerticesForElementTypeMSH(type))) return 0; } int indices[30]; for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &indices[j]); std::vector<MVertex*> vertices; if(vertexVector.size()){ if(!getVertices(numVertices, indices, vertexVector, vertices)) return 0; } else{ if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0; } createElementMSH(this, num, type, physical, elementary, partition, vertices, points, elements, physicals); if(progress && (i % progress == progress - 1)) Msg(PROGRESS, "Read %d elements", i + 1); } } else{ int numElementsPartial = 0; while(numElementsPartial < numElements){ int header[3]; if(fread(header, sizeof(int), 3, fp) != 3) return 0; if(swap) swapBytes((char*)header, sizeof(int), 3); int type = header[0]; int numElms = header[1]; int numTags = header[2]; int numVertices = getNumVerticesForElementTypeMSH(type); unsigned int n = 1 + numTags + numVertices; int *data = new int[n]; for(int i = 0; i < numElms; i++) { if(fread(data, sizeof(int), n, fp) != n) return 0; if(swap) swapBytes((char*)data, sizeof(int), n); int num = data[0]; int physical = (numTags > 0) ? data[4 - numTags] : 0; int elementary = (numTags > 1) ? data[4 - numTags + 1] : 0; int partition = (numTags > 2) ? data[4 - numTags + 2] : 0; int *indices = &data[numTags + 1]; std::vector<MVertex*> vertices; if(vertexVector.size()){ if(!getVertices(numVertices, indices, vertexVector, vertices)) return 0; } else{ if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0; } createElementMSH(this, num, type, physical, elementary, partition, vertices, points, elements, physicals); if(progress && ((numElementsPartial + i) % progress == progress - 1)) Msg(PROGRESS, "Read %d elements", i + 1); } delete [] data; numElementsPartial += numElms; } } if(progress) 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 < (int)(sizeof(elements)/sizeof(elements[0])); i++) storeElementsInEntities(this, elements[i]); // treat points separately { std::map<int, std::vector<MVertex*> >::const_iterator it = points.begin(); for(; it != points.end(); ++it){ GVertex *v = vertexByTag(it->first); if(!v){ v = new gmshVertex(this, it->first); add(v); } for(unsigned int i = 0; i < it->second.size(); i++) v->mesh_vertices.push_back(it->second[i]); } } // 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 if(vertexVector.size()) storeVerticesInEntities(vertexVector); else storeVerticesInEntities(vertexMap); // store the physical tags for(int i = 0; i < 4; i++) storePhysicalTagsInEntities(this, i, physicals[i]); fclose(fp); return 1; } static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elements, int t1, int t2=0, int t3=0) { if(!binary) return; int numTags = 3; int data[3]; if(elements.count(t1)){ data[0] = t1; data[1] = elements[t1]; data[2] = numTags; fwrite(data, sizeof(int), 3, fp); } else if(t2 && elements.count(t2)){ data[0] = t2; data[1] = elements[t2]; data[2] = numTags; fwrite(data, sizeof(int), 3, fp); } else if(t3 && elements.count(t3)){ data[0] = t3; data[1] = elements[t3]; data[2] = numTags; fwrite(data, sizeof(int), 3, fp); } } template<class T> static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, int saveAll, double version, bool binary, int &num, int elementary, std::vector<int> &physicals) { for(unsigned int i = 0; i < ele.size(); i++) if(saveAll) ele[i]->writeMSH(fp, version, binary, ++num, elementary, 0); else for(unsigned int j = 0; j < physicals.size(); j++) ele[i]->writeMSH(fp, version, binary, ++num, elementary, physicals[j]); } int GModel::writeMSH(const std::string &name, double version, bool binary, bool saveAll, double scalingFactor) { FILE *fp = fopen(name.c_str(), binary ? "wb" : "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } // binary format exists only in version 2 if(binary) version = 2.0; // if there are no physicals we save all the elements if(noPhysicalGroups()) saveAll = true; // get the number of vertices and renumber the vertices in a // continuous sequence int numVertices = renumberMeshVertices(); // get the number of elements (we assume that all the elements in a // list have the same type, i.e., they are all of the same // polynomial order) std::map<int,int> elements; for(viter it = firstVertex(); it != lastVertex(); ++it){ int p = (saveAll ? 1 : (*it)->physicals.size()); int n = p * (*it)->mesh_vertices.size(); if(n) elements[PNT_1] += n; } for(eiter it = firstEdge(); it != lastEdge(); ++it){ int p = (saveAll ? 1 : (*it)->physicals.size()); int n = p * (*it)->lines.size(); if(n) elements[(*it)->lines[0]->getTypeForMSH()] += n; } for(fiter it = firstFace(); it != lastFace(); ++it){ int p = (saveAll ? 1 : (*it)->physicals.size()); int n = p * (*it)->triangles.size(); if(n) elements[(*it)->triangles[0]->getTypeForMSH()] += n; n = p * (*it)->quadrangles.size(); if(n) elements[(*it)->quadrangles[0]->getTypeForMSH()] += n; } for(riter it = firstRegion(); it != lastRegion(); ++it){ int p = (saveAll ? 1 : (*it)->physicals.size()); int n = p * (*it)->tetrahedra.size(); if(n) elements[(*it)->tetrahedra[0]->getTypeForMSH()] += n; n = p * (*it)->hexahedra.size(); if(n) elements[(*it)->hexahedra[0]->getTypeForMSH()] += n; n = p * (*it)->prisms.size(); if(n) elements[(*it)->prisms[0]->getTypeForMSH()] += n; n = p * (*it)->pyramids.size(); if(n) elements[(*it)->pyramids[0]->getTypeForMSH()] += n; } int numElements = 0; std::map<int,int>::const_iterator it = elements.begin(); for(; it != elements.end(); ++it) numElements += it->second; if(version >= 2.0){ fprintf(fp, "$MeshFormat\n"); fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double)); if(binary){ int one = 1; fwrite(&one, sizeof(int), 1, fp); fprintf(fp, "\n"); } 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, binary, 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, binary, 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, binary, 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, binary, scalingFactor); if(binary) fprintf(fp, "\n"); 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; writeElementHeaderMSH(binary, fp, elements, PNT_1); for(viter it = firstVertex(); it != lastVertex(); ++it) writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, LIN_2, LIN_3); for(eiter it = firstEdge(); it != lastEdge(); ++it) writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, TRI_3, TRI_6); for(fiter it = firstFace(); it != lastFace(); ++it) writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, QUA_4, QUA_9, QUA_8); for(fiter it = firstFace(); it != lastFace(); ++it) writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, TET_4, TET_10); for(riter it = firstRegion(); it != lastRegion(); ++it) writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, HEX_8, HEX_27, HEX_20); for(riter it = firstRegion(); it != lastRegion(); ++it) writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, PRI_6, PRI_18, PRI_15); for(riter it = firstRegion(); it != lastRegion(); ++it) writeElementsMSH(fp, (*it)->prisms, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); writeElementHeaderMSH(binary, fp, elements, PYR_5, PYR_14, PYR_13); for(riter it = firstRegion(); it != lastRegion(); ++it) writeElementsMSH(fp, (*it)->pyramids, saveAll, version, binary, num, (*it)->tag(), (*it)->physicals); if(binary) fprintf(fp, "\n"); if(version >= 2.0){ fprintf(fp, "$EndElements\n"); } else{ fprintf(fp, "$ENDELM\n"); } fclose(fp); 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, (*it)->tag()); for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++) (*it)->hexahedra[i]->writePOS(fp, scalingFactor, (*it)->tag()); for(unsigned int i = 0; i < (*it)->prisms.size(); i++) (*it)->prisms[i]->writePOS(fp, scalingFactor, (*it)->tag()); for(unsigned int i = 0; i < (*it)->pyramids.size(); i++) (*it)->pyramids[i]->writePOS(fp, scalingFactor, (*it)->tag()); } 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, (*it)->tag()); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writePOS(fp, scalingFactor, (*it)->tag()); } 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, (*it)->tag()); } fprintf(fp, "};\n"); } fclose(fp); return 1; } 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 Msg(INFO, "Mesh is in binary format"); 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 double lc = norm(SVector3(bbox.max(), bbox.min())); 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])); } fclose(fp); return 1; } int GModel::writeSTL(const std::string &name, bool binary, double scalingFactor) { FILE *fp = fopen(name.c_str(), binary ? "wb" : "w"); if(!fp){ Msg(GERROR, "Unable to open file '%s'", name.c_str()); return 0; } if(!binary) fprintf(fp, "solid Created by Gmsh\n"); else{ char header[80]; strncpy(header, "Created by Gmsh", 80); fwrite(header, sizeof(char), 80, fp); unsigned int nfacets = 0; for(fiter it = firstFace(); it != lastFace(); ++it) nfacets += (*it)->triangles.size() + 2 * (*it)->quadrangles.size(); fwrite(&nfacets, sizeof(unsigned int), 1, fp); } for(fiter it = firstFace(); it != lastFace(); ++it) { for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeSTL(fp, binary, scalingFactor); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writeSTL(fp, binary, scalingFactor); } if(!binary) fprintf(fp, "endsolid Created by Gmsh\n"); fclose(fp); return 1; } static int skipUntil(FILE *fp, char *key) { char str[256]; while(fscanf(fp, "%s", str)){ if(!strcmp(str, key)) return 1; } return 0; } static int readVerticesVRML(FILE *fp, std::vector<MVertex*> &vertexVector, std::vector<MVertex*> &allVertexVector) { double x, y, z; if(fscanf(fp, " [ %lf %lf %lf", &x, &y, &z) != 3) return 0; vertexVector.push_back(new MVertex(x, y, z)); while(fscanf(fp, " , %lf %lf %lf", &x, &y, &z) == 3) vertexVector.push_back(new MVertex(x, y, z)); for(unsigned int i = 0; i < vertexVector.size(); i++) allVertexVector.push_back(vertexVector[i]); Msg(INFO, "%d vertices", vertexVector.size()); return 1; } static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int region, std::map<int, std::vector<MElement*> > elements[3], bool strips=false) { int i; std::vector<int> idx; if(fscanf(fp, " [ %d", &i) != 1) return 0; idx.push_back(i); while(fscanf(fp, " , %d", &i) == 1){ if(i != -1){ idx.push_back(i); } else{ std::vector<MVertex*> vertices; if(!getVertices(idx.size(), &idx[0], vertexVector, vertices)) return 0; idx.clear(); if(vertices.size() < 2){ Msg(INFO, "Skipping %d-vertex element", (int)vertices.size()); } else if(vertices.size() == 2){ elements[0][region].push_back(new MLine(vertices)); } else if(vertices.size() == 3){ elements[1][region].push_back(new MTriangle(vertices)); } else if(!strips && vertices.size() == 4){ elements[2][region].push_back(new MQuadrangle(vertices)); } else if(strips){ // triangle strip for(unsigned int j = 2; j < vertices.size(); j++){ if(j % 2) elements[1][region].push_back (new MTriangle(vertices[j], vertices[j - 1], vertices[j - 2])); else elements[1][region].push_back (new MTriangle(vertices[j - 2], vertices[j - 1], vertices[j])); } } else{ // import polygon as triangle fan for(unsigned int j = 2; j < vertices.size(); j++){ elements[1][region].push_back (new MTriangle(vertices[0], vertices[j-1], vertices[j])); } } } } if(idx.size()){ Msg(GERROR, "Prematured end of VRML file"); return 0; } Msg(INFO, "%d elements", elements[0][region].size() + elements[1][region].size() + elements[2][region].size()); return 1; } int GModel::readVRML(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; } // This is by NO means a complete VRML/Inventor parser (but it's // sufficient for reading simple Inventor files... which is all I // need) std::vector<MVertex*> vertexVector, allVertexVector; std::map<int, std::vector<MElement*> > elements[3]; int region = 0; char buffer[256], str[256]; while(!feof(fp)) { if(!fgets(buffer, sizeof(buffer), fp)) break; if(buffer[0] != '#'){ // skip comments sscanf(buffer, "%s", str); if(!strcmp(str, "IndexedTriangleStripSet")){ region++; vertexVector.clear(); if(!skipUntil(fp, "vertex")) break; if(!readVerticesVRML(fp, vertexVector, allVertexVector)) break; if(!skipUntil(fp, "coordIndex")) break; if(!readElementsVRML(fp, vertexVector, region, elements, true)) break; } else if(!strcmp(str, "Coordinate3")){ vertexVector.clear(); if(!skipUntil(fp, "point")) break; if(!readVerticesVRML(fp, vertexVector, allVertexVector)) break; } else if(!strcmp(str, "IndexedFaceSet") || !strcmp(str, "IndexedLineSet")){ region++; if(!skipUntil(fp, "coordIndex")) break; if(!readElementsVRML(fp, vertexVector, region, elements)) break; } else if(!strcmp(str, "DEF")){ char str1[256], str2[256]; if(!sscanf(buffer, "%s %s %s", str1, str2, str)) break; if(!strcmp(str, "IndexedFaceSet") || !strcmp(str, "IndexedLineSet")){ region++; if(!skipUntil(fp, "coordIndex")) break; if(!readElementsVRML(fp, vertexVector, region, elements)) break; } } } } for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) storeElementsInEntities(this, elements[i]); associateEntityWithVertices(this); storeVerticesInEntities(allVertexVector); 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; } int GModel::readUNV(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]; std::map<int, MVertex*> vertexMap; std::map<int, std::vector<MElement*> > elements[7]; std::map<int, std::map<int, std::string> > physicals[4]; while(!feof(fp)) { if(!fgets(buffer, sizeof(buffer), fp)) break; if(!strncmp(buffer, " -1", 6)){ if(!fgets(buffer, sizeof(buffer), fp)) break; if(!strncmp(buffer, " -1", 6)) if(!fgets(buffer, sizeof(buffer), fp)) break; int record = 0; sscanf(buffer, "%d", &record); if(record == 2411){ // nodes while(fgets(buffer, sizeof(buffer), fp)){ if(!strncmp(buffer, " -1", 6)) break; int num, dum; if(sscanf(buffer, "%d %d %d %d", &num, &dum, &dum, &dum) != 4) break; if(!fgets(buffer, sizeof(buffer), fp)) break; double x, y, z; for(unsigned int i = 0; i < strlen(buffer); i++) if(buffer[i] == 'D') buffer[i] = 'E'; if(sscanf(buffer, "%lf %lf %lf", &x, &y, &z) != 3) break; vertexMap[num] = new MVertex(x, y, z); } } else if(record == 2412){ // elements while(fgets(buffer, sizeof(buffer), fp)){ if(strlen(buffer) < 3) continue; // possible line ending after last fscanf if(!strncmp(buffer, " -1", 6)) break; int num, type, elementary, physical, color, numNodes; if(sscanf(buffer, "%d %d %d %d %d %d", &num, &type, &elementary, &physical, &color, &numNodes) != 6) break; if(elementary < 0) elementary = 1; if(physical < 0) physical = 0; if(type >= 21 && type <= 24){ // beam elements if(!fgets(buffer, sizeof(buffer), fp)) break; int dum; if(sscanf(buffer, "%d %d %d", &dum, &dum, &dum) != 3) break; } int n[30]; for(int i = 0; i < numNodes; i++) if(!fscanf(fp, "%d", &n[i])) return 0; std::vector<MVertex*> vertices; if(!getVertices(numNodes, n, vertexMap, vertices)) return 0; int dim = -1; switch(type){ case 11: case 21: case 22: case 31: elements[0][elementary].push_back(new MLine(vertices, num)); dim = 1; break; case 23: case 24: case 32: elements[0][elementary].push_back (new MLine3(vertices[0], vertices[2], vertices[1], num)); dim = 1; break; case 41: case 51: case 61: case 74: case 81: case 91: elements[1][elementary].push_back(new MTriangle(vertices, num)); dim = 2; break; case 42: case 52: case 62: case 72: case 82: case 92: elements[1][elementary].push_back (new MTriangle6(vertices[0], vertices[2], vertices[4], vertices[1], vertices[3], vertices[5], num)); dim = 2; break; case 44: case 54: case 64: case 71: case 84: case 94: elements[2][elementary].push_back(new MQuadrangle(vertices, num)); dim = 2; break; case 45: case 55: case 65: case 75: case 85: case 95: elements[2][elementary].push_back (new MQuadrangle8(vertices[0], vertices[2], vertices[4], vertices[6], vertices[1], vertices[3], vertices[5], vertices[7], num)); dim = 2; break; case 111: elements[3][elementary].push_back(new MTetrahedron(vertices, num)); dim = 3; break; case 118: elements[3][elementary].push_back (new MTetrahedron10(vertices[0], vertices[2], vertices[4], vertices[9], vertices[1], vertices[3], vertices[5], vertices[8], vertices[6], vertices[7], num)); dim = 3; break; case 104: case 115: elements[4][elementary].push_back(new MHexahedron(vertices, num)); dim = 3; break; case 105: case 116: elements[4][elementary].push_back (new MHexahedron20(vertices[0], vertices[2], vertices[4], vertices[6], vertices[12], vertices[14], vertices[16], vertices[18], vertices[1], vertices[7], vertices[8], vertices[3], vertices[9], vertices[5], vertices[10], vertices[11], vertices[13], vertices[19], vertices[15], vertices[17], num)); dim = 3; break; case 101: case 112: elements[5][elementary].push_back(new MPrism(vertices, num)); dim = 3; break; case 102: case 113: elements[5][elementary].push_back (new MPrism15(vertices[0], vertices[2], vertices[4], vertices[9], vertices[11], vertices[13], vertices[1], vertices[5], vertices[6], vertices[3], vertices[7], vertices[8], vertices[10], vertices[14], vertices[12], num)); dim = 3; break; } if(dim >= 0 && physical && (!physicals[dim].count(elementary) || !physicals[dim][elementary].count(physical))) physicals[dim][elementary][physical] = "unnamed"; } } } } for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) storeElementsInEntities(this, elements[i]); associateEntityWithVertices(this); storeVerticesInEntities(vertexMap); for(int i = 0; i < 4; i++) storePhysicalTagsInEntities(this, i, physicals[i]); fclose(fp); return 1; } template<class T> static void writeElementsUNV(FILE *fp, const std::vector<T*> &ele, int saveAll, int &num, int elementary, std::vector<int> &physicals) { for(unsigned int i = 0; i < ele.size(); i++) if(saveAll) ele[i]->writeUNV(fp, ++num, elementary, 0); else for(unsigned int j = 0; j < physicals.size(); j++) ele[i]->writeUNV(fp, ++num, elementary, physicals[j]); } int GModel::writeUNV(const std::string &name, 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(noPhysicalGroups()) saveAll = true; renumberMeshVertices(); // nodes fprintf(fp, "%6d\n", -1); fprintf(fp, "%6d\n", 2411); 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); // elements fprintf(fp, "%6d\n", -1); fprintf(fp, "%6d\n", 2412); int num = 0; for(eiter it = firstEdge(); it != lastEdge(); ++it){ writeElementsUNV(fp, (*it)->lines, saveAll, num, (*it)->tag(), (*it)->physicals); } for(fiter it = firstFace(); it != lastFace(); ++it){ writeElementsUNV(fp, (*it)->triangles, saveAll, num, (*it)->tag(), (*it)->physicals); writeElementsUNV(fp, (*it)->quadrangles, saveAll, num, (*it)->tag(), (*it)->physicals); } for(riter it = firstRegion(); it != lastRegion(); ++it){ writeElementsUNV(fp, (*it)->tetrahedra, saveAll, num, (*it)->tag(), (*it)->physicals); writeElementsUNV(fp, (*it)->hexahedra, saveAll, num, (*it)->tag(), (*it)->physicals); writeElementsUNV(fp, (*it)->prisms, saveAll, num, (*it)->tag(), (*it)->physicals); } 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, "Medit mesh import only available for ASCII files"); return 0; } std::vector<MVertex*> vertexVector; std::map<int, std::vector<MElement*> > elements[3]; while(!feof(fp)) { if(!fgets(buffer, sizeof(buffer), fp)) break; if(buffer[0] != '#'){ // skip comments sscanf(buffer, "%s", str); if(!strcmp(str, "Dimension")){ if(!fgets(buffer, sizeof(buffer), fp)) break; } else if(!strcmp(str, "Vertices")){ if(!fgets(buffer, sizeof(buffer), fp)) break; int nbv; sscanf(buffer, "%d", &nbv); Msg(INFO, "%d vertices", nbv); vertexVector.resize(nbv); for(int i = 0; i < nbv; i++) { if(!fgets(buffer, sizeof(buffer), fp)) break; int dum; double x, y, z; sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &dum); vertexVector[i] = new MVertex(x, y, z); } } else if(!strcmp(str, "Triangles")){ if(!fgets(buffer, sizeof(buffer), fp)) break; int nbe; sscanf(buffer, "%d", &nbe); Msg(INFO, "%d triangles", nbe); for(int i = 0; i < nbe; i++) { if(!fgets(buffer, sizeof(buffer), fp)) break; int n[3], cl; sscanf(buffer, "%d %d %d %d", &n[0], &n[1], &n[2], &cl); for(int j = 0; j < 3; j++) n[j]--; std::vector<MVertex*> vertices; if(!getVertices(3, n, vertexVector, vertices)) return 0; elements[0][cl].push_back(new MTriangle(vertices)); } } else if(!strcmp(str, "Quadrilaterals")) { if(!fgets(buffer, sizeof(buffer), fp)) break; int nbe; sscanf(buffer, "%d", &nbe); Msg(INFO, "%d quadrangles", nbe); for(int i = 0; i < nbe; i++) { if(!fgets(buffer, sizeof(buffer), fp)) break; int n[4], cl; sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl); for(int j = 0; j < 4; j++) n[j]--; std::vector<MVertex*> vertices; if(!getVertices(4, n, vertexVector, vertices)) return 0; elements[1][cl].push_back(new MQuadrangle(vertices)); } } else if(!strcmp(str, "Tetrahedra")) { if(!fgets(buffer, sizeof(buffer), fp)) break; int nbe; sscanf(buffer, "%d", &nbe); Msg(INFO, "%d tetrahedra", nbe); for(int i = 0; i < nbe; i++) { if(!fgets(buffer, sizeof(buffer), fp)) break; int n[4], cl; sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl); for(int j = 0; j < 4; j++) n[j]--; std::vector<MVertex*> vertices; if(!getVertices(4, n, vertexVector, vertices)) return 0; elements[2][cl].push_back(new MTetrahedron(vertices)); } } } } for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) storeElementsInEntities(this, elements[i]); associateEntityWithVertices(this); storeVerticesInEntities(vertexVector); 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, numTetrahedra = 0; for(fiter it = firstFace(); it != lastFace(); ++it){ numTriangles += (*it)->triangles.size(); numQuadrangles += (*it)->quadrangles.size(); } for(riter it = firstRegion(); it != lastRegion(); ++it){ numTetrahedra += (*it)->tetrahedra.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()); } if(numTetrahedra){ fprintf(fp, " Tetrahedra\n"); fprintf(fp, " %d\n", numTetrahedra); for(riter it = firstRegion(); it != lastRegion(); ++it) for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++) (*it)->tetrahedra[i]->writeMESH(fp, (*it)->tag()); } fprintf(fp, " End\n"); fclose(fp); return 1; } static int getFormatBDF(char *buffer, int &keySize) { if(buffer[keySize] == '*'){ keySize++; return 2; } // long fields for(unsigned int i = 0; i < strlen(buffer); i++) if(buffer[i] == ',') return 0; // free fields return 1; // small fields; } static double atofBDF(char *str) { int len = strlen(str); // classic numbers with E or D exponent notation for(int i = 0; i < len; i++){ if(str[i] == 'E' || str[i] == 'e') { return atof(str); } else if(str[i] == 'D' || str[i] == 'd'){ str[i] = 'E'; return atof(str); } } // special Nastran floating point format (e.g. "-7.-1" instead of // "-7.E-01" or "2.3+2" instead of "2.3E+02") char tmp[32]; int j = 0, leading_minus = 1; for(int i = 0; i < len; i++){ if(leading_minus && str[i] != ' ' && str[i] != '-') leading_minus = 0; if(!leading_minus && str[i] == '-') tmp[j++] = 'E'; if(str[i] == '+') tmp[j++] = 'E'; tmp[j++] = str[i]; } tmp[j] = '\0'; return atof(tmp); } static int readVertexBDF(FILE *fp, char *buffer, int keySize, int *num, double *x, double *y, double *z) { char tmp[5][32]; int j = keySize; switch(getFormatBDF(buffer, keySize)){ case 0: // free field for(int i = 0; i < 5; i++){ tmp[i][16] = '\0'; strncpy(tmp[i], &buffer[j + 1], 16); for(int k = 0; k < 16; k++){ if(tmp[i][k] == ',') tmp[i][k] = '\0'; } j++; while(j < (int)strlen(buffer) && buffer[j] != ',') j++; } break; case 1: // small field for(int i = 0; i < 5; i++) tmp[i][8] = '\0'; strncpy(tmp[0], &buffer[8], 8); strncpy(tmp[2], &buffer[24], 8); strncpy(tmp[3], &buffer[32], 8); strncpy(tmp[4], &buffer[40], 8); break; case 2: // long field for(int i = 0; i < 5; i++) tmp[i][16] = '\0'; strncpy(tmp[0], &buffer[8], 16); strncpy(tmp[2], &buffer[40], 16); strncpy(tmp[3], &buffer[56], 16); char buffer2[256]; for(unsigned int i = 0; i < sizeof(buffer2); i++) buffer2[i] = '\0'; if(!fgets(buffer2, sizeof(buffer2), fp)) return 0; strncpy(tmp[4], &buffer2[8], 16); break; } *num = atoi(tmp[0]); *x = atofBDF(tmp[2]); *y = atofBDF(tmp[3]); *z = atofBDF(tmp[4]); return 1; } static bool emptyFieldBDF(char *field, int length) { for(int i = 0; i < length; i++) if(field[i] != '\0' && field[i] != ' ' && field[i] != '\n' && field[i] != '\r') return false; return true; } static void readLineBDF(char *buffer, int format, std::vector<char*> &fields) { int cmax = (format == 2) ? 16 : 8; // max char per (center) field int nmax = (format == 2) ? 4 : 8; // max num of (center) fields per line if(format == 0){ // free fields for(unsigned int i = 0; i < strlen(buffer); i++){ if(buffer[i] == ',') fields.push_back(&buffer[i + 1]); } } else{ // small or long fields for(int i = 0; i < nmax + 1; i++){ if(!emptyFieldBDF(&buffer[8 + cmax * i], cmax)) fields.push_back(&buffer[8 + cmax * i]); } } } static int readElementBDF(FILE *fp, char *buffer, int keySize, int numVertices, int &num, int ®ion, std::vector<MVertex*> &vertices, std::map<int, MVertex*> &vertexMap) { char buffer2[256], buffer3[256]; std::vector<char*> fields; int format = getFormatBDF(buffer, keySize); for(unsigned int i = 0; i < sizeof(buffer2); i++) buffer2[i] = buffer3[i] = '\0'; readLineBDF(buffer, format, fields); if((fields.size() - 2 < abs(numVertices)) || (numVertices < 0 && (fields.size() == 9))){ if(fields.size() == 9) fields.pop_back(); if(!fgets(buffer2, sizeof(buffer2), fp)) return 0; readLineBDF(buffer2, format, fields); } if((fields.size() - 2 < abs(numVertices)) || (numVertices < 0 && (fields.size() == 17))){ if(fields.size() == 17) fields.pop_back(); if(!fgets(buffer3, sizeof(buffer3), fp)) return 0; readLineBDF(buffer3, format, fields); } // negative 'numVertices' gives the minimum required number of vertices if(fields.size() - 2 < abs(numVertices)){ Msg(GERROR, "Wrong number of vertices %d for element", fields.size() - 2); return 0; } int n[30], cmax = (format == 2) ? 16 : 8; // max char per (center) field char tmp[32]; tmp[cmax] = '\0'; strncpy(tmp, fields[0], cmax); num = atoi(tmp); strncpy(tmp, fields[1], cmax); region = atoi(tmp); for(unsigned int i = 2; i < fields.size(); i++){ strncpy(tmp, fields[i], cmax); n[i - 2] = atoi(tmp); } // ignore the extra fields when we know how many vertices we need int numCheck = (numVertices > 0) ? numVertices : fields.size() - 2; if(!getVertices(numCheck, n, vertexMap, vertices)) return 0; return 1; } int GModel::readBDF(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]; std::map<int, MVertex*> vertexMap; std::map<int, std::vector<MElement*> > elements[7]; // nodes can be defined after elements, so parse the file twice while(!feof(fp)) { for(unsigned int i = 0; i < sizeof(buffer); i++) buffer[i] = '\0'; if(!fgets(buffer, sizeof(buffer), fp)) break; if(buffer[0] != '$'){ // skip comments if(!strncmp(buffer, "GRID", 4)){ int num; double x, y, z; if(!readVertexBDF(fp, buffer, 4, &num, &x, &y, &z)) break; vertexMap[num] = new MVertex(x, y, z); } } } Msg(INFO, "%d vertices", vertexMap.size()); rewind(fp); while(!feof(fp)) { for(unsigned int i = 0; i < sizeof(buffer); i++) buffer[i] = '\0'; if(!fgets(buffer, sizeof(buffer), fp)) break; if(buffer[0] != '$'){ // skip comments int num, region; std::vector<MVertex*> vertices; if(!strncmp(buffer, "CBAR", 4)){ if(readElementBDF(fp, buffer, 4, 2, num, region, vertices, vertexMap)) elements[0][region].push_back(new MLine(vertices, num)); } else if(!strncmp(buffer, "CROD", 4)){ if(readElementBDF(fp, buffer, 4, 2, num, region, vertices, vertexMap)) elements[0][region].push_back(new MLine(vertices, num)); } else if(!strncmp(buffer, "CBEAM", 5)){ if(readElementBDF(fp, buffer, 5, 2, num, region, vertices, vertexMap)) elements[0][region].push_back(new MLine(vertices, num)); } else if(!strncmp(buffer, "CTRIA3", 6)){ if(readElementBDF(fp, buffer, 6, 3, num, region, vertices, vertexMap)) elements[1][region].push_back(new MTriangle(vertices, num)); } else if(!strncmp(buffer, "CTRIA6", 6)){ if(readElementBDF(fp, buffer, 6, 6, num, region, vertices, vertexMap)) elements[1][region].push_back(new MTriangle6(vertices, num)); } else if(!strncmp(buffer, "CQUAD4", 6)){ if(readElementBDF(fp, buffer, 6, 4, num, region, vertices, vertexMap)) elements[2][region].push_back(new MQuadrangle(vertices, num)); } else if(!strncmp(buffer, "CQUAD8", 6)){ if(readElementBDF(fp, buffer, 6, 8, num, region, vertices, vertexMap)) elements[2][region].push_back(new MQuadrangle8(vertices, num)); } else if(!strncmp(buffer, "CQUAD", 5)){ if(readElementBDF(fp, buffer, 5, -4, num, region, vertices, vertexMap)) if(vertices.size() == 9) elements[2][region].push_back(new MQuadrangle9(vertices, num)); else if(vertices.size() == 8) elements[2][region].push_back(new MQuadrangle8(vertices, num)); else elements[2][region].push_back(new MQuadrangle(vertices, num)); } else if(!strncmp(buffer, "CTETRA", 6)){ if(readElementBDF(fp, buffer, 6, -4, num, region, vertices, vertexMap)) if(vertices.size() == 10) elements[3][region].push_back(new MTetrahedron10(vertices, num)); else elements[3][region].push_back(new MTetrahedron(vertices, num)); } else if(!strncmp(buffer, "CHEXA", 5)){ if(readElementBDF(fp, buffer, 5, -8, num, region, vertices, vertexMap)) if(vertices.size() == 20) elements[4][region].push_back(new MHexahedron20(vertices, num)); else elements[4][region].push_back(new MHexahedron(vertices, num)); } else if(!strncmp(buffer, "CPENTA", 6)){ if(readElementBDF(fp, buffer, 6, -6, num, region, vertices, vertexMap)) if(vertices.size() == 15) elements[5][region].push_back(new MPrism15(vertices, num)); else elements[5][region].push_back(new MPrism(vertices, num)); } else if(!strncmp(buffer, "CPYRAM", 6)){ if(readElementBDF(fp, buffer, 6, 5, num, region, vertices, vertexMap)) elements[6][region].push_back(new MPyramid(vertices, num)); } } } for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) storeElementsInEntities(this, elements[i]); associateEntityWithVertices(this); storeVerticesInEntities(vertexMap); fclose(fp); return 1; } int GModel::writeBDF(const std::string &name, int format, 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; } renumberMeshVertices(); fprintf(fp, "$ Created by Gmsh\n"); for(viter it = firstVertex(); it != lastVertex(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor); for(eiter it = firstEdge(); it != lastEdge(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor); for(fiter it = firstFace(); it != lastFace(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor); for(riter it = firstRegion(); it != lastRegion(); ++it) for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor); for(eiter it = firstEdge(); it != lastEdge(); ++it){ for(unsigned int i = 0; i < (*it)->lines.size(); i++) (*it)->lines[i]->writeBDF(fp, format, (*it)->tag()); } for(fiter it = firstFace(); it != lastFace(); ++it){ for(unsigned int i = 0; i < (*it)->triangles.size(); i++) (*it)->triangles[i]->writeBDF(fp, format, (*it)->tag()); for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++) (*it)->quadrangles[i]->writeBDF(fp, format, (*it)->tag()); } for(riter it = firstRegion(); it != lastRegion(); ++it){ for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++) (*it)->tetrahedra[i]->writeBDF(fp, format, (*it)->tag()); for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++) (*it)->hexahedra[i]->writeBDF(fp, format, (*it)->tag()); for(unsigned int i = 0; i < (*it)->prisms.size(); i++) (*it)->prisms[i]->writeBDF(fp, format, (*it)->tag()); for(unsigned int i = 0; i < (*it)->pyramids.size(); i++) (*it)->pyramids[i]->writeBDF(fp, format, (*it)->tag()); } fprintf(fp, "ENDDATA\n"); fclose(fp); return 1; }