1 #include <mystdlib.h>
2 #include "meshing.hpp"
3
4 #ifdef WIN32
5 #define COMMASIGN ':'
6 #else
7 #define COMMASIGN ','
8 #endif
9
10
11 namespace netgen
12 {
13
14 extern const char * tetrules[];
15
LoadVMatrixLine(istream & ist,DenseMatrix & m,int line)16 void LoadVMatrixLine (istream & ist, DenseMatrix & m, int line)
17 {
18 char ch;
19 int pnum;
20 float f;
21
22 ist >> ch;
23 while (ch != '}')
24 {
25 ist.putback (ch);
26 ist >> f;
27 ist >> ch;
28 ist >> pnum;
29
30 if (ch == 'x' || ch == 'X')
31 m.Elem(line, 3 * pnum - 2) = f;
32 if (ch == 'y' || ch == 'Y')
33 m.Elem(line, 3 * pnum - 1) = f;
34 if (ch == 'z' || ch == 'Z')
35 m.Elem(line, 3 * pnum ) = f;
36
37 if (ch == 'p' || ch == 'P')
38 {
39 m.Elem(line , 3 * pnum-2) = f;
40 m.Elem(line+1, 3 * pnum-1) = f;
41 m.Elem(line+2, 3 * pnum ) = f;
42 }
43
44 ist >> ch;
45 if (ch == COMMASIGN)
46 ist >> ch;
47 }
48 }
49
50
51
52
53
NeighbourTrianglePoint(const threeint & t1,const threeint & t2) const54 int vnetrule :: NeighbourTrianglePoint (const threeint & t1, const threeint & t2) const
55 {
56 NgArray<int> tr1(3);
57 NgArray<int> tr2(3);
58 tr1.Elem(1)=t1.i1;
59 tr1.Elem(2)=t1.i2;
60 tr1.Elem(3)=t1.i3;
61 tr2.Elem(1)=t2.i1;
62 tr2.Elem(2)=t2.i2;
63 tr2.Elem(3)=t2.i3;
64
65
66 int ret=0;
67
68 for (int i=1; i<=3; i++)
69 {
70 for (int j=1; j<=3; j++)
71 {
72 if ((tr1.Get(i)==tr2.Get(j) && tr1.Get((i%3)+1)==tr2.Get((j%3)+1)) ||
73 (tr1.Get(i)==tr2.Get((j%3)+1) && tr1.Get((i%3)+1)==tr2.Get(j)))
74 {ret = tr2.Get((j+1)%3+1);}
75 }
76 }
77
78 return ret;
79
80 }
81
LoadRule(istream & ist)82 void vnetrule :: LoadRule (istream & ist)
83 {
84 char buf[256];
85 char ch, ok;
86 Point3d p;
87 Element2d face(TRIG);
88 int i, j, i1, i2, i3, fs, ii, ii1, ii2, ii3;
89 twoint edge;
90 DenseMatrix tempoldutonewu(30, 20),
91 tempoldutofreezone(30, 20),
92 tempoldutofreezonelimit(30, 20),
93 tfz(20, 20),
94 tfzl(20, 20);
95
96 tempoldutonewu = 0;
97 tempoldutofreezone = 0;
98 tfz = 0;
99 tfzl = 0;
100
101
102 noldp = 0;
103 noldf = 0;
104
105 ist.get (buf, sizeof(buf), '"');
106 ist.get (ch);
107 ist.get (buf, sizeof(buf), '"');
108 ist.get (ch);
109
110 delete [] name;
111 name = new char[strlen (buf) + 1];
112 strcpy (name, buf);
113 // (*mycout) << "Rule " << name << " found." << endl;
114
115 do
116 {
117 ist >> buf;
118
119 if (strcmp (buf, "quality") == 0)
120
121 {
122 ist >> quality;
123 }
124
125 else if (strcmp (buf, "flags") == 0)
126 {
127 ist >> ch;
128 while (ch != ';')
129 {
130 flags.Append (ch);
131 ist >> ch;
132 }
133 }
134
135 else if (strcmp (buf, "mappoints") == 0)
136 {
137 ist >> ch;
138
139 while (ch == '(')
140 {
141 ist >> p.X();
142 ist >> ch; // ','
143 ist >> p.Y();
144 ist >> ch; // ','
145 ist >> p.Z();
146 ist >> ch; // ')'
147
148 points.Append (p);
149 noldp++;
150
151 tolerances.SetSize (noldp);
152 tolerances.Elem(noldp) = 1;
153
154 ist >> ch;
155 while (ch != ';')
156 {
157 if (ch == '{')
158 {
159 ist >> tolerances.Elem(noldp);
160 ist >> ch; // '}'
161 }
162
163 ist >> ch;
164 }
165
166 ist >> ch;
167 }
168
169 ist.putback (ch);
170 }
171
172
173 else if (strcmp (buf, "mapfaces") == 0)
174 {
175 ist >> ch;
176
177 while (ch == '(')
178 {
179 face.SetType(TRIG);
180 ist >> face.PNum(1);
181 ist >> ch; // ','
182 ist >> face.PNum(2);
183 ist >> ch; // ','
184 ist >> face.PNum(3);
185 ist >> ch; // ')' or ','
186 if (ch == COMMASIGN)
187 {
188 face.SetType(QUAD);
189 ist >> face.PNum(4);
190 ist >> ch; // ')'
191 }
192 faces.Append (face);
193 noldf++;
194
195 ist >> ch;
196 while (ch != ';')
197 {
198 if (ch == 'd')
199 {
200 delfaces.Append (noldf);
201 ist >> ch; // 'e'
202 ist >> ch; // 'l'
203 }
204
205 ist >> ch;
206 }
207
208 ist >> ch;
209 }
210
211 ist.putback (ch);
212 }
213
214 else if (strcmp (buf, "mapedges") == 0)
215 {
216 ist >> ch;
217
218 while (ch == '(')
219 {
220 ist >> edge.i1;
221 ist >> ch; // ','
222 ist >> edge.i2;
223 ist >> ch; // ')'
224
225 edges.Append (edge);
226
227 ist >> ch;
228 while (ch != ';')
229 {
230 ist >> ch;
231 }
232
233 ist >> ch;
234 }
235
236 ist.putback (ch);
237 }
238
239
240 else if (strcmp (buf, "newpoints") == 0)
241 {
242 ist >> ch;
243
244 while (ch == '(')
245 {
246 ist >> p.X();
247 ist >> ch; // ','
248 ist >> p.Y();
249 ist >> ch; // ','
250 ist >> p.Z();
251 ist >> ch; // ')'
252
253 points.Append (p);
254
255 ist >> ch;
256 while (ch != ';')
257 {
258 if (ch == '{')
259 {
260 LoadVMatrixLine (ist, tempoldutonewu,
261 3 * (points.Size()-noldp) - 2);
262
263 ist >> ch; // '{'
264 LoadVMatrixLine (ist, tempoldutonewu,
265 3 * (points.Size()-noldp) - 1);
266
267 ist >> ch; // '{'
268 LoadVMatrixLine (ist, tempoldutonewu,
269 3 * (points.Size()-noldp) );
270 }
271
272 ist >> ch;
273 }
274
275 ist >> ch;
276 }
277
278 ist.putback (ch);
279 }
280
281 else if (strcmp (buf, "newfaces") == 0)
282 {
283 ist >> ch;
284
285 while (ch == '(')
286 {
287 face.SetType(TRIG);
288 ist >> face.PNum(1);
289 ist >> ch; // ','
290 ist >> face.PNum(2);
291 ist >> ch; // ','
292 ist >> face.PNum(3);
293 ist >> ch; // ')' or ','
294 if (ch == COMMASIGN)
295 {
296 face.SetType(QUAD);
297 ist >> face.PNum(4);
298 ist >> ch; // ')'
299 }
300 faces.Append (face);
301
302 ist >> ch;
303 while (ch != ';')
304 {
305 ist >> ch;
306 }
307
308 ist >> ch;
309 }
310
311 ist.putback (ch);
312 }
313
314 else if (strcmp (buf, "freezone") == 0)
315 {
316 ist >> ch;
317
318 while (ch == '(')
319 {
320 ist >> p.X();
321 ist >> ch; // ','
322 ist >> p.Y();
323 ist >> ch; // ','
324 ist >> p.Z();
325 ist >> ch; // ')'
326
327 freezone.Append (p);
328
329 ist >> ch;
330 while (ch != ';')
331 {
332 if (ch == '{')
333 {
334 LoadVMatrixLine (ist, tempoldutofreezone,
335 3 * freezone.Size() - 2);
336
337 ist >> ch; // '{'
338 LoadVMatrixLine (ist, tempoldutofreezone,
339 3 * freezone.Size() - 1);
340
341 ist >> ch; // '{'
342 LoadVMatrixLine (ist, tempoldutofreezone,
343 3 * freezone.Size() );
344 }
345
346 ist >> ch;
347 }
348
349 ist >> ch;
350 }
351
352 ist.putback (ch);
353 }
354 else if (strcmp (buf, "freezone2") == 0)
355 {
356 int k, nfp;
357
358 nfp = 0;
359 ist >> ch;
360
361 DenseMatrix hm1(3, 50), hm2(50, 50), hm3(50, 50);
362 hm3 = 0;
363
364 while (ch == '{')
365 {
366 hm1 = 0;
367 nfp++;
368 LoadVMatrixLine (ist, hm1, 1);
369
370 for (i = 1; i <= points.Size(); i++)
371 tfz.Elem(nfp, i) = hm1.Get(1, 3*i-2);
372
373
374 p.X() = p.Y() = p.Z() = 0;
375 for (i = 1; i <= points.Size(); i++)
376 {
377 p.X() += hm1.Get(1, 3*i-2) * points.Get(i).X();
378 p.Y() += hm1.Get(1, 3*i-2) * points.Get(i).Y();
379 p.Z() += hm1.Get(1, 3*i-2) * points.Get(i).Z();
380 }
381 freezone.Append (p);
382 freezonelimit.Append (p);
383
384 hm2 = 0;
385 for (i = 1; i <= 3 * noldp; i++)
386 hm2.Elem(i, i) = 1;
387 for (i = 1; i <= 3 * noldp; i++)
388 for (j = 1; j <= 3 * (points.Size() - noldp); j++)
389 hm2.Elem(j + 3 * noldp, i) = tempoldutonewu.Get(j, i);
390
391 for (i = 1; i <= 3; i++)
392 for (j = 1; j <= 3 * noldp; j++)
393 {
394 double sum = 0;
395 for (k = 1; k <= 3 * points.Size(); k++)
396 sum += hm1.Get(i, k) * hm2.Get(k, j);
397
398 hm3.Elem(i + 3 * (nfp-1), j) = sum;
399 }
400
401 // (*testout) << "freepoint: " << p << endl;
402
403 while (ch != ';')
404 ist >> ch;
405
406 ist >> ch;
407 }
408
409 tfzl = tfz;
410
411 tempoldutofreezone = hm3;
412 tempoldutofreezonelimit = hm3;
413 ist.putback(ch);
414 }
415
416 else if (strcmp (buf, "freezonelimit") == 0)
417 {
418 int k, nfp;
419 nfp = 0;
420 ist >> ch;
421
422 DenseMatrix hm1(3, 50), hm2(50, 50), hm3(50, 50);
423 hm3 = 0;
424
425 while (ch == '{')
426 {
427 hm1 = 0;
428 nfp++;
429 LoadVMatrixLine (ist, hm1, 1);
430
431 for (i = 1; i <= points.Size(); i++)
432 tfzl.Elem(nfp, i) = hm1.Get(1, 3*i-2);
433
434
435 p.X() = p.Y() = p.Z() = 0;
436 for (i = 1; i <= points.Size(); i++)
437 {
438 p.X() += hm1.Get(1, 3*i-2) * points.Get(i).X();
439 p.Y() += hm1.Get(1, 3*i-2) * points.Get(i).Y();
440 p.Z() += hm1.Get(1, 3*i-2) * points.Get(i).Z();
441 }
442 freezonelimit.Elem(nfp) = p;
443
444 hm2 = 0;
445 for (i = 1; i <= 3 * noldp; i++)
446 hm2.Elem(i, i) = 1;
447 for (i = 1; i <= 3 * noldp; i++)
448 for (j = 1; j <= 3 * (points.Size() - noldp); j++)
449 hm2.Elem(j + 3 * noldp, i) = tempoldutonewu.Get(j, i);
450
451 for (i = 1; i <= 3; i++)
452 for (j = 1; j <= 3 * noldp; j++)
453 {
454 double sum = 0;
455 for (k = 1; k <= 3 * points.Size(); k++)
456 sum += hm1.Get(i, k) * hm2.Get(k, j);
457
458 hm3.Elem(i + 3 * (nfp-1), j) = sum;
459 }
460
461 // (*testout) << "freepoint: " << p << endl;
462
463 while (ch != ';')
464 ist >> ch;
465
466 ist >> ch;
467 }
468
469 tempoldutofreezonelimit = hm3;
470 ist.putback(ch);
471 }
472
473 else if (strcmp (buf, "freeset") == 0)
474 {
475 freesets.Append (new NgArray<int>);
476
477 ist >> ch;
478
479 while (ch != ';')
480 {
481 ist.putback (ch);
482 ist >> i;
483 freesets.Last()->Append(i);
484 ist >> ch;
485 }
486 }
487
488 else if (strcmp (buf, "elements") == 0)
489 {
490 ist >> ch;
491
492 while (ch == '(')
493 {
494 elements.Append (Element(TET));
495
496 // elements.Last().SetNP(1);
497 ist >> elements.Last().PNum(1);
498 ist >> ch; // ','
499
500 if (ch == COMMASIGN)
501 {
502 // elements.Last().SetNP(2);
503 ist >> elements.Last().PNum(2);
504 ist >> ch; // ','
505 }
506 if (ch == COMMASIGN)
507 {
508 // elements.Last().SetNP(3);
509 ist >> elements.Last().PNum(3);
510 ist >> ch; // ','
511 }
512 if (ch == COMMASIGN)
513 {
514 // elements.Last().SetNP(4);
515 elements.Last().SetType(TET);
516 ist >> elements.Last().PNum(4);
517 ist >> ch; // ','
518 }
519 if (ch == COMMASIGN)
520 {
521 // elements.Last().SetNP(5);
522 elements.Last().SetType(PYRAMID);
523 ist >> elements.Last().PNum(5);
524 ist >> ch; // ','
525 }
526 if (ch == COMMASIGN)
527 {
528 // elements.Last().SetNP(6);
529 elements.Last().SetType(PRISM);
530 ist >> elements.Last().PNum(6);
531 ist >> ch; // ','
532 }
533
534 if (ch == COMMASIGN)
535 {
536 // elements.Last().SetNP(6);
537 elements.Last().SetType(HEX);
538 ist >> elements.Last().PNum(7);
539 ist >> ch; // ','
540 }
541 if (ch == COMMASIGN)
542 {
543 // elements.Last().SetNP(6);
544 elements.Last().SetType(HEX);
545 ist >> elements.Last().PNum(8);
546 ist >> ch; // ','
547 }
548
549 /*
550 orientations.Append (fourint());
551 orientations.Last().i1 = elements.Last().PNum(1);
552 orientations.Last().i2 = elements.Last().PNum(2);
553 orientations.Last().i3 = elements.Last().PNum(3);
554 orientations.Last().i4 = elements.Last().PNum(4);
555 */
556
557 ist >> ch;
558 while (ch != ';')
559 {
560 ist >> ch;
561 }
562
563 ist >> ch;
564 }
565
566 ist.putback (ch);
567 }
568
569 else if (strcmp (buf, "orientations") == 0)
570
571 {
572 ist >> ch;
573
574 while (ch == '(')
575 {
576 // fourint a = fourint();
577 orientations.Append (fourint());
578
579 ist >> orientations.Last().i1;
580 ist >> ch; // ','
581 ist >> orientations.Last().i2;
582 ist >> ch; // ','
583 ist >> orientations.Last().i3;
584 ist >> ch; // ','
585 ist >> orientations.Last().i4;
586 ist >> ch; // ','
587
588
589 ist >> ch;
590 while (ch != ';')
591 {
592 ist >> ch;
593 }
594
595 ist >> ch;
596 }
597
598 ist.putback (ch);
599 }
600
601
602 else if (strcmp (buf, "endrule") != 0)
603 {
604 PrintSysError ("Parser3d, unknown token " , buf);
605 }
606 }
607 while (!ist.eof() && strcmp (buf, "endrule") != 0);
608
609
610 // (*testout) << endl;
611 // (*testout) << Name() << endl;
612 // (*testout) << "no1 = " << GetNO() << endl;
613
614 oldutonewu.SetSize (3 * (points.Size() - noldp), 3 * noldp);
615 oldutonewu = 0;
616
617 for (i = 1; i <= oldutonewu.Height(); i++)
618 for (j = 1; j <= oldutonewu.Width(); j++)
619 oldutonewu.Elem(i, j) = tempoldutonewu.Elem(i, j);
620
621
622 /*
623 oldutofreezone = new SparseMatrixFlex (3 * freezone.Size(), 3 * noldp);
624 oldutofreezonelimit = new SparseMatrixFlex (3 * freezone.Size(), 3 * noldp);
625
626 oldutofreezone -> SetSymmetric(0);
627 oldutofreezonelimit -> SetSymmetric(0);
628 */
629
630 /*
631 oldutofreezone = new DenseMatrix (3 * freezone.Size(), 3 * noldp);
632 oldutofreezonelimit = new DenseMatrix (3 * freezone.Size(), 3 * noldp);
633
634 for (i = 1; i <= oldutofreezone->Height(); i++)
635 for (j = 1; j <= oldutofreezone->Width(); j++)
636 // if (j == 4 || j >= 7)
637 {
638 if (tempoldutofreezone.Elem(i, j))
639 (*oldutofreezone)(i, j) = tempoldutofreezone(i, j);
640 if (tempoldutofreezonelimit.Elem(i, j))
641 (*oldutofreezonelimit)(i, j) = tempoldutofreezonelimit(i, j);
642 }
643 */
644
645
646
647
648 oldutofreezone = new DenseMatrix (freezone.Size(), points.Size());
649 oldutofreezonelimit = new DenseMatrix (freezone.Size(), points.Size());
650 // oldutofreezone = new SparseMatrixFlex (freezone.Size(), points.Size());
651 // oldutofreezonelimit = new SparseMatrixFlex (freezone.Size(), points.Size());
652
653 for (i = 1; i <= freezone.Size(); i++)
654 for (j = 1; j <= points.Size(); j++)
655 {
656 if (tfz.Elem(i, j))
657 (*oldutofreezone).Elem(i, j) = tfz.Elem(i, j);
658 if (tfzl.Elem(i, j))
659 (*oldutofreezonelimit).Elem(i, j) = tfzl.Elem(i, j);
660 }
661
662 /*
663 (*testout) << "Rule " << Name() << endl;
664 (*testout) << "oldutofreezone = " << (*oldutofreezone) << endl;
665 (*testout) << "oldutofreezonelimit = " << (*oldutofreezonelimit) << endl;
666 */
667
668 freezonepi.SetSize (freezone.Size());
669 for (i = 1; i <= freezonepi.Size(); i++)
670 freezonepi.Elem(i) = 0;
671 for (i = 1; i <= freezone.Size(); i++)
672 for (j = 1; j <= noldp; j++)
673 if (Dist (freezone.Get(i), points.Get(j)) < 1e-8)
674 freezonepi.Elem(i) = j;
675
676
677
678
679 for (i = 1; i <= elements.Size(); i++)
680 {
681 if (elements.Elem(i).GetNP() == 4)
682 {
683 orientations.Append (fourint());
684 orientations.Last().i1 = elements.Get(i).PNum(1);
685 orientations.Last().i2 = elements.Get(i).PNum(2);
686 orientations.Last().i3 = elements.Get(i).PNum(3);
687 orientations.Last().i4 = elements.Get(i).PNum(4);
688 }
689 if (elements.Elem(i).GetNP() == 5)
690 {
691 orientations.Append (fourint());
692 orientations.Last().i1 = elements.Get(i).PNum(1);
693 orientations.Last().i2 = elements.Get(i).PNum(2);
694 orientations.Last().i3 = elements.Get(i).PNum(3);
695 orientations.Last().i4 = elements.Get(i).PNum(5);
696
697 orientations.Append (fourint());
698 orientations.Last().i1 = elements.Get(i).PNum(1);
699 orientations.Last().i2 = elements.Get(i).PNum(3);
700 orientations.Last().i3 = elements.Get(i).PNum(4);
701 orientations.Last().i4 = elements.Get(i).PNum(5);
702 }
703 }
704
705
706
707 if (freesets.Size() == 0)
708 {
709 freesets.Append (new NgArray<int>);
710 for (i = 1; i <= freezone.Size(); i++)
711 freesets.Elem(1)->Append(i);
712 }
713
714
715 // testout << "Freezone: " << endl;
716
717 // for (i = 1; i <= freezone.Size(); i++)
718 // (*testout) << "freepoint: " << freezone.Get(i) << endl;
719 Vector vp(points.Size()), vfp(freezone.Size());
720
721
722 if (quality < 100)
723 {
724 for (int i = 1; i <= 3; i++)
725 {
726 for (int j = 1; j <= points.Size(); j++)
727 vp(j-1) = points.Get(j).X(i);
728 oldutofreezone->Mult(vp, vfp);
729 for (int j = 1; j <= freezone.Size(); j++)
730 freezone.Elem(j).X(i) = vfp(j-1);
731 }
732 // for (i = 1; i <= freezone.Size(); i++)
733 // (*testout) << "freepoint: " << freezone.Get(i) << endl;
734 }
735
736
737 for (fs = 1; fs <= freesets.Size(); fs++)
738 {
739 freefaces.Append (new NgArray<threeint>);
740
741 NgArray<int> & freeset = *freesets.Elem(fs);
742 NgArray<threeint> & freesetfaces = *freefaces.Last();
743
744 for (ii1 = 1; ii1 <= freeset.Size(); ii1++)
745 for (ii2 = 1; ii2 <= freeset.Size(); ii2++)
746 for (ii3 = 1; ii3 <= freeset.Size(); ii3++)
747 if (ii1 < ii2 && ii1 < ii3 && ii2 != ii3)
748 {
749 i1 = freeset.Get(ii1);
750 i2 = freeset.Get(ii2);
751 i3 = freeset.Get(ii3);
752
753 Vec3d v1, v2, n;
754
755 v1 = freezone.Get(i3) - freezone.Get(i1);
756 v2 = freezone.Get(i2) - freezone.Get(i1);
757 n = Cross (v1, v2);
758 n /= n.Length();
759 // (*testout) << "i1,2,3 = " << i1 << ", " << i2 << ", " << i3 << endl;
760 // (*testout) << "v1 = " << v1 << " v2 = " << v2 << " n = " << n << endl;
761 ok = 1;
762 for (ii = 1; ii <= freeset.Size(); ii++)
763 {
764 i = freeset.Get(ii);
765 // (*testout) << "i = " << i << endl;
766 if (i != i1 && i != i2 && i != i3)
767 if ( (freezone.Get(i) - freezone.Get(i1)) * n < 0 ) ok = 0;
768 }
769
770 if (ok)
771 {
772 freesetfaces.Append (threeint());
773 freesetfaces.Last().i1 = i1;
774 freesetfaces.Last().i2 = i2;
775 freesetfaces.Last().i3 = i3;
776 }
777 }
778 }
779
780 for (fs = 1; fs <= freesets.Size(); fs++)
781 {
782 freefaceinequ.Append (new DenseMatrix (freefaces.Get(fs)->Size(), 4));
783 }
784
785
786 {
787 int minn;
788 // NgArray<int> pnearness (noldp);
789 pnearness.SetSize (noldp);
790
791 for (i = 1; i <= pnearness.Size(); i++)
792 pnearness.Elem(i) = INT_MAX/10;
793
794 for (j = 1; j <= GetNP(1); j++)
795 pnearness.Elem(GetPointNr (1, j)) = 0;
796
797 do
798 {
799 ok = 1;
800
801 for (i = 1; i <= noldf; i++)
802 {
803 minn = INT_MAX/10;
804 for (j = 1; j <= GetNP(i); j++)
805 minn = min2 (minn, pnearness.Get(GetPointNr (i, j)));
806
807 for (j = 1; j <= GetNP(i); j++)
808 if (pnearness.Get(GetPointNr (i, j)) > minn+1)
809 {
810 ok = 0;
811 pnearness.Elem(GetPointNr (i, j)) = minn+1;
812 }
813 }
814
815 for (i = 1; i <= edges.Size(); i++)
816 {
817 int pi1 = edges.Get(i).i1;
818 int pi2 = edges.Get(i).i2;
819
820 if (pnearness.Get(pi1) > pnearness.Get(pi2)+1)
821 {
822 ok = 0;
823 pnearness.Elem(pi1) = pnearness.Get(pi2)+1;
824 }
825 if (pnearness.Get(pi2) > pnearness.Get(pi1)+1)
826 {
827 ok = 0;
828 pnearness.Elem(pi2) = pnearness.Get(pi1)+1;
829 }
830 }
831
832
833 for (i = 1; i <= elements.Size(); i++)
834 if (elements.Get(i).GetNP() == 6) // prism rule
835 {
836 for (j = 1; j <= 3; j++)
837 {
838 int pi1 = elements.Get(i).PNum(j);
839 int pi2 = elements.Get(i).PNum(j+3);
840
841 if (pnearness.Get(pi1) > pnearness.Get(pi2)+1)
842 {
843 ok = 0;
844 pnearness.Elem(pi1) = pnearness.Get(pi2)+1;
845 }
846 if (pnearness.Get(pi2) > pnearness.Get(pi1)+1)
847 {
848 ok = 0;
849 pnearness.Elem(pi2) = pnearness.Get(pi1)+1;
850 }
851 }
852 }
853 }
854 while (!ok);
855
856 maxpnearness = 0;
857 for (i = 1; i <= pnearness.Size(); i++)
858 maxpnearness = max2 (maxpnearness, pnearness.Get(i));
859
860
861 fnearness.SetSize (noldf);
862
863 for (i = 1; i <= noldf; i++)
864 {
865 fnearness.Elem(i) = 0;
866 for (j = 1; j <= GetNP(i); j++)
867 fnearness.Elem(i) += pnearness.Get(GetPointNr (i, j));
868 }
869
870 // (*testout) << "rule " << name << ", pnear = " << pnearness << endl;
871 }
872
873
874 //Table of edges:
875 for (fs = 1; fs <= freesets.Size(); fs++)
876 {
877 freeedges.Append (new NgArray<twoint>);
878
879 // NgArray<int> & freeset = *freesets.Get(fs);
880 NgArray<twoint> & freesetedges = *freeedges.Last();
881 NgArray<threeint> & freesetfaces = *freefaces.Get(fs);
882 int k,l;
883 INDEX ind;
884
885 for (k = 1; k <= freesetfaces.Size(); k++)
886 {
887 // threeint tr = freesetfaces.Get(k);
888
889 for (l = k+1; l <= freesetfaces.Size(); l++)
890 {
891 ind = NeighbourTrianglePoint(freesetfaces.Get(k), freesetfaces.Get(l));
892 if (!ind) continue;
893
894 INDEX_3 f1(freesetfaces.Get(k).i1,
895 freesetfaces.Get(k).i2,
896 freesetfaces.Get(k).i3);
897 INDEX_3 f2(freesetfaces.Get(l).i1,
898 freesetfaces.Get(l).i2,
899 freesetfaces.Get(l).i3);
900 INDEX_2 ed(0, 0);
901 for (int f11 = 1; f11 <= 3; f11++)
902 for (int f12 = 1; f12 <= 3; f12++)
903 if (f11 != f12)
904 for (int f21 = 1; f21 <= 3; f21++)
905 for (int f22 = 1; f22 <= 3; f22++)
906 if (f1.I(f11) == f2.I(f21) && f1.I(f12) == f2.I(f22))
907 {
908 ed.I(1) = f1.I(f11);
909 ed.I(2) = f1.I(f12);
910 }
911 // (*testout) << "ed = " << ed.I(1) << "-" << ed.I(2) << endl;
912 // (*testout) << "ind = " << ind << " ed = " << ed << endl;
913 for (int eli = 1; eli <= GetNOldF(); eli++)
914 {
915 if (GetNP(eli) == 4)
916 {
917 for (int elr = 1; elr <= 4; elr++)
918 {
919 if (GetPointNrMod (eli, elr) == ed.I(1) &&
920 GetPointNrMod (eli, elr+2) == ed.I(2))
921 {
922 /*
923 (*testout) << "ed is diagonal of rectangle" << endl;
924 (*testout) << "ed = " << ed.I(1) << "-" << ed.I(2) << endl;
925 (*testout) << "ind = " << ind << endl;
926 */
927 ind = 0;
928 }
929
930 }
931 }
932 }
933
934 if (ind)
935 {
936 /*
937 (*testout) << "new edge from face " << k
938 << " = (" << freesetfaces.Get(k).i1
939 << ", " << freesetfaces.Get(k).i2
940 << ", " << freesetfaces.Get(k).i3
941 << "), point " << ind << endl;
942 */
943 freesetedges.Append(twoint(k,ind));
944 }
945 }
946 }
947 }
948
949 }
950
951
952
953
954
LoadRules(const char * filename,const char ** prules)955 void Meshing3 :: LoadRules (const char * filename, const char ** prules)
956 {
957 char buf[256];
958 istream * ist;
959 char *tr1 = NULL;
960
961 if (filename)
962 {
963 PrintMessage (3, "rule-filename = ", filename);
964 ist = new ifstream (filename);
965 }
966 else
967 {
968 /* connect tetrules to one string */
969 PrintMessage (3, "Use internal rules");
970 if (!prules) prules = tetrules;
971
972 const char ** hcp = prules;
973 size_t len = 0;
974 while (*hcp)
975 {
976 len += strlen (*hcp);
977 hcp++;
978 }
979 tr1 = new char[len+1];
980 tr1[0] = 0;
981 hcp = prules; // tetrules;
982
983
984 char * tt1 = tr1;
985 while (*hcp)
986 {
987 strcat (tt1, *hcp);
988 tt1 += strlen (*hcp);
989 hcp++;
990 }
991
992
993 #ifdef WIN32
994 // VC++ 2005 workaround
995 for(size_t i=0; i<len; i++)
996 if(tr1[i] == ',')
997 tr1[i] = ':';
998 #endif
999
1000 ist = new istringstream (tr1);
1001 }
1002
1003 if (!ist->good())
1004 {
1005 cerr << "Rule description file " << filename << " not found" << endl;
1006 delete ist;
1007 exit (1);
1008 }
1009
1010 while (!ist->eof())
1011 {
1012 buf[0] = 0;
1013 (*ist) >> buf;
1014
1015 if (strcmp (buf, "rule") == 0)
1016 {
1017 vnetrule * rule = new vnetrule;
1018 rule -> LoadRule(*ist);
1019 rules.Append (rule);
1020 if (!rule->TestOk())
1021 {
1022 PrintSysError ("Parser3d: Rule ", rules.Size(), " not ok");
1023 exit (1);
1024 }
1025 }
1026 else if (strcmp (buf, "tolfak") == 0)
1027 {
1028 (*ist) >> tolfak;
1029 }
1030 }
1031 delete ist;
1032 delete [] tr1;
1033 }
1034 }
1035