1 #ifdef OCCGEOMETRY 2 3 #include <mystdlib.h> 4 #include <occgeom.hpp> 5 #include <meshing.hpp> 6 7 8 namespace netgen 9 { 10 11 #include "occmeshsurf.hpp" 12 13 #define TCL_OK 0 14 #define TCL_ERROR 1 15 16 #define DIVIDEEDGESECTIONS 1000 17 #define IGNORECURVELENGTH 1e-4 18 #define VSMALL 1e-10 19 20 21 bool merge_solids = 1; 22 23 24 // can you please explain what you intend to compute here (JS) !!! Dist(Line l)25 double Line :: Dist (Line l) 26 { 27 Vec<3> n = p1-p0; 28 Vec<3> q = l.p1-l.p0; 29 double nq = n*q; 30 31 Point<3> p = p0 + 0.5*n; 32 double lambda = (p-l.p0)*n / (nq + VSMALL); 33 34 if (lambda >= 0 && lambda <= 1) 35 { 36 double d = (p-l.p0-lambda*q).Length(); 37 // if (d < 1e-3) d = 1e99; 38 return d; 39 } 40 else 41 return 1e99; 42 } 43 44 45 Length()46 double Line :: Length () 47 { 48 return (p1-p0).Length(); 49 } 50 51 52 occ2ng(const gp_Pnt & p)53 inline Point<3> occ2ng (const gp_Pnt & p) 54 { 55 return Point<3> (p.X(), p.Y(), p.Z()); 56 } 57 58 59 ComputeH(double kappa)60 double ComputeH (double kappa) 61 { 62 double hret; 63 kappa *= mparam.curvaturesafety; 64 65 if (mparam.maxh * kappa < 1) 66 hret = mparam.maxh; 67 else 68 hret = 1 / (kappa + VSMALL); 69 70 if (mparam.maxh < hret) 71 hret = mparam.maxh; 72 73 return (hret); 74 } 75 76 77 78 RestrictHTriangle(gp_Pnt2d & par0,gp_Pnt2d & par1,gp_Pnt2d & par2,BRepLProp_SLProps * prop,Mesh & mesh,int depth,double h=0)79 void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, 80 BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0) 81 { 82 int ls = -1; 83 84 gp_Pnt pnt0,pnt1,pnt2; 85 86 prop->SetParameters (par0.X(), par0.Y()); 87 pnt0 = prop->Value(); 88 89 prop->SetParameters (par1.X(), par1.Y()); 90 pnt1 = prop->Value(); 91 92 prop->SetParameters (par2.X(), par2.Y()); 93 pnt2 = prop->Value(); 94 95 double aux; 96 double maxside = pnt0.Distance(pnt1); 97 ls = 2; 98 aux = pnt1.Distance(pnt2); 99 if(aux > maxside) 100 { 101 maxside = aux; 102 ls = 0; 103 } 104 aux = pnt2.Distance(pnt0); 105 if(aux > maxside) 106 { 107 maxside = aux; 108 ls = 1; 109 } 110 111 112 113 gp_Pnt2d parmid; 114 115 parmid.SetX( (par0.X()+par1.X()+par2.X()) / 3 ); 116 parmid.SetY( (par0.Y()+par1.Y()+par2.Y()) / 3 ); 117 118 if (depth%3 == 0) 119 { 120 double curvature = 0; 121 122 prop->SetParameters (parmid.X(), parmid.Y()); 123 if (!prop->IsCurvatureDefined()) 124 { 125 (*testout) << "curvature not defined!" << endl; 126 return; 127 } 128 curvature = max(fabs(prop->MinCurvature()), 129 fabs(prop->MaxCurvature())); 130 131 prop->SetParameters (par0.X(), par0.Y()); 132 if (!prop->IsCurvatureDefined()) 133 { 134 (*testout) << "curvature not defined!" << endl; 135 return; 136 } 137 curvature = max(curvature,max(fabs(prop->MinCurvature()), 138 fabs(prop->MaxCurvature()))); 139 140 prop->SetParameters (par1.X(), par1.Y()); 141 if (!prop->IsCurvatureDefined()) 142 { 143 (*testout) << "curvature not defined!" << endl; 144 return; 145 } 146 curvature = max(curvature,max(fabs(prop->MinCurvature()), 147 fabs(prop->MaxCurvature()))); 148 149 prop->SetParameters (par2.X(), par2.Y()); 150 if (!prop->IsCurvatureDefined()) 151 { 152 (*testout) << "curvature not defined!" << endl; 153 return; 154 } 155 curvature = max(curvature,max(fabs(prop->MinCurvature()), 156 fabs(prop->MaxCurvature()))); 157 158 //(*testout) << "curvature " << curvature << endl; 159 160 if (curvature < 1e-3) 161 { 162 //(*testout) << "curvature too small (" << curvature << ")!" << endl; 163 return; 164 // return war bis 10.2.05 auskommentiert 165 } 166 167 168 169 h = ComputeH (curvature+1e-10); 170 171 if(h < 1e-4*maxside) 172 return; 173 174 175 if (h > 30) return; 176 } 177 178 if (h < maxside && depth < 10) 179 { 180 //cout << "\r h " << h << flush; 181 gp_Pnt2d pm; 182 183 //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; 184 //cout << "par0 " << par0.X() << " " << par0.Y() 185 //<< " par1 " << par1.X() << " " << par1.Y() 186 // << " par2 " << par2.X() << " " << par2.Y()<< endl; 187 188 if(ls == 0) 189 { 190 pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); 191 RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); 192 RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); 193 } 194 else if(ls == 1) 195 { 196 pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); 197 RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); 198 RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); 199 } 200 else if(ls == 2) 201 { 202 pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); 203 RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); 204 RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); 205 } 206 207 } 208 else 209 { 210 gp_Pnt pnt; 211 Point3d p3d; 212 213 prop->SetParameters (parmid.X(), parmid.Y()); 214 pnt = prop->Value(); 215 p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); 216 mesh.RestrictLocalH (p3d, h); 217 218 p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z()); 219 mesh.RestrictLocalH (p3d, h); 220 221 p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z()); 222 mesh.RestrictLocalH (p3d, h); 223 224 p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z()); 225 mesh.RestrictLocalH (p3d, h); 226 227 //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl; 228 229 } 230 } 231 232 233 DivideEdge(TopoDS_Edge & edge,Array<MeshPoint> & ps,Array<double> & params,Mesh & mesh)234 void DivideEdge (TopoDS_Edge & edge, Array<MeshPoint> & ps, 235 Array<double> & params, Mesh & mesh) 236 { 237 double s0, s1; 238 double maxh = mparam.maxh; 239 int nsubedges = 1; 240 gp_Pnt pnt, oldpnt; 241 double svalue[DIVIDEEDGESECTIONS]; 242 243 GProp_GProps system; 244 BRepGProp::LinearProperties(edge, system); 245 double L = system.Mass(); 246 247 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); 248 249 double hvalue[DIVIDEEDGESECTIONS+1]; 250 hvalue[0] = 0; 251 pnt = c->Value(s0); 252 253 double olddist = 0; 254 double dist = 0; 255 256 int tmpVal = (int)(DIVIDEEDGESECTIONS); 257 258 for (int i = 1; i <= tmpVal; i++) 259 { 260 oldpnt = pnt; 261 pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); 262 hvalue[i] = hvalue[i-1] + 263 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* 264 pnt.Distance(oldpnt); 265 266 //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) 267 // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; 268 269 olddist = dist; 270 dist = pnt.Distance(oldpnt); 271 } 272 273 // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS])); 274 nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5))); 275 276 ps.SetSize(nsubedges-1); 277 params.SetSize(nsubedges+1); 278 279 int i = 1; 280 int i1 = 0; 281 do 282 { 283 if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i) 284 { 285 params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0); 286 pnt = c->Value(params[i]); 287 ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z())); 288 i++; 289 } 290 i1++; 291 if (i1 > DIVIDEEDGESECTIONS) 292 { 293 nsubedges = i; 294 ps.SetSize(nsubedges-1); 295 params.SetSize(nsubedges+1); 296 cout << "divide edge: local h too small" << endl; 297 } 298 } while (i < nsubedges); 299 300 params[0] = s0; 301 params[nsubedges] = s1; 302 303 if (params[nsubedges] <= params[nsubedges-1]) 304 { 305 cout << "CORRECTED" << endl; 306 ps.SetSize (nsubedges-2); 307 params.SetSize (nsubedges); 308 params[nsubedges] = s1; 309 } 310 } 311 312 313 314 OCCFindEdges(OCCGeometry & geom,Mesh & mesh)315 void OCCFindEdges (OCCGeometry & geom, Mesh & mesh) 316 { 317 const char * savetask = multithread.task; 318 multithread.task = "Edge meshing"; 319 320 (*testout) << "edge meshing" << endl; 321 322 int nvertices = geom.vmap.Extent(); 323 int nedges = geom.emap.Extent(); 324 325 (*testout) << "nvertices = " << nvertices << endl; 326 (*testout) << "nedges = " << nedges << endl; 327 328 double eps = 1e-6 * geom.GetBoundingBox().Diam(); 329 330 for (int i = 1; i <= nvertices; i++) 331 { 332 gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i))); 333 MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) ); 334 335 bool exists = 0; 336 if (merge_solids) 337 for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++) 338 if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps) 339 { 340 exists = 1; 341 break; 342 } 343 344 if (!exists) 345 mesh.AddPoint (mp); 346 } 347 348 (*testout) << "different vertices = " << mesh.GetNP() << endl; 349 350 351 int first_ep = mesh.GetNP()+1; 352 353 Array<int> face2solid[2]; 354 for (int i = 0; i<2; i++) 355 { 356 face2solid[i].SetSize (geom.fmap.Extent()); 357 face2solid[i] = 0; 358 } 359 360 int solidnr = 0; 361 for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next()) 362 { 363 solidnr++; 364 for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next()) 365 { 366 TopoDS_Face face = TopoDS::Face(exp1.Current()); 367 int facenr = geom.fmap.FindIndex(face); 368 369 if (face2solid[0][facenr-1] == 0) 370 face2solid[0][facenr-1] = solidnr; 371 else 372 face2solid[1][facenr-1] = solidnr; 373 } 374 } 375 376 377 int total = 0; 378 for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) 379 for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next()) 380 for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next()) 381 total++; 382 383 384 int facenr = 0; 385 int edgenr = 0; 386 387 388 (*testout) << "faces = " << geom.fmap.Extent() << endl; 389 int curr = 0; 390 391 for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) 392 { 393 TopoDS_Face face = TopoDS::Face(geom.fmap(i3)); 394 facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS 395 396 int solidnr0 = face2solid[0][i3-1]; 397 int solidnr1 = face2solid[1][i3-1]; 398 399 /* auskommentiert am 3.3.05 von robert 400 for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next()) 401 { 402 TopoDS_Face face2 = TopoDS::Face(exp2.Current()); 403 if (geom.fmap.FindIndex(face2) == facenr) 404 { 405 // if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1); 406 } 407 } 408 */ 409 410 mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); 411 412 // Philippose - 06/07/2009 413 // Add the face colour to the mesh data 414 Quantity_Color face_colour; 415 416 if(!(geom.face_colours.IsNull()) 417 && (geom.face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour))) 418 { 419 mesh.GetFaceDescriptor(facenr).SetSurfColour(Vec3d(face_colour.Red(),face_colour.Green(),face_colour.Blue())); 420 } 421 else 422 { 423 mesh.GetFaceDescriptor(facenr).SetSurfColour(Vec3d(0.0,1.0,0.0)); 424 } 425 // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG) 426 427 428 Handle(Geom_Surface) occface = BRep_Tool::Surface(face); 429 430 for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next()) 431 { 432 TopoDS_Shape wire = exp2.Current(); 433 434 for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next()) 435 { 436 curr++; 437 (*testout) << "edge nr " << curr << endl; 438 439 multithread.percent = 100 * curr / double (total); 440 if (multithread.terminate) return; 441 442 TopoDS_Edge edge = TopoDS::Edge (exp3.Current()); 443 if (BRep_Tool::Degenerated(edge)) 444 { 445 //(*testout) << "ignoring degenerated edge" << endl; 446 continue; 447 } 448 449 if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == 450 geom.vmap.FindIndex(TopExp::LastVertex (edge))) 451 { 452 GProp_GProps system; 453 BRepGProp::LinearProperties(edge, system); 454 455 if (system.Mass() < eps) 456 { 457 cout << "ignoring edge " << geom.emap.FindIndex (edge) 458 << ". closed edge with length < " << eps << endl; 459 continue; 460 } 461 } 462 463 464 Handle(Geom2d_Curve) cof; 465 double s0, s1; 466 cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); 467 468 int geomedgenr = geom.emap.FindIndex(edge); 469 470 Array <MeshPoint> mp; 471 Array <double> params; 472 473 DivideEdge (edge, mp, params, mesh); 474 475 Array <int> pnums; 476 pnums.SetSize (mp.Size()+2); 477 478 if (!merge_solids) 479 { 480 pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)); 481 pnums[pnums.Size()-1] = geom.vmap.FindIndex (TopExp::LastVertex (edge)); 482 } 483 else 484 { 485 Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); 486 Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); 487 488 pnums[0] = -1; 489 pnums.Last() = -1; 490 for (PointIndex pi = 1; pi < first_ep; pi++) 491 { 492 if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi; 493 if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; 494 } 495 } 496 497 498 for (int i = 1; i <= mp.Size(); i++) 499 { 500 bool exists = 0; 501 int j; 502 for (j = first_ep; j <= mesh.GetNP(); j++) 503 if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) 504 { 505 exists = 1; 506 break; 507 } 508 509 if (exists) 510 pnums[i] = j; 511 else 512 { 513 mesh.AddPoint (mp[i-1]); 514 (*testout) << "add meshpoint " << mp[i-1] << endl; 515 pnums[i] = mesh.GetNP(); 516 } 517 } 518 (*testout) << "NP = " << mesh.GetNP() << endl; 519 520 //(*testout) << pnums[pnums.Size()-1] << endl; 521 522 for (int i = 1; i <= mp.Size()+1; i++) 523 { 524 edgenr++; 525 Segment seg; 526 527 seg[0] = pnums[i-1]; 528 seg[1] = pnums[i]; 529 seg.edgenr = edgenr; 530 seg.si = facenr; 531 seg.epgeominfo[0].dist = params[i-1]; 532 seg.epgeominfo[1].dist = params[i]; 533 seg.epgeominfo[0].edgenr = geomedgenr; 534 seg.epgeominfo[1].edgenr = geomedgenr; 535 536 gp_Pnt2d p2d; 537 p2d = cof->Value(params[i-1]); 538 // if (i == 1) p2d = cof->Value(s0); 539 seg.epgeominfo[0].u = p2d.X(); 540 seg.epgeominfo[0].v = p2d.Y(); 541 p2d = cof->Value(params[i]); 542 // if (i == mp.Size()+1) p2d = cof -> Value(s1); 543 seg.epgeominfo[1].u = p2d.X(); 544 seg.epgeominfo[1].v = p2d.Y(); 545 546 /* 547 if (occface->IsUPeriodic()) 548 { 549 cout << "U Periodic" << endl; 550 if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > 551 fabs(seg.epgeominfo[1].u- 552 (seg.epgeominfo[0].u-occface->UPeriod()))) 553 seg.epgeominfo[0].u = p2d.X()+occface->UPeriod(); 554 555 if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > 556 fabs(seg.epgeominfo[1].u- 557 (seg.epgeominfo[0].u+occface->UPeriod()))) 558 seg.epgeominfo[0].u = p2d.X()-occface->UPeriod(); 559 } 560 561 if (occface->IsVPeriodic()) 562 { 563 cout << "V Periodic" << endl; 564 if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > 565 fabs(seg.epgeominfo[1].v- 566 (seg.epgeominfo[0].v-occface->VPeriod()))) 567 seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod(); 568 569 if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > 570 fabs(seg.epgeominfo[1].v- 571 (seg.epgeominfo[0].v+occface->VPeriod()))) 572 seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod(); 573 } 574 */ 575 576 if (edge.Orientation() == TopAbs_REVERSED) 577 { 578 swap (seg[0], seg[1]); 579 swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist); 580 swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); 581 swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); 582 } 583 584 mesh.AddSegment (seg); 585 586 //edgesegments[geomedgenr-1]->Append(mesh.GetNSeg()); 587 588 } 589 } 590 } 591 } 592 593 // for(i=1; i<=mesh.GetNSeg(); i++) 594 // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si 595 // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; 596 // exit(10); 597 598 mesh.CalcSurfacesOfNode(); 599 multithread.task = savetask; 600 } 601 602 603 604 OCCMeshSurface(OCCGeometry & geom,Mesh & mesh,int perfstepsend)605 void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend) 606 { 607 int i, j, k; 608 int changed; 609 610 const char * savetask = multithread.task; 611 multithread.task = "Surface meshing"; 612 613 geom.facemeshstatus = 0; 614 615 int noldp = mesh.GetNP(); 616 617 double starttime = GetTime(); 618 619 Array<int> glob2loc(noldp); 620 621 //int projecttype = PARAMETERSPACE; 622 623 int projecttype = PARAMETERSPACE; 624 625 int notrys = 1; 626 627 int surfmesherror = 0; 628 629 for (k = 1; k <= mesh.GetNFD(); k++) 630 { 631 if(1==0 && !geom.fvispar[k-1].IsDrawable()) 632 { 633 (*testout) << "ignoring face " << k << endl; 634 cout << "ignoring face " << k << endl; 635 continue; 636 } 637 638 (*testout) << "mesh face " << k << endl; 639 multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); 640 geom.facemeshstatus[k-1] = -1; 641 642 643 /* 644 if (k != 42) 645 { 646 cout << "skipped" << endl; 647 continue; 648 } 649 */ 650 651 652 FaceDescriptor & fd = mesh.GetFaceDescriptor(k); 653 654 int oldnf = mesh.GetNSE(); 655 656 Box<3> bb = geom.GetBoundingBox(); 657 658 // int projecttype = PLANESPACE; 659 660 Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype); 661 662 if (meshing.GetProjectionType() == PLANESPACE) 663 PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); 664 else 665 PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)"); 666 667 if (surfmesherror) 668 cout << "Surface meshing error occured before (in " << surfmesherror << " faces)" << endl; 669 670 // Meshing2OCCSurfaces meshing(f2, bb); 671 meshing.SetStartTime (starttime); 672 673 //(*testout) << "Face " << k << endl << endl; 674 675 676 if (meshing.GetProjectionType() == PLANESPACE) 677 { 678 int cntp = 0; 679 glob2loc = 0; 680 for (i = 1; i <= mesh.GetNSeg(); i++) 681 { 682 Segment & seg = mesh.LineSegment(i); 683 if (seg.si == k) 684 { 685 for (j = 1; j <= 2; j++) 686 { 687 int pi = (j == 1) ? seg[0] : seg[1]; 688 if (!glob2loc.Get(pi)) 689 { 690 meshing.AddPoint (mesh.Point(pi), pi); 691 cntp++; 692 glob2loc.Elem(pi) = cntp; 693 } 694 } 695 } 696 } 697 698 for (i = 1; i <= mesh.GetNSeg(); i++) 699 { 700 Segment & seg = mesh.LineSegment(i); 701 if (seg.si == k) 702 { 703 PointGeomInfo gi0, gi1; 704 gi0.trignum = gi1.trignum = k; 705 gi0.u = seg.epgeominfo[0].u; 706 gi0.v = seg.epgeominfo[0].v; 707 gi1.u = seg.epgeominfo[1].u; 708 gi1.v = seg.epgeominfo[1].v; 709 710 meshing.AddBoundaryElement (glob2loc.Get(seg[0]), glob2loc.Get(seg[1]), gi0, gi1); 711 //(*testout) << gi0.u << " " << gi0.v << endl; 712 //(*testout) << gi1.u << " " << gi1.v << endl; 713 } 714 } 715 } 716 else 717 { 718 int cntp = 0; 719 720 for (i = 1; i <= mesh.GetNSeg(); i++) 721 if (mesh.LineSegment(i).si == k) 722 cntp+=2; 723 724 725 Array< PointGeomInfo > gis; 726 727 gis.SetAllocSize (cntp); 728 gis.SetSize (0); 729 730 for (i = 1; i <= mesh.GetNSeg(); i++) 731 { 732 Segment & seg = mesh.LineSegment(i); 733 if (seg.si == k) 734 { 735 PointGeomInfo gi0, gi1; 736 gi0.trignum = gi1.trignum = k; 737 gi0.u = seg.epgeominfo[0].u; 738 gi0.v = seg.epgeominfo[0].v; 739 gi1.u = seg.epgeominfo[1].u; 740 gi1.v = seg.epgeominfo[1].v; 741 742 int locpnum[2] = {0, 0}; 743 744 for (j = 0; j < 2; j++) 745 { 746 PointGeomInfo gi = (j == 0) ? gi0 : gi1; 747 748 int l; 749 for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) 750 { 751 double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); 752 753 if (dist < 1e-10) 754 locpnum[j] = l+1; 755 } 756 757 if (locpnum[j] == 0) 758 { 759 int pi = (j == 0) ? seg[0] : seg[1]; 760 meshing.AddPoint (mesh.Point(pi), pi); 761 762 gis.SetSize (gis.Size()+1); 763 gis[l] = gi; 764 locpnum[j] = l+1; 765 } 766 } 767 768 meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); 769 //(*testout) << gi0.u << " " << gi0.v << endl; 770 //(*testout) << gi1.u << " " << gi1.v << endl; 771 772 } 773 } 774 } 775 776 777 778 779 780 // Philippose - 15/01/2009 781 double maxh = geom.face_maxh[k-1]; 782 //double maxh = mparam.maxh; 783 mparam.checkoverlap = 0; 784 // int noldpoints = mesh->GetNP(); 785 int noldsurfel = mesh.GetNSE(); 786 787 GProp_GProps sprops; 788 BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops); 789 meshing.SetMaxArea(2.*sprops.Mass()); 790 791 MESHING2_RESULT res; 792 793 try { 794 res = meshing.GenerateMesh (mesh, mparam, maxh, k); 795 } 796 797 catch (SingularMatrixException) 798 { 799 (*myerr) << "Singular Matrix" << endl; 800 res = MESHING2_GIVEUP; 801 } 802 803 catch (UVBoundsException) 804 { 805 (*myerr) << "UV bounds exceeded" << endl; 806 res = MESHING2_GIVEUP; 807 } 808 809 projecttype = PARAMETERSPACE; 810 811 if (res != MESHING2_OK) 812 { 813 if (notrys == 1) 814 { 815 for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++) 816 mesh.DeleteSurfaceElement (i); 817 818 mesh.Compress(); 819 820 cout << "retry Surface " << k << endl; 821 822 k--; 823 projecttype*=-1; 824 notrys++; 825 continue; 826 } 827 else 828 { 829 geom.facemeshstatus[k-1] = -1; 830 PrintError ("Problem in Surface mesh generation"); 831 surfmesherror++; 832 // throw NgException ("Problem in Surface mesh generation"); 833 } 834 } 835 else 836 { 837 geom.facemeshstatus[k-1] = 1; 838 } 839 840 notrys = 1; 841 842 for (i = oldnf+1; i <= mesh.GetNSE(); i++) 843 mesh.SurfaceElement(i).SetIndex (k); 844 845 } 846 847 // ofstream problemfile("occmesh.rep"); 848 849 // problemfile << "SURFACEMESHING" << endl << endl; 850 851 if (surfmesherror) 852 { 853 cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl; 854 cout << "SURFACE MESHING ERROR OCCURED IN " << surfmesherror << " FACES:" << endl; 855 for (int i = 1; i <= geom.fmap.Extent(); i++) 856 if (geom.facemeshstatus[i-1] == -1) 857 { 858 cout << "Face " << i << endl; 859 // problemfile << "problem with face " << i << endl; 860 // problemfile << "vertices: " << endl; 861 TopExp_Explorer exp0,exp1,exp2; 862 for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() ) 863 { 864 TopoDS_Wire wire = TopoDS::Wire(exp0.Current()); 865 for ( exp1.Init(wire,TopAbs_EDGE); exp1.More(); exp1.Next() ) 866 { 867 TopoDS_Edge edge = TopoDS::Edge(exp1.Current()); 868 for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() ) 869 { 870 TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current()); 871 gp_Pnt point = BRep_Tool::Pnt(vertex); 872 // problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; 873 } 874 } 875 } 876 // problemfile << endl; 877 878 } 879 cout << endl << endl; 880 cout << "for more information open IGES/STEP Topology Explorer" << endl; 881 // problemfile.close(); 882 throw NgException ("Problem in Surface mesh generation"); 883 } 884 else 885 { 886 // problemfile << "OK" << endl << endl; 887 // problemfile.close(); 888 } 889 890 891 892 893 if (multithread.terminate || perfstepsend < MESHCONST_OPTSURFACE) 894 return; 895 896 multithread.task = "Optimizing surface"; 897 898 static int timer_opt2d = NgProfiler::CreateTimer ("Optimization 2D"); 899 NgProfiler::StartTimer (timer_opt2d); 900 901 for (k = 1; k <= mesh.GetNFD(); k++) 902 { 903 // if (k != 42) continue; 904 // if (k != 36) continue; 905 906 // (*testout) << "optimize face " << k << endl; 907 multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); 908 909 FaceDescriptor & fd = mesh.GetFaceDescriptor(k); 910 911 PrintMessage (1, "Optimize Surface ", k); 912 for (i = 1; i <= mparam.optsteps2d; i++) 913 { 914 // (*testout) << "optstep " << i << endl; 915 if (multithread.terminate) return; 916 917 { 918 MeshOptimize2dOCCSurfaces meshopt(geom); 919 meshopt.SetFaceIndex (k); 920 meshopt.SetImproveEdges (0); 921 meshopt.SetMetricWeight (mparam.elsizeweight); 922 //meshopt.SetMetricWeight (0.2); 923 meshopt.SetWriteStatus (0); 924 925 // (*testout) << "EdgeSwapping (mesh, (i > mparam.optsteps2d/2))" << endl; 926 meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2)); 927 } 928 929 if (multithread.terminate) return; 930 { 931 MeshOptimize2dOCCSurfaces meshopt(geom); 932 meshopt.SetFaceIndex (k); 933 meshopt.SetImproveEdges (0); 934 //meshopt.SetMetricWeight (0.2); 935 meshopt.SetMetricWeight (mparam.elsizeweight); 936 meshopt.SetWriteStatus (0); 937 938 // (*testout) << "ImproveMesh (mesh)" << endl; 939 meshopt.ImproveMesh (mesh, mparam); 940 } 941 942 { 943 MeshOptimize2dOCCSurfaces meshopt(geom); 944 meshopt.SetFaceIndex (k); 945 meshopt.SetImproveEdges (0); 946 //meshopt.SetMetricWeight (0.2); 947 meshopt.SetMetricWeight (mparam.elsizeweight); 948 meshopt.SetWriteStatus (0); 949 950 // (*testout) << "CombineImprove (mesh)" << endl; 951 meshopt.CombineImprove (mesh); 952 } 953 954 if (multithread.terminate) return; 955 { 956 MeshOptimize2dOCCSurfaces meshopt(geom); 957 meshopt.SetFaceIndex (k); 958 meshopt.SetImproveEdges (0); 959 //meshopt.SetMetricWeight (0.2); 960 meshopt.SetMetricWeight (mparam.elsizeweight); 961 meshopt.SetWriteStatus (0); 962 963 // (*testout) << "ImproveMesh (mesh)" << endl; 964 meshopt.ImproveMesh (mesh, mparam); 965 } 966 } 967 968 } 969 970 971 mesh.CalcSurfacesOfNode(); 972 mesh.Compress(); 973 974 NgProfiler::StopTimer (timer_opt2d); 975 976 multithread.task = savetask; 977 } 978 979 980 OCCSetLocalMeshSize(OCCGeometry & geom,Mesh & mesh)981 void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh) 982 { 983 mesh.SetGlobalH (mparam.maxh); 984 mesh.SetMinimalH (mparam.minh); 985 986 Array<double> maxhdom; 987 maxhdom.SetSize (geom.NrSolids()); 988 maxhdom = mparam.maxh; 989 990 mesh.SetMaxHDomain (maxhdom); 991 992 Box<3> bb = geom.GetBoundingBox(); 993 bb.Increase (bb.Diam()/10); 994 995 mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); 996 997 if (mparam.uselocalh) 998 { 999 const char * savetask = multithread.task; 1000 multithread.percent = 0; 1001 1002 mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); 1003 1004 int nedges = geom.emap.Extent(); 1005 1006 double maxedgelen = 0; 1007 double minedgelen = 1e99; 1008 1009 multithread.task = "Setting local mesh size (elements per edge)"; 1010 1011 // setting elements per edge 1012 1013 for (int i = 1; i <= nedges && !multithread.terminate; i++) 1014 { 1015 TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); 1016 multithread.percent = 100 * (i-1)/double(nedges); 1017 if (BRep_Tool::Degenerated(e)) continue; 1018 1019 GProp_GProps system; 1020 BRepGProp::LinearProperties(e, system); 1021 double len = system.Mass(); 1022 1023 if (len < IGNORECURVELENGTH) 1024 { 1025 (*testout) << "ignored" << endl; 1026 continue; 1027 } 1028 1029 double localh = len/mparam.segmentsperedge; 1030 double s0, s1; 1031 1032 // Philippose - 23/01/2009 1033 // Find all the parent faces of a given edge 1034 // and limit the mesh size of the edge based on the 1035 // mesh size limit of the face 1036 TopTools_IndexedDataMapOfShapeListOfShape edge_face_map; 1037 edge_face_map.Clear(); 1038 1039 TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map); 1040 const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e); 1041 1042 TopTools_ListIteratorOfListOfShape parent_face_list; 1043 1044 for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next()) 1045 { 1046 TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); 1047 1048 int face_index = geom.fmap.FindIndex(parent_face); 1049 1050 if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); 1051 } 1052 1053 Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); 1054 1055 maxedgelen = max (maxedgelen, len); 1056 minedgelen = min (minedgelen, len); 1057 1058 // Philippose - 23/01/2009 1059 // Modified the calculation of maxj, because the 1060 // method used so far always results in maxj = 2, 1061 // which causes the localh to be set only at the 1062 // starting, mid and end of the edge. 1063 // Old Algorithm: 1064 // int maxj = 2 * (int) ceil (localh/len); 1065 int maxj = max((int) ceil(len/localh), 2); 1066 1067 for (int j = 0; j <= maxj; j++) 1068 { 1069 gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); 1070 mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh); 1071 } 1072 } 1073 1074 multithread.task = "Setting local mesh size (edge curvature)"; 1075 1076 // setting edge curvature 1077 1078 int nsections = 20; 1079 1080 for (int i = 1; i <= nedges && !multithread.terminate; i++) 1081 { 1082 double maxcur = 0; 1083 multithread.percent = 100 * (i-1)/double(nedges); 1084 TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); 1085 if (BRep_Tool::Degenerated(edge)) continue; 1086 double s0, s1; 1087 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); 1088 BRepAdaptor_Curve brepc(edge); 1089 BRepLProp_CLProps prop(brepc, 2, 1e-5); 1090 1091 for (int j = 1; j <= nsections; j++) 1092 { 1093 double s = s0 + j/(double) nsections * (s1-s0); 1094 prop.SetParameter (s); 1095 double curvature = prop.Curvature(); 1096 if(curvature> maxcur) maxcur = curvature; 1097 1098 if (curvature >= 1e99) 1099 continue; 1100 1101 gp_Pnt pnt = c->Value (s); 1102 1103 mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature))); 1104 } 1105 // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl; 1106 } 1107 1108 multithread.task = "Setting local mesh size (face curvature)"; 1109 1110 // setting face curvature 1111 1112 int nfaces = geom.fmap.Extent(); 1113 1114 for (int i = 1; i <= nfaces && !multithread.terminate; i++) 1115 { 1116 multithread.percent = 100 * (i-1)/double(nfaces); 1117 TopoDS_Face face = TopoDS::Face(geom.fmap(i)); 1118 TopLoc_Location loc; 1119 Handle(Geom_Surface) surf = BRep_Tool::Surface (face); 1120 Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); 1121 1122 if (triangulation.IsNull()) continue; 1123 1124 BRepAdaptor_Surface sf(face, Standard_True); 1125 BRepLProp_SLProps prop(sf, 2, 1e-5); 1126 1127 int ntriangles = triangulation -> NbTriangles(); 1128 for (int j = 1; j <= ntriangles; j++) 1129 { 1130 gp_Pnt p[3]; 1131 gp_Pnt2d par[3]; 1132 1133 for (int k = 1; k <=3; k++) 1134 { 1135 int n = triangulation->Triangles()(j)(k); 1136 p[k-1] = triangulation->Nodes()(n).Transformed(loc); 1137 par[k-1] = triangulation->UVNodes()(n); 1138 } 1139 1140 //double maxside = 0; 1141 //maxside = max (maxside, p[0].Distance(p[1])); 1142 //maxside = max (maxside, p[0].Distance(p[2])); 1143 //maxside = max (maxside, p[1].Distance(p[2])); 1144 //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; 1145 1146 RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0); 1147 //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; 1148 } 1149 } 1150 1151 // setting close edges 1152 1153 if (occparam.resthcloseedgeenable) 1154 { 1155 multithread.task = "Setting local mesh size (close edges)"; 1156 1157 int sections = 100; 1158 1159 Array<Line> lines(sections*nedges); 1160 1161 Box3dTree* searchtree = 1162 new Box3dTree (bb.PMin(), bb.PMax()); 1163 1164 int nlines = 0; 1165 for (int i = 1; i <= nedges && !multithread.terminate; i++) 1166 { 1167 TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); 1168 if (BRep_Tool::Degenerated(edge)) continue; 1169 1170 double s0, s1; 1171 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); 1172 BRepAdaptor_Curve brepc(edge); 1173 BRepLProp_CLProps prop(brepc, 1, 1e-5); 1174 prop.SetParameter (s0); 1175 1176 gp_Vec d0 = prop.D1().Normalized(); 1177 double s_start = s0; 1178 int count = 0; 1179 for (int j = 1; j <= sections; j++) 1180 { 1181 double s = s0 + (s1-s0)*(double)j/(double)sections; 1182 prop.SetParameter (s); 1183 gp_Vec d1 = prop.D1().Normalized(); 1184 double cosalpha = fabs(d0*d1); 1185 if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI))) 1186 { 1187 count++; 1188 gp_Pnt p0 = c->Value (s_start); 1189 gp_Pnt p1 = c->Value (s); 1190 lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); 1191 lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); 1192 1193 Box3d box; 1194 box.SetPoint (Point3d(lines[nlines].p0)); 1195 box.AddPoint (Point3d(lines[nlines].p1)); 1196 1197 searchtree->Insert (box.PMin(), box.PMax(), nlines+1); 1198 nlines++; 1199 1200 s_start = s; 1201 d0 = d1; 1202 } 1203 } 1204 } 1205 1206 Array<int> linenums; 1207 1208 for (int i = 0; i < nlines; i++) 1209 { 1210 multithread.percent = (100*i)/double(nlines); 1211 Line & line = lines[i]; 1212 1213 Box3d box; 1214 box.SetPoint (Point3d(line.p0)); 1215 box.AddPoint (Point3d(line.p1)); 1216 double maxhline = max (mesh.GetH(box.PMin()), 1217 mesh.GetH(box.PMax())); 1218 box.Increase(maxhline); 1219 1220 double mindist = 1e99; 1221 linenums.SetSize(0); 1222 searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums); 1223 1224 for (int j = 0; j < linenums.Size(); j++) 1225 { 1226 int num = linenums[j]-1; 1227 if (i == num) continue; 1228 if ((line.p0-lines[num].p0).Length2() < 1e-15) continue; 1229 if ((line.p0-lines[num].p1).Length2() < 1e-15) continue; 1230 if ((line.p1-lines[num].p0).Length2() < 1e-15) continue; 1231 if ((line.p1-lines[num].p1).Length2() < 1e-15) continue; 1232 mindist = min (mindist, line.Dist(lines[num])); 1233 } 1234 1235 mindist /= (occparam.resthcloseedgefac + VSMALL); 1236 1237 if (mindist < 1e-3) 1238 { 1239 (*testout) << "extremely small local h: " << mindist 1240 << " --> setting to 1e-3" << endl; 1241 (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl; 1242 mindist = 1e-3; 1243 } 1244 1245 mesh.RestrictLocalHLine(line.p0, line.p1, mindist); 1246 } 1247 } 1248 1249 multithread.task = savetask; 1250 1251 } 1252 1253 // Philippose - 09/03/2009 1254 // Added the capability to load the mesh size from a 1255 // file also for OpenCascade Geometry 1256 // Note: 1257 // ** If the "uselocalh" option is ticked in 1258 // the "mesh options...insider" menu, the mesh 1259 // size will be further modified by the topology 1260 // analysis routines. 1261 // ** To use the mesh size file as the sole source 1262 // for defining the mesh size, uncheck the "uselocalh" 1263 // option. 1264 mesh.LoadLocalMeshSize (mparam.meshsizefilename); 1265 } 1266 1267 1268 OCCGenerateMesh(OCCGeometry & geom,Mesh * & mesh,MeshingParameters & mparam,int perfstepsstart,int perfstepsend)1269 int OCCGenerateMesh (OCCGeometry & geom, Mesh *& mesh, MeshingParameters & mparam, 1270 int perfstepsstart, int perfstepsend) 1271 { 1272 multithread.percent = 0; 1273 1274 if (perfstepsstart <= MESHCONST_ANALYSE) 1275 { 1276 delete mesh; 1277 mesh = new Mesh(); 1278 mesh->geomtype = Mesh::GEOM_OCC; 1279 1280 OCCSetLocalMeshSize(geom,*mesh); 1281 } 1282 1283 if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE) 1284 return TCL_OK; 1285 1286 if (perfstepsstart <= MESHCONST_MESHEDGES) 1287 { 1288 OCCFindEdges (geom, *mesh); 1289 1290 /* 1291 cout << "Removing redundant points" << endl; 1292 1293 int i, j; 1294 int np = mesh->GetNP(); 1295 Array<int> equalto; 1296 1297 equalto.SetSize (np); 1298 equalto = 0; 1299 1300 for (i = 1; i <= np; i++) 1301 { 1302 for (j = i+1; j <= np; j++) 1303 { 1304 if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12)) 1305 equalto[j-1] = i; 1306 } 1307 } 1308 1309 for (i = 1; i <= np; i++) 1310 if (equalto[i-1]) 1311 { 1312 cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl; 1313 for (j = 1; j <= mesh->GetNSeg(); j++) 1314 { 1315 Segment & seg = mesh->LineSegment(j); 1316 if (seg[0] == i) seg[0] = equalto[i-1]; 1317 if (seg[1] == i) seg[1] = equalto[i-1]; 1318 } 1319 } 1320 1321 cout << "Removing degenerated segments" << endl; 1322 for (j = 1; j <= mesh->GetNSeg(); j++) 1323 { 1324 Segment & seg = mesh->LineSegment(j); 1325 if (seg[0] == seg[1]) 1326 { 1327 mesh->DeleteSegment(j); 1328 cout << "Deleting Segment " << j << endl; 1329 } 1330 } 1331 1332 mesh->Compress(); 1333 */ 1334 1335 /* 1336 for (int i = 1; i <= geom.fmap.Extent(); i++) 1337 { 1338 Handle(Geom_Surface) hf1 = 1339 BRep_Tool::Surface(TopoDS::Face(geom.fmap(i))); 1340 for (int j = i+1; j <= geom.fmap.Extent(); j++) 1341 { 1342 Handle(Geom_Surface) hf2 = 1343 BRep_Tool::Surface(TopoDS::Face(geom.fmap(j))); 1344 if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl; 1345 } 1346 } 1347 */ 1348 1349 #ifdef LOG_STREAM 1350 (*logout) << "Edges meshed" << endl 1351 << "time = " << GetTime() << " sec" << endl 1352 << "points: " << mesh->GetNP() << endl; 1353 #endif 1354 } 1355 1356 if (multithread.terminate || perfstepsend <= MESHCONST_MESHEDGES) 1357 return TCL_OK; 1358 1359 if (perfstepsstart <= MESHCONST_MESHSURFACE) 1360 { 1361 OCCMeshSurface (geom, *mesh, perfstepsend); 1362 if (multithread.terminate) return TCL_OK; 1363 1364 #ifdef LOG_STREAM 1365 (*logout) << "Surfaces meshed" << endl 1366 << "time = " << GetTime() << " sec" << endl 1367 << "points: " << mesh->GetNP() << endl; 1368 #endif 1369 1370 #ifdef STAT_STREAM 1371 (*statout) << mesh->GetNSeg() << " & " 1372 << mesh->GetNSE() << " & - &" 1373 << GetTime() << " & " << endl; 1374 #endif 1375 1376 // MeshQuality2d (*mesh); 1377 mesh->CalcSurfacesOfNode(); 1378 } 1379 1380 if (multithread.terminate || perfstepsend <= MESHCONST_OPTSURFACE) 1381 return TCL_OK; 1382 1383 if (perfstepsstart <= MESHCONST_MESHVOLUME) 1384 { 1385 multithread.task = "Volume meshing"; 1386 1387 MESHING3_RESULT res = MeshVolume (mparam, *mesh); 1388 1389 /* 1390 ofstream problemfile("occmesh.rep",ios_base::app); 1391 1392 problemfile << "VOLUMEMESHING" << endl << endl; 1393 if(res != MESHING3_OK) 1394 problemfile << "ERROR" << endl << endl; 1395 else 1396 problemfile << "OK" << endl 1397 << mesh->GetNE() << " elements" << endl << endl; 1398 1399 problemfile.close(); 1400 */ 1401 1402 if (res != MESHING3_OK) return TCL_ERROR; 1403 1404 if (multithread.terminate) return TCL_OK; 1405 1406 RemoveIllegalElements (*mesh); 1407 if (multithread.terminate) return TCL_OK; 1408 1409 MeshQuality3d (*mesh); 1410 1411 #ifdef STAT_STREAM 1412 (*statout) << GetTime() << " & "; 1413 #endif 1414 1415 #ifdef LOG_STREAM 1416 (*logout) << "Volume meshed" << endl 1417 << "time = " << GetTime() << " sec" << endl 1418 << "points: " << mesh->GetNP() << endl; 1419 #endif 1420 } 1421 1422 if (multithread.terminate || perfstepsend <= MESHCONST_MESHVOLUME) 1423 return TCL_OK; 1424 1425 if (perfstepsstart <= MESHCONST_OPTVOLUME) 1426 { 1427 multithread.task = "Volume optimization"; 1428 1429 OptimizeVolume (mparam, *mesh); 1430 if (multithread.terminate) return TCL_OK; 1431 1432 #ifdef STAT_STREAM 1433 (*statout) << GetTime() << " & " 1434 << mesh->GetNE() << " & " 1435 << mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl; 1436 #endif 1437 1438 #ifdef LOG_STREAM 1439 (*logout) << "Volume optimized" << endl 1440 << "time = " << GetTime() << " sec" << endl 1441 << "points: " << mesh->GetNP() << endl; 1442 #endif 1443 1444 // cout << "Optimization complete" << endl; 1445 1446 } 1447 1448 (*testout) << "NP: " << mesh->GetNP() << endl; 1449 for (int i = 1; i <= mesh->GetNP(); i++) 1450 (*testout) << mesh->Point(i) << endl; 1451 1452 (*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl; 1453 for (int i = 1; i <= mesh->GetNSeg(); i++) 1454 (*testout) << mesh->LineSegment(i) << endl; 1455 1456 return TCL_OK; 1457 } 1458 } 1459 1460 #endif 1461