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