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