diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp index f4449d24c1dc312fea32afc2a9f3164b6e2c07b3..917723d14d9c876d48adce8dd61b84e0401417a9 100644 --- a/Mesh/BDS.cpp +++ b/Mesh/BDS.cpp @@ -18,7 +18,10 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace *gf) { FILE *f = Fopen(iii, "w"); - if(!f) return; + if(!f){ + Msg::Error("Could not open file '%s'", iii); + return; + } fprintf(f, "View \"scalar\" {\n"); std::list<BDS_Face*>::iterator tit = t.begin(); std::list<BDS_Face*>::iterator tite = t.end(); @@ -323,12 +326,12 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal, if (selfIntersection) return 0; -// if(ix > 300){ -// Msg::Warning("edge %d %d cannot be recovered after %d iterations, trying again", -// num1, num2, ix); -// ix = 0; -// } -// printf("%d %d\n",intersected.size(),ix); + // if(ix > 300){ + // Msg::Warning("edge %d %d cannot be recovered after %d iterations, trying again", + // num1, num2, ix); + // ix = 0; + // } + // printf("%d %d\n",intersected.size(),ix); if(!intersected.size() || ix > 1000){ BDS_Edge *eee = find_edge(num1, num2); diff --git a/Mesh/BGMBase.cpp b/Mesh/BGMBase.cpp index 0f0f6fe95d0f559bc7c96ef83f906b394e782eb1..232102ba2dc7873873a60e5fdd157d340be777ac 100644 --- a/Mesh/BGMBase.cpp +++ b/Mesh/BGMBase.cpp @@ -15,7 +15,12 @@ void BGMBase::export_scalar(const std::string &filename, const DoubleStorageType &_whatToPrint) const { - FILE *f = Fopen (filename.c_str(),"w"); + FILE *f = Fopen(filename.c_str(), "w"); + if(!f){ + Msg::Error("Could not open file '%s'", filename.c_str()); + return; + } + fprintf(f,"View \"Background Mesh\"{\n"); const MElement *elem; @@ -64,7 +69,12 @@ void BGMBase::export_scalar(const std::string &filename, void BGMBase::export_vector(const std::string &filename, const VectorStorageType &_whatToPrint) const { - FILE *f = Fopen (filename.c_str(),"w"); + FILE *f = Fopen(filename.c_str(), "w"); + if(!f){ + Msg::Error("Could not open file '%s'", filename.c_str()); + return; + } + fprintf(f,"View \"Background Mesh\"{\n"); const MElement *elem; @@ -117,7 +127,12 @@ void BGMBase::export_vector(const std::string &filename, void BGMBase::export_tensor_as_vectors(const std::string &filename, const TensorStorageType &_whatToPrint) const { - FILE *f = Fopen (filename.c_str(),"w"); + FILE *f = Fopen(filename.c_str(), "w"); + if(!f){ + Msg::Error("Could not open file '%s'", filename.c_str()); + return; + } + fprintf(f,"View \"Background Mesh\"{\n"); TensorStorageType::const_iterator it = _whatToPrint.begin(); diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp index 8cbb926e45a7abb73cdf41dd665e6fad96daadfb..e943c4d25caac079f1b45086244c509960dfd108 100644 --- a/Mesh/BackgroundMesh.cpp +++ b/Mesh/BackgroundMesh.cpp @@ -648,9 +648,9 @@ double backgroundMesh::getAngle(double u, double v, double w) const void backgroundMesh::print(const std::string &filename, GFace *gf, const std::map<MVertex*,double> &_whatToPrint, int smooth) { - FILE *f = Fopen (filename.c_str(),"w"); + FILE *f = Fopen(filename.c_str(), "w"); if(!f){ - Msg::Warning("Could not open file '%s'", filename.c_str()); + Msg::Error("Could not open file '%s'", filename.c_str()); return; } fprintf(f, "View \"Background Mesh\"{\n"); diff --git a/Mesh/BackgroundMesh2D.cpp b/Mesh/BackgroundMesh2D.cpp index 0c0ce3414d54c4aa30c8b035ab3dad3eba993abb..c023bc0362eba2b0cd628c9a28e9774cc2f17d66 100644 --- a/Mesh/BackgroundMesh2D.cpp +++ b/Mesh/BackgroundMesh2D.cpp @@ -571,7 +571,11 @@ void frameFieldBackgroundMesh2D::computeSmoothness() void frameFieldBackgroundMesh2D::exportCrossField(const std::string &filename) { - FILE *f = Fopen (filename.c_str(),"w"); + FILE *f = Fopen(filename.c_str(), "w"); + if(!f){ + Msg::Error("Could not open file '%s'", filename.c_str()); + return; + } fprintf(f,"View \"Cross Field\"{\n"); std::vector<double> deltas(2); deltas[0] = 0.; diff --git a/Mesh/BackgroundMesh3D.cpp b/Mesh/BackgroundMesh3D.cpp index 8ab79c42c29c9781797105a5d213446ff38d33c3..9ff8296e71ea566f58eeb5a20ad545b12b9e79fa 100644 --- a/Mesh/BackgroundMesh3D.cpp +++ b/Mesh/BackgroundMesh3D.cpp @@ -1139,21 +1139,26 @@ void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filename) { - FILE *f = Fopen (filename.c_str(),"w"); - fprintf(f,"View \"Background Mesh\"{\n"); + FILE *f = Fopen(filename.c_str(), "w"); + if(!f){ + Msg::Error("Could not open file '%s'", filename.c_str()); + return; + } + + fprintf(f, "View \"Background Mesh\"{\n"); set<const MVertex*> done; - for (unsigned int ie=0;ie<getNumMeshElements();ie++){// for all elements + for (unsigned int ie = 0; ie < getNumMeshElements(); ie++){// for all elements const MElement *e = getElement(ie); - for (int iv = 0;iv<e->getNumVertices();iv++){ + for (int iv = 0; iv < e->getNumVertices(); iv++){ const MVertex *v = e->getVertex(iv); set<const MVertex*>::iterator itfind = done.find(v); if (itfind!=done.end()) continue; done.insert(v); STensor3 cf; eval_approximate_crossfield(v,cf); - for (int i=0;i<3;i++){ + for (int i = 0; i < 3; i++){ double vs = get_vectorial_smoothness(i,v); fprintf(f,"VP(%g,%g,%g){%g,%g,%g};\n", v->x(), v->y(), v->z(), cf(0,i)*vs,cf(1,i)*vs,cf(2,i)*vs); diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp index 10a483ae7a832bfa31521066971460e8a7514190..3a743493e53c98ad948cf19a37981a8338c5688b 100644 --- a/Mesh/CenterlineField.cpp +++ b/Mesh/CenterlineField.cpp @@ -575,8 +575,9 @@ void Centerline::computeRadii() void Centerline::buildKdTree() { - FILE * f = Fopen("myPOINTS.pos","w"); - fprintf(f, "View \"\"{\n"); + FILE *f = Fopen("myPOINTS.pos","w"); + + if(f) fprintf(f, "View \"\"{\n"); int nbPL = 3; //10 points per line //int nbNodes = (lines.size()+1) + (nbPL*lines.size()); @@ -586,35 +587,37 @@ void Centerline::buildKdTree() int ind = 0; std::map<MVertex*, int>::iterator itp = colorp.begin(); while (itp != colorp.end()){ - MVertex *v = itp->first; - nodes[ind][0] = v->x(); - nodes[ind][1] = v->y(); - nodes[ind][2] = v->z(); - itp++; ind++; + MVertex *v = itp->first; + nodes[ind][0] = v->x(); + nodes[ind][1] = v->y(); + nodes[ind][2] = v->z(); + itp++; ind++; } for(unsigned int k = 0; k < lines.size(); ++k){ - MVertex *v0 = lines[k]->getVertex(0); - MVertex *v1 = lines[k]->getVertex(1); - SVector3 P0(v0->x(),v0->y(), v0->z()); - SVector3 P1(v1->x(),v1->y(), v1->z()); - for (int j= 1; j < nbPL+1; j++){ - double inc = (double)j/(double)(nbPL+1); - SVector3 Pj = P0+inc*(P1-P0); - nodes[ind][0] = Pj.x(); - nodes[ind][1] = Pj.y(); - nodes[ind][2] = Pj.z(); - ind++; - } - } + MVertex *v0 = lines[k]->getVertex(0); + MVertex *v1 = lines[k]->getVertex(1); + SVector3 P0(v0->x(),v0->y(), v0->z()); + SVector3 P1(v1->x(),v1->y(), v1->z()); + for (int j= 1; j < nbPL+1; j++){ + double inc = (double)j/(double)(nbPL+1); + SVector3 Pj = P0+inc*(P1-P0); + nodes[ind][0] = Pj.x(); + nodes[ind][1] = Pj.y(); + nodes[ind][2] = Pj.z(); + ind++; + } + } - kdtree = new ANNkd_tree(nodes, nbNodes, 3); + kdtree = new ANNkd_tree(nodes, nbNodes, 3); - for(int i = 0; i < nbNodes; ++i){ - fprintf(f, "SP(%g,%g,%g){%g};\n", - nodes[i][0], nodes[i][1],nodes[i][2],1.0); - } - fprintf(f,"};\n"); - fclose(f); + if(f){ + for(int i = 0; i < nbNodes; ++i){ + fprintf(f, "SP(%g,%g,%g){%g};\n", + nodes[i][0], nodes[i][1],nodes[i][2],1.0); + } + fprintf(f,"};\n"); + fclose(f); + } } void Centerline::createSplitCompounds() @@ -1033,9 +1036,10 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad) // char name[256]; // sprintf(name, "myCUT-%d.pos", step); // FILE * f2 = Fopen(name,"w"); - // fprintf(f2, "View \"\"{\n"); - // std::set<MEdge,Less_Edge>::iterator itp = newCut.begin(); - // while (itp != newCut.end()){ + // if(f2){ + // fprintf(f2, "View \"\"{\n"); + // std::set<MEdge,Less_Edge>::iterator itp = newCut.begin(); + // while (itp != newCut.end()){ // MEdge l = *itp; // fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n", // l.getVertex(0)->x(), l.getVertex(0)->y(), l.getVertex(0)->z(), @@ -1045,6 +1049,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad) // } // fprintf(f2,"};\n"); // fclose(f2); + // } } } @@ -1282,46 +1287,51 @@ void Centerline::computeCrossField(double x,double y,double z, void Centerline::printSplit() const { FILE * f = Fopen("mySPLIT.pos","w"); - fprintf(f, "View \"\"{\n"); - for(unsigned int i = 0; i < edges.size(); ++i){ - std::vector<MLine*> lines = edges[i].lines; - for(unsigned int k = 0; k < lines.size(); ++k){ - MLine *l = lines[k]; - fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n", - l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(), - l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(), - (double)edges[i].tag, (double)edges[i].tag); + if(f){ + fprintf(f, "View \"\"{\n"); + for(unsigned int i = 0; i < edges.size(); ++i){ + std::vector<MLine*> lines = edges[i].lines; + for(unsigned int k = 0; k < lines.size(); ++k){ + MLine *l = lines[k]; + fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n", + l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(), + l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(), + (double)edges[i].tag, (double)edges[i].tag); + } } + fprintf(f,"};\n"); + fclose(f); } - fprintf(f,"};\n"); - fclose(f); // FILE * f3 = Fopen("myJUNCTIONS.pos","w"); - // fprintf(f3, "View \"\"{\n"); - // std::set<MVertex*>::const_iterator itj = junctions.begin(); - // while (itj != junctions.end()){ - // MVertex *v = *itj; - // fprintf(f3, "SP(%g,%g,%g){%g};\n", + // if(f3){ + // fprintf(f3, "View \"\"{\n"); + // std::set<MVertex*>::const_iterator itj = junctions.begin(); + // while (itj != junctions.end()){ + // MVertex *v = *itj; + // fprintf(f3, "SP(%g,%g,%g){%g};\n", // v->x(), v->y(), v->z(), // (double)v->getNum()); - // itj++; + // itj++; + // } + // fprintf(f3,"};\n"); + // fclose(f3); // } - // fprintf(f3,"};\n"); - // fclose(f3); FILE * f4 = Fopen("myRADII.pos","w"); - fprintf(f4, "View \"\"{\n"); - for(unsigned int i = 0; i < lines.size(); ++i){ - MLine *l = lines[i]; - std::map<MLine*,double>::const_iterator itc = radiusl.find(l); - fprintf(f4, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n", - l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(), - l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(), - itc->second,itc->second); + if(f4){ + fprintf(f4, "View \"\"{\n"); + for(unsigned int i = 0; i < lines.size(); ++i){ + MLine *l = lines[i]; + std::map<MLine*,double>::const_iterator itc = radiusl.find(l); + fprintf(f4, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n", + l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(), + l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(), + itc->second,itc->second); + } + fprintf(f4,"};\n"); + fclose(f4); } - fprintf(f4,"};\n"); - fclose(f4); - } #endif diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp index 751e0951112769086c2a0b6023e0d2fd4b96527f..e0c69e4c040fe90ab277ef505575094ddceb65b7 100644 --- a/Mesh/DivideAndConquer.cpp +++ b/Mesh/DivideAndConquer.cpp @@ -571,8 +571,12 @@ void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const void DocRecord::makePosView(std::string fileName, GFace *gf) { FILE *f = Fopen(fileName.c_str(),"w"); - if (_adjacencies){ - fprintf(f,"View \"voronoi\" {\n"); + if(!f){ + Msg::Error("Could not open file '%s'", fileName.c_str()); + return; + } + if (_adjacencies){ + fprintf(f, "View \"voronoi\" {\n"); for(PointNumero i = 0; i < numPoints; i++) { std::vector<SPoint2> pts; double pc[2] = {(double)points[i].where.h, (double)points[i].where.v}; @@ -607,7 +611,11 @@ void DocRecord::makePosView(std::string fileName, GFace *gf) void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf, GEdge *ge) { FILE *f = Fopen(fileName.c_str(),"w"); - if (_adjacencies){ + if(!f){ + Msg::Error("Could not open file '%s'", fileName.c_str()); + return; + } + if (_adjacencies){ fprintf(f,"View \"medial axis\" {\n"); for(PointNumero i = 0; i < numPoints; i++) { std::vector<SPoint2> pts; @@ -625,8 +633,8 @@ void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf GPoint p1(pp1.x(), pp1.y(), 0.0); GPoint p2(pp2.x(), pp2.y(), 0.0); if (gf) { - p1 = gf->point(p1.x(), p1.y()); - p2 = gf->point(p2.x(), p2.y()); + p1 = gf->point(p1.x(), p1.y()); + p2 = gf->point(p2.x(), p2.y()); } double P1[3] = {p1.x(), p1.y(), p1.z()}; double P2[3] = {p2.x(), p2.y(), p2.z()}; @@ -642,14 +650,13 @@ void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf p1.x(), p1.y(), p1.z(), p2.x(), p2.y(), p2.z(), (double)i, (double)i); - } + } } - } + } } fprintf(f,"};\n"); } fclose(f); - } void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle, diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp index 99b7150cd59253b9efb82c0107f1971c60e65d37..243be7e5561028fa8589ed8164f03c0a5381cde4 100644 --- a/Mesh/Generator.cpp +++ b/Mesh/Generator.cpp @@ -267,7 +267,13 @@ static void PrintMesh2dStatistics(GModel *m) statreport = Fopen(CTX::instance()->meshStatReportFileName.c_str(), "w"); else if(CTX::instance()->createAppendMeshStatReport == 2) statreport = Fopen(CTX::instance()->meshStatReportFileName.c_str(), "a"); - else return; + else + return; + + if(statreport){ + Msg::Error("Could not open file '%s'", CTX::instance()->meshStatReportFileName.c_str()); + return; + } double worst = 1, best = 0, avg = 0; double e_long = 0, e_short = 1.e22, e_avg = 0; diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp index 3a21106316dbf9237c64a7857ef5480c26796c51..ce39d1f7b02d1302ad66ae1f8a6a6596a1f98792 100644 --- a/Mesh/directions3D.cpp +++ b/Mesh/directions3D.cpp @@ -33,27 +33,27 @@ void Frame_field::init_region(GRegion* gr){ GFace* gf; std::list<GFace*> faces; std::list<GFace*>::iterator it; - + Nearest_point::init_region(gr); - + faces = gr->faces(); - + field.clear(); labels.clear(); - + for(it=faces.begin();it!=faces.end();it++){ gf = *it; init_face(gf); } - + ANNpointArray duplicate = annAllocPts(field.size(),3); - + for(i=0;i<field.size();i++){ duplicate[i][0] = field[i].first.x(); duplicate[i][1] = field[i].first.y(); duplicate[i][2] = field[i].first.z(); } - + kd_tree = new ANNkd_tree(duplicate,field.size(),3); #endif } @@ -67,12 +67,12 @@ void Frame_field::init_face(GFace* gf){ SVector3 v2; SVector3 v3; STensor3 m(1.0); - + for(i=0;i<gf->storage1.size();i++){ point = gf->storage1[i]; v1 = gf->storage2[i]; v2 = gf->storage3[i]; - + v1.normalize(); v2.normalize(); v3 = crossprod(v1,v2); @@ -98,7 +98,7 @@ STensor3 Frame_field::search(double x,double y,double z){ double distance1; double distance2; double e2; - + index1 = 0; index2 = 0; distance1 = 1000000.0; @@ -108,31 +108,31 @@ STensor3 Frame_field::search(double x,double y,double z){ ANNpoint query; ANNidxArray indices; ANNdistArray distances; - + if(field.size()<=1){ return STensor3(1.0); } - + query = annAllocPt(3); query[0] = x; query[1] = y; query[2] = z; - + indices = new ANNidx[2]; distances = new ANNdist[2]; - + double e = 0.0; kd_tree->annkSearch(query,2,indices,distances,e); index1 = indices[0]; index2 = indices[1]; distance1 = distances[0]; distance2 = distances[1]; - + annDeallocPt(query); delete[] indices; delete[] distances; #endif - + if(fabs(sqrt(distance2)-sqrt(distance1))<e2){ if(labels[index2]<labels[index1]){ return field[index2].second; @@ -155,21 +155,21 @@ STensor3 Frame_field::combine(double x,double y,double z){ SVector3 vec1,vec2,vec3; SVector3 final1,final2; STensor3 m(1.0),m2(1.0); - + m = search(x,y,z); m2 = m; ok = Nearest_point::search(x,y,z,vec); vec.normalize(); - + if(ok){ vec1 = SVector3(m.get_m11(),m.get_m21(),m.get_m31()); vec2 = SVector3(m.get_m12(),m.get_m22(),m.get_m32()); vec3 = SVector3(m.get_m13(),m.get_m23(),m.get_m33()); - + val1 = fabs(dot(vec,vec1)); val2 = fabs(dot(vec,vec2)); val3 = fabs(dot(vec,vec3)); - + if(val1<=val2 && val1<=val3){ other = vec1; } @@ -179,12 +179,12 @@ STensor3 Frame_field::combine(double x,double y,double z){ else{ other = vec3; } - + final1 = crossprod(vec,other); final1.normalize(); final2 = crossprod(vec,final1); final2.normalize(); - + m2.set_m11(vec.x()); m2.set_m21(vec.y()); m2.set_m31(vec.z()); @@ -195,7 +195,7 @@ STensor3 Frame_field::combine(double x,double y,double z){ m2.set_m23(final2.y()); m2.set_m33(final2.z()); } - + return m2; } @@ -216,17 +216,17 @@ void Frame_field::print_field1(){ SPoint3 point; SPoint3 p1,p2,p3,p4,p5,p6; STensor3 m(1.0); - + k = 0.05; std::ofstream file("frame1.pos"); file << "View \"cross field\" {\n"; color1 = 10.0; color2 = 20.0; - + for(i=0;i<field.size();i++){ point = field[i].first; m = field[i].second; - + p1 = SPoint3(point.x() + k*m.get_m11(), point.y() + k*m.get_m21(), point.z() + k*m.get_m31()); @@ -245,7 +245,7 @@ void Frame_field::print_field1(){ p6 = SPoint3(point.x() - k*m.get_m13(), point.y() - k*m.get_m23(), point.z() - k*m.get_m33()); - + print_segment(point,p1,color1,color2,file); print_segment(point,p2,color1,color2,file); print_segment(point,p3,color1,color2,file); @@ -253,7 +253,7 @@ void Frame_field::print_field1(){ print_segment(point,p5,color1,color2,file); print_segment(point,p6,color1,color2,file); } - + file << "};\n"; } @@ -269,13 +269,13 @@ void Frame_field::print_field2(GRegion* gr){ MVertex* vertex; MElement* element; STensor3 m(1.0); - + k = 0.05; std::ofstream file("frame2.pos"); file << "View \"cross field\" {\n"; color1 = 10.0; color2 = 20.0; - + for(i=0;i<gr->getNumMeshElements();i++){ element = gr->getMeshElement(i); for(j=0;j<element->getNumVertices();j++){ @@ -283,7 +283,7 @@ void Frame_field::print_field2(GRegion* gr){ if(vertex->onWhat()->dim()>2){ point = SPoint3(vertex->x(),vertex->y(),vertex->z()); m = search(vertex->x(),vertex->y(),vertex->z()); - + p1 = SPoint3(point.x() + k*m.get_m11(), point.y() + k*m.get_m21(), point.z() + k*m.get_m31()); @@ -302,7 +302,7 @@ void Frame_field::print_field2(GRegion* gr){ p6 = SPoint3(point.x() - k*m.get_m13(), point.y() - k*m.get_m23(), point.z() - k*m.get_m33()); - + print_segment(point,p1,color1,color2,file); print_segment(point,p2,color1,color2,file); print_segment(point,p3,color1,color2,file); @@ -312,7 +312,7 @@ void Frame_field::print_field2(GRegion* gr){ } } } - + file << "};\n"; } @@ -857,8 +857,8 @@ void Frame_field::recur_connect_vert(FILE *fi, int count, MVertex *v, STensor3 &cross, std::multimap<MVertex*,MVertex*> &v2v, - std::set<MVertex*> &touched){ - + std::set<MVertex*> &touched) +{ if (touched.find(v) != touched.end()) return; touched.insert(v); @@ -870,70 +870,71 @@ void Frame_field::recur_connect_vert(FILE *fi, int count, MVertex *nextV = it->second; if (touched.find(nextV) == touched.end()){ - //compute dot product (N0,R0,A0) dot (Ni,Ri,Ai)^T - //where N,R,A are the 3 directions - std::map<MVertex*, STensor3>::iterator iter = crossField.find(nextV); - STensor3 nextCross = iter->second; - STensor3 nextCrossT = nextCross.transpose(); - STensor3 prod = cross.operator*=(nextCrossT); - fullMatrix<double> mat(3,3); prod.getMat(mat); - - //find biggest dot product - fullVector<int> Id(3); - Id(0) = Id(1) = Id(2) = 0; - for (int j = 0; j < 3; j++){ - double maxVal = 0.0; - for (int i = 0; i < 3; i++){ - double val = fabs(mat(i,j)); - if( val > maxVal ){ - maxVal = val; - Id(j) = i; - } + //compute dot product (N0,R0,A0) dot (Ni,Ri,Ai)^T + //where N,R,A are the 3 directions + std::map<MVertex*, STensor3>::iterator iter = crossField.find(nextV); + STensor3 nextCross = iter->second; + STensor3 nextCrossT = nextCross.transpose(); + STensor3 prod = cross.operator*=(nextCrossT); + fullMatrix<double> mat(3,3); prod.getMat(mat); + + //find biggest dot product + fullVector<int> Id(3); + Id(0) = Id(1) = Id(2) = 0; + for (int j = 0; j < 3; j++){ + double maxVal = 0.0; + for (int i = 0; i < 3; i++){ + double val = fabs(mat(i,j)); + if( val > maxVal ){ + maxVal = val; + Id(j) = i; + } + } } - } - //check - if (Id(0) +Id(1)+ Id(2) != 3 || (Id(0) == 1 && Id(1)==1 && Id(2)==1)) { - std::cout << "This should not happen: sum should be 0+1+2" << std::endl; - printf("Id =%d %d %d \n", Id(0), Id(1), Id(2)); - return; - } + //check + if (Id(0) +Id(1)+ Id(2) != 3 || (Id(0) == 1 && Id(1)==1 && Id(2)==1)) { + std::cout << "This should not happen: sum should be 0+1+2" << std::endl; + printf("Id =%d %d %d \n", Id(0), Id(1), Id(2)); + return; + } - //create new cross - fullMatrix<double> newmat(3,3); - for (int i = 0; i < 3; i++){ - for (int j = 0; j < 3; j++){ - newmat(i,j) = nextCross(Id(i),j) ; - } - } + //create new cross + fullMatrix<double> newmat(3,3); + for (int i = 0; i < 3; i++){ + for (int j = 0; j < 3; j++){ + newmat(i,j) = nextCross(Id(i),j) ; + } + } - STensor3 newcross(0.0); - newcross.setMat(newmat); - crossField[iter->first] = newcross; - - //print info - printf("************** COUNT = %d \n", count); - if (Id(0) == 0 && Id(1) == 1 && Id(2) == 2) - printf("orient OK Id=%d %d %d\n", Id(0), Id(1), Id(2)); - else{ - printf("change orientation Id=%d %d %d \n", Id(0), Id(1), Id(2)); - cross.print("cross "); - nextCross.print("nextCross "); - prod.print("product"); - newcross.print("newcross"); - } + STensor3 newcross(0.0); + newcross.setMat(newmat); + crossField[iter->first] = newcross; + + //print info + printf("************** COUNT = %d \n", count); + if (Id(0) == 0 && Id(1) == 1 && Id(2) == 2) + printf("orient OK Id=%d %d %d\n", Id(0), Id(1), Id(2)); + else{ + printf("change orientation Id=%d %d %d \n", Id(0), Id(1), Id(2)); + cross.print("cross "); + nextCross.print("nextCross "); + prod.print("product"); + newcross.print("newcross"); + } - fprintf(fi,"SP(%g,%g,%g) {%g};\n",nextV->x(),nextV->y(),nextV->z(), (double)count); + if(fi) + fprintf(fi,"SP(%g,%g,%g) {%g};\n",nextV->x(),nextV->y(),nextV->z(), (double)count); - //continue recursion - recur_connect_vert (fi, count, nextV, newcross, v2v,touched); + //continue recursion + recur_connect_vert (fi, count, nextV, newcross, v2v,touched); } } - } -void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){ +void Frame_field::continuousCrossField(GRegion *gr, GFace *gf) +{ printf("continuous cross field \n"); //start from a vertex of a face @@ -969,14 +970,16 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){ STensor3 bCross = iter->second; FILE *fi = Fopen ("cross_recur.pos","w"); - fprintf(fi,"View \"\"{\n"); - fprintf(fi,"SP(%g,%g,%g) {%g};\n",beginV->x(),beginV->y(),beginV->z(), 0.0); + if(fi){ + fprintf(fi,"View \"\"{\n"); + fprintf(fi,"SP(%g,%g,%g) {%g};\n",beginV->x(),beginV->y(),beginV->z(), 0.0); + } int count = 0; - - recur_connect_vert (fi, count, beginV,bCross,v2v,touched); - - fprintf(fi,"};\n"); - fclose (fi); + recur_connect_vert(fi, count, beginV,bCross,v2v,touched); + if(fi){ + fprintf(fi,"};\n"); + fclose (fi); + } //printf("touched =%d vert =%d \n", touched.size(), vertex_to_vertices.size()); } @@ -1141,9 +1144,9 @@ void Size_field::init_region(GRegion* gr){ ANNpoint query; ANNidxArray indices; ANNdistArray distances; - + faces = gr->faces(); - + field.clear(); for(it=faces.begin();it!=faces.end();it++){ @@ -1158,47 +1161,47 @@ void Size_field::init_region(GRegion* gr){ } ANNpointArray duplicate = annAllocPts(field.size(),3); - + for(i=0;i<field.size();i++){ duplicate[i][0] = field[i].first.x(); duplicate[i][1] = field[i].first.y(); duplicate[i][2] = field[i].first.z(); } - + kd_tree = new ANNkd_tree(duplicate,field.size(),3); boundary.clear(); - + query = annAllocPt(3); indices = new ANNidx[1]; distances = new ANNdist[1]; - + index = 0; e = 0.0; - + for(it=faces.begin();it!=faces.end();it++){ gf = *it; - + for(i=0;i<gf->getNumMeshElements();i++){ - element = gf->getMeshElement(i); - + element = gf->getMeshElement(i); + for(j=0;j<element->getNumVertices();j++){ vertex = element->getVertex(j); - + query[0] = vertex->x(); query[1] = vertex->y(); query[2] = vertex->z(); - + kd_tree->annkSearch(query,1,indices,distances,e); index = indices[0]; - + boundary.insert(std::pair<MVertex*,double>(vertex,field[index].second)); } } } - + octree = new MElementOctree(model); - + annDeallocPt(query); delete[] indices; delete[] distances; @@ -1416,7 +1419,7 @@ void Size_field::clear(){ delete kd_tree; annClose(); #endif - + } /****************class Nearest_point****************/ diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp index 33bd7df668c527054b3b51e826dc2e546bb8a724..8ca50b2b48580a0db25be165edc862447843fb33 100644 --- a/Mesh/meshGEdge.cpp +++ b/Mesh/meshGEdge.cpp @@ -318,7 +318,8 @@ static void printFandPrimitive(int tag, std::vector<IntPoint> &Points) */ // new algo for recombining + splitting -static int increaseN (int N){ +static int increaseN (int N) +{ if (((N+1)/2 - 1) % 2 != 0) return N+2; return N; } diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index 9494fe7ff0684fea42dba1a16552308ef2c1c4c2..2e700563277b98622850e55692f857bd07bf9ca4 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -779,7 +779,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) edges.insert(edges.begin(), embedded_edges.begin(),embedded_edges.end()); std::list<GEdge*>::iterator ite = edges.begin(); FILE *ff2 = Fopen ("tato.pos","w"); - fprintf(ff2,"View \" \"{\n"); + if(ff2) fprintf(ff2,"View \" \"{\n"); std::set<MVertex*> verts; while(ite != edges.end()){ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){ @@ -812,13 +812,14 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) MQuadrangle *qq = new MQuadrangle(v11,v21,v22,v12); myCol.push_back(qq); blQuads.push_back(qq); - fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n", - v11->x(),v11->y(),v11->z(), - v12->x(),v12->y(),v12->z(), - v22->x(),v22->y(),v22->z(), - v21->x(),v21->y(),v21->z()); + if(ff2) + fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n", + v11->x(),v11->y(),v11->z(), + v12->x(),v12->y(),v12->z(), + v22->x(),v22->y(),v22->z(), + v21->x(),v21->y(),v21->z()); } - // int M = std::max(c1._column.size(),c2._column.size()); + // int M = std::max(c1._column.size(),c2._column.size()); for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0]; _columns->_elemColumns[myCol[0]] = myCol; } @@ -852,20 +853,22 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) MQuadrangle *qq = new MQuadrangle(v11,v12,v22,v21); myCol.push_back(qq); blQuads.push_back(qq); - fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n", - v11->x(),v11->y(),v11->z(), - v12->x(),v12->y(),v12->z(), - v22->x(),v22->y(),v22->z(), - v21->x(),v21->y(),v21->z()); + if(ff2) + fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n", + v11->x(),v11->y(),v11->z(), + v12->x(),v12->y(),v12->z(), + v22->x(),v22->y(),v22->z(), + v21->x(),v21->y(),v21->z()); } else { MTriangle *qq = new MTriangle(v,v22,v21); myCol.push_back(qq); blTris.push_back(qq); - fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n", - v->x(),v->y(),v->z(), - v22->x(),v22->y(),v22->z(), - v21->x(),v21->y(),v21->z()); + if(ff2) + fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n", + v->x(),v->y(),v->z(), + v22->x(),v22->y(),v22->z(), + v21->x(),v21->y(),v21->z()); } } } @@ -873,8 +876,10 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) _columns->_elemColumns[myCol[0]] = myCol; } - fprintf(ff2,"};\n"); - fclose(ff2); + if(ff2){ + fprintf(ff2,"};\n"); + fclose(ff2); + } // filterOverlappingElements (blTris,blQuads,_columns->_elemColumns,_columns->_toFirst); @@ -900,16 +905,18 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) std::set<MEdge,Less_Edge>::iterator it = bedges.begin(); FILE *ff = Fopen ("toto.pos","w"); - fprintf(ff,"View \" \"{\n"); + if(ff) fprintf(ff,"View \" \"{\n"); for (; it != bedges.end(); ++it){ ne.lines.push_back(new MLine (it->getVertex(0),it->getVertex(1))); - fprintf(ff,"SL (%g,%g,%g,%g,%g,%g){1,1};\n", - it->getVertex(0)->x(),it->getVertex(0)->y(),it->getVertex(0)->z(), - it->getVertex(1)->x(),it->getVertex(1)->y(),it->getVertex(1)->z()); + if(ff) + fprintf(ff,"SL (%g,%g,%g,%g,%g,%g){1,1};\n", + it->getVertex(0)->x(),it->getVertex(0)->y(),it->getVertex(0)->z(), + it->getVertex(1)->x(),it->getVertex(1)->y(),it->getVertex(1)->z()); + } + if(ff){ + fprintf(ff,"};\n"); + fclose(ff); } - fprintf(ff,"};\n"); - fclose(ff); - hop.push_back(&ne); deMeshGFace kil_; @@ -1359,7 +1366,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, } - // effectively recover the medge + // effectively recover the medge ite = edges.begin(); while(ite != edges.end()){ if(!(*ite)->isMeshDegenerated()){ diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp index b0fa51697012160d1a3c3a4b36aa05dcd8c3552a..14cad8fad2278fd4bce86292988023faba3f9cdb 100644 --- a/Mesh/meshGFaceDelaunayInsertion.cpp +++ b/Mesh/meshGFaceDelaunayInsertion.cpp @@ -9,11 +9,9 @@ #include "GmshMessage.h" #include "OS.h" #include "robustPredicates.h" -// PEB MODIF #include "BackgroundMesh.h" #include "surfaceFiller.h" #include "pointInsertion.h" -// END PEB MODIF #include "meshGFaceDelaunayInsertion.h" #include "meshGFaceOptimize.h" #include "meshGFace.h" @@ -29,28 +27,20 @@ #include "intersectCurveSurface.h" #include "HilbertCurve.h" -double LIMIT_ = 0.5 * sqrt(2.0) * 1; -int N_GLOBAL_SEARCH; -int N_SEARCH; -double DT_INSERT_VERTEX; +static double LIMIT_ = 0.5 * sqrt(2.0) * 1; +static int N_GLOBAL_SEARCH; +static int N_SEARCH; +static double DT_INSERT_VERTEX; int MTri3::radiusNorm = 2; -/* -static bool isBoundary(MTri3 *t, double limit_, int &active) -{ - if (t->isDeleted()) return false; - for (active = 0; active < 3; active++){ - MTri3 *neigh = t->getNeigh(active); - if (!neigh) return true; - } - return false; -} -*/ - template <class ITERATOR> void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData * data) { FILE *ff = Fopen (name,"w"); + if(!ff){ + Msg::Error("Could not open file '%s'", name); + return; + } fprintf(ff,"View\"test\"{\n"); while ( it != end ){ MTri3 *worst = *it; @@ -91,7 +81,6 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData * data) fclose (ff); } - static bool isActive(MTri3 *t, double limit_, int &active) { if (t->isDeleted()) return false; @@ -1535,9 +1524,9 @@ void buildBackGroundMesh (GFace *gf, std::map<MVertex* , MVertex*>* equivalence, std::map<MVertex*, SPoint2> * parametricCoordinates) { - // TODO PEB : - // this is now done in the new backgroundMesh !!! - // on le vire ? On insère ici les operations sur le new backgmesh ? + // TODO PEB : + // this is now done in the new backgroundMesh !!! + // on le vire ? On insère ici les operations sur le new backgmesh ? // parce que les opérations qui suivent sur l'ancien BGM sont inutiles maintenant... ! if (!old_algo_hexa()) return; @@ -1969,26 +1958,19 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad, double t2 = Cpu(); Msg::Debug("Delaunay 2D done for %d points : CPU = %g, %d global searches, AVG walk size %g , AVG cavity size %g", v.size(), t2-t1,NB_GLOBAL_SEARCH,1.+AVG_ITER/v.size(),AVG_CAVSIZE/v.size()); - // printf("%g %g %g %g --> %g(%g)\n",ta,tb,tc,td,t2-t1,ta+tb+tc+td); if (edgesToRecover)recoverEdges (t,*edgesToRecover); - // FILE *f = fopen ("tri.pos","w"); - // fprintf(f,"View \"\"{\n"); for (size_t i = 0;i<t.size();i++){ if (t[i]->isDeleted() || (removeBox && triOnBox (t[i]->tri(),box))) delete t[i]->tri(); else { result.push_back(t[i]->tri()); - // t[i]->tri()->writePOS (f, false,false,true,false,false,false); } delete t[i]; } if (removeBox){for (int i=0;i<4;i++)delete box[i];} else {for (int i=0;i<4;i++)v.push_back(box[i]);} - - // fprintf(f,"};\n"); - // fclose(f); } bool swapedge (MVertex *v1 ,MVertex *v2, MVertex *v3, MVertex *v4, MTri3* t1, int iLocalEdge){ diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp index c0c5ee2de62878390bef065493cd169406091116..99652442405669c271f7a3a78927f70547ba79f0 100644 --- a/Mesh/meshGFaceElliptic.cpp +++ b/Mesh/meshGFaceElliptic.cpp @@ -34,13 +34,18 @@ (1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4) static void printQuads(GFace *gf, fullMatrix<SPoint2> uv, - std::vector<MQuadrangle*> quads, int iter){ - + std::vector<MQuadrangle*> quads, int iter) +{ if(!CTX::instance()->mesh.saveAll) return; char name[234]; - sprintf(name,"quadUV_%d_%d.pos", gf->tag(), iter); - FILE *f = Fopen(name,"w"); + sprintf(name, "quadUV_%d_%d.pos", gf->tag(), iter); + FILE *f = Fopen(name, "w"); + if(!f){ + Msg::Error("Could not open file '%s'", name); + return; + } + fprintf(f,"View \"%s\" {\n",name); for (int i = 1; i < uv.size1()-1; i++) @@ -51,15 +56,19 @@ static void printQuads(GFace *gf, fullMatrix<SPoint2> uv, fclose(f); char name3[234]; - sprintf(name3,"quadXYZ_%d_%d.pos", gf->tag(), iter); + sprintf(name3, "quadXYZ_%d_%d.pos", gf->tag(), iter); FILE *f3 = Fopen(name3,"w"); + if(!f3){ + Msg::Error("Could not open file '%s'", name3); + return; + } + fprintf(f3,"View \"%s\" {\n",name3); for (unsigned int i = 0; i < quads.size(); i++){ quads[i]->writePOS(f3,true,false,false,false,false,false); } fprintf(f3,"};\n"); fclose(f3); - } static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<MVertex*> vert2, @@ -83,10 +92,14 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M p2.push_back(pj); } - char name[234]; sprintf(name,"paramGrid_%d.pos", gf->tag()); FILE *f = Fopen(name,"w"); + if(!f){ + Msg::Error("Could not open file '%s'", name); + return; + } + fprintf(f,"View \"%s\" {\n",name); // for (unsigned int i = 0; i < p1.size(); i++) @@ -107,8 +120,12 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M char name2[234]; sprintf(name2,"paramEdges_%d.pos", gf->tag()); FILE *f2 = Fopen(name2,"w"); - fprintf(f2,"View \"%s\" {\n",name2); + if(!f2){ + Msg::Error("Could not open file '%s'", name2); + return; + } + fprintf(f2,"View \"%s\" {\n",name2); for (unsigned int i = 0; i < e01.size(); i++){ SPoint2 pi; reparamMeshVertexOnFace(e01[i], gf, pi); fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 1); @@ -144,21 +161,23 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M char name3[234]; sprintf(name3,"quadXYZ_%d.pos", gf->tag()); FILE *f3 = Fopen(name3,"w"); + if(!f3){ + Msg::Error("Could not open file '%s'", name2); + return; + } + fprintf(f3,"View \"%s\" {\n",name3); for (unsigned int i = 0; i < quads.size(); i++){ quads[i]->writePOS(f3,true,false,false,false,false,false); } fprintf(f3,"};\n"); fclose(f3); - - return; - } static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv, std::vector<std::vector<MVertex*> > &tab, - std::vector<MQuadrangle*> &newq, std::vector<MVertex*> &newv){ - + std::vector<MQuadrangle*> &newq, std::vector<MVertex*> &newv) +{ newq.clear(); newv.clear(); @@ -193,8 +212,8 @@ static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv, } static std::vector<MVertex*> saturateEdgeRegular (GFace *gf, SPoint2 p1, SPoint2 p2, - int M, std::vector<SPoint2> &pe){ - + int M, std::vector<SPoint2> &pe) +{ std::vector<MVertex*> pts; for (int i=1;i<M;i++){ double s = ((double)i/((double)(M))); @@ -210,10 +229,11 @@ static std::vector<MVertex*> saturateEdgeRegular (GFace *gf, SPoint2 p1, SPoint2 } return pts; } + static std::vector<MVertex*> saturateEdgeHarmonic (GFace *gf, SPoint2 p1, SPoint2 p2, double H, double L, - int M, std::vector<SPoint2> &pe){ - + int M, std::vector<SPoint2> &pe) +{ std::vector<MVertex*> pts; for (int i=1;i<M;i++){ double y = ((double)(i))*H/M; @@ -239,13 +259,13 @@ static void transfiniteSmoother(GFace* gf, std::vector<std::vector<MVertex*> > &tab, std::vector<MQuadrangle*> &newq, std::vector<MVertex*> &newv, - bool isPeriodic=false){ - - int M = uv.size1(); - int N = uv.size2(); - int jStart = isPeriodic ? 0 : 1; + bool isPeriodic=false) +{ + int M = uv.size1(); + int N = uv.size2(); + int jStart = isPeriodic ? 0 : 1; - int numSmooth = 150; + int numSmooth = 150; fullMatrix<SPoint2> uvold = uv; for(int k = 0; k < numSmooth; k++){ double norm = 0.0; @@ -293,7 +313,6 @@ static void transfiniteSmoother(GFace* gf, //final print createQuadsFromUV(gf, uv, tab, newq, newv); printQuads(gf, uv, newq, numSmooth); - } //elliptic surface grid generator @@ -303,8 +322,8 @@ static void ellipticSmoother( GFace* gf, std::vector<std::vector<MVertex*> > &tab, std::vector<MQuadrangle*> &newq, std::vector<MVertex*> &newv, - bool isPeriodic=false){ - + bool isPeriodic=false) +{ printQuads(gf, uv, newq, 0); int nbSmooth = 70; @@ -395,9 +414,6 @@ static void ellipticSmoother( GFace* gf, createQuadsFromUV(gf, uv, tab, newq, newv); printQuads(gf, uv, newq, nbSmooth); - - //exit(1); - } //create initial grid points MxN using transfinite interpolation @@ -427,7 +443,6 @@ static void createRegularGrid (GFace *gf, fullMatrix<SPoint2> &uv, std::vector<std::vector<MVertex*> > &tab) { - int M = e23.size(); int N = e12.size(); @@ -528,7 +543,8 @@ static void createRegularGridPeriodic (GFace *gf,int sign2, char name3[234]; sprintf(name3,"quadParam_%d.pos", gf->tag()); FILE *f3 = Fopen(name3,"w"); - fprintf(f3,"View \"%s\" {\n",name3); + + if(f3) fprintf(f3,"View \"%s\" {\n",name3); tab.resize(M+2); for(int i = 0; i < M+2; i++) tab[i].resize(N+2); @@ -577,24 +593,24 @@ static void createRegularGridPeriodic (GFace *gf,int sign2, } //print - for (int i=0;i<N+2;i++) - for (int j=0;j<M+2;j++) - fprintf(f3,"SP(%g,%g,%g) {%d};\n", uv(j,i).x(), uv(j,i).y(), 0.0, j); - - fprintf(f3,"};\n"); - fclose(f3); + if(f3){ + for (int i=0;i<N+2;i++) + for (int j=0;j<M+2;j++) + fprintf(f3,"SP(%g,%g,%g) {%d};\n", uv(j,i).x(), uv(j,i).y(), 0.0, j); + fprintf(f3,"};\n"); + fclose(f3); + } } -static void updateFaceQuads(GFace *gf, std::vector<MQuadrangle*> &quads, std::vector<MVertex*> &newv){ - +static void updateFaceQuads(GFace *gf, std::vector<MQuadrangle*> &quads, std::vector<MVertex*> &newv) +{ for (unsigned int i = 0; i < quads.size(); i++){ gf->quadrangles.push_back(quads[i]); } for(unsigned int i = 0; i < newv.size(); i++){ gf->mesh_vertices.push_back(newv[i]); } - } static bool computeRingVertices(GFace *gf, Centerline *center, diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp index a061d38082c9a4d9a2d4feec9936788032ef646a..a74eb44f4b71333c5f07ede956f1d3fa93744a6d 100644 --- a/Mesh/meshGFaceOptimize.cpp +++ b/Mesh/meshGFaceOptimize.cpp @@ -1020,10 +1020,10 @@ void createRegularMesh (GFace *gf, int N = e12.size(); int M = e23.size(); - char name3[234]; + char name3[234]; sprintf(name3,"quadParam_%d.pos", gf->tag()); FILE *f3 = Fopen(name3,"w"); - fprintf(f3,"View \"%s\" {\n",name3); + if(f3) fprintf(f3,"View \"%s\" {\n",name3); //create points using transfinite interpolation @@ -1059,7 +1059,7 @@ void createRegularMesh (GFace *gf, c1.x(),c2.x(),c3.x(),c4.x(),u,v); double Vp = TRAN_QUA(p12.y(), p23.y(), p34.y(), p41.y(), c1.y(),c2.y(),c3.y(),c4.y(),u,v); - fprintf(f3,"SP(%g,%g,%g) {%d};\n", Up, Vp, 0.0, 1); + if(f3) fprintf(f3,"SP(%g,%g,%g) {%d};\n", Up, Vp, 0.0, 1); // printf("v1=%d v2=%d v3=%d v4=%d \n", v1->getNum(), v2->getNum(), v3->getNum(), v4->getNum()); // printf("c1=%g %g, c2=%g %g, c3=%g %g, c4=%g,%g \n", c1.x(),c1.y(),c2.x(),c2.y(),c3.x(),c3.y(),c4.x(),c4.y()); @@ -1069,7 +1069,7 @@ void createRegularMesh (GFace *gf, (p34.x() && p34.y() == -1.0) || (p41.x() && p41.y() == -1.0)) { Msg::Error("Wrong param -1"); - fclose(f3); + if(f3) fclose(f3); return; } @@ -1085,9 +1085,10 @@ void createRegularMesh (GFace *gf, } } - fprintf(f3,"};\n"); - fclose(f3); - + if(f3){ + fprintf(f3,"};\n"); + fclose(f3); + } } void updateQuadCavity (GFace *gf, @@ -1351,18 +1352,20 @@ struct quadBlob { char name[234]; sprintf(name,"blob%d_%d_%d.pos", iBlob, iter, method); FILE *f = Fopen(name,"w"); - fprintf(f,"View \"%s\" {\n",name); - for (unsigned int i = 0; i < quads.size(); i++){ - quads[i]->writePOS(f,true,false,false,false,false,false); - } - for (unsigned int i = 0; i < bnodes.size(); i++){ - fprintf(f,"SP(%g,%g,%g) {%d};\n", - bnodes[i]->x(), - bnodes[i]->y(), - bnodes[i]->z(),topologicalAngle(bnodes[i])); + if(f){ + fprintf(f,"View \"%s\" {\n",name); + for (unsigned int i = 0; i < quads.size(); i++){ + quads[i]->writePOS(f,true,false,false,false,false,false); + } + for (unsigned int i = 0; i < bnodes.size(); i++){ + fprintf(f,"SP(%g,%g,%g) {%d};\n", + bnodes[i]->x(), + bnodes[i]->y(), + bnodes[i]->z(),topologicalAngle(bnodes[i])); + } + fprintf(f,"};\n"); + fclose(f); } - fprintf(f,"};\n"); - fclose(f); } void smooth (int); bool meshable (int iter) @@ -2026,15 +2029,18 @@ struct opti_data_vertex_relocation { } } - void print_cavity (char *name){ + void print_cavity (char *name) + { FILE *f = Fopen(name,"w"); - fprintf(f,"View \"\"{\n"); - for (unsigned int i=0;i<e.size();++i){ - MElement *el = e[i]; - el->writePOS(f,false,false,false,true,false,false); + if(f){ + fprintf(f,"View \"\"{\n"); + for (unsigned int i=0;i<e.size();++i){ + MElement *el = e[i]; + el->writePOS(f,false,false,false,true,false,false); + } + fprintf(f,"};"); + fclose (f); } - fprintf(f,"};"); - fclose (f); } /// quality measure for a quad @@ -3641,7 +3647,7 @@ void recombineIntoQuads(GFace *gf, // optistatus[0] = (ITERB == 1) ?splitFlatQuads(gf, .01, prioritory) : 0; //optistatus[1] = removeTwoQuadsNodes(gf); - //optistatus[4] = + //optistatus[4] = // _defectsRemovalBunin(gf,maxCavitySize); //optistatus[2] = removeDiamonds(gf) ; diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp index fec89fa4b35c69f0deefe04b70e40d494a2a555f..447d0c6aac64b6c49aab0944b5399b43794884f0 100644 --- a/Mesh/meshGRegion.cpp +++ b/Mesh/meshGRegion.cpp @@ -162,10 +162,9 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates) } //print voronoi poles - FILE *outfile; //smooth_normals *snorm = gr->model()->normals; - outfile = Fopen("nodes.pos", "w"); - fprintf(outfile, "View \"Voronoi poles\" {\n"); + FILE *outfile = Fopen("nodes.pos", "w"); + if(outfile) fprintf(outfile, "View \"Voronoi poles\" {\n"); std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.begin(); for(; itmap != node2Tet.end(); itmap++){ //MVertex *v = itmap->first; @@ -176,35 +175,36 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates) for(; it != allTets.end(); it++){ double radius = (*it)->getCircumRadius(); if (radius > maxRadius){ - maxRadius = radius; - poleTet = *it; + maxRadius = radius; + poleTet = *it; } } //if (v->onWhat()->dim() == 2 ){ - SPoint3 pc = poleTet->circumcenter(); - //double nx,ny,nz; - // SVector3 vN = snorm->get(v->x(), v->y(), v->z(), nx,ny,nz); - // SVector3 pcv(pc.x()-nx, pc.y()-ny,pc.z()-nz); - // printf("nx=%g ny=%g nz=%g dot=%g \n", nx,ny,nz, dot(vN, pcv)); - // if ( dot(vN, pcv) > 0.0 ) - double x[3] = {pc.x(), pc.y(), pc.z()}; - double uvw[3]; - poleTet->xyz2uvw(x, uvw); - //bool inside = poleTet->isInside(uvw[0], uvw[1], uvw[2]); - //if (inside){ - fprintf(outfile,"SP(%g,%g,%g) {%g};\n", - pc.x(), pc.y(), pc.z(), maxRadius); - candidates.insert(pc); - //} + SPoint3 pc = poleTet->circumcenter(); + //double nx,ny,nz; + // SVector3 vN = snorm->get(v->x(), v->y(), v->z(), nx,ny,nz); + // SVector3 pcv(pc.x()-nx, pc.y()-ny,pc.z()-nz); + // printf("nx=%g ny=%g nz=%g dot=%g \n", nx,ny,nz, dot(vN, pcv)); + // if ( dot(vN, pcv) > 0.0 ) + double x[3] = {pc.x(), pc.y(), pc.z()}; + double uvw[3]; + poleTet->xyz2uvw(x, uvw); + //bool inside = poleTet->isInside(uvw[0], uvw[1], uvw[2]); + //if (inside){ + if(outfile) + fprintf(outfile,"SP(%g,%g,%g) {%g};\n", pc.x(), pc.y(), pc.z(), maxRadius); + candidates.insert(pc); //} + //} + } + if(outfile){ + fprintf(outfile,"};\n"); + fclose(outfile); } - fprintf(outfile,"};\n"); - fclose(outfile); - //print scalar lines - FILE *outfile2; - outfile2 = Fopen("edges.pos", "w"); - fprintf(outfile2, "View \"Voronoi edges\" {\n"); + // print scalar lines + FILE *outfile2 = Fopen("edges.pos", "w"); + if(outfile2) fprintf(outfile2, "View \"Voronoi edges\" {\n"); std::map<MFace, std::vector<MTetrahedron*> , Less_Face >::iterator itmap2 = face2Tet.begin(); for(; itmap2 != face2Tet.end(); itmap2++){ std::vector<MTetrahedron*> allTets = itmap2->second; @@ -214,13 +214,15 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates) //std::set<SPoint3>::const_iterator it1 = candidates.find(pc1); //std::set<SPoint3>::const_iterator it2 = candidates.find(pc2); //if( it1 != candidates.end() || it2 != candidates.end()) + if(outfile2) fprintf(outfile2,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n", - pc1.x(), pc1.y(), pc1.z(), pc2.x(), pc2.y(), pc2.z(), - allTets[0]->getCircumRadius(),allTets[1]->getCircumRadius()); - } - fprintf(outfile2,"};\n"); - fclose(outfile2); - + pc1.x(), pc1.y(), pc1.z(), pc2.x(), pc2.y(), pc2.z(), + allTets[0]->getCircumRadius(),allTets[1]->getCircumRadius()); + } + if(outfile2){ + fprintf(outfile2,"};\n"); + fclose(outfile2); + } } @@ -978,15 +980,17 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr, spl // ------------------------------------------------------------------------------------ { FILE *ff2 = fopen ("tato3D.pos","w"); - fprintf(ff2,"View \" \"{\n"); - for (unsigned int i = 0; i < blPrisms.size();i++){ - blPrisms[i]->writePOS(ff2,1,0,0,0,0,0); - } - for (unsigned int i = 0; i < blHexes.size();i++){ - blHexes[i]->writePOS(ff2,1,0,0,0,0,0); + if(ff2){ + fprintf(ff2,"View \" \"{\n"); + for (unsigned int i = 0; i < blPrisms.size();i++){ + blPrisms[i]->writePOS(ff2,1,0,0,0,0,0); + } + for (unsigned int i = 0; i < blHexes.size();i++){ + blHexes[i]->writePOS(ff2,1,0,0,0,0,0); + } + fprintf(ff2,"};\n"); + fclose(ff2); } - fprintf(ff2,"};\n"); - fclose(ff2); } for (unsigned int i = 0; i < blPrisms.size();i++){ @@ -1008,14 +1012,15 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr, spl std::map<MFace,MElement*,Less_Face>::iterator it = bfaces.begin(); FILE *ff = fopen ("toto3D.pos","w"); - fprintf(ff,"View \" \"{\n"); + if(ff) fprintf(ff,"View \" \"{\n"); for (; it != bfaces.end(); ++it){ if (it->first.getNumVertices() == 3){ nf->triangles.push_back(new MTriangle (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2))); - fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n", - it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), - it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), - it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z()); + if(ff) + fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n", + it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), + it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), + it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z()); } else { @@ -1061,35 +1066,38 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr, spl nf->triangles.push_back(new MTriangle (it->first.getVertex(1),it->first.getVertex(2),newv)); nf->triangles.push_back(new MTriangle (it->first.getVertex(2),it->first.getVertex(3),newv)); nf->triangles.push_back(new MTriangle (it->first.getVertex(3),it->first.getVertex(0),newv)); - - fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", - it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), - it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), - newv->x(),newv->y(),newv->z()); - - fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", - it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), - it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(), - newv->x(),newv->y(),newv->z()); - - fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", - it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(), - it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(), - newv->x(),newv->y(),newv->z()); - fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", - it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(), - it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), - newv->x(),newv->y(),newv->z()); - - fprintf(ff,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){3,3,3,3};\n", - it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), - it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), - it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(), - it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z()); + if(ff){ + fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", + it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), + it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), + newv->x(),newv->y(),newv->z()); + + fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", + it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), + it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(), + newv->x(),newv->y(),newv->z()); + + fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", + it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(), + it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(), + newv->x(),newv->y(),newv->z()); + fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n", + it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(), + it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), + newv->x(),newv->y(),newv->z()); + + fprintf(ff,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){3,3,3,3};\n", + it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(), + it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(), + it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(), + it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z()); + } } } - fprintf(ff,"};\n"); - fclose(ff); + if(ff){ + fprintf(ff,"};\n"); + fclose(ff); + } printf("discrete face with %d %d elems\n", (int)nf->triangles.size(), (int)nf->quadrangles.size()); @@ -1652,12 +1660,14 @@ void meshNormalsPointOutOfTheRegion(GRegion *gr) } // FILE *fp = Fopen("debug.pos", "w"); - // fprintf(fp, "View \"debug\" {\n"); - // for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++) - // for(unsigned int i = 0; i < (*it)->triangles.size(); i++) - // (*it)->triangles[i]->writePOS(fp, 1., (*it)->tag()); - // fprintf(fp, "};\n"); - // fclose(fp); + // if(fp){ + // fprintf(fp, "View \"debug\" {\n"); + // for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++) + // for(unsigned int i = 0; i < (*it)->triangles.size(); i++) + // (*it)->triangles[i]->writePOS(fp, 1., (*it)->tag()); + // fprintf(fp, "};\n"); + // fclose(fp); + // } } void meshGRegion::operator() (GRegion *gr) diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp index 9c210689eeca626f29df3b375d7d1bed56d19da3..18b57cd21ae4db07594f6d8242ce48caae14d1c7 100644 --- a/Mesh/meshGRegionDelaunayInsertion.cpp +++ b/Mesh/meshGRegionDelaunayInsertion.cpp @@ -408,19 +408,21 @@ void recurFindCavity(std::vector<faceXtet> & shell, void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false ) { FILE *f = Fopen (fn,"w"); - fprintf(f,"View \"\"{\n"); - std::list<MTet4*>::iterator ittet = cavity.begin(); - std::list<MTet4*>::iterator ittete = cavity.end(); - while (ittet != ittete){ - MTet4 *tet = *ittet; - if (force || !tet->isDeleted()){ - MTetrahedron *t = tet->tet(); - t->writePOS (f, false,false,false,true,false,false); + if(f){ + fprintf(f,"View \"\"{\n"); + std::list<MTet4*>::iterator ittet = cavity.begin(); + std::list<MTet4*>::iterator ittete = cavity.end(); + while (ittet != ittete){ + MTet4 *tet = *ittet; + if (force || !tet->isDeleted()){ + MTetrahedron *t = tet->tet(); + t->writePOS (f, false,false,false,true,false,false); + } + ittet++; } - ittet++; + fprintf(f,"};\n"); + fclose(f); } - fprintf(f,"};\n"); - fclose(f); } diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp index dd3b2a99f9e4a9c252bb39b6f6965c5e0cef4f9a..9973dd9f412b4d02fc471a7d96a775d0ba29a6c1 100644 --- a/Mesh/multiscalePartition.cpp +++ b/Mesh/multiscalePartition.cpp @@ -275,27 +275,29 @@ static void printLevel(std::vector<MElement *> &elements, int recur, int region) bool binary = false; FILE *fp = Fopen (fn, "w"); - fprintf(fp, "$MeshFormat\n"); - fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double)); - fprintf(fp, "$EndMeshFormat\n"); - - fprintf(fp,"$Nodes\n%d\n", (int)vs.size()); - std::set<MVertex*> :: iterator it = vs.begin(); - int index = 1; - for (; it != vs.end() ; ++it){ - (*it)->setIndex(index++); - fprintf(fp,"%d %g %g %g\n",(*it)->getIndex(), - (*it)->x(),(*it)->y(),(*it)->z()); - } - fprintf(fp,"$EndNodes\n"); + if(fp){ + fprintf(fp, "$MeshFormat\n"); + fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double)); + fprintf(fp, "$EndMeshFormat\n"); + + fprintf(fp,"$Nodes\n%d\n", (int)vs.size()); + std::set<MVertex*> :: iterator it = vs.begin(); + int index = 1; + for (; it != vs.end() ; ++it){ + (*it)->setIndex(index++); + fprintf(fp,"%d %g %g %g\n",(*it)->getIndex(), + (*it)->x(),(*it)->y(),(*it)->z()); + } + fprintf(fp,"$EndNodes\n"); - fprintf(fp,"$Elements\n%d\n", (int)elements.size()); - for (unsigned int i = 0; i < elements.size(); i++){ - elements[i]->writeMSH(fp, version); - } - fprintf(fp,"$EndElements\n%d\n", (int)elements.size()); + fprintf(fp,"$Elements\n%d\n", (int)elements.size()); + for (unsigned int i = 0; i < elements.size(); i++){ + elements[i]->writeMSH(fp, version); + } + fprintf(fp,"$EndElements\n%d\n", (int)elements.size()); - fclose(fp); + fclose(fp); + } } */ diff --git a/Mesh/pointInsertion.cpp b/Mesh/pointInsertion.cpp index b6945fb3bb78b8358ff28be2ccd17575692b9d84..c548c5874500bcc85076e5bdf42afb167069f5e8 100644 --- a/Mesh/pointInsertion.cpp +++ b/Mesh/pointInsertion.cpp @@ -412,19 +412,21 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, // surface mesh char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag()); FILE *f = Fopen(ccc,"w"); - fprintf(f,"View \"\"{\n"); - for (unsigned int i=0;i<vertices.size();i++){ - vertices[i]->print(f,i); - if(vertices[i]->_v->onWhat() == gf) { - packed.push_back(vertices[i]->_v); - metrics.push_back(vertices[i]->_meshMetric); - SPoint2 midpoint; - reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint); + if(f){ + fprintf(f,"View \"\"{\n"); + for (unsigned int i=0;i<vertices.size();i++){ + vertices[i]->print(f,i); + if(vertices[i]->_v->onWhat() == gf) { + packed.push_back(vertices[i]->_v); + metrics.push_back(vertices[i]->_meshMetric); + SPoint2 midpoint; + reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint); + } + delete vertices[i]; } - delete vertices[i]; + fprintf(f,"};"); + fclose(f); } - fprintf(f,"};"); - fclose(f); } bool Filler3D::treat_region(GRegion *gr) diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp index aae6d520d04f11b9f695069d13d6c4e0f42d93ee..e205176e72be10e2360e5807538cbfb0452d9059 100644 --- a/Mesh/surfaceFiller.cpp +++ b/Mesh/surfaceFiller.cpp @@ -801,11 +801,11 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed std::set<MVertex*>::iterator it = bnd_vertices.begin() ; char NAME[345]; sprintf(NAME,"crossReal%d.pos",gf->tag()); - FILE *crossf=NULL; + FILE *crossf = NULL; if (debug){ crossf = Fopen (NAME,"w"); } - if (crossf)fprintf(crossf,"View \"\"{\n"); + if (crossf) fprintf(crossf,"View \"\"{\n"); for (; it != bnd_vertices.end() ; ++it){ SPoint2 midpoint; //compute4neighbors_RK2 (gf, *it, midpoint, goNonLinear, newp, metricField,crossf); @@ -887,9 +887,9 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed // surface mesh char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag()); FILE *f = Fopen(ccc,"w"); - fprintf(f,"View \"\"{\n"); + if(f) fprintf(f, "View \"\"{\n"); for (unsigned int i=0;i<vertices.size();i++){ - vertices[i]->print(f,i); + if(f) vertices[i]->print(f,i); if(vertices[i]->_v->onWhat() == gf) { packed.push_back(vertices[i]->_v); metrics.push_back(vertices[i]->_meshMetric); @@ -898,8 +898,10 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed } delete vertices[i]; } - fprintf(f,"};"); - fclose(f); + if(f){ + fprintf(f,"};"); + fclose(f); + } } @@ -946,7 +948,7 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec char NAME[345]; sprintf(NAME,"crossReal%d.pos",gf->tag()); FILE *crossf = Fopen (NAME,"w"); - if (crossf)fprintf(crossf,"View \"\"{\n"); + if (crossf) fprintf(crossf,"View \"\"{\n"); for (; it != bnd_vertices.end() ; ++it){ SPoint2 midpoint; compute4neighbors (gf, *it, midpoint, goNonLinear, newp, metricField,crossf); @@ -1009,10 +1011,10 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec // surface mesh char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag()); FILE *f = Fopen(ccc,"w"); - fprintf(f,"View \"\"{\n"); + if(f) fprintf(f,"View \"\"{\n"); for (unsigned int i=0;i<vertices.size();i++){ // if(vertices[i]->_v->onWhat() != gf) - vertices[i]->print(f,i); + if(f) vertices[i]->print(f,i); if(vertices[i]->_v->onWhat() == gf) { packed.push_back(vertices[i]->_v); metrics.push_back(vertices[i]->_meshMetric); @@ -1026,8 +1028,10 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec } delete vertices[i]; } - fprintf(f,"};"); - fclose(f); + if(f){ + fprintf(f,"};"); + fclose(f); + } // printf("packed.size = %d\n",packed.size()); // delete rtree; }