1 #include <mystdlib.h>
2 
3 
4 #include <myadt.hpp>
5 
6 #include <linalg.hpp>
7 #include <csg.hpp>
8 #include <meshing.hpp>
9 
10 
11 namespace netgen
12 {
13 
14   NgArray<SpecialPoint> global_specpoints;  // for visualization
15   //static NgArray<MeshPoint> spoints;
16 
17 #define TCL_OK 0
18 #define TCL_ERROR 1
19 
20 
21 
FindPoints(CSGeometry & geom,NgArray<SpecialPoint> & specpoints,NgArray<MeshPoint> & spoints,Mesh & mesh)22   static void FindPoints (CSGeometry & geom,
23                           NgArray<SpecialPoint> &  specpoints,
24                           NgArray<MeshPoint> & spoints,
25                           Mesh & mesh)
26   {
27     PrintMessage (1, "Start Findpoints");
28 
29     const char * savetask = multithread.task;
30     multithread.task = "Find points";
31 
32     mesh.pointelements.SetSize(0);
33     for (int i = 0; i < geom.GetNUserPoints(); i++)
34       {
35         auto up = geom.GetUserPoint(i);
36 	auto pnum = mesh.AddPoint(up);
37 	mesh.Points().Last().Singularity (geom.GetUserPointRefFactor(i));
38 	mesh.AddLockedPoint (PointIndex (i+1));
39         int index = up.GetIndex();
40         if (index == -1)
41           index = mesh.AddCD3Name (up.GetName())+1;
42         // cout << "adding 0d element, pnum = " << pnum << ", material index = " << index << endl;
43         mesh.pointelements.Append (Element0d(pnum, index));
44       }
45 
46     SpecialPointCalculation spc;
47 
48     spc.SetIdEps(geom.GetIdEps());
49 
50     if (spoints.Size() == 0)
51       spc.CalcSpecialPoints (geom, spoints);
52 
53     PrintMessage (2, "Analyze spec points");
54     spc.AnalyzeSpecialPoints (geom, spoints, specpoints);
55 
56     {
57       static mutex mut;
58       lock_guard<mutex> guard(mut);
59       global_specpoints = specpoints;
60     }
61 
62     PrintMessage (5, "done");
63 
64     (*testout) << specpoints.Size() << " special points:" << endl;
65     for (int i = 0; i < specpoints.Size(); i++)
66       specpoints[i].Print (*testout);
67 
68     /*
69       for (int i = 1; i <= geom.identifications.Size(); i++)
70       geom.identifications.Elem(i)->IdentifySpecialPoints (specpoints);
71     */
72     multithread.task = savetask;
73   }
74 
75 
76 
77 
78 
79 
FindEdges(CSGeometry & geom,Mesh & mesh,NgArray<SpecialPoint> & specpoints,NgArray<MeshPoint> & spoints,MeshingParameters & mparam,const bool setmeshsize=false)80   static void FindEdges (CSGeometry & geom, Mesh & mesh,
81                          NgArray<SpecialPoint> &  specpoints,
82                          NgArray<MeshPoint> & spoints,
83                          MeshingParameters & mparam,
84                          const bool setmeshsize = false)
85   {
86     EdgeCalculation ec (geom, specpoints, mparam);
87     ec.SetIdEps(geom.GetIdEps());
88     ec.Calc (mparam.maxh, mesh);
89 
90     for (int i = 0; i < geom.singedges.Size(); i++)
91       {
92 	geom.singedges[i]->FindPointsOnEdge (mesh);
93 	if(setmeshsize)
94 	  geom.singedges[i]->SetMeshSize(mesh,10.*geom.BoundingBox().Diam());
95       }
96     for (int i = 0; i < geom.singpoints.Size(); i++)
97       geom.singpoints[i]->FindPoints (mesh);
98 
99     for (int i = 1; i <= mesh.GetNSeg(); i++)
100       {
101 	//(*testout) << "segment " << mesh.LineSegment(i) << endl;
102 	int ok = 0;
103 	for (int k = 1; k <= mesh.GetNFD(); k++)
104 	  if (mesh.GetFaceDescriptor(k).SegmentFits (mesh.LineSegment(i)))
105 	    {
106 	      ok = k;
107 	      //(*testout) << "fits to " << k << endl;
108 	    }
109 
110 	if (!ok)
111 	  {
112 	    ok = mesh.AddFaceDescriptor (FaceDescriptor (mesh.LineSegment(i)));
113 	    //(*testout) << "did not find, now " << ok << endl;
114 	  }
115 
116 	//(*testout) << "change from " << mesh.LineSegment(i).si;
117 	mesh.LineSegment(i).si = ok;
118 	//(*testout) << " to " << mesh.LineSegment(i).si << endl;
119       }
120 
121     for(int k = 1; k<=mesh.GetNFD(); k++)
122       {
123 	*testout << "face: " << k << endl
124 		 << "FD: " << mesh.GetFaceDescriptor(k) << endl;
125       }
126 
127     if (geom.identifications.Size())
128       {
129 	PrintMessage (3, "Find Identifications");
130 	for (int i = 0; i < geom.identifications.Size(); i++)
131 	  {
132 	    geom.identifications[i]->IdentifyPoints (mesh);
133 	    //(*testout) << "identification " << i << " is "
134 	    //	       << *geom.identifications[i] << endl;
135 
136 	  }
137 	for (int i = 0; i < geom.identifications.Size(); i++)
138 	  geom.identifications[i]->IdentifyFaces (mesh);
139       }
140 
141 
142     // find intersecting segments
143     PrintMessage (3, "Check intersecting edges");
144 
145     Point3d pmin, pmax;
146     mesh.GetBox (pmin, pmax);
147     BoxTree<3> segtree (pmin, pmax);
148 
149     for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
150       {
151 	if (mesh[si].seginfo)
152 	  {
153 	    Box<3> hbox;
154 	    hbox.Set (mesh[mesh[si][0]]);
155 	    hbox.Add (mesh[mesh[si][1]]);
156 	    segtree.Insert (hbox.PMin(), hbox.PMax(), si);
157 	  }
158       }
159 
160     NgArray<int> loc;
161     if (!ec.point_on_edge_problem)
162       for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
163 	{
164 	  if (!mesh[si].seginfo) continue;
165 
166 	  Box<3> hbox;
167 	  hbox.Set (mesh[mesh[si][0]]);
168 	  hbox.Add (mesh[mesh[si][1]]);
169 	  hbox.Increase (1e-6);
170 	  segtree.GetIntersecting (hbox.PMin(), hbox.PMax(), loc);
171 
172 	  // for (SegmentIndex sj = 0; sj < si; sj++)
173 	  for (int j = 0; j < loc.Size(); j++)
174 	    {
175 	      SegmentIndex sj = loc[j];
176 	      if (sj >= si) continue;
177 	      if (!mesh[si].seginfo || !mesh[sj].seginfo) continue;
178 	      if (mesh[mesh[si][0]].GetLayer() != mesh[mesh[sj][1]].GetLayer()) continue;
179 
180 	      Point<3> pi1 = mesh[mesh[si][0]];
181 	      Point<3> pi2 = mesh[mesh[si][1]];
182 	      Point<3> pj1 = mesh[mesh[sj][0]];
183 	      Point<3> pj2 = mesh[mesh[sj][1]];
184 	      Vec<3> vi = pi2 - pi1;
185 	      Vec<3> vj = pj2 - pj1;
186 
187 	      if (sqr (vi * vj) > (1.-1e-6) * Abs2 (vi) * Abs2 (vj)) continue;
188 
189 	      // pi1 + vi t = pj1 + vj s
190 	      Mat<3,2> mat;
191 	      Vec<3> rhs;
192 	      Vec<2> sol;
193 
194 	      for (int jj = 0; jj < 3; jj++)
195 		{
196 		  mat(jj,0) = vi(jj);
197 		  mat(jj,1) = -vj(jj);
198 		  rhs(jj) = pj1(jj)-pi1(jj);
199 		}
200 
201 	      mat.Solve (rhs, sol);
202 
203 	      //(*testout) << "mat " << mat << endl << "rhs " << rhs << endl << "sol " << sol << endl;
204 
205 	      if (sol(0) > 1e-6 && sol(0) < 1-1e-6 &&
206 		  sol(1) > 1e-6 && sol(1) < 1-1e-6 &&
207 		  Abs (rhs - mat*sol) < 1e-6)
208 		{
209 		  Point<3> ip = pi1 + sol(0) * vi;
210 
211 		  //(*testout) << "ip " << ip << endl;
212 
213 		  Point<3> pip = ip;
214 		  ProjectToEdge (geom.GetSurface (mesh[si].surfnr1),
215 				 geom.GetSurface (mesh[si].surfnr2), pip);
216 
217 		  //(*testout) << "Dist (ip, pip_si) " << Dist (ip, pip) << endl;
218 		  if (Dist (ip, pip) > 1e-6*geom.MaxSize()) continue;
219 		  pip = ip;
220 		  ProjectToEdge (geom.GetSurface (mesh[sj].surfnr1),
221 				 geom.GetSurface (mesh[sj].surfnr2), pip);
222 
223 		  //(*testout) << "Dist (ip, pip_sj) " << Dist (ip, pip) << endl;
224 		  if (Dist (ip, pip) > 1e-6*geom.MaxSize()) continue;
225 
226 
227 
228 		  cout << "Intersection at " << ip << endl;
229 
230 		  geom.AddUserPoint (ip);
231 		  spoints.Append (MeshPoint (ip, mesh[mesh[si][0]].GetLayer()));
232 		  mesh.AddPoint (ip);
233 
234 		  (*testout) << "found intersection at " << ip << endl;
235 		  (*testout) << "sol = " << sol << endl;
236 		  (*testout) << "res = " << (rhs - mat*sol) << endl;
237 		  (*testout) << "segs = " << pi1 << " - " << pi2 << endl;
238 		  (*testout) << "and = " << pj1 << " - " << pj2 << endl << endl;
239 		}
240 	    }
241 	}
242   }
243 
244 
245 
246 
247 
248 
MeshSurface(CSGeometry & geom,Mesh & mesh,MeshingParameters & mparam)249   static void MeshSurface (CSGeometry & geom, Mesh & mesh, MeshingParameters & mparam)
250   {
251     const char * savetask = multithread.task;
252     multithread.task = "Surface meshing";
253 
254     NgArray<Segment> segments;
255     int noldp = mesh.GetNP();
256 
257     double starttime = GetTime();
258 
259     // find master faces from identified
260     NgArray<int> masterface(mesh.GetNFD());
261     for (int i = 1; i <= mesh.GetNFD(); i++)
262       masterface.Elem(i) = i;
263 
264     NgArray<INDEX_2> fpairs;
265     bool changed;
266     do
267       {
268 	changed = 0;
269 	for (int i = 0; i < geom.identifications.Size(); i++)
270 	  {
271 	    geom.identifications[i]->GetIdentifiedFaces (fpairs);
272 
273 	    for (int j = 0; j < fpairs.Size(); j++)
274 	      {
275 		if (masterface.Get(fpairs[j].I1()) <
276 		    masterface.Get(fpairs[j].I2()))
277 		  {
278 		    changed = 1;
279 		    masterface.Elem(fpairs[j].I2()) =
280 		      masterface.Elem(fpairs[j].I1());
281 		  }
282 		if (masterface.Get(fpairs[j].I2()) <
283 		    masterface.Get(fpairs[j].I1()))
284 		  {
285 		    changed = 1;
286 		    masterface.Elem(fpairs[j].I1()) =
287 		      masterface.Elem(fpairs[j].I2());
288 		  }
289 	      }
290 	  }
291       }
292     while (changed);
293 
294 
295     int bccnt=0;
296     for (int k = 0; k < geom.GetNSurf(); k++)
297       bccnt = max2 (bccnt, geom.GetSurface(k)->GetBCProperty());
298 
299     for (int k = 1; k <= mesh.GetNFD(); k++)
300       {
301 	bool increased = false;
302 
303 	FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
304 	const Surface * surf = geom.GetSurface(fd.SurfNr());
305 
306 	if (fd.TLOSurface() &&
307 	    geom.GetTopLevelObject(fd.TLOSurface()-1) -> GetBCProp() > 0)
308 	  fd.SetBCProperty (geom.GetTopLevelObject(fd.TLOSurface()-1) -> GetBCProp());
309 	else if (surf -> GetBCProperty() != -1)
310 	  fd.SetBCProperty (surf->GetBCProperty());
311 	else
312 	  {
313 	    bccnt++;
314 	    fd.SetBCProperty (bccnt);
315 	    increased = true;
316 	  }
317 
318 	for (int l = 0; l < geom.bcmodifications.Size(); l++)
319 	  {
320 	    if (geom.GetSurfaceClassRepresentant (fd.SurfNr()) ==
321 		geom.GetSurfaceClassRepresentant (geom.bcmodifications[l].si) &&
322 		(fd.DomainIn() == geom.bcmodifications[l].tlonr+1 ||
323 		 fd.DomainOut() == geom.bcmodifications[l].tlonr+1))
324 	      {
325 		if(geom.bcmodifications[l].bcname == NULL)
326 		  fd.SetBCProperty (geom.bcmodifications[l].bcnr);
327 		else
328 		  {
329 		    if(!increased)
330 		      {
331 			bccnt++;
332 			fd.SetBCProperty (bccnt);
333 			increased = true;
334 		      }
335 		  }
336 	      }
337 	  }
338       }
339 
340     mesh.SetNBCNames( bccnt );
341 
342     for (int k = 1; k <= mesh.GetNFD(); k++)
343       {
344 	FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
345 	const Surface * surf = geom.GetSurface(fd.SurfNr());
346 	if (fd.TLOSurface() )
347 	  {
348 	    int bcp = fd.BCProperty();
349 	    string nextbcname = geom.GetTopLevelObject(fd.TLOSurface()-1) -> GetBCName();
350 	    if ( nextbcname != "default" )
351 	      mesh.SetBCName ( bcp - 1 , nextbcname );
352 	  }
353 	else // if (surf -> GetBCProperty() != -1)
354 	  {
355 	    int bcp = fd.BCProperty();
356 	    string nextbcname = surf->GetBCName();
357 	    if ( nextbcname != "default" )
358 	      mesh.SetBCName ( bcp - 1, nextbcname );
359 	  }
360       }
361 
362     for (int k = 1; k <= mesh.GetNFD(); k++)
363       {
364 	FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
365 	fd.SetBCName ( mesh.GetBCNamePtr ( fd.BCProperty() - 1 ) );
366       }
367 
368     //!!
369 
370     for (int k = 1; k <= mesh.GetNFD(); k++)
371       {
372 	FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
373 	//const Surface * surf = geom.GetSurface(fd.SurfNr());
374 
375 	for (int l = 0; l < geom.bcmodifications.Size(); l++)
376 	  {
377 	    if (geom.GetSurfaceClassRepresentant (fd.SurfNr()) ==
378 		geom.GetSurfaceClassRepresentant (geom.bcmodifications[l].si) &&
379 		(fd.DomainIn() == geom.bcmodifications[l].tlonr+1 ||
380 		 fd.DomainOut() == geom.bcmodifications[l].tlonr+1) &&
381 		geom.bcmodifications[l].bcname != NULL
382 		)
383 	      {
384 		int bcp = fd.BCProperty();
385 		mesh.SetBCName ( bcp - 1, *(geom.bcmodifications[l].bcname) );
386 		fd.SetBCName ( mesh.GetBCNamePtr ( bcp - 1) );
387 	      }
388 	  }
389       }
390 
391     for(int k = 0; k<geom.bcmodifications.Size(); k++)
392       {
393 	delete geom.bcmodifications[k].bcname;
394 	geom.bcmodifications[k].bcname = NULL;
395       }
396 
397     //!!
398 
399 
400     for (int j = 0; j < geom.singfaces.Size(); j++)
401       {
402 	NgArray<int> surfs;
403 	geom.GetIndependentSurfaceIndices (geom.singfaces[j]->GetSolid(),
404 					   geom.BoundingBox(), surfs);
405 	for (int k = 1; k <= mesh.GetNFD(); k++)
406 	  {
407 	    FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
408 	    for (int l = 0; l < surfs.Size(); l++)
409 	      if (surfs[l] == fd.SurfNr())
410 		{
411 		  if (geom.singfaces[j]->GetDomainNr() == fd.DomainIn())
412 		    fd.SetDomainInSingular (1);
413 		  if (geom.singfaces[j]->GetDomainNr() == fd.DomainOut())
414 		    fd.SetDomainOutSingular (1);
415 		}
416 	  }
417       }
418 
419 
420     // assemble edge hash-table
421     mesh.CalcSurfacesOfNode();
422 
423     for (int k = 1; k <= mesh.GetNFD(); k++)
424       {
425 	multithread.percent = 100.0 * k / (mesh.GetNFD()+1e-10);
426 
427 	if (masterface.Get(k) != k)
428 	  continue;
429 
430 	FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
431 
432 	(*testout) << "Surface " << k << endl;
433 	(*testout) << "Face Descriptor: " << fd << endl;
434 	PrintMessage (1, "Surface ", k, " / ", mesh.GetNFD());
435 
436 	int oldnf = mesh.GetNSE();
437 
438 	const Surface * surf =
439 	  geom.GetSurface((mesh.GetFaceDescriptor(k).SurfNr()));
440 
441 
442 	Meshing2Surfaces meshing(geom, *surf, mparam, geom.BoundingBox());
443 	meshing.SetStartTime (starttime);
444 
445         double eps = 1e-8 * geom.MaxSize();
446 	for (PointIndex pi = PointIndex::BASE; pi < noldp+PointIndex::BASE; pi++)
447 	  {
448 	    // if(surf->PointOnSurface(mesh[pi]))
449 	    meshing.AddPoint (mesh[pi], pi, NULL,
450 			      (surf->PointOnSurface(mesh[pi], eps) != 0));
451 	  }
452 
453 	segments.SetSize (0);
454 
455 	for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
456 	  if (mesh[si].si == k)
457 	    {
458 	      segments.Append (mesh[si]);
459 	      (*testout) << "appending segment " << mesh[si] << endl;
460 	      //<< " from " << mesh[mesh[si][0]]
461 	      //	 << " to " <<mesh[mesh[si][1]]<< endl;
462 	    }
463 
464 	(*testout) << "num-segments " << segments.Size() << endl;
465 
466 	for (int i = 1; i <= geom.identifications.Size(); i++)
467 	  {
468 	    geom.identifications.Get(i)->
469 	      BuildSurfaceElements(segments, mesh, surf);
470 	  }
471 
472 	for (int si = 0; si < segments.Size(); si++)
473 	  {
474 	    PointGeomInfo gi;
475 	    gi.trignum = k;
476 	    meshing.AddBoundaryElement (segments[si][0] + 1 - PointIndex::BASE,
477 					segments[si][1] + 1 - PointIndex::BASE,
478 					gi, gi);
479 	  }
480 
481 	double maxh = mparam.maxh;
482 	if (fd.DomainIn() != 0)
483 	  {
484 	    const Solid * s1 =
485 	      geom.GetTopLevelObject(fd.DomainIn()-1) -> GetSolid();
486 	    if (s1->GetMaxH() < maxh)
487 	      maxh = s1->GetMaxH();
488 	    maxh = min2(maxh, geom.GetTopLevelObject(fd.DomainIn()-1)->GetMaxH());
489 	  }
490 	if (fd.DomainOut() != 0)
491 	  {
492 	    const Solid * s1 =
493 	      geom.GetTopLevelObject(fd.DomainOut()-1) -> GetSolid();
494 	    if (s1->GetMaxH() < maxh)
495 	      maxh = s1->GetMaxH();
496 	    maxh = min2(maxh, geom.GetTopLevelObject(fd.DomainOut()-1)->GetMaxH());
497 	  }
498 	if (fd.TLOSurface() != 0)
499 	  {
500 	    double hi = geom.GetTopLevelObject(fd.TLOSurface()-1) -> GetMaxH();
501 	    if (hi < maxh) maxh = hi;
502 	  }
503 
504 	(*testout) << "domin = " << fd.DomainIn() << ", domout = " << fd.DomainOut()
505 		   << ", tlo-surf = " << fd.TLOSurface()
506 		   << " mpram.maxh = " << mparam.maxh << ", maxh = " << maxh << endl;
507 
508 	mparam.checkoverlap = 0;
509 
510 	MESHING2_RESULT res =
511 	  meshing.GenerateMesh (mesh, mparam, maxh, k);
512 
513 	if (res != MESHING2_OK)
514 	  {
515 	    PrintError ("Problem in Surface mesh generation");
516 	    throw NgException ("Problem in Surface mesh generation");
517 	  }
518 
519 	if (multithread.terminate) return;
520 
521 	for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
522 	  mesh[sei].SetIndex (k);
523 
524         auto n_illegal_trigs = mesh.FindIllegalTrigs();
525         PrintMessage (3, n_illegal_trigs, " illegal triangles");
526 
527 	//      mesh.CalcSurfacesOfNode();
528 
529 	if (segments.Size())
530 	  {
531 	    // surface was meshed, not copied
532 
533 	    static int timer = NgProfiler::CreateTimer ("total surface mesh optimization");
534 	    NgProfiler::RegionTimer reg (timer);
535 
536 
537 	    PrintMessage (2, "Optimize Surface");
538 	    for (int i = 1; i <= mparam.optsteps2d; i++)
539 	      {
540 		if (multithread.terminate) return;
541 
542 		{
543 		  MeshOptimize2d meshopt(mesh);
544 		  meshopt.SetFaceIndex (k);
545 		  meshopt.SetImproveEdges (0);
546 		  meshopt.SetMetricWeight (mparam.elsizeweight);
547 		  meshopt.SetWriteStatus (0);
548 
549 		  meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
550 		}
551 
552 		if (multithread.terminate) return;
553 		{
554 		  //		mesh.CalcSurfacesOfNode();
555 
556 		  MeshOptimize2d meshopt(mesh);
557 		  meshopt.SetFaceIndex (k);
558 		  meshopt.SetImproveEdges (0);
559 		  meshopt.SetMetricWeight (mparam.elsizeweight);
560 		  meshopt.SetWriteStatus (0);
561 
562 		  meshopt.ImproveMesh(mparam);
563 		}
564 
565 		{
566 		  MeshOptimize2d meshopt(mesh);
567 		  meshopt.SetFaceIndex (k);
568 		  meshopt.SetImproveEdges (0);
569 		  meshopt.SetMetricWeight (mparam.elsizeweight);
570 		  meshopt.SetWriteStatus (0);
571 
572 		  meshopt.CombineImprove();
573 		  //		mesh.CalcSurfacesOfNode();
574 		}
575 
576 		if (multithread.terminate) return;
577 		{
578 		  MeshOptimize2d meshopt(mesh);
579 		  meshopt.SetFaceIndex (k);
580 		  meshopt.SetImproveEdges (0);
581 		  meshopt.SetMetricWeight (mparam.elsizeweight);
582 		  meshopt.SetWriteStatus (0);
583 
584 		  meshopt.ImproveMesh(mparam);
585 		}
586 	      }
587 	  }
588 
589 
590 	PrintMessage (3, (mesh.GetNSE() - oldnf), " elements, ", mesh.GetNP(), " points");
591 
592 	mparam.Render();
593       }
594 
595     mesh.Compress();
596 
597     do
598       {
599 	changed = 0;
600 	for (int k = 1; k <= mesh.GetNFD(); k++)
601 	  {
602 	    multithread.percent = 100.0 * k / (mesh.GetNFD()+1e-10);
603 
604 	    if (masterface.Get(k) == k)
605 	      continue;
606 
607 	    FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
608 
609 	    (*testout) << "Surface " << k << endl;
610 	    (*testout) << "Face Descriptor: " << fd << endl;
611 	    PrintMessage (2, "Surface ", k);
612 
613 	    int oldnf = mesh.GetNSE();
614 
615 	    const Surface * surf =
616 	      geom.GetSurface((mesh.GetFaceDescriptor(k).SurfNr()));
617 
618 	    /*
619 	      if (surf -> GetBCProperty() != -1)
620 	      fd.SetBCProperty (surf->GetBCProperty());
621 	      else
622 	      {
623 	      bccnt++;
624 	      fd.SetBCProperty (bccnt);
625 	      }
626 	    */
627 
628 	    segments.SetSize (0);
629 	    for (int i = 1; i <= mesh.GetNSeg(); i++)
630 	      {
631 		Segment * seg = &mesh.LineSegment(i);
632 		if (seg->si == k)
633 		  segments.Append (*seg);
634 	      }
635 
636 	    for (int i = 1; i <= geom.identifications.Size(); i++)
637 	      {
638 		geom.identifications.Elem(i)->GetIdentifiedFaces (fpairs);
639 		int found = 0;
640 		for (int j = 1; j <= fpairs.Size(); j++)
641 		  if (fpairs.Get(j).I1() == k || fpairs.Get(j).I2() == k)
642 		    found = 1;
643 
644 		if (!found)
645 		  continue;
646 
647 		geom.identifications.Get(i)->
648 		  BuildSurfaceElements(segments, mesh, surf);
649 		if (!segments.Size())
650 		  break;
651 	      }
652 
653 
654 	    if (multithread.terminate) return;
655 
656 	    for (SurfaceElementIndex  sei = oldnf; sei < mesh.GetNSE(); sei++)
657 	      mesh[sei].SetIndex (k);
658 
659 
660 	    if (!segments.Size())
661 	      {
662 		masterface.Elem(k) = k;
663 		changed = 1;
664 	      }
665 
666 	    PrintMessage (3, (mesh.GetNSE() - oldnf), " elements, ", mesh.GetNP(), " points");
667 	  }
668 
669         mparam.Render();
670       }
671     while (changed);
672 
673 
674     mesh.SplitSeparatedFaces();
675     mesh.CalcSurfacesOfNode();
676 
677     multithread.task = savetask;
678   }
679 
680 
681 
CSGGenerateMesh(CSGeometry & geom,shared_ptr<Mesh> & mesh,MeshingParameters & mparam)682   int CSGGenerateMesh (CSGeometry & geom,
683 		       shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
684   {
685     NgArray<SpecialPoint> specpoints;
686     NgArray<MeshPoint> spoints;
687 
688 
689     if (mesh && mesh->GetNSE() &&
690 	!geom.GetNSolids())
691       {
692 	if (mparam.perfstepsstart < MESHCONST_MESHVOLUME)
693 	  mparam.perfstepsstart = MESHCONST_MESHVOLUME;
694       }
695 
696     if (mparam.perfstepsstart <= MESHCONST_ANALYSE)
697       {
698         if (mesh)
699           mesh -> DeleteMesh();
700         else
701           mesh = make_shared<Mesh>();
702 
703 	mesh->SetGlobalH (mparam.maxh);
704 	mesh->SetMinimalH (mparam.minh);
705 
706 	NgArray<double> maxhdom(geom.GetNTopLevelObjects());
707 	for (int i = 0; i < maxhdom.Size(); i++)
708 	  maxhdom[i] = geom.GetTopLevelObject(i)->GetMaxH();
709 
710 	mesh->SetMaxHDomain (maxhdom);
711 
712 	if (mparam.uselocalh)
713 	  {
714 	    double maxsize = geom.MaxSize();
715 	    mesh->SetLocalH (Point<3>(-maxsize, -maxsize, -maxsize),
716 			     Point<3>(maxsize, maxsize, maxsize),
717 			     mparam.grading);
718 
719 	    mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
720             for (auto mspnt : mparam.meshsize_points)
721               mesh -> RestrictLocalH (mspnt.pnt, mspnt.h);
722 	  }
723 
724 	spoints.SetSize(0);
725 	FindPoints (geom, specpoints, spoints, *mesh);
726 
727 	PrintMessage (5, "find points done");
728 
729 #ifdef LOG_STREAM
730 	(*logout) << "Special points found" << endl
731 		  << "time = " << GetTime() << " sec" << endl
732 		  << "points: " << mesh->GetNP() << endl << endl;
733 #endif
734       }
735 
736 
737     if (multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE)
738       return TCL_OK;
739 
740 
741     if (mparam.perfstepsstart <= MESHCONST_MESHEDGES)
742       {
743 	FindEdges (geom, *mesh, specpoints, spoints, mparam, true);
744 	if (multithread.terminate) return TCL_OK;
745 #ifdef LOG_STREAM
746 	(*logout) << "Edges meshed" << endl
747 		  << "time = " << GetTime() << " sec" << endl
748 		  << "points: " << mesh->GetNP() << endl;
749 #endif
750 
751 
752 	if (multithread.terminate)
753 	  return TCL_OK;
754 
755 	if (mparam.uselocalh)
756 	  {
757 	    mesh->CalcLocalH(mparam.grading);
758 	    mesh->DeleteMesh();
759 
760 	    FindPoints (geom, specpoints, spoints, *mesh);
761 	    if (multithread.terminate) return TCL_OK;
762 	    FindEdges (geom, *mesh, specpoints, spoints, mparam, true);
763 	    if (multithread.terminate) return TCL_OK;
764 
765 	    mesh->DeleteMesh();
766 
767 	    FindPoints (geom, specpoints, spoints, *mesh);
768 	    if (multithread.terminate) return TCL_OK;
769 	    FindEdges (geom, *mesh, specpoints, spoints, mparam);
770 	    if (multithread.terminate) return TCL_OK;
771 	  }
772       }
773 
774     if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
775       return TCL_OK;
776 
777 
778     if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
779       {
780 	MeshSurface (geom, *mesh, mparam);
781 	if (multithread.terminate) return TCL_OK;
782 
783 #ifdef LOG_STREAM
784 	(*logout) << "Surfaces meshed" << endl
785 		  << "time = " << GetTime() << " sec" << endl
786 		  << "points: " << mesh->GetNP() << endl;
787 #endif
788 
789         /*
790 	if (mparam.uselocalh)
791 	  {
792 	    mesh->CalcLocalH(mparam.grading);
793 	    mesh->DeleteMesh();
794 
795 	    FindPoints (geom, *mesh);
796 	    if (multithread.terminate) return TCL_OK;
797 	    FindEdges (geom, *mesh, mparam);
798 	    if (multithread.terminate) return TCL_OK;
799 
800 	    MeshSurface (geom, *mesh, mparam);
801 	    if (multithread.terminate) return TCL_OK;
802 	  }
803         */
804 
805 #ifdef LOG_STREAM
806 	(*logout) << "Surfaces remeshed" << endl
807 		  << "time = " << GetTime() << " sec" << endl
808 		  << "points: " << mesh->GetNP() << endl;
809 #endif
810 
811 #ifdef STAT_STREAM
812 	(*statout) << mesh->GetNSeg() << " & "
813 		   << mesh->GetNSE() << " & - &"
814 		   << GetTime() << " & " << endl;
815 #endif
816 
817 	MeshQuality2d (*mesh);
818 	mesh->CalcSurfacesOfNode();
819       }
820 
821     if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
822       return TCL_OK;
823 
824 
825     if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
826       {
827 	multithread.task = "Volume meshing";
828 
829 	MESHING3_RESULT res =
830 	  MeshVolume (mparam, *mesh);
831 
832 	if (res != MESHING3_OK) return TCL_ERROR;
833 
834 	if (multithread.terminate) return TCL_OK;
835 
836 	RemoveIllegalElements (*mesh);
837 	if (multithread.terminate) return TCL_OK;
838 
839 	MeshQuality3d (*mesh);
840 
841 	for (int i = 0; i < geom.GetNTopLevelObjects(); i++)
842 	  mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str());
843 
844 
845 #ifdef STAT_STREAM
846 	(*statout) << GetTime() << " & ";
847 #endif
848 
849 #ifdef LOG_STREAM
850 	(*logout) << "Volume meshed" << endl
851 		  << "time = " << GetTime() << " sec" << endl
852 		  << "points: " << mesh->GetNP() << endl;
853 #endif
854       }
855 
856     if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME)
857       return TCL_OK;
858 
859 
860     if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME)
861       {
862 	multithread.task = "Volume optimization";
863 
864 	OptimizeVolume (mparam, *mesh);
865 	if (multithread.terminate) return TCL_OK;
866 
867 #ifdef STAT_STREAM
868 	(*statout) << GetTime() << " & "
869 		   << mesh->GetNE() << " & "
870 		   << mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl;
871 #endif
872 
873 #ifdef LOG_STREAM
874 	(*logout) << "Volume optimized" << endl
875 		  << "time = " << GetTime() << " sec" << endl
876 		  << "points: " << mesh->GetNP() << endl;
877 #endif
878       }
879 
880     mesh -> OrderElements();
881     return TCL_OK;
882   }
883 }
884