1 #include <mystdlib.h>
2 #include "meshing.hpp"
3
4
5 namespace netgen
6 {
7 extern double minother;
8 extern double minwithoutother;
9
10
CalcElementBadness(const Array<Point3d> & points,const Element & elem)11 static double CalcElementBadness (const Array<Point3d> & points,
12 const Element & elem)
13 {
14 double vol, l, l4, l5, l6;
15 if (elem.GetNP() != 4)
16 {
17 if (elem.GetNP() == 5)
18 {
19 double z = points.Get(elem.PNum(5)).Z();
20 if (z > -1e-8) return 1e8;
21 return (-1 / z) - z; // - 2;
22 }
23 return 0;
24 }
25
26 Vec3d v1 = points.Get(elem.PNum(2)) - points.Get(elem.PNum(1));
27 Vec3d v2 = points.Get(elem.PNum(3)) - points.Get(elem.PNum(1));
28 Vec3d v3 = points.Get(elem.PNum(4)) - points.Get(elem.PNum(1));
29
30 vol = - (Cross (v1, v2) * v3);
31 l4 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(3)));
32 l5 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(4)));
33 l6 = Dist (points.Get(elem.PNum(3)), points.Get(elem.PNum(4)));
34
35 l = v1.Length() + v2.Length() + v3.Length() + l4 + l5 + l6;
36
37 // testout << "vol = " << vol << " l = " << l << endl;
38 if (vol < 1e-8) return 1e10;
39 // (*testout) << "l^3/vol = " << (l*l*l / vol) << endl;
40
41 double err = pow (l*l*l/vol, 1.0/3.0) / 12;
42 return err;
43 }
44
45
46
47
48
49
ApplyRules(Array<Point3d> & lpoints,Array<int> & allowpoint,Array<MiniElement2d> & lfaces,INDEX lfacesplit,INDEX_2_HASHTABLE<int> & connectedpairs,Array<Element> & elements,Array<INDEX> & delfaces,int tolerance,double sloppy,int rotind1,float & retminerr)50 int Meshing3 :: ApplyRules
51 (
52 Array<Point3d> & lpoints, // in: local points, out: old+new local points
53 Array<int> & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means
54 Array<MiniElement2d> & lfaces, // in: local faces, out: old+new local faces
55 INDEX lfacesplit, // for local faces in outer radius
56 INDEX_2_HASHTABLE<int> & connectedpairs, // connected pairs for prism-meshing
57 Array<Element> & elements, // out: new elements
58 Array<INDEX> & delfaces, // out: face indices of faces to delete
59 int tolerance, // quality class: 1 best
60 double sloppy, // quality strength
61 int rotind1, // how to rotate base element
62 float & retminerr // element error
63 )
64
65 {
66 NgProfiler::RegionTimer regtot(97);
67
68 int i, j, k, ri, nfok, npok, incnpok, refpi, locpi, locfi, locfr;
69 float hf, err, minerr, teterr, minteterr;
70 char ok, found, hc;
71 vnetrule * rule;
72 Vector oldu, newu, newu1, newu2, allp;
73 Vec3d ui;
74 Point3d np;
75 int oldnp, noldlp, noldlf;
76 const MiniElement2d * locface = NULL;
77 int loktestmode;
78
79
80 Array<int> pused; // point is already mapped
81 Array<char> fused; // face is already mapped
82 Array<int> pmap; // map of reference point to local point
83 Array<char> pfixed; // point mapped by face-map
84 Array<int> fmapi; // face in reference is mapped to face nr ...
85 Array<int> fmapr; // face in reference is rotated to map
86 Array<Point3d> transfreezone; // transformed free-zone
87 INDEX_2_CLOSED_HASHTABLE<int> ledges(100); // edges in local environment
88
89 Array<Point3d> tempnewpoints;
90 Array<MiniElement2d> tempnewfaces;
91 Array<int> tempdelfaces;
92 Array<Element> tempelements;
93 Array<Box3d> triboxes; // bounding boxes of local faces
94
95 Array<int, PointIndex::BASE> pnearness;
96 Array<int> fnearness;
97
98 static int cnt = 0;
99 cnt++;
100
101 delfaces.SetSize (0);
102 elements.SetSize (0);
103
104 // determine topological distance of faces and points to
105 // base element
106
107 pnearness.SetSize (lpoints.Size());
108 fnearness.SetSize (lfacesplit);
109
110 pnearness = INT_MAX/10;
111 for (j = 0; j < lfaces[0].GetNP(); j++)
112 pnearness[lfaces[0][j]] = 0;
113
114 NgProfiler::RegionTimer reg2(98);
115
116 NgProfiler::StartTimer (90);
117
118 for (int loop = 0; loop < 2; loop++)
119 {
120
121 for (i = 0; i < lfacesplit; i++)
122 {
123 const MiniElement2d & hface = lfaces[i];
124
125 int minn = INT_MAX-1;
126 for (j = 0; j < hface.GetNP(); j++)
127 {
128 int hi = pnearness[hface[j]];
129 if (hi < minn) minn = hi;
130 }
131 if (minn < INT_MAX/10)
132 for (j = 0; j < hface.GetNP(); j++)
133 if (pnearness[hface[j]] > minn+1)
134 pnearness[hface[j]] = minn+1;
135 }
136
137 for (i = 1; i <= connectedpairs.GetNBags(); i++)
138 for (j = 1; j <= connectedpairs.GetBagSize(i); j++)
139 {
140 INDEX_2 edge;
141 int val;
142 connectedpairs.GetData (i, j, edge, val);
143
144 if (pnearness[edge.I1()] > pnearness[edge.I2()] + 1)
145 pnearness[edge.I1()] = pnearness[edge.I2()] + 1;
146
147 if (pnearness[edge.I2()] > pnearness[edge.I1()] + 1)
148 pnearness[edge.I2()] = pnearness[edge.I1()] + 1;
149 }
150
151 }
152
153 for (i = 0; i < fnearness.Size(); i++)
154 {
155 int sum = 0;
156 for (j = 0; j < lfaces[i].GetNP(); j++)
157 sum += pnearness[lfaces[i][j]];
158 fnearness[i] = sum;
159 }
160
161
162 NgProfiler::StopTimer (90);
163 NgProfiler::StartTimer (91);
164
165 // find bounding boxes of faces
166
167 triboxes.SetSize (lfaces.Size());
168 for (i = 0; i < lfaces.Size(); i++)
169 {
170 const MiniElement2d & face = lfaces[i];
171 triboxes[i].SetPoint (lpoints.Get(face[0]));
172 for (j = 1; j < face.GetNP(); j++)
173 triboxes[i].AddPoint (lpoints.Get(face[j]));
174 }
175
176 NgProfiler::StopTimer (91);
177 NgProfiler::StartTimer (92);
178
179
180 bool useedges = 0;
181 for (ri = 0; ri < rules.Size(); ri++)
182 if (rules[ri]->GetNEd()) useedges = 1;
183
184 if (useedges)
185 {
186 ledges.SetSize (5 * lfacesplit);
187
188 for (j = 0; j < lfacesplit; j++)
189 // if (fnearness[j] <= 5)
190 {
191 const MiniElement2d & face = lfaces[j];
192 int newp, oldp;
193
194 newp = face[face.GetNP()-1];
195 for (k = 0; k < face.GetNP(); k++)
196 {
197 oldp = newp;
198 newp = face[k];
199 ledges.Set (INDEX_2::Sort(oldp, newp), 1);
200 }
201 }
202 }
203
204 NgProfiler::StopTimer (92);
205
206 NgProfiler::RegionTimer reg3(99);
207
208 pused.SetSize (lpoints.Size());
209 fused.SetSize (lfaces.Size());
210
211 found = 0;
212 minerr = tolfak * tolerance * tolerance;
213 minteterr = sloppy * tolerance;
214
215 if (testmode)
216 (*testout) << "cnt = " << cnt << " class = " << tolerance << endl;
217
218
219
220 // impossible, if no rule can be applied at any tolerance class
221 bool impossible = 1;
222
223
224 // check each rule:
225
226 for (ri = 1; ri <= rules.Size(); ri++)
227 {
228 int base = (lfaces[0].GetNP() == 3) ? 100 : 200;
229 NgProfiler::RegionTimer regx1(base);
230 NgProfiler::RegionTimer regx(base+ri);
231
232 sprintf (problems.Elem(ri), "");
233
234 rule = rules.Get(ri);
235
236 if (rule->GetNP(1) != lfaces[0].GetNP())
237 continue;
238
239 if (rule->GetQuality() > tolerance)
240 {
241 if (rule->GetQuality() < 100) impossible = 0;
242
243 if (testmode)
244 sprintf (problems.Elem(ri), "Quality not ok");
245 continue;
246 }
247
248 if (testmode)
249 sprintf (problems.Elem(ri), "no mapping found");
250
251 loktestmode = testmode || rule->TestFlag ('t') || tolerance > 5;
252
253 if (loktestmode)
254 (*testout) << "Rule " << ri << " = " << rule->Name() << endl;
255
256 pmap.SetSize (rule->GetNP());
257 fmapi.SetSize (rule->GetNF());
258 fmapr.SetSize (rule->GetNF());
259
260 fused = 0;
261 pused = 0;
262 pmap = 0;
263 fmapi = 0;
264 for (i = 1; i <= fmapr.Size(); i++)
265 fmapr.Set(i, rule->GetNP(i));
266
267 fused[0] = 1;
268 fmapi[0] = 1;
269 fmapr[0] = rotind1;
270
271
272 for (j = 1; j <= lfaces.Get(1).GetNP(); j++)
273 {
274 locpi = lfaces[0].PNumMod (j+rotind1);
275 pmap.Set (rule->GetPointNr (1, j), locpi);
276 pused.Elem(locpi)++;
277 }
278
279 /*
280 map all faces
281 nfok .. first nfok-1 faces are mapped properly
282 */
283
284 nfok = 2;
285 NgProfiler::RegionTimer regfa(300);
286 NgProfiler::RegionTimer regx2(base+50+ri);
287 while (nfok >= 2)
288 {
289
290 if (nfok <= rule->GetNOldF())
291 {
292 // not all faces mapped
293
294 ok = 0;
295 locfi = fmapi.Get(nfok);
296 locfr = fmapr.Get(nfok);
297
298 int actfnp = rule->GetNP(nfok);
299
300 while (!ok)
301 {
302 locfr++;
303 if (locfr == actfnp + 1)
304 {
305 locfr = 1;
306 locfi++;
307 if (locfi > lfacesplit) break;
308 }
309
310
311 if (fnearness.Get(locfi) > rule->GetFNearness (nfok) ||
312 fused.Get(locfi) ||
313 actfnp != lfaces.Get(locfi).GetNP() )
314 {
315 // face not feasible in any rotation
316
317 locfr = actfnp;
318 }
319 else
320 {
321
322 ok = 1;
323
324 locface = &lfaces.Get(locfi);
325
326
327 // reference point already mapped differently ?
328 for (j = 1; j <= actfnp && ok; j++)
329 {
330 locpi = pmap.Get(rule->GetPointNr (nfok, j));
331
332 if (locpi && locpi != locface->PNumMod(j+locfr))
333 ok = 0;
334 }
335
336 // local point already used or point outside tolerance ?
337 for (j = 1; j <= actfnp && ok; j++)
338 {
339 refpi = rule->GetPointNr (nfok, j);
340
341 if (pmap.Get(refpi) == 0)
342 {
343 locpi = locface->PNumMod (j + locfr);
344
345 if (pused.Get(locpi))
346 ok = 0;
347 else
348 {
349 const Point3d & lp = lpoints.Get(locpi);
350 const Point3d & rp = rule->GetPoint(refpi);
351
352 if ( Dist2 (lp, rp) * rule->PointDistFactor(refpi) > minerr)
353 {
354 impossible = 0;
355 ok = 0;
356 }
357 }
358 }
359 }
360 }
361 }
362
363
364 if (ok)
365 {
366 // map face nfok
367
368 fmapi.Set (nfok, locfi);
369 fmapr.Set (nfok, locfr);
370 fused.Set (locfi, 1);
371
372 for (j = 1; j <= rule->GetNP (nfok); j++)
373 {
374 locpi = locface->PNumMod(j+locfr);
375
376 if (rule->GetPointNr (nfok, j) <= 3 &&
377 pmap.Get(rule->GetPointNr(nfok, j)) != locpi)
378 (*testout) << "change face1 point, mark1" << endl;
379
380 pmap.Set(rule->GetPointNr (nfok, j), locpi);
381 pused.Elem(locpi)++;
382 }
383
384 nfok++;
385 }
386 else
387 {
388 // backtrack one face
389 fmapi.Set (nfok, 0);
390 fmapr.Set (nfok, rule->GetNP(nfok));
391 nfok--;
392
393 fused.Set (fmapi.Get(nfok), 0);
394 for (j = 1; j <= rule->GetNP (nfok); j++)
395 {
396 refpi = rule->GetPointNr (nfok, j);
397 pused.Elem(pmap.Get(refpi))--;
398
399 if (pused.Get(pmap.Get(refpi)) == 0)
400 {
401 pmap.Set(refpi, 0);
402 }
403 }
404 }
405 }
406
407 else
408
409 {
410 NgProfiler::RegionTimer regfb(301);
411
412 // all faces are mapped
413 // now map all isolated points:
414
415 if (loktestmode)
416 {
417 (*testout) << "Faces Ok" << endl;
418 sprintf (problems.Elem(ri), "Faces Ok");
419 }
420
421 npok = 1;
422 incnpok = 1;
423
424 pfixed.SetSize (pmap.Size());
425 for (i = 1; i <= pmap.Size(); i++)
426 pfixed.Set(i, (pmap.Get(i) != 0) );
427
428 while (npok >= 1)
429 {
430
431 if (npok <= rule->GetNOldP())
432 {
433
434 if (pfixed.Get(npok))
435
436 {
437 if (incnpok)
438 npok++;
439 else
440 npok--;
441 }
442
443 else
444
445 {
446 locpi = pmap.Elem(npok);
447 ok = 0;
448
449 if (locpi)
450 pused.Elem(locpi)--;
451
452 while (!ok && locpi < lpoints.Size())
453 {
454 ok = 1;
455 locpi++;
456
457 if (pused.Get(locpi) ||
458 pnearness.Get(locpi) > rule->GetPNearness(npok))
459 {
460 ok = 0;
461 }
462 else if (allowpoint.Get(locpi) != 2)
463 {
464 ok = 0;
465 if (allowpoint.Get(locpi) == 1)
466 impossible = 0;
467 }
468 else
469 {
470 const Point3d & lp = lpoints.Get(locpi);
471 const Point3d & rp = rule->GetPoint(npok);
472
473 if ( Dist2 (lp, rp) * rule->PointDistFactor(npok) > minerr)
474 {
475 ok = 0;
476 impossible = 0;
477 }
478 }
479 }
480
481
482 if (ok)
483 {
484 pmap.Set (npok, locpi);
485
486 if (npok <= 3)
487 (*testout) << "set face1 point, mark3" << endl;
488
489 pused.Elem(locpi)++;
490 npok++;
491 incnpok = 1;
492 }
493
494 else
495
496 {
497 pmap.Set (npok, 0);
498
499 if (npok <= 3)
500 (*testout) << "set face1 point, mark4" << endl;
501
502 npok--;
503 incnpok = 0;
504 }
505 }
506 }
507
508 else
509
510 {
511 NgProfiler::RegionTimer regfa2(302);
512
513 // all points are mapped
514
515 if (loktestmode)
516 {
517 (*testout) << "Mapping found!!: Rule " << rule->Name() << endl;
518 for (i = 1; i <= pmap.Size(); i++)
519 (*testout) << pmap.Get(i) << " ";
520 (*testout) << endl;
521 sprintf (problems.Elem(ri), "mapping found");
522 (*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl;
523 }
524
525 ok = 1;
526
527
528 // check mapedges:
529 for (i = 1; i <= rule->GetNEd(); i++)
530 {
531 INDEX_2 in2(pmap.Get(rule->GetEdge(i).i1),
532 pmap.Get(rule->GetEdge(i).i2));
533 in2.Sort();
534 if (!ledges.Used (in2)) ok = 0;
535 }
536
537
538 // check prism edges:
539 for (i = 1; i <= rule->GetNE(); i++)
540 {
541 const Element & el = rule->GetElement (i);
542 if (el.GetType() == PRISM)
543 {
544 for (j = 1; j <= 3; j++)
545 {
546 INDEX_2 in2(pmap.Get(el.PNum(j)),
547 pmap.Get(el.PNum(j+3)));
548 in2.Sort();
549 if (!connectedpairs.Used (in2)) ok = 0;
550 }
551 }
552 if (el.GetType() == PYRAMID)
553 {
554 if (loktestmode)
555 (*testout) << "map pyramid, rule = " << rule->Name() << endl;
556 for (j = 1; j <= 2; j++)
557 {
558 INDEX_2 in2;
559 if (j == 1)
560 {
561 in2.I1() = pmap.Get(el.PNum(2));
562 in2.I2() = pmap.Get(el.PNum(3));
563 }
564 else
565 {
566 in2.I1() = pmap.Get(el.PNum(1));
567 in2.I2() = pmap.Get(el.PNum(4));
568 }
569 in2.Sort();
570 if (!connectedpairs.Used (in2))
571 {
572 ok = 0;
573 if (loktestmode)
574 (*testout) << "no pair" << endl;
575 }
576 }
577 }
578
579 }
580
581
582
583 for (i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++)
584 fmapi.Set(i, 0);
585
586
587 if (ok)
588 {
589 foundmap.Elem(ri)++;
590 }
591
592
593
594
595 // deviation of existing points
596
597 oldu.SetSize (3 * rule->GetNOldP());
598 newu.SetSize (3 * (rule->GetNP() - rule->GetNOldP()));
599 allp.SetSize (3 * rule->GetNP());
600
601 for (i = 1; i <= rule->GetNOldP(); i++)
602 {
603 const Point3d & lp = lpoints.Get(pmap.Get(i));
604 const Point3d & rp = rule->GetPoint(i);
605 oldu (3*i-3) = lp.X()-rp.X();
606 oldu (3*i-2) = lp.Y()-rp.Y();
607 oldu (3*i-1) = lp.Z()-rp.Z();
608
609 allp (3*i-3) = lp.X();
610 allp (3*i-2) = lp.Y();
611 allp (3*i-1) = lp.Z();
612 }
613
614 if (rule->GetNP() > rule->GetNOldP())
615 {
616 newu.SetSize (rule->GetOldUToNewU().Height());
617 rule->GetOldUToNewU().Mult (oldu, newu);
618 }
619
620 // int idiff = 3 * (rule->GetNP()-rule->GetNOldP());
621 int idiff = 3 * rule->GetNOldP();
622 for (i = rule->GetNOldP()+1; i <= rule->GetNP(); i++)
623 {
624 const Point3d & rp = rule->GetPoint(i);
625 allp (3*i-3) = rp.X() + newu(3*i-3 - idiff);
626 allp (3*i-2) = rp.Y() + newu(3*i-2 - idiff);
627 allp (3*i-1) = rp.Z() + newu(3*i-1 - idiff);
628 }
629
630 rule->SetFreeZoneTransformation (allp,
631 tolerance + int(sloppy));
632
633 if (!rule->ConvexFreeZone())
634 {
635 ok = 0;
636 sprintf (problems.Elem(ri), "Freezone not convex");
637
638 if (loktestmode)
639 (*testout) << "Freezone not convex" << endl;
640 }
641
642 if (loktestmode)
643 {
644 const Array<Point3d> & fz = rule->GetTransFreeZone();
645 (*testout) << "Freezone: " << endl;
646 for (i = 1; i <= fz.Size(); i++)
647 (*testout) << fz.Get(i) << endl;
648 }
649
650
651 // check freezone:
652
653 for (i = 1; i <= lpoints.Size(); i++)
654 {
655 if ( !pused.Get(i) )
656 {
657 const Point3d & lp = lpoints.Get(i);
658
659 if (rule->fzbox.IsIn (lp))
660 {
661 if (rule->IsInFreeZone(lp))
662 {
663 if (loktestmode)
664 {
665 (*testout) << "Point " << i
666 << " in Freezone" << endl;
667 sprintf (problems.Elem(ri),
668 "locpoint %d in Freezone", i);
669 }
670 ok = 0;
671 break;
672 }
673 }
674 }
675 }
676
677 for (i = 1; i <= lfaces.Size() && ok; i++)
678 {
679 static Array<int> lpi(4);
680
681 if (!fused.Get(i))
682 {
683 int triin;
684 const MiniElement2d & lfacei = lfaces.Get(i);
685
686 if (!triboxes.Elem(i).Intersect (rule->fzbox))
687 triin = 0;
688 else
689 {
690 int li, lj;
691 for (li = 1; li <= lfacei.GetNP(); li++)
692 {
693 int lpii = 0;
694 int pi = lfacei.PNum(li);
695 for (lj = 1; lj <= rule->GetNOldP(); lj++)
696 if (pmap.Get(lj) == pi)
697 lpii = lj;
698 lpi.Elem(li) = lpii;
699 }
700
701
702 if (lfacei.GetNP() == 3)
703 {
704 triin = rule->IsTriangleInFreeZone
705 (
706 lpoints.Get(lfacei.PNum(1)),
707 lpoints.Get(lfacei.PNum(2)),
708 lpoints.Get(lfacei.PNum(3)), lpi, 1
709 );
710 }
711 else
712 {
713 triin = rule->IsQuadInFreeZone
714 (
715 lpoints.Get(lfacei.PNum(1)),
716 lpoints.Get(lfacei.PNum(2)),
717 lpoints.Get(lfacei.PNum(3)),
718 lpoints.Get(lfacei.PNum(4)),
719 lpi, 1
720 );
721 }
722 }
723
724
725 if (triin == -1)
726 {
727 ok = 0;
728 }
729
730 if (triin == 1)
731 {
732 #ifdef TEST_JS
733 ok = 0;
734
735 if (loktestmode)
736 {
737 (*testout) << "El with " << lfaces.Get(i).GetNP() << " points in freezone: "
738 << lfaces.Get(i).PNum(1) << " - "
739 << lfaces.Get(i).PNum(2) << " - "
740 << lfaces.Get(i).PNum(3) << " - "
741 << lfaces.Get(i).PNum(4) << endl;
742 for (int lj = 1; lj <= lfaces.Get(i).GetNP(); lj++)
743 (*testout) << lpoints.Get(lfaces.Get(i).PNum(lj)) << " ";
744
745 (*testout) << endl;
746
747 sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
748 lfaces.Get(i).PNum(1), lfaces.Get(i).PNum(2),
749 lfaces.Get(i).PNum(3));
750 }
751 #else
752 if (loktestmode)
753 {
754 if (lfacei.GetNP() == 3)
755 {
756 (*testout) << "Triangle in freezone: "
757 << lfacei.PNum(1) << " - "
758 << lfacei.PNum(2) << " - "
759 << lfacei.PNum(3)
760 << ", or "
761 << lpoints.Get(lfacei.PNum(1)) << " - "
762 << lpoints.Get(lfacei.PNum(2)) << " - "
763 << lpoints.Get(lfacei.PNum(3))
764 << endl;
765 (*testout) << "lpi = " << lpi.Get(1) << ", "
766 << lpi.Get(2) << ", " << lpi.Get(3) << endl;
767 }
768 else
769 (*testout) << "Quad in freezone: "
770 << lfacei.PNum(1) << " - "
771 << lfacei.PNum(2) << " - "
772 << lfacei.PNum(3) << " - "
773 << lfacei.PNum(4)
774 << ", or "
775 << lpoints.Get(lfacei.PNum(1)) << " - "
776 << lpoints.Get(lfacei.PNum(2)) << " - "
777 << lpoints.Get(lfacei.PNum(3)) << " - "
778 << lpoints.Get(lfacei.PNum(4))
779 << endl;
780
781 sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
782 int(lfaces.Get(i).PNum(1)),
783 int(lfaces.Get(i).PNum(2)),
784 int(lfaces.Get(i).PNum(3)));
785 }
786
787 hc = 0;
788 for (k = rule->GetNOldF() + 1; k <= rule->GetNF(); k++)
789 {
790 if (rule->GetPointNr(k, 1) <= rule->GetNOldP() &&
791 rule->GetPointNr(k, 2) <= rule->GetNOldP() &&
792 rule->GetPointNr(k, 3) <= rule->GetNOldP())
793 {
794 for (j = 1; j <= 3; j++)
795 if (lfaces.Get(i).PNumMod(j ) == pmap.Get(rule->GetPointNr(k, 1)) &&
796 lfaces.Get(i).PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) &&
797 lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2)))
798 {
799 fmapi.Elem(k) = i;
800 hc = 1;
801
802
803 // (*testout) << "found from other side: "
804 // << rule->Name()
805 // << " ( " << pmap.Get (rule->GetPointNr(k, 1))
806 // << " - " << pmap.Get (rule->GetPointNr(k, 2))
807 // << " - " << pmap.Get (rule->GetPointNr(k, 3)) << " ) "
808 // << endl;
809
810 strcpy (problems.Elem(ri), "other");
811 }
812 }
813 }
814
815 if (!hc)
816 {
817 if (loktestmode)
818 {
819 (*testout) << "Triangle in freezone: "
820 << lfaces.Get(i).PNum(1) << " - "
821 << lfaces.Get(i).PNum(2) << " - "
822 << lfaces.Get(i).PNum(3) << endl;
823
824 sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
825 int (lfaces.Get(i).PNum(1)),
826 int (lfaces.Get(i).PNum(2)),
827 int (lfaces.Get(i).PNum(3)));
828 }
829 ok = 0;
830 }
831 #endif
832 }
833 }
834
835 }
836
837
838 if (ok)
839 {
840 err = 0;
841 for (i = 1; i <= rule->GetNOldP(); i++)
842 {
843 hf = rule->CalcPointDist (i, lpoints.Get(pmap.Get(i)));
844 if (hf > err) err = hf;
845 }
846
847
848 if (loktestmode)
849 {
850 (*testout) << "Rule ok" << endl;
851 sprintf (problems.Elem(ri), "Rule ok, err = %f", err);
852 }
853
854
855 // newu = rule->GetOldUToNewU() * oldu;
856
857 // set new points:
858
859 oldnp = rule->GetNOldP();
860 noldlp = lpoints.Size();
861 noldlf = lfaces.Size();
862
863
864 for (i = oldnp + 1; i <= rule->GetNP(); i++)
865 {
866 np = rule->GetPoint(i);
867 np.X() += newu (3 * (i-oldnp) - 3);
868 np.Y() += newu (3 * (i-oldnp) - 2);
869 np.Z() += newu (3 * (i-oldnp) - 1);
870
871 pmap.Elem(i) = lpoints.Append (np);
872 }
873
874 // Set new Faces:
875
876 for (i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++)
877 if (!fmapi.Get(i))
878 {
879 MiniElement2d nface(rule->GetNP(i));
880 for (j = 1; j <= nface.GetNP(); j++)
881 nface.PNum(j) = pmap.Get(rule->GetPointNr (i, j));
882
883 lfaces.Append (nface);
884 }
885
886
887 // Delete old Faces:
888
889 for (i = 1; i <= rule->GetNDelF(); i++)
890 delfaces.Append (fmapi.Get(rule->GetDelFace(i)));
891 for (i = rule->GetNOldF()+1; i <= rule->GetNF(); i++)
892 if (fmapi.Get(i))
893 {
894 delfaces.Append (fmapi.Get(i));
895 fmapi.Elem(i) = 0;
896 }
897
898
899 // check orientation
900 for (i = 1; i <= rule->GetNO() && ok; i++)
901 {
902 const fourint * fouri;
903
904 fouri = &rule->GetOrientation(i);
905 Vec3d v1 (lpoints.Get(pmap.Get(fouri->i1)),
906 lpoints.Get(pmap.Get(fouri->i2)));
907 Vec3d v2 (lpoints.Get(pmap.Get(fouri->i1)),
908 lpoints.Get(pmap.Get(fouri->i3)));
909 Vec3d v3 (lpoints.Get(pmap.Get(fouri->i1)),
910 lpoints.Get(pmap.Get(fouri->i4)));
911
912 Vec3d n;
913 Cross (v1, v2, n);
914 //if (n * v3 >= -1e-7*n.Length()*v3.Length()) // OR -1e-7???
915 if (n * v3 >= -1e-9)
916 {
917 if (loktestmode)
918 {
919 sprintf (problems.Elem(ri), "Orientation wrong");
920 (*testout) << "Orientation wrong ("<< n*v3 << ")" << endl;
921 }
922 ok = 0;
923 }
924 }
925
926
927
928 // new points in free-zone ?
929 for (i = rule->GetNOldP() + 1; i <= rule->GetNP() && ok; i++)
930 if (!rule->IsInFreeZone (lpoints.Get(pmap.Get(i))))
931 {
932 if (loktestmode)
933 {
934 (*testout) << "Newpoint " << lpoints.Get(pmap.Get(i))
935 << " outside convex hull" << endl;
936 sprintf (problems.Elem(ri), "newpoint outside convex hull");
937 }
938 ok = 0;
939
940 }
941
942 // insert new elements
943
944 for (i = 1; i <= rule->GetNE(); i++)
945 {
946 elements.Append (rule->GetElement(i));
947 for (j = 1; j <= elements.Get(i).NP(); j++)
948 elements.Elem(i).PNum(j) = pmap.Get(elements.Get(i).PNum(j));
949 }
950
951
952 // Calculate Element badness
953
954 teterr = 0;
955 for (i = 1; i <= elements.Size(); i++)
956 {
957 hf = CalcElementBadness (lpoints, elements.Get(i));
958 if (hf > teterr) teterr = hf;
959 }
960
961 /*
962 // keine gute Erfahrung am 25.1.2000, js
963 if (ok && teterr < 100 &&
964 (rule->TestFlag('b') || tolerance > 10) )
965 {
966 (*mycout) << "Reset teterr "
967 << rule->Name()
968 << " err = " << teterr
969 << endl;
970 teterr = 1;
971 }
972 */
973
974 // compare edgelength
975 if (rule->TestFlag('l'))
976 {
977 double oldlen = 0;
978 double newlen = 0;
979
980 for (i = 1; i <= rule->GetNDelF(); i++)
981 {
982 const Element2d & face =
983 rule->GetFace (rule->GetDelFace(i));
984 for (j = 1; j <= 3; j++)
985 {
986 const Point3d & p1 =
987 lpoints.Get(pmap.Get(face.PNumMod(j)));
988 const Point3d & p2 =
989 lpoints.Get(pmap.Get(face.PNumMod(j+1)));
990 oldlen += Dist(p1, p2);
991 }
992 }
993
994 for (i = rule->GetNOldF()+1; i <= rule->GetNF(); i++)
995 {
996 const Element2d & face = rule->GetFace (i);
997 for (j = 1; j <= 3; j++)
998 {
999 const Point3d & p1 =
1000 lpoints.Get(pmap.Get(face.PNumMod(j)));
1001 const Point3d & p2 =
1002 lpoints.Get(pmap.Get(face.PNumMod(j+1)));
1003 newlen += Dist(p1, p2);
1004 }
1005 }
1006
1007 if (oldlen < newlen)
1008 {
1009 ok = 0;
1010 if (loktestmode)
1011 sprintf (problems.Elem(ri), "oldlen < newlen");
1012 }
1013 }
1014
1015
1016 if (loktestmode)
1017 (*testout) << "ok = " << int(ok)
1018 << "teterr = " << teterr
1019 << "minteterr = " << minteterr << endl;
1020
1021
1022 if (ok && teterr < tolerance)
1023 {
1024 canuse.Elem(ri) ++;
1025 /*
1026 (*testout) << "can use rule " << rule->Name()
1027 << ", err = " << teterr << endl;
1028 for (i = 1; i <= pmap.Size(); i++)
1029 (*testout) << pmap.Get(i) << " ";
1030 (*testout) << endl;
1031 */
1032
1033 if (strcmp (problems.Elem(ri), "other") == 0)
1034 {
1035 if (teterr < minother)
1036 minother = teterr;
1037 }
1038 else
1039 {
1040 if (teterr < minwithoutother)
1041 minwithoutother = teterr;
1042 }
1043 }
1044
1045
1046 if (teterr > minteterr) impossible = 0;
1047
1048 if (ok && teterr < minteterr)
1049 {
1050
1051 if (loktestmode)
1052 (*testout) << "use rule" << endl;
1053
1054 found = ri;
1055 minteterr = teterr;
1056
1057 if (testmode)
1058 {
1059 for (i = 1; i <= rule->GetNOldP(); i++)
1060 {
1061 (*testout) << "P" << i << ": Ref: "
1062 << rule->GetPoint (i) << " is: "
1063 << lpoints.Get(pmap.Get(i)) << endl;
1064 }
1065 }
1066
1067 tempnewpoints.SetSize (0);
1068 for (i = noldlp+1; i <= lpoints.Size(); i++)
1069 tempnewpoints.Append (lpoints.Get(i));
1070
1071 tempnewfaces.SetSize (0);
1072 for (i = noldlf+1; i <= lfaces.Size(); i++)
1073 tempnewfaces.Append (lfaces.Get(i));
1074
1075 tempdelfaces.SetSize (0);
1076 for (i = 1; i <= delfaces.Size(); i++)
1077 tempdelfaces.Append (delfaces.Get(i));
1078
1079 tempelements.SetSize (0);
1080 for (i = 1; i <= elements.Size(); i++)
1081 tempelements.Append (elements.Get(i));
1082 }
1083
1084
1085 lpoints.SetSize (noldlp);
1086 lfaces.SetSize (noldlf);
1087 delfaces.SetSize (0);
1088 elements.SetSize (0);
1089 }
1090
1091 npok = rule->GetNOldP();
1092 incnpok = 0;
1093 }
1094 }
1095
1096 nfok = rule->GetNOldF();
1097
1098 for (j = 1; j <= rule->GetNP (nfok); j++)
1099 {
1100 refpi = rule->GetPointNr (nfok, j);
1101 pused.Elem(pmap.Get(refpi))--;
1102
1103 if (pused.Get(pmap.Get(refpi)) == 0)
1104 {
1105 pmap.Set(refpi, 0);
1106 }
1107 }
1108
1109 }
1110 }
1111 if (loktestmode)
1112 (*testout) << "end rule" << endl;
1113 }
1114
1115 if (found)
1116 {
1117 for (i = 1; i <= tempnewpoints.Size(); i++)
1118 lpoints.Append (tempnewpoints.Get(i));
1119 for (i = 1; i <= tempnewfaces.Size(); i++)
1120 if (tempnewfaces.Get(i).PNum(1))
1121 lfaces.Append (tempnewfaces.Get(i));
1122 for (i = 1; i <= tempdelfaces.Size(); i++)
1123 delfaces.Append (tempdelfaces.Get(i));
1124 for (i = 1; i <= tempelements.Size(); i++)
1125 elements.Append (tempelements.Get(i));
1126 }
1127
1128 retminerr = minerr;
1129
1130
1131 if (impossible && found == 0)
1132 return -1;
1133
1134 return found;
1135 }
1136 }
1137