// $Id: GModelIO_Geo.cpp,v 1.17 2008-03-18 08:41:21 remacle Exp $ // // Copyright (C) 1997-2008 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 "GModel.h" #include "Geo.h" #include "OpenFile.h" #include "Tools.h" #include "Numeric.h" #include "Message.h" #include "gmshVertex.h" #include "gmshFace.h" #include "gmshEdge.h" #include "gmshRegion.h" #include "Parser.h" // for Symbol_T int GModel::readGEO(const std::string &name) { ParseFile((char*)name.c_str(), 1); return importGEOInternals(); } void GModel::createGEOInternals() { _geo_internals = new GEO_Internals; } void GModel::deleteGEOInternals() { delete _geo_internals; _geo_internals = 0; } int GModel::importGEOInternals() { if(Tree_Nbr(_geo_internals->Points)) { List_T *points = Tree2List(_geo_internals->Points); for(int i = 0; i < List_Nbr(points); i++){ Vertex *p; List_Read(points, i, &p); GVertex *v = getVertexByTag(p->Num); if(!v){ v = new gmshVertex(this, p); add(v); } if(!p->Visible) v->setVisibility(0); } List_Delete(points); } if(Tree_Nbr(_geo_internals->Curves)) { List_T *curves = Tree2List(_geo_internals->Curves); for(int i = 0; i < List_Nbr(curves); i++){ Curve *c; List_Read(curves, i, &c); if(c->Num >= 0 && c->beg && c->end){ GEdge *e = getEdgeByTag(c->Num); if(!e){ e = new gmshEdge(this, c, getVertexByTag(c->beg->Num), getVertexByTag(c->end->Num)); add(e); } else e->resetMeshAttributes(); if(!c->Visible) e->setVisibility(0); if(c->Color.type) e->setColor(c->Color.mesh); } } List_Delete(curves); } if(Tree_Nbr(_geo_internals->Surfaces)) { List_T *surfaces = Tree2List(_geo_internals->Surfaces); for(int i = 0; i < List_Nbr(surfaces); i++){ Surface *s; List_Read(surfaces, i, &s); GFace *f = getFaceByTag(s->Num); if(!f){ f = new gmshFace(this, s); add(f); } else f->resetMeshAttributes(); if(!s->Visible) f->setVisibility(0); if(s->Color.type) f->setColor(s->Color.mesh); } List_Delete(surfaces); } if(Tree_Nbr(_geo_internals->Volumes)) { List_T *volumes = Tree2List(_geo_internals->Volumes); for(int i = 0; i < List_Nbr(volumes); i++){ Volume *v; List_Read(volumes, i, &v); GRegion *r = getRegionByTag(v->Num); if(!r){ r = new gmshRegion(this, v); add(r); } else r->resetMeshAttributes(); if(!v->Visible) r->setVisibility(0); if(v->Color.type) r->setColor(v->Color.mesh); } List_Delete(volumes); } for(int i = 0; i < List_Nbr(_geo_internals->PhysicalGroups); i++){ PhysicalGroup *p; List_Read(_geo_internals->PhysicalGroups, i, &p); for(int j = 0; j < List_Nbr(p->Entities); j++){ int num; List_Read(p->Entities, j, &num); GEntity *ge = 0; switch(p->Typ){ case MSH_PHYSICAL_POINT: ge = getVertexByTag(abs(num)); break; case MSH_PHYSICAL_LINE: ge = getEdgeByTag(abs(num)); break; case MSH_PHYSICAL_SURFACE: ge = getFaceByTag(abs(num)); break; case MSH_PHYSICAL_VOLUME: ge = getRegionByTag(abs(num)); break; } int pnum = sign(num) * p->Num; if(ge && std::find(ge->physicals.begin(), ge->physicals.end(), pnum) == ge->physicals.end()) ge->physicals.push_back(pnum); } } Msg(DEBUG, "Gmsh model imported:"); Msg(DEBUG, "%d Vertices", vertices.size()); Msg(DEBUG, "%d Edges", edges.size()); Msg(DEBUG, "%d Faces", faces.size()); Msg(DEBUG, "%d Regions", regions.size()); return 1; } class writeFieldOptionGEO{ private : FILE *geo; Field *field; public : writeFieldOptionGEO(FILE *fp,Field *_field) { geo = fp ? fp : stdout; field=_field;} void operator() (std::pair<const char *,FieldOption *> it) { std::string v; it.second->get_text_representation(v); fprintf(geo,"Field[%i].%s = %s;\n",field->id,it.first,v.c_str()); } }; class writeFieldGEO{ private : FILE *geo; public : writeFieldGEO(FILE *fp) { geo = fp ? fp : stdout; } void operator() (std::pair<int, Field *> it) { fprintf(geo,"Field[%i] = %s;\n",it.first,it.second->get_name()); std::for_each(it.second->options.begin(),it.second->options.end(),writeFieldOptionGEO(geo,it.second)); } }; class writeGVertexGEO { private : FILE *geo; public : writeGVertexGEO(FILE *fp) { geo = fp ? fp : stdout; } void operator() (GVertex *gv) { if(gv->getNativeType() == GEntity::GmshModel){ Vertex *v = (Vertex*)gv->getNativePtr(); if(!v) return; fprintf(geo, "Point (%d) = {%.16g, %.16g, %.16g, %.16g};\n", v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, v->lc); } else{ fprintf(geo, "Point (%d) = {%.16g, %.16g, %.16g, %.16g};\n", gv->tag(), gv->x(), gv->y(), gv->z(), gv->prescribedMeshSizeAtVertex()); } } }; class writeGEdgeGEO { private : FILE *geo; public : writeGEdgeGEO(FILE *fp) { geo = fp ? fp : stdout; } void operator () (GEdge *ge) { if(ge->geomType() == GEntity::DiscreteCurve) return; if(ge->getNativeType() == GEntity::GmshModel){ Curve *c = (Curve *)ge->getNativePtr(); if(!c || c->Num < 0) return; switch (c->Typ) { case MSH_SEGM_LINE: fprintf(geo, "Line (%d) = ", c->Num); break; case MSH_SEGM_CIRC: case MSH_SEGM_CIRC_INV: fprintf(geo, "Circle (%d) = ", c->Num); break; case MSH_SEGM_ELLI: case MSH_SEGM_ELLI_INV: fprintf(geo, "Ellipse (%d) = ", c->Num); break; case MSH_SEGM_NURBS: fprintf(geo, "Nurbs (%d) = {", c->Num); for(int i = 0; i < List_Nbr(c->Control_Points); i++) { Vertex *v; List_Read(c->Control_Points, i, &v); if(!i) fprintf(geo, "%d", v->Num); else fprintf(geo, ", %d", v->Num); if(i % 8 == 7 && i != List_Nbr(c->Control_Points) - 1) fprintf(geo, "\n"); } fprintf(geo, "}\n"); fprintf(geo, " Knots {"); for(int j = 0; j < List_Nbr(c->Control_Points) + c->degre + 1; j++) { if(!j) fprintf(geo, "%.16g", c->k[j]); else fprintf(geo, ", %.16g", c->k[j]); if(j % 5 == 4 && j != List_Nbr(c->Control_Points) + c->degre) fprintf(geo, "\n "); } fprintf(geo, "}\n"); fprintf(geo, " Order %d;\n", c->degre); return; case MSH_SEGM_SPLN: fprintf(geo, "CatmullRom (%d) = ", c->Num); break; case MSH_SEGM_BSPLN: fprintf(geo, "BSpline (%d) = ", c->Num); break; case MSH_SEGM_BEZIER: fprintf(geo, "Bezier (%d) = ", c->Num); break; default: Msg(GERROR, "Unknown curve type %d", c->Typ); return; } for(int i = 0; i < List_Nbr(c->Control_Points); i++) { Vertex *v; List_Read(c->Control_Points, i, &v); if(i) fprintf(geo, ", %d", v->Num); else fprintf(geo, "{%d", v->Num); if(i % 6 == 7) fprintf(geo, "\n"); } fprintf(geo, "};\n"); } else{ if(ge->getBeginVertex() && ge->getEndVertex()){ if(ge->geomType() == GEntity::Line){ fprintf(geo, "Line (%d) = {%d, %d};\n", ge->tag(), ge->getBeginVertex()->tag(), ge->getEndVertex()->tag()); } else{ // approximate all other curves by splines Range<double> bounds = ge->parBounds(0); double umin = bounds.low(); double umax = bounds.high(); fprintf(geo, "p%d = newp;\n", ge->tag()); for(int i = 1; i < ge->minimumDrawSegments(); i++){ double u = umin + (double)i / ge->minimumDrawSegments() * (umax - umin); GPoint p = ge->point(u); fprintf(geo, "Point (p%d + %d) = {%.16g, %.16g, %.16g, 1.e+22};\n", ge->tag(), i, p.x(), p.y(), p.z()); } fprintf(geo, "CatmullRom (%d) = {%d", ge->tag(), ge->getBeginVertex()->tag()); for(int i = 1; i < ge->minimumDrawSegments(); i++) fprintf(geo, ", p%d + %d", ge->tag(), i); fprintf(geo, ", %d};\n", ge->getEndVertex()->tag()); } } } } }; class writeGFaceGEO { private : FILE *geo; public : writeGFaceGEO(FILE *fp) { geo = fp ? fp : stdout; } void operator () (GFace *gf) { if(gf->geomType() == GEntity::DiscreteSurface) return; std::list<GEdge*> edges = gf->edges(); std::list<int> orientations = gf->orientations(); if(edges.size() && orientations.size() == edges.size()){ std::vector<int> num, ori; for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++) num.push_back((*it)->tag()); for(std::list<int>::iterator it = orientations.begin(); it != orientations.end(); it++) ori.push_back((*it) > 0 ? 1 : -1); fprintf(geo, "Line Loop (%d) = ", gf->tag()); for(unsigned int i = 0; i < num.size(); i++){ if(i) fprintf(geo, ", %d", num[i] * ori[i]); else fprintf(geo, "{%d", num[i] * ori[i]); } fprintf(geo, "};\n"); if(gf->geomType() == GEntity::Plane){ fprintf(geo, "Plane Surface (%d) = {%d};\n", gf->tag(), gf->tag()); } else if(edges.size() == 3 || edges.size() == 4){ fprintf(geo, "Ruled Surface (%d) = {%d};\n", gf->tag(), gf->tag()); } else{ Msg(GERROR, "Skipping surface %d in export", gf->tag()); } } } }; class writeGRegionGEO { private : FILE *geo; public : writeGRegionGEO(FILE *fp) { geo = fp ? fp : stdout; } void operator () (GRegion *gr) { if(gr->geomType() == GEntity::DiscreteVolume) return; std::list<GFace*> faces = gr->faces(); if(faces.size()){ fprintf(geo, "Surface Loop (%d) = ", gr->tag()); for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++) { if(it != faces.begin()) fprintf(geo, ", %d", (*it)->tag()); else fprintf(geo, "{%d", (*it)->tag()); } fprintf(geo, "};\n"); fprintf(geo, "Volume (%d) = {%d};\n", gr->tag(), gr->tag()); } } }; class writePhysicalGroupGEO { private : FILE *geo; int dim; bool printLabels; std::map<int, std::string> &oldLabels, &newLabels; public : writePhysicalGroupGEO(FILE *fp, int i, bool labels, std::map<int, std::string> &o, std::map<int, std::string> &n) : dim(i), printLabels(labels), oldLabels(o), newLabels(n) { geo = fp ? fp : stdout; } void operator () (std::pair<const int, std::vector<GEntity *> > &g) { std::string oldName, newName; if(printLabels){ if(oldLabels.count(g.first)) { oldName = oldLabels[g.first]; fprintf(geo, "%s = %d;\n", oldName.c_str(), g.first); } else if(newLabels.count(g.first)) { newName = newLabels[g.first]; } } switch (dim) { case 0: fprintf(geo, "Physical Point"); break; case 1: fprintf(geo, "Physical Line"); break; case 2: fprintf(geo, "Physical Surface"); break; case 3: fprintf(geo, "Physical Volume"); break; } if(oldName.size()) fprintf(geo, " (%s) = {", oldName.c_str()); else if(newName.size()) fprintf(geo, " (\"%s\") = {", newName.c_str()); else fprintf(geo, " (%d) = {", g.first); for(unsigned int i = 0; i < g.second.size(); i++) { if(i) fprintf(geo, ", "); fprintf(geo, "%d", g.second[i]->tag()); } fprintf(geo, "};\n"); } }; int GModel::writeGEO(const std::string &name, bool printLabels) { FILE *fp = fopen(name.c_str(), "w"); std::for_each(firstVertex(), lastVertex(), writeGVertexGEO(fp)); std::for_each(firstEdge(), lastEdge(), writeGEdgeGEO(fp)); std::for_each(firstFace(), lastFace(), writeGFaceGEO(fp)); std::for_each(firstRegion(), lastRegion(), writeGRegionGEO(fp)); // get "old-style" labels from parser std::map<int, std::string> labels; List_T *old = Tree2List(Symbol_T); for(int i = 0; i < List_Nbr(old); i++) { Symbol *s = (Symbol *)List_Pointer(old, i); for(int j = 0; j < List_Nbr(s->val); j++) { double tag; List_Read(s->val, j, &tag); labels[(int)tag] = std::string(s->Name); } } List_Delete(old); std::map<int, std::vector<GEntity*> > groups[4]; getPhysicalGroups(groups); for(int i = 0; i < 4; i++) std::for_each(groups[i].begin(), groups[i].end(), writePhysicalGroupGEO(fp, i, printLabels, labels, physicalNames)); std::for_each(fields.begin(),fields.end(), writeFieldGEO(fp)); if(fields.background_field>0)fprintf(fp,"Background Field = %i;\n",fields.background_field); if(fp) fclose(fp); return 1; }