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