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