1 #include <mystdlib.h>
2 #include "meshing.hpp"
3 
4 namespace netgen
5 {
6 
7 double minother;
8 double minwithoutother;
9 
10 
11 
12 
13 
MeshingStat3d()14 MeshingStat3d :: MeshingStat3d ()
15 {
16   cntsucc = cnttrials = cntelem = qualclass = 0;
17   vol0 = h = 1;
18   problemindex = 1;
19 }
20 
21 
Meshing3(const string & rulefilename)22 Meshing3 :: Meshing3 (const string & rulefilename)
23 {
24   tolfak = 1;
25 
26   LoadRules (rulefilename.c_str(), NULL);
27   adfront = new AdFront3;
28 
29   problems.SetSize (rules.Size());
30   foundmap.SetSize (rules.Size());
31   canuse.SetSize (rules.Size());
32   ruleused.SetSize (rules.Size());
33 
34   for (int i = 1; i <= rules.Size(); i++)
35     {
36       problems.Elem(i) = new char[255];
37       foundmap.Elem(i) = 0;
38       canuse.Elem(i) = 0;
39       ruleused.Elem(i) = 0;
40     }
41 }
42 
43 
Meshing3(const char ** rulep)44 Meshing3 :: Meshing3 (const char ** rulep)
45 {
46   tolfak = 1;
47 
48   LoadRules (NULL, rulep);
49   adfront = new AdFront3;
50 
51   problems.SetSize (rules.Size());
52   foundmap.SetSize (rules.Size());
53   canuse.SetSize (rules.Size());
54   ruleused.SetSize (rules.Size());
55 
56   for (int i = 0; i < rules.Size(); i++)
57     {
58       problems[i] = new char[255];
59       foundmap[i] = 0;
60       canuse[i]   = 0;
61       ruleused[i] = 0;
62     }
63 }
64 
~Meshing3()65 Meshing3 :: ~Meshing3 ()
66 {
67   delete adfront;
68   for (int i = 0; i < rules.Size(); i++)
69     {
70       delete [] problems[i];
71       delete rules[i];
72     }
73 }
74 
75 
76 
77 /*
78   // was war das ????
79 static double CalcLocH (const NgArray<Point3d> & locpoints,
80 			const NgArray<MiniElement2d> & locfaces,
81 			double h)
82 {
83   return h;
84 
85 
86   int i, j;
87   double hi, h1, d, dn, sum, weight, wi;
88   Point3d p0, pc;
89   Vec3d n, v1, v2;
90 
91   p0.X() = p0.Y() = p0.Z() = 0;
92   for (j = 1; j <= 3; j++)
93     {
94       p0.X() += locpoints.Get(locfaces.Get(1).PNum(j)).X();
95       p0.Y() += locpoints.Get(locfaces.Get(1).PNum(j)).Y();
96       p0.Z() += locpoints.Get(locfaces.Get(1).PNum(j)).Z();
97     }
98   p0.X() /= 3; p0.Y() /= 3; p0.Z() /= 3;
99 
100   v1 = locpoints.Get(locfaces.Get(1).PNum(2)) -
101     locpoints.Get(locfaces.Get(1).PNum(1));
102   v2 = locpoints.Get(locfaces.Get(1).PNum(3)) -
103     locpoints.Get(locfaces.Get(1).PNum(1));
104 
105   h1 = v1.Length();
106   n = Cross (v2, v1);
107   n /= n.Length();
108 
109   sum = 0;
110   weight = 0;
111 
112   for(int i = 1; i <= locfaces.Size(); i++)
113     {
114       pc.X() = pc.Y() = pc.Z() = 0;
115       for (j = 1; j <= 3; j++)
116 	{
117 	  pc.X() += locpoints.Get(locfaces.Get(i).PNum(j)).X();
118 	  pc.Y() += locpoints.Get(locfaces.Get(i).PNum(j)).Y();
119 	  pc.Z() += locpoints.Get(locfaces.Get(i).PNum(j)).Z();
120 	}
121       pc.X() /= 3; pc.Y() /= 3; pc.Z() /= 3;
122 
123       d = Dist (p0, pc);
124       dn = n * (pc - p0);
125       hi = Dist (locpoints.Get(locfaces.Get(i).PNum(1)),
126 		 locpoints.Get(locfaces.Get(i).PNum(2)));
127 
128       if (dn > -0.2 * h1)
129 	{
130 	  wi = 1 / (h1 + d);
131 	  wi *= wi;
132 	}
133       else
134 	wi = 0;
135 
136       sum += hi * wi;
137       weight += wi;
138     }
139 
140   return sum/weight;
141 }
142 */
143 
AddPoint(const Point3d & p,PointIndex globind)144 PointIndex Meshing3 :: AddPoint (const Point3d & p, PointIndex globind)
145 {
146   return adfront -> AddPoint (p, globind);
147 }
148 
AddBoundaryElement(const Element2d & elem)149 void Meshing3 :: AddBoundaryElement (const Element2d & elem)
150 {
151   MiniElement2d mini(elem.GetNP());
152   for (int j = 0; j < elem.GetNP(); j++)
153     mini[j] = elem[j];
154   adfront -> AddFace(mini);
155 }
156 
157 
AddBoundaryElement(const MiniElement2d & elem)158 void Meshing3 :: AddBoundaryElement (const MiniElement2d & elem)
159 {
160   adfront -> AddFace(elem);
161 }
162 
AddConnectedPair(const INDEX_2 & apair)163 int Meshing3 :: AddConnectedPair (const INDEX_2 & apair)
164 {
165   return adfront -> AddConnectedPair (apair);
166 }
167 
168 MESHING3_RESULT Meshing3 ::
GenerateMesh(Mesh & mesh,const MeshingParameters & mp)169 GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
170 {
171   static Timer t("Meshing3::GenerateMesh"); RegionTimer reg(t);
172   // static Timer meshing3_timer_a("Meshing3::GenerateMesh a", 2);
173   // static Timer meshing3_timer_b("Meshing3::GenerateMesh b", 2);
174   // static Timer meshing3_timer_c("Meshing3::GenerateMesh c", 1);
175   // static Timer meshing3_timer_d("Meshing3::GenerateMesh d", 2);
176   // static int meshing3_timer = NgProfiler::CreateTimer ("Meshing3::GenerateMesh");
177   // static int meshing3_timer_a = NgProfiler::CreateTimer ("Meshing3::GenerateMesh a");
178   // static int meshing3_timer_b = NgProfiler::CreateTimer ("Meshing3::GenerateMesh b");
179   // static int meshing3_timer_c = NgProfiler::CreateTimer ("Meshing3::GenerateMesh c");
180   // static int meshing3_timer_d = NgProfiler::CreateTimer ("Meshing3::GenerateMesh d");
181   // NgProfiler::RegionTimer reg (meshing3_timer);
182 
183 
184   NgArray<Point3d, PointIndex::BASE> locpoints;      // local points
185   NgArray<MiniElement2d> locfaces;                   // local faces
186   NgArray<PointIndex, PointIndex::BASE> pindex;      // mapping from local to front point numbering
187   NgArray<int, PointIndex::BASE> allowpoint;         // point is allowed ?
188   NgArray<INDEX> findex;                             // mapping from local to front face numbering
189   //INDEX_2_HASHTABLE<int> connectedpairs(100);    // connecgted pairs for prism meshing
190 
191   NgArray<Point3d, PointIndex::BASE> plainpoints;    // points in reference coordinates
192   NgArray<int> delpoints, delfaces;   // points and lines to be deleted
193   NgArray<Element> locelements;       // new generated elements
194 
195   int j, oldnp, oldnf;
196   int found;
197   referencetransform trans;
198   int rotind;
199   Point3d inp;
200   float err;
201 
202   INDEX locfacesplit;             //index for faces in outer area
203 
204   bool loktestmode = false;
205 
206   int uselocalh = mp.uselocalh;
207 
208   // int giveuptol = mp.giveuptol; //
209   MeshingStat3d stat;      // statistics
210   int plotstat_oldne = -1;
211 
212 
213   // for star-shaped domain meshing
214   NgArray<MeshPoint, PointIndex::BASE> grouppoints;
215   NgArray<MiniElement2d> groupfaces;
216   NgArray<PointIndex, PointIndex::BASE> grouppindex;
217   NgArray<INDEX> groupfindex;
218 
219 
220   float minerr;
221   int hasfound;
222   double tetvol;
223   // int giveup = 0;
224 
225 
226   NgArray<Point3d> tempnewpoints;
227   NgArray<MiniElement2d> tempnewfaces;
228   NgArray<int> tempdelfaces;
229   NgArray<Element> templocelements;
230 
231 
232   stat.h = mp.maxh;
233 
234   adfront->SetStartFront (mp.baseelnp);
235 
236 
237   found = 0;
238   stat.vol0 = adfront -> Volume();
239   tetvol = 0;
240 
241   stat.qualclass = 1;
242 
243   while (1)
244     {
245       if (multithread.terminate)
246 	throw NgException ("Meshing stopped");
247 
248       // break if advancing front is empty
249       if (!mp.baseelnp && adfront->Empty())
250 	break;
251 
252       // break, if advancing front has no elements with
253       // mp.baseelnp nodes
254       if (mp.baseelnp && adfront->Empty (mp.baseelnp))
255 	break;
256 
257       locpoints.SetSize(0);
258       locfaces.SetSize(0);
259       locelements.SetSize(0);
260       pindex.SetSize(0);
261       findex.SetSize(0);
262 
263       INDEX_2_HASHTABLE<int> connectedpairs(100);  // connected pairs for prism meshing
264 
265       // select base-element (will be locface[1])
266       // and get local environment of radius (safety * h)
267 
268 
269       int baseelem = adfront -> SelectBaseElement ();
270       if (mp.baseelnp && adfront->GetFace (baseelem).GetNP() != mp.baseelnp)
271 	{
272 	  adfront->IncrementClass (baseelem);
273 	  continue;
274 	}
275 
276       const MiniElement2d & bel = adfront->GetFace (baseelem);
277       const Point<3> p1 = adfront->GetPoint (bel[0]);
278       const Point<3> p2 = adfront->GetPoint (bel[1]);
279       const Point<3> p3 = adfront->GetPoint (bel[2]);
280 
281 
282       Point<3> pmid = Center (p1, p2, p3);
283 
284       double his = (Dist (p1, p2) + Dist(p1, p3) + Dist(p2, p3)) / 3;
285       double hshould = mesh.GetH (pmid);
286 
287       if (adfront->GetFace (baseelem).GetNP() == 4)
288 	hshould = max2 (his, hshould);
289 
290       double hmax = (his > hshould) ? his : hshould;
291 
292       // qualclass should come from baseelem !!!!!
293       double hinner = hmax * (1 + stat.qualclass);
294       double houter = hmax * (1 + 2 * stat.qualclass);
295 
296       // meshing3_timer_a.Start();
297       stat.qualclass =
298         adfront -> GetLocals (baseelem, locpoints, locfaces,
299 			      pindex, findex, connectedpairs,
300 			      houter, hinner,
301 			      locfacesplit);
302       // meshing3_timer_a.Stop();
303 
304       // (*testout) << "locfaces = " << endl << locfaces << endl;
305 
306 
307       //loktestmode = 1;
308       testmode = loktestmode;  //changed
309       // loktestmode = testmode =  (adfront->GetFace (baseelem).GetNP() == 4) && (rules.Size() == 5);
310 
311       loktestmode = stat.qualclass > 5;
312 
313 
314       if (loktestmode)
315 	{
316 	  (*testout) << "baseel = " << baseelem << ", ind = " << findex.Get(1) << endl;
317           int pi1 = pindex[locfaces[0].PNum(1)];
318           int pi2 = pindex[locfaces[0].PNum(2)];
319           int pi3 = pindex[locfaces[0].PNum(3)];
320 	  (*testout) << "pi = " << pi1 << ", " << pi2 << ", " << pi3 << endl;
321 	}
322 
323 
324       if (testmode)
325 	{
326 	  (*testout) << "baseelem = " << baseelem << " qualclass = " << stat.qualclass << endl;
327 	  (*testout) << "locpoints = " << endl << locpoints << endl;
328 	  (*testout) << "connected = " << endl << connectedpairs << endl;
329 	}
330 
331 
332 
333       // loch = CalcLocH (locpoints, locfaces, h);
334 
335       stat.nff = adfront->GetNF();
336       stat.vol = adfront->Volume();
337       if (stat.vol < 0) break;
338 
339       oldnp = locpoints.Size();
340       oldnf = locfaces.Size();
341 
342 
343       allowpoint.SetSize(locpoints.Size());
344       if (uselocalh && stat.qualclass <= 3)
345 	for(int i = 1; i <= allowpoint.Size(); i++)
346 	  {
347 	    allowpoint.Elem(i) =
348 	      (mesh.GetH (locpoints.Get(i)) > 0.4 * hshould / mp.sloppy) ? 2 : 1;
349 	  }
350       else
351 	allowpoint = 2;
352 
353 
354 
355       if (stat.qualclass >= mp.starshapeclass &&
356 	  mp.baseelnp != 4)
357 	{
358 	  // NgProfiler::RegionTimer reg1 (meshing3_timer_b);
359 	  // star-shaped domain removing
360 
361 	  grouppoints.SetSize (0);
362 	  groupfaces.SetSize (0);
363 	  grouppindex.SetSize (0);
364 	  groupfindex.SetSize (0);
365 
366 	  adfront -> GetGroup (findex[0], grouppoints, groupfaces,
367 			       grouppindex, groupfindex);
368 
369 	  bool onlytri = 1;
370 	  for (auto i : groupfaces.Range())
371 	    if (groupfaces[i].GetNP() != 3)
372 	      onlytri = 0;
373 
374 	  if (onlytri && groupfaces.Size() <= 20 + 2*stat.qualclass &&
375 	      FindInnerPoint (grouppoints, groupfaces, inp))
376 	    {
377 	      (*testout) << "inner point found" << endl;
378 
379 	      for(int i = 1; i <= groupfaces.Size(); i++)
380 		adfront -> DeleteFace (groupfindex.Get(i));
381 
382 	      for(int i = 1; i <= groupfaces.Size(); i++)
383 		for (j = 1; j <= locfaces.Size(); j++)
384 		  if (findex.Get(j) == groupfindex.Get(i))
385 		    delfaces.Append (j);
386 
387 
388 	      delfaces.SetSize (0);
389 
390 	      INDEX npi;
391 	      Element newel(TET);
392 
393 	      npi = mesh.AddPoint (inp);
394 	      newel.SetNP(4);
395 	      newel.PNum(4) = npi;
396 
397 	      for(int i = 1; i <= groupfaces.Size(); i++)
398 		{
399 		  for (j = 1; j <= 3; j++)
400 		    {
401 		      newel.PNum(j) =
402 			adfront->GetGlobalIndex
403 			(grouppindex.Get(groupfaces.Get(i).PNum(j)));
404 		    }
405 		  mesh.AddVolumeElement (newel);
406 		}
407 	      continue;
408 	    }
409 	}
410 
411       found = 0;
412       hasfound = 0;
413       minerr = 1e6;
414 
415       //      int optother = 0;
416 
417       /*
418       for(int i = 1; i <= locfaces.Size(); i++)
419 	{
420 	  (*testout) << "Face " << i << ": ";
421 	  for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
422 	    (*testout) << pindex.Get(locfaces.Get(i).PNum(j)) << " ";
423 	  (*testout) << endl;
424 	}
425       for(int i = 1; i <= locpoints.Size(); i++)
426 	{
427 	  (*testout) << "p" << i
428 		     << ", gi = " << pindex.Get(i)
429 		     << " = " << locpoints.Get(i) << endl;
430 	}
431 	*/
432 
433       minother = 1e10;
434       minwithoutother = 1e10;
435 
436       bool impossible = 1;
437 
438       for (rotind = 1; rotind <= locfaces[0].GetNP(); rotind++)
439 	{
440 	  // set transformatino to reference coordinates
441 
442 	  if (locfaces[0].GetNP() == 3)
443 	    {
444 	      trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)],
445 			 locpoints[locfaces[0].PNumMod(2+rotind)],
446 			 locpoints[locfaces[0].PNumMod(3+rotind)], hshould);
447 	    }
448 	  else
449 	    {
450 	      trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)],
451 			 locpoints[locfaces[0].PNumMod(2+rotind)],
452 			 locpoints[locfaces[0].PNumMod(4+rotind)], hshould);
453 	    }
454 
455 	  // trans.ToPlain (locpoints, plainpoints);
456 
457           plainpoints.SetSize (locpoints.Size());
458           for (auto i : locpoints.Range())
459             trans.ToPlain (locpoints[i], plainpoints[i]);
460 
461           for (auto i : allowpoint.Range())
462             if (plainpoints[i].Z() > 0)
463               allowpoint[i] = false;
464 
465 	  stat.cnttrials++;
466 
467 
468 	  if (stat.cnttrials % 100 == 0)
469 	    {
470 	      (*testout) << "\n";
471 	      for(int i = 1; i <= canuse.Size(); i++)
472 	      {
473 		(*testout) << foundmap.Get(i) << "/"
474 			   << canuse.Get(i) << "/"
475 			   << ruleused.Get(i) << " map/can/use rule " << rules.Get(i)->Name() << "\n";
476 	      }
477 	      (*testout) << endl;
478 	    }
479 
480 	  // NgProfiler::StartTimer (meshing3_timer_c);
481           // meshing3_timer_c.Start();
482 
483 	  found = ApplyRules (plainpoints, allowpoint,
484 			      locfaces, locfacesplit, connectedpairs,
485 			      locelements, delfaces,
486 			      stat.qualclass, mp.sloppy, rotind, err);
487 
488 	  if (found >= 0) impossible = 0;
489 	  if (found < 0) found = 0;
490 
491           // meshing3_timer_c.Stop();
492 	  // NgProfiler::StopTimer (meshing3_timer_c);
493 
494 	  if (!found) loktestmode = 0;
495 
496 	  // NgProfiler::RegionTimer reg2 (meshing3_timer_d);
497 
498 	  if (loktestmode)
499 	    {
500 	      (*testout) << "plainpoints = " << endl << plainpoints << endl;
501 	      (*testout) << "Applyrules found " << found << endl;
502 	    }
503 
504 	  if (found) stat.cntsucc++;
505 
506 	  locpoints.SetSize (plainpoints.Size());
507 	  for (int i = oldnp+1; i <= plainpoints.Size(); i++)
508 	    trans.FromPlain (plainpoints.Elem(i), locpoints.Elem(i));
509 
510 
511 
512 	  // avoid meshing from large to small mesh-size
513 	  if (uselocalh && found && stat.qualclass <= 3)
514 	    {
515 	      for (int i = 1; i <= locelements.Size(); i++)
516 		{
517 		  Point3d pmin = locpoints[locelements.Get(i).PNum(1)];
518 		  Point3d pmax = pmin;
519 		  for (j = 2; j <= 4; j++)
520 		    {
521 		      const Point3d & hp = locpoints[locelements.Get(i).PNum(j)];
522 		      pmin.SetToMin (hp);
523 		      pmax.SetToMax (hp);
524 		    }
525 
526 		  if (mesh.GetMinH (pmin, pmax) < 0.4 * hshould / mp.sloppy)
527 		    found = 0;
528 		}
529 	    }
530 	  if (found)
531 	    {
532 	      for (int i = 1; i <= locelements.Size(); i++)
533 		for (int j = 1; j <= 4; j++)
534 		  {
535 		    const Point3d & hp = locpoints[locelements.Get(i).PNum(j)];
536 		    if (Dist (hp, pmid) > hinner)
537 		      found = 0;
538 		  }
539 	    }
540 
541 
542 	  if (found)
543 	    ruleused.Elem(found)++;
544 
545 
546 	  // plotstat->Plot(stat);
547 
548 	  if (stat.cntelem != plotstat_oldne)
549 	    {
550 	      plotstat_oldne = stat.cntelem;
551 
552 	      PrintMessageCR (5, "El: ", stat.cntelem,
553 			      //	    << " trials: " << stat.cnttrials
554 			      " faces: ", stat.nff,
555 			      " vol = ", float(100 * stat.vol / stat.vol0));
556 
557 	      multithread.percent = 100 -  100.0 * stat.vol / stat.vol0;
558 	    }
559 
560 
561 	  if (found && (!hasfound || err < minerr) )
562 	    {
563 
564 	      if (testmode)
565 		{
566 		  (*testout) << "found is active, 3" << endl;
567 		  for(int i = 1; i <= plainpoints.Size(); i++)
568 		    {
569 		      (*testout) << "p";
570 		      if (i <= pindex.Size())
571 			(*testout) << pindex.Get(i) << ": ";
572 		      else
573 			(*testout) << "new: ";
574 		      (*testout) << plainpoints.Get(i) << endl;
575 		    }
576 		}
577 
578 
579 
580 	      hasfound = found;
581 	      minerr = err;
582 
583 	      tempnewpoints.SetSize (0);
584 	      for(int i = oldnp+1; i <= locpoints.Size(); i++)
585 		tempnewpoints.Append (locpoints.Get(i));
586 
587 	      tempnewfaces.SetSize (0);
588 	      for(int i = oldnf+1; i <= locfaces.Size(); i++)
589 		tempnewfaces.Append (locfaces.Get(i));
590 
591 	      tempdelfaces.SetSize (0);
592 	      for(int i = 1; i <= delfaces.Size(); i++)
593 		tempdelfaces.Append (delfaces.Get(i));
594 
595 	      templocelements.SetSize (0);
596 	      for(int i = 1; i <= locelements.Size(); i++)
597 		templocelements.Append (locelements.Get(i));
598 
599 	      /*
600 	      optother =
601 		strcmp (problems[found], "other") == 0;
602 	      */
603 	    }
604 
605 	  locpoints.SetSize (oldnp);
606 	  locfaces.SetSize (oldnf);
607 	  delfaces.SetSize (0);
608 	  locelements.SetSize (0);
609 	}
610 
611 
612 
613       if (hasfound)
614 	{
615 
616 	  /*
617 	  if (optother)
618 	    (*testout) << "Other is optimal" << endl;
619 
620 	  if (minother < minwithoutother)
621 	    {
622 	      (*testout) << "Other is better, " << minother << " less " << minwithoutother << endl;
623 	    }
624 	    */
625 
626 	  for(int i = 1; i <= tempnewpoints.Size(); i++)
627 	    locpoints.Append (tempnewpoints.Get(i));
628 	  for(int i = 1; i <= tempnewfaces.Size(); i++)
629 	    locfaces.Append (tempnewfaces.Get(i));
630 	  for(int i = 1; i <= tempdelfaces.Size(); i++)
631 	    delfaces.Append (tempdelfaces.Get(i));
632 	  for(int i = 1; i <= templocelements.Size(); i++)
633 	    locelements.Append (templocelements.Get(i));
634 
635 
636 	  if (loktestmode)
637 	    {
638 	      (*testout) << "apply rule" << endl;
639 	      for(int i = 1; i <= locpoints.Size(); i++)
640 		{
641 		  (*testout) << "p";
642 		  if (i <= pindex.Size())
643 		    (*testout) << pindex.Get(i) << ": ";
644 		  else
645 		    (*testout) << "new: ";
646 		  (*testout) << locpoints.Get(i) << endl;
647 		}
648 	    }
649 
650 
651 
652 	  pindex.SetSize(locpoints.Size());
653 
654 	  for (int i = oldnp+1; i <= locpoints.Size(); i++)
655 	    {
656 	      PointIndex globind = mesh.AddPoint (locpoints.Get(i));
657 	      pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind);
658 	    }
659 
660 	  for (int i = 1; i <= locelements.Size(); i++)
661 	    {
662 	      Point3d * hp1, * hp2, * hp3, * hp4;
663 	      hp1 = &locpoints[locelements.Get(i).PNum(1)];
664 	      hp2 = &locpoints[locelements.Get(i).PNum(2)];
665 	      hp3 = &locpoints[locelements.Get(i).PNum(3)];
666 	      hp4 = &locpoints[locelements.Get(i).PNum(4)];
667 
668 	      tetvol += (1.0 / 6.0) * ( Cross ( *hp2 - *hp1, *hp3 - *hp1) * (*hp4 - *hp1) );
669 
670 	      for (j = 1; j <= locelements.Get(i).NP(); j++)
671 		locelements.Elem(i).PNum(j) =
672 		  adfront -> GetGlobalIndex (pindex[locelements.Get(i).PNum(j)]);
673 
674 	      mesh.AddVolumeElement (locelements.Get(i));
675 	      stat.cntelem++;
676 	    }
677 
678 	  for(int i = oldnf+1; i <= locfaces.Size(); i++)
679 	    {
680 	      for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
681 		locfaces.Elem(i).PNum(j) =
682 		  pindex[locfaces.Get(i).PNum(j)];
683 	      // (*testout) << "add face " << locfaces.Get(i) << endl;
684 	      adfront->AddFace (locfaces.Get(i));
685 	    }
686 
687 	  for(int i = 1; i <= delfaces.Size(); i++)
688 	    adfront->DeleteFace (findex.Get(delfaces.Get(i)));
689 	}
690       else
691 	{
692 	  adfront->IncrementClass (findex.Get(1));
693 	  if (impossible && mp.check_impossible)
694 	    {
695 	      (*testout) << "skip face since it is impossible" << endl;
696 	      for (j = 0; j < 100; j++)
697 		adfront->IncrementClass (findex.Get(1));
698 	    }
699 	}
700 
701       locelements.SetSize (0);
702       delpoints.SetSize(0);
703       delfaces.SetSize(0);
704 
705       if (stat.qualclass >= mp.giveuptol)
706 	break;
707     }
708 
709   PrintMessage (5, "");  // line feed after statistics
710 
711   for(int i = 1; i <= ruleused.Size(); i++)
712     (*testout) << setw(4) << ruleused.Get(i)
713 	       << " times used rule " << rules.Get(i) -> Name() << endl;
714 
715 
716   if (!mp.baseelnp && adfront->Empty())
717     return MESHING3_OK;
718 
719   if (mp.baseelnp && adfront->Empty (mp.baseelnp))
720     return MESHING3_OK;
721 
722   if (stat.vol < -1e-15)
723     return MESHING3_NEGVOL;
724 
725   return MESHING3_NEGVOL;
726 }
727 
728 
729 
730 
731 enum blocktyp { BLOCKUNDEF, BLOCKINNER, BLOCKBOUND, BLOCKOUTER };
732 
BlockFill(Mesh & mesh,double gh)733 void Meshing3 :: BlockFill (Mesh & mesh, double gh)
734 {
735   static Timer t("Mesing3::BlockFill"); RegionTimer reg(t);
736 
737   PrintMessage (3, "Block-filling called (obsolete) ");
738 
739   int i, j(0), i1, i2, i3, j1, j2, j3;
740   int n1, n2, n3, n, min1, min2, min3, max1, max2, max3;
741   int changed, filled;
742   double xmin(0), xmax(0), ymin(0), ymax(0), zmin(0), zmax(0);
743   double xminb, xmaxb, yminb, ymaxb, zminb, zmaxb;
744   //double rad = 0.7 * gh;
745 
746   for(int i = 1; i <= adfront->GetNP(); i++)
747     {
748       const Point3d & p = adfront->GetPoint(PointIndex(i));
749       if (i == 1)
750 	{
751 	  xmin = xmax = p.X();
752 	  ymin = ymax = p.Y();
753 	  zmin = zmax = p.Z();
754 	}
755       else
756 	{
757 	  if (p.X() < xmin) xmin = p.X();
758 	  if (p.X() > xmax) xmax = p.X();
759 	  if (p.Y() < ymin) ymin = p.Y();
760 	  if (p.Y() > ymax) ymax = p.Y();
761 	  if (p.Z() < zmin) zmin = p.Z();
762 	  if (p.Z() > zmax) zmax = p.Z();
763 	}
764     }
765 
766   xmin -= 5 * gh;
767   ymin -= 5 * gh;
768   zmin -= 5 * gh;
769 
770   n1 = int ((xmax-xmin) / gh + 5);
771   n2 = int ((ymax-ymin) / gh + 5);
772   n3 = int ((zmax-zmin) / gh + 5);
773   n = n1 * n2 * n3;
774 
775   PrintMessage (5, "n1 = ", n1, " n2 = ", n2, " n3 = ", n3);
776 
777   NgArray<blocktyp> inner(n);
778   NgArray<PointIndex> pointnr(n);
779   NgArray<int> frontpointnr(n);
780 
781 
782   // initialize inner to 1
783 
784   for(int i = 1; i <= n; i++)
785     inner.Elem(i) = BLOCKUNDEF;
786 
787 
788   // set blocks cutting surfaces to 0
789 
790   for(int i = 1; i <= adfront->GetNF(); i++)
791     {
792       const MiniElement2d & el = adfront->GetFace(i);
793       xminb = xmax; xmaxb = xmin;
794       yminb = ymax; ymaxb = ymin;
795       zminb = zmax; zmaxb = zmin;
796 
797       for (j = 1; j <= 3; j++)
798 	{
799 	  const Point3d & p = adfront->GetPoint (el.PNum(j));
800 	  if (p.X() < xminb) xminb = p.X();
801 	  if (p.X() > xmaxb) xmaxb = p.X();
802 	  if (p.Y() < yminb) yminb = p.Y();
803 	  if (p.Y() > ymaxb) ymaxb = p.Y();
804 	  if (p.Z() < zminb) zminb = p.Z();
805 	  if (p.Z() > zmaxb) zmaxb = p.Z();
806 	}
807 
808 
809 
810       double filldist = 0.2; // globflags.GetNumFlag ("filldist", 0.4);
811       xminb -= filldist * gh;
812       xmaxb += filldist * gh;
813       yminb -= filldist * gh;
814       ymaxb += filldist * gh;
815       zminb -= filldist * gh;
816       zmaxb += filldist * gh;
817 
818       min1 = int ((xminb - xmin) / gh) + 1;
819       max1 = int ((xmaxb - xmin) / gh) + 1;
820       min2 = int ((yminb - ymin) / gh) + 1;
821       max2 = int ((ymaxb - ymin) / gh) + 1;
822       min3 = int ((zminb - zmin) / gh) + 1;
823       max3 = int ((zmaxb - zmin) / gh) + 1;
824 
825 
826       for (i1 = min1; i1 <= max1; i1++)
827 	for (i2 = min2; i2 <= max2; i2++)
828 	  for (i3 = min3; i3 <= max3; i3++)
829 	    inner.Elem(i3 + (i2-1) * n3 + (i1-1) * n2 * n3) = BLOCKBOUND;
830     }
831 
832 
833 
834 
835   while (1)
836     {
837       int undefi = 0;
838       Point3d undefp;
839 
840       for (i1 = 1; i1 <= n1 && !undefi; i1++)
841 	for (i2 = 1; i2 <= n2 && !undefi; i2++)
842 	  for (i3 = 1; i3 <= n3 && !undefi; i3++)
843 	    {
844 	      i = i3 + (i2-1) * n3 + (i1-1) * n2 * n3;
845 	      if (inner.Elem(i) == BLOCKUNDEF)
846 		{
847 		  undefi = i;
848 		  undefp.X() = xmin + (i1-0.5) * gh;
849 		  undefp.Y() = ymin + (i2-0.5) * gh;
850 		  undefp.Z() = zmin + (i3-0.5) * gh;
851 		}
852 	    }
853 
854       if (!undefi)
855 	break;
856 
857       //      PrintMessage (5, "Test point: ", undefp);
858 
859       if (adfront -> Inside (undefp))
860 	{
861 	  //	  (*mycout) << "inner" << endl;
862 	  inner.Elem(undefi) = BLOCKINNER;
863 	}
864       else
865 	{
866 	  //	  (*mycout) << "outer" << endl;
867 	  inner.Elem(undefi) = BLOCKOUTER;
868 	}
869 
870       do
871 	{
872 	  changed = 0;
873 	  for (i1 = 1; i1 <= n1; i1++)
874 	    for (i2 = 1; i2 <= n2; i2++)
875 	      for (i3 = 1; i3 <= n3; i3++)
876 		{
877 		  i = i3 + (i2-1) * n3 + (i1-1) * n2 * n3;
878 
879 		  for (int k = 1; k <= 3; k++)
880 		    {
881 		      switch (k)
882 			{
883 			case 1: j = i + n2 * n3; break;
884 			case 2: j = i + n3; break;
885 			case 3: j = i + 1; break;
886 			}
887 
888 		      if (j > n1 * n2 * n3) continue;
889 
890 		      if (inner.Elem(i) == BLOCKOUTER && inner.Elem(j) == BLOCKUNDEF)
891 			{
892 			  changed = 1;
893 			  inner.Elem(j) = BLOCKOUTER;
894 			}
895 		      if (inner.Elem(j) == BLOCKOUTER && inner.Elem(i) == BLOCKUNDEF)
896 			{
897 			  changed = 1;
898 			  inner.Elem(i) = BLOCKOUTER;
899 			}
900 		      if (inner.Elem(i) == BLOCKINNER && inner.Elem(j) == BLOCKUNDEF)
901 			{
902 			  changed = 1;
903 			  inner.Elem(j) = BLOCKINNER;
904 			}
905 		      if (inner.Elem(j) == BLOCKINNER && inner.Elem(i) == BLOCKUNDEF)
906 			{
907 			  changed = 1;
908 			  inner.Elem(i) = BLOCKINNER;
909 			}
910 		    }
911 		}
912 	}
913       while (changed);
914 
915     }
916 
917 
918 
919   filled = 0;
920   for(int i = 1; i <= n; i++)
921     if (inner.Elem(i) == BLOCKINNER)
922       {
923 	filled++;
924       }
925   PrintMessage (5, "Filled blocks: ", filled);
926 
927   for(int i = 1; i <= n; i++)
928     {
929       pointnr.Elem(i) = 0;
930       frontpointnr.Elem(i) = 0;
931     }
932 
933   for (i1 = 1; i1 <= n1-1; i1++)
934     for (i2 = 1; i2 <= n2-1; i2++)
935       for (i3 = 1; i3 <= n3-1; i3++)
936 	{
937 	  i = i3 + (i2-1) * n3 + (i1-1) * n2 * n3;
938 	  if (inner.Elem(i) == BLOCKINNER)
939 	    {
940 	      for (j1 = i1; j1 <= i1+1; j1++)
941 		for (j2 = i2; j2 <= i2+1; j2++)
942 		  for (j3 = i3; j3 <= i3+1; j3++)
943 		    {
944 		      j = j3 + (j2-1) * n3 + (j1-1) * n2 * n3;
945 		      if (pointnr.Get(j) == 0)
946 			{
947 			  Point3d hp(xmin + (j1-1) * gh,
948 				     ymin + (j2-1) * gh,
949 				     zmin + (j3-1) * gh);
950 			  pointnr.Elem(j) = mesh.AddPoint (hp);
951 			  frontpointnr.Elem(j) =
952 			    AddPoint (hp, pointnr.Elem(j));
953 
954 			}
955 		    }
956 	    }
957 	}
958 
959 
960   for (i1 = 2; i1 <= n1-1; i1++)
961     for (i2 = 2; i2 <= n2-1; i2++)
962       for (i3 = 2; i3 <= n3-1; i3++)
963 	{
964 	  i = i3 + (i2-1) * n3 + (i1-1) * n2 * n3;
965 	  if (inner.Elem(i) == BLOCKINNER)
966 	    {
967 	      int pn[9];
968 	      pn[1] = pointnr.Get(i);
969 	      pn[2] = pointnr.Get(i+1);
970 	      pn[3] = pointnr.Get(i+n3);
971 	      pn[4] = pointnr.Get(i+n3+1);
972 	      pn[5] = pointnr.Get(i+n2*n3);
973 	      pn[6] = pointnr.Get(i+n2*n3+1);
974 	      pn[7] = pointnr.Get(i+n2*n3+n3);
975 	      pn[8] = pointnr.Get(i+n2*n3+n3+1);
976 	      static int elind[][4] =
977 	      {
978 		{ 1, 8, 2, 4 },
979 		{ 1, 8, 4, 3 },
980 		{ 1, 8, 3, 7 },
981 		{ 1, 8, 7, 5 },
982 		{ 1, 8, 5, 6 },
983 		{ 1, 8, 6, 2 }
984 	      };
985 	      for (j = 1; j <= 6; j++)
986 		{
987 		  Element el(4);
988 		  for (int k = 1; k <= 4;  k++)
989 		    el.PNum(k) = pn[elind[j-1][k-1]];
990 
991 		  mesh.AddVolumeElement (el);
992 		}
993 	    }
994 	}
995 
996 
997 
998   for (i1 = 2; i1 <= n1-1; i1++)
999     for (i2 = 2; i2 <= n2-1; i2++)
1000       for (i3 = 2; i3 <= n3-1; i3++)
1001 	{
1002 	  i = i3 + (i2-1) * n3 + (i1-1) * n2 * n3;
1003 	  if (inner.Elem(i) == BLOCKINNER)
1004 	    {
1005 	      int pi1(0), pi2(0), pi3(0), pi4(0);
1006 
1007 	      int pn1 = frontpointnr.Get(i);
1008 	      int pn2 = frontpointnr.Get(i+1);
1009 	      int pn3 = frontpointnr.Get(i+n3);
1010 	      int pn4 = frontpointnr.Get(i+n3+1);
1011 	      int pn5 = frontpointnr.Get(i+n2*n3);
1012 	      int pn6 = frontpointnr.Get(i+n2*n3+1);
1013 	      int pn7 = frontpointnr.Get(i+n2*n3+n3);
1014 	      int pn8 = frontpointnr.Get(i+n2*n3+n3+1);
1015 
1016 	      for (int k = 1; k <= 6; k++)
1017 		{
1018 		  switch (k)
1019 		    {
1020 		    case 1: // j3 = i3+1
1021 		      j = i + 1;
1022 		      pi1 = pn2;
1023 		      pi2 = pn6;
1024 		      pi3 = pn4;
1025 		      pi4 = pn8;
1026 		      break;
1027 		    case 2: // j3 = i3-1
1028 		      j = i - 1;
1029 		      pi1 = pn1;
1030 		      pi2 = pn3;
1031 		      pi3 = pn5;
1032 		      pi4 = pn7;
1033 		      break;
1034 		    case 3: // j2 = i2+1
1035 		      j = i + n3;
1036 		      pi1 = pn3;
1037 		      pi2 = pn4;
1038 		      pi3 = pn7;
1039 		      pi4 = pn8;
1040 		      break;
1041 		    case 4: // j2 = i2-1
1042 		      j = i - n3;
1043 		      pi1 = pn1;
1044 		      pi2 = pn5;
1045 		      pi3 = pn2;
1046 		      pi4 = pn6;
1047 		      break;
1048 		    case 5: // j1 = i1+1
1049 		      j = i + n3*n2;
1050 		      pi1 = pn5;
1051 		      pi2 = pn7;
1052 		      pi3 = pn6;
1053 		      pi4 = pn8;
1054 		      break;
1055 		    case 6: // j1 = i1-1
1056 		      j = i - n3*n2;
1057 		      pi1 = pn1;
1058 		      pi2 = pn2;
1059 		      pi3 = pn3;
1060 		      pi4 = pn4;
1061 		      break;
1062 		    }
1063 
1064 		  if (inner.Get(j) == BLOCKBOUND)
1065 		    {
1066 		      MiniElement2d face;
1067 		      face.PNum(1) = pi4;
1068 		      face.PNum(2) = pi1;
1069 		      face.PNum(3) = pi3;
1070 		      AddBoundaryElement (face);
1071 
1072 		      face.PNum(1) = pi1;
1073 		      face.PNum(2) = pi4;
1074 		      face.PNum(3) = pi2;
1075 		      AddBoundaryElement (face);
1076 
1077 		    }
1078 		}
1079 	    }
1080 	}
1081 }
1082 
1083 
1084 /*
1085 static const AdFront3 * locadfront;
1086 static int TestInner (const Point3d & p)
1087 {
1088   return locadfront->Inside (p);
1089 }
1090 static int TestSameSide (const Point3d & p1, const Point3d & p2)
1091 {
1092   return locadfront->SameSide (p1, p2);
1093 }
1094 */
1095 
1096 
1097 
BlockFillLocalH(Mesh & mesh,const MeshingParameters & mp)1098 void Meshing3 :: BlockFillLocalH (Mesh & mesh,
1099 				  const MeshingParameters & mp)
1100 {
1101   static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
1102 
1103   double filldist = mp.filldist;
1104 
1105   // (*testout) << "blockfill local h" << endl;
1106   // (*testout) << "rel filldist = " << filldist << endl;
1107   PrintMessage (3, "blockfill local h");
1108 
1109 
1110   NgArray<Point<3> > npoints;
1111 
1112   adfront -> CreateTrees();
1113 
1114   Box<3> bbox ( Box<3>::EMPTY_BOX );
1115   double maxh = 0;
1116 
1117   for (int i = 1; i <= adfront->GetNF(); i++)
1118     {
1119       const MiniElement2d & el = adfront->GetFace(i);
1120       for (int j = 1; j <= 3; j++)
1121 	{
1122 	  const Point3d & p1 = adfront->GetPoint (el.PNumMod(j));
1123 	  const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1));
1124 
1125 	  double hi = Dist (p1, p2);
1126 	  if (hi > maxh) maxh = hi;
1127 
1128 	  bbox.Add (p1);
1129 	}
1130     }
1131 
1132 
1133   Point3d mpmin = bbox.PMin();
1134   Point3d mpmax = bbox.PMax();
1135   Point3d mpc = Center (mpmin, mpmax);
1136   double d = max3(mpmax.X()-mpmin.X(),
1137 		  mpmax.Y()-mpmin.Y(),
1138 		  mpmax.Z()-mpmin.Z()) / 2;
1139   mpmin = mpc - Vec3d (d, d, d);
1140   mpmax = mpc + Vec3d (d, d, d);
1141   Box3d meshbox (mpmin, mpmax);
1142 
1143   LocalH loch2 (mpmin, mpmax, 1);
1144 
1145   if (mp.maxh < maxh) maxh = mp.maxh;
1146 
1147   auto loch_ptr = mesh.LocalHFunction().Copy(bbox);
1148   auto & loch = *loch_ptr;
1149 
1150   bool changed;
1151   static Timer t1("loop1");
1152   t1.Start();
1153   do
1154     {
1155       loch.ClearFlags();
1156 
1157       static Timer tbox("adfront-bbox");
1158       tbox.Start();
1159       for (int i = 1; i <= adfront->GetNF(); i++)
1160 	{
1161 	  const MiniElement2d & el = adfront->GetFace(i);
1162 
1163 	  Box<3> bbox (adfront->GetPoint (el[0]));
1164 	  bbox.Add (adfront->GetPoint (el[1]));
1165 	  bbox.Add (adfront->GetPoint (el[2]));
1166 
1167 
1168 	  double filld = filldist * bbox.Diam();
1169 	  bbox.Increase (filld);
1170 
1171       	  loch.CutBoundary (bbox); // .PMin(), bbox.PMax());
1172 	}
1173       tbox.Stop();
1174 
1175       //      locadfront = adfront;
1176       loch.FindInnerBoxes (adfront, NULL);
1177 
1178       npoints.SetSize(0);
1179       loch.GetInnerPoints (npoints);
1180 
1181       changed = false;
1182       for (int i = 1; i <= npoints.Size(); i++)
1183 	{
1184 	  if (loch.GetH(npoints.Get(i)) > 1.5 * maxh)
1185 	    {
1186 	      loch.SetH (npoints.Get(i), maxh);
1187 	      changed = true;
1188 	    }
1189 	}
1190     }
1191   while (changed);
1192   t1.Stop();
1193 
1194   if (debugparam.slowchecks)
1195     (*testout) << "Blockfill with points: " << endl;
1196   for (int i = 1; i <= npoints.Size(); i++)
1197     {
1198       if (meshbox.IsIn (npoints.Get(i)))
1199 	{
1200 	  PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
1201 	  adfront->AddPoint (npoints.Get(i), gpnum);
1202 
1203 	  if (debugparam.slowchecks)
1204 	    {
1205 	      (*testout) << npoints.Get(i) << endl;
1206 	      if (!adfront->Inside(npoints.Get(i)))
1207 		{
1208 		  cout << "add outside point" << endl;
1209 		  (*testout) << "outside" << endl;
1210 		}
1211 	    }
1212 
1213 	}
1214     }
1215 
1216 
1217 
1218   // find outer points
1219 
1220   static Timer tloch2("build loch2");
1221   tloch2.Start();
1222   loch2.ClearFlags();
1223 
1224   for (int i = 1; i <= adfront->GetNF(); i++)
1225     {
1226       const MiniElement2d & el = adfront->GetFace(i);
1227       Point3d pmin = adfront->GetPoint (el.PNum(1));
1228       Point3d pmax = pmin;
1229 
1230       for (int j = 2; j <= 3; j++)
1231 	{
1232 	  const Point3d & p = adfront->GetPoint (el.PNum(j));
1233 	  pmin.SetToMin (p);
1234 	  pmax.SetToMax (p);
1235 	}
1236 
1237       loch2.SetH (Center (pmin, pmax), Dist (pmin, pmax));
1238     }
1239 
1240   for (int i = 1; i <= adfront->GetNF(); i++)
1241     {
1242       const MiniElement2d & el = adfront->GetFace(i);
1243       Point3d pmin = adfront->GetPoint (el.PNum(1));
1244       Point3d pmax = pmin;
1245 
1246       for (int j = 2; j <= 3; j++)
1247 	{
1248 	  const Point3d & p = adfront->GetPoint (el.PNum(j));
1249 	  pmin.SetToMin (p);
1250 	  pmax.SetToMax (p);
1251 	}
1252 
1253       double filld = filldist * Dist (pmin, pmax);
1254       pmin = pmin - Vec3d (filld, filld, filld);
1255       pmax = pmax + Vec3d (filld, filld, filld);
1256       // loch2.CutBoundary (pmin, pmax);
1257       loch2.CutBoundary (Box<3> (pmin, pmax)); // pmin, pmax);
1258     }
1259   tloch2.Stop();
1260 
1261   // locadfront = adfront;
1262   loch2.FindInnerBoxes (adfront, NULL);
1263 
1264   npoints.SetSize(0);
1265   loch2.GetOuterPoints (npoints);
1266 
1267   for (int i = 1; i <= npoints.Size(); i++)
1268     {
1269       if (meshbox.IsIn (npoints.Get(i)))
1270 	{
1271 	  PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
1272 	  adfront->AddPoint (npoints.Get(i), gpnum);
1273 	}
1274     }
1275 }
1276 
1277 }
1278