1 #include <mystdlib.h>
2 #include <myadt.hpp>
3
4 #include <linalg.hpp>
5 #include <csg.hpp>
6 #include <meshing.hpp>
7
8
9 namespace netgen
10 {
Identification(int anr,const CSGeometry & ageom)11 Identification :: Identification (int anr, const CSGeometry & ageom)
12 : geom(ageom), identfaces(10)
13 {
14 nr = anr;
15 }
16
~Identification()17 Identification :: ~Identification ()
18 {
19 ;
20 }
21
22
operator <<(ostream & ost,Identification & ident)23 ostream & operator<< (ostream & ost, Identification & ident)
24 {
25 ident.Print (ost);
26 return ost;
27 }
28
29
30 /*
31 void Identification :: IdentifySpecialPoints (ARRAY<class SpecialPoint> & points)
32 {
33 ;
34 }
35 */
36
37
38 int Identification ::
Identifyable(const SpecialPoint & sp1,const SpecialPoint & sp2,const TABLE<int> & specpoint2solid,const TABLE<int> & specpoint2surface) const39 Identifyable (const SpecialPoint & sp1, const SpecialPoint & sp2,
40 const TABLE<int> & specpoint2solid,
41 const TABLE<int> & specpoint2surface) const
42 {
43 cout << "Identification::Identifyable called for base-class" << endl;
44 return 0;
45 }
46
47 int Identification ::
Identifyable(const Point<3> & p1,const Point<3> & sp2) const48 Identifyable (const Point<3> & p1, const Point<3> & sp2) const
49 {
50 cout << "Identification::Identifyable called for base-class" << endl;
51 return 0;
52 }
53
54
55 int Identification ::
IdentifyableCandidate(const SpecialPoint & sp1) const56 IdentifyableCandidate (const SpecialPoint & sp1) const
57 {
58 return 1;
59 }
60
61
62 int Identification ::
ShortEdge(const SpecialPoint & sp1,const SpecialPoint & sp2) const63 ShortEdge (const SpecialPoint & sp1, const SpecialPoint & sp2) const
64 {
65 return 0;
66 }
67
GetIdentifiedPoint(class Mesh & mesh,int pi)68 int Identification :: GetIdentifiedPoint (class Mesh & mesh, int pi)
69 {
70 cout << "Identification::GetIdentifiedPoint called for base-class" << endl;
71 return -1;
72 }
73
IdentifyPoints(Mesh & mesh)74 void Identification :: IdentifyPoints (Mesh & mesh)
75 {
76 cout << "Identification::IdentifyPoints called for base-class" << endl;
77 ;
78 }
79
IdentifyFaces(class Mesh & mesh)80 void Identification :: IdentifyFaces (class Mesh & mesh)
81 {
82 cout << "Identification::IdentifyFaces called for base-class" << endl;
83 ;
84 }
85
86 void Identification ::
BuildSurfaceElements(ARRAY<Segment> & segs,Mesh & mesh,const Surface * surf)87 BuildSurfaceElements (ARRAY<Segment> & segs,
88 Mesh & mesh, const Surface * surf)
89 {
90 cout << "Identification::BuildSurfaceElements called for base-class" << endl;
91 ;
92 }
93
94
95 void Identification ::
BuildVolumeElements(ARRAY<class Element2d> & surfels,class Mesh & mesh)96 BuildVolumeElements (ARRAY<class Element2d> & surfels,
97 class Mesh & mesh)
98 {
99 ;
100 }
101
102 void Identification ::
GetIdentifiedFaces(ARRAY<INDEX_2> & idfaces) const103 GetIdentifiedFaces (ARRAY<INDEX_2> & idfaces) const
104 {
105 idfaces.SetSize(0);
106 for (int i = 1; i <= identfaces.GetNBags(); i++)
107 for (int j = 1; j <= identfaces.GetBagSize(i); j++)
108 {
109 INDEX_2 i2;
110 int val;
111 identfaces.GetData (i, j, i2, val);
112 idfaces.Append (i2);
113 }
114 }
115
116
117
118
119 PeriodicIdentification ::
PeriodicIdentification(int anr,const CSGeometry & ageom,const Surface * as1,const Surface * as2)120 PeriodicIdentification (int anr,
121 const CSGeometry & ageom,
122 const Surface * as1,
123 const Surface * as2)
124 : Identification(anr, ageom)
125 {
126 s1 = as1;
127 s2 = as2;
128 }
129
~PeriodicIdentification()130 PeriodicIdentification :: ~PeriodicIdentification ()
131 {
132 ;
133 }
134
135 /*
136 void PeriodicIdentification :: IdentifySpecialPoints
137 (ARRAY<class SpecialPoint> & points)
138 {
139 int i, j;
140 int bestj;
141 double bestval, val;
142
143 for (i = 1; i <= points.Size(); i++)
144 {
145 Point<3> p1 = points.Get(i).p;
146 Point<3> hp1 = p1;
147 s1->Project (hp1);
148 if (Dist (p1, hp1) > 1e-6) continue;
149
150 Vec<3> n1;
151 s1->GetNormalVector (p1, n1);
152 n1 /= n1.Length();
153 if ( fabs(n1 * points.Get(i).v) > 1e-3)
154 continue;
155
156 bestval = 1e8;
157 bestj = 1;
158 for (j = 1; j <= points.Size(); j++)
159 {
160 Point<3> p2= points.Get(j).p;
161 Point<3> hp2 = p2;
162 s2->Project (hp2);
163 if (Dist (p2, hp2) > 1e-6) continue;
164
165 Vec<3> n2;
166 s2->GetNormalVector (p2, n2);
167 n2 /= n2.Length();
168 if ( fabs(n2 * points.Get(j).v) > 1e-3)
169 continue;
170
171
172 Vec<3> v(p1, p2);
173 double vl = v.Length();
174 double cl = fabs (v*n1);
175
176 val = 1 - cl*cl/(vl*vl);
177
178 val += (points.Get(i).v - points.Get(j).v).Length();
179
180 if (val < bestval)
181 {
182 bestj = j;
183 bestval = val;
184 }
185 }
186
187 (*testout) << "Identify Periodic special points: pi = "
188 << points.Get(i).p << ", vi = " << points.Get(i).v
189 << " pj = " << points.Get(bestj).p
190 << ", vj = " << points.Get(bestj).v
191 << " bestval = " << bestval << endl;
192 }
193 }
194 */
195
196 int PeriodicIdentification ::
Identifyable(const SpecialPoint & sp1,const SpecialPoint & sp2,const TABLE<int> & specpoint2solid,const TABLE<int> & specpoint2surface) const197 Identifyable (const SpecialPoint & sp1, const SpecialPoint & sp2,
198 const TABLE<int> & specpoint2solid,
199 const TABLE<int> & specpoint2surface) const
200 {
201 int i;
202 double val;
203
204 SpecialPoint hsp1 = sp1;
205 SpecialPoint hsp2 = sp2;
206
207 for (i = 1; i <= 1; i++)
208 {
209 // Swap (hsp1, hsp2);
210
211 if (!s1->PointOnSurface (hsp1.p))
212 continue;
213
214 Vec<3> n1;
215 n1 = s1->GetNormalVector (hsp1.p);
216 n1 /= n1.Length();
217 if ( fabs(n1 * hsp1.v) > 1e-3)
218 continue;
219
220
221 if (!s2->PointOnSurface(hsp2.p))
222 continue;
223
224 Vec<3> n2;
225 n2 = s2->GetNormalVector (hsp2.p);
226 n2 /= n2.Length();
227 if ( fabs(n2 * hsp2.v) > 1e-3)
228 continue;
229
230
231 Vec<3> v = hsp2.p - hsp1.p;
232 double vl = v.Length();
233 double cl = fabs (v*n1);
234
235
236 val = 1 - cl*cl/(vl*vl);
237 val += (hsp1.v - hsp2.v).Length();
238
239 if (val < 1e-6)
240 return 1;
241 }
242
243 return 0;
244 }
245
246 int PeriodicIdentification ::
Identifyable(const Point<3> & p1,const Point<3> & p2) const247 Identifyable (const Point<3> & p1, const Point<3> & p2) const
248 {
249 return (s1->PointOnSurface (p1) &&
250 s2->PointOnSurface (p2));
251 }
252
253
254
255
256 int PeriodicIdentification ::
GetIdentifiedPoint(class Mesh & mesh,int pi)257 GetIdentifiedPoint (class Mesh & mesh, int pi)
258 {
259 const Surface * sold, *snew;
260 const Point<3> & p = mesh.Point (pi);
261
262 if (s1->PointOnSurface (p))
263 {
264 snew = s2;
265 }
266 else
267 {
268 if (s2->PointOnSurface (p))
269 {
270 snew = s1;
271 }
272 else
273 {
274 cerr << "GetIdenfifiedPoint: Not possible" << endl;
275 exit (1);
276 }
277 }
278
279 // project to other surface
280 Point<3> hp = p;
281 snew->Project (hp);
282
283 int i;
284 int newpi = 0;
285 for (i = 1; i <= mesh.GetNP(); i++)
286 if (Dist2 (mesh.Point(i), hp) < 1e-12)
287 {
288 newpi = i;
289 break;
290 }
291 if (!newpi)
292 newpi = mesh.AddPoint (hp);
293
294 if (snew == s2)
295 mesh.GetIdentifications().Add (pi, newpi, nr);
296 else
297 mesh.GetIdentifications().Add (newpi, pi, nr);
298
299 mesh.GetIdentifications().SetType(nr,Identifications::PERIODIC);
300
301 /*
302 (*testout) << "Identify points(periodic), nr = " << nr << ": " << mesh.Point(pi)
303 << " and " << mesh.Point(newpi)
304 << ((snew == s2) ? "" : " inverse")
305 << endl;
306 */
307 return newpi;
308 }
309
310
IdentifyPoints(class Mesh & mesh)311 void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh)
312 {
313 int i, j;
314 for (i = 1; i <= mesh.GetNP(); i++)
315 {
316 Point<3> p = mesh.Point(i);
317 if (s1->PointOnSurface (p))
318 {
319 Point<3> pp = p;
320 s2->Project (pp);
321 for (j = 1; j <= mesh.GetNP(); j++)
322 if (Dist2(mesh.Point(j), pp) < 1e-6)
323 {
324 mesh.GetIdentifications().Add (i, j, nr);
325 /*
326 (*testout) << "Identify points(periodic:), nr = " << nr << ": "
327 << mesh.Point(i) << " - " << mesh.Point(j) << endl;
328 */
329 }
330 }
331 }
332
333 mesh.GetIdentifications().SetType(nr,Identifications::PERIODIC);
334 }
335
336
IdentifyFaces(class Mesh & mesh)337 void PeriodicIdentification :: IdentifyFaces (class Mesh & mesh)
338 {
339 int i, j, k, l;
340 int fi1, fi2, side;
341 for (i = 1; i <= mesh.GetNFD(); i++)
342 for (j = 1; j <= mesh.GetNFD(); j++)
343 {
344 int surfi = mesh.GetFaceDescriptor(i).SurfNr();
345 int surfj = mesh.GetFaceDescriptor(j).SurfNr();
346 if (surfi == surfj)
347 continue;
348
349 if (geom.GetSurface (surfi) != s1 ||
350 geom.GetSurface (surfj) != s2)
351 continue;
352
353 int idok = 1;
354
355
356 // (*testout) << "check faces " << i << " and " << j << endl;
357 for (side = 1; side <= 2 && idok; side++)
358 {
359 if (side == 1)
360 {
361 fi1 = i;
362 fi2 = j;
363 }
364 else
365 {
366 fi1 = j;
367 fi2 = i;
368 }
369
370 for (k = 1; k <= mesh.GetNSeg(); k++)
371 {
372 const Segment & seg1 = mesh.LineSegment(k);
373 if (seg1.si != fi1)
374 continue;
375
376 int foundother = 0;
377 for (l = 1; l <= mesh.GetNSeg(); l++)
378 {
379 const Segment & seg2 = mesh.LineSegment(l);
380 if (seg2.si != fi2)
381 continue;
382
383 // (*testout) << "seg1 = " << seg1.p1 << "-" << seg1.p2 << ", seg2 = " << seg2.p1 << "-" << seg2.p2;
384
385 if (side == 1)
386 {
387 if (mesh.GetIdentifications().Get (seg1.p1, seg2.p1) &&
388 mesh.GetIdentifications().Get (seg1.p2, seg2.p2))
389 {
390 foundother = 1;
391 break;
392 }
393
394 if (mesh.GetIdentifications().Get (seg1.p1, seg2.p2) &&
395 mesh.GetIdentifications().Get (seg1.p2, seg2.p1))
396 {
397 foundother = 1;
398 break;
399 }
400 }
401 else
402 {
403 if (mesh.GetIdentifications().Get (seg2.p1, seg1.p1) &&
404 mesh.GetIdentifications().Get (seg2.p2, seg1.p2))
405 {
406 foundother = 1;
407 break;
408 }
409
410 if (mesh.GetIdentifications().Get (seg2.p1, seg1.p2) &&
411 mesh.GetIdentifications().Get (seg2.p2, seg1.p1))
412 {
413 foundother = 1;
414 break;
415 }
416 }
417 }
418
419 if (!foundother)
420 {
421 idok = 0;
422 break;
423 }
424 }
425 }
426
427
428 if (idok)
429 {
430 // (*testout) << "Identify faces " << i << " and " << j << endl;
431 INDEX_2 fpair(i,j);
432 fpair.Sort();
433 identfaces.Set (fpair, 1);
434 }
435 }
436 }
437
438
439
440 void PeriodicIdentification ::
BuildSurfaceElements(ARRAY<Segment> & segs,Mesh & mesh,const Surface * surf)441 BuildSurfaceElements (ARRAY<Segment> & segs,
442 Mesh & mesh, const Surface * surf)
443 {
444 int found = 0;
445 int fother;
446
447 int facei = segs.Get(1).si;
448 int surfnr = mesh.GetFaceDescriptor(facei).SurfNr();
449
450 if (geom.GetSurface(surfnr) == s1 ||
451 geom.GetSurface(surfnr) == s2)
452 {
453 ARRAY<int> copy_points;
454
455 for (int i = 1; i <= mesh.GetNSE(); i++)
456 {
457 const Element2d & sel = mesh.SurfaceElement(i);
458 INDEX_2 fpair (facei, sel.GetIndex());
459 fpair.Sort();
460 if (identfaces.Used (fpair))
461 {
462 for (int k = 0; k < sel.GetNP(); k++)
463 if (!copy_points.Contains (sel[k]))
464 copy_points.Append (sel[k]);
465 }
466 }
467 BubbleSort (copy_points);
468 for (int k = 0; k < copy_points.Size(); k++)
469 GetIdentifiedPoint (mesh, copy_points[k]);
470
471
472
473 for (int i = 1; i <= mesh.GetNSE(); i++)
474 {
475 const Element2d & sel = mesh.SurfaceElement(i);
476 INDEX_2 fpair (facei, sel.GetIndex());
477 fpair.Sort();
478 if (identfaces.Used (fpair))
479 {
480 found = 1;
481 fother = sel.GetIndex();
482
483 // copy element
484 Element2d newel(sel.GetType());
485 newel.SetIndex (facei);
486 for (int k = 0; k < sel.GetNP(); k++)
487 newel[k] = GetIdentifiedPoint (mesh, sel[k]);
488
489 Vec<3> nt = Cross (Point<3> (mesh[newel[1]])- Point<3> (mesh[newel[0]]),
490 Point<3> (mesh[newel[2]])- Point<3> (mesh[newel[0]]));
491
492 Vec<3> nsurf = geom.GetSurface (surfnr)->GetNormalVector (mesh[newel[0]]);
493 if (nsurf * nt < 0)
494 Swap (newel[0], newel[2]);
495
496 mesh.AddSurfaceElement (newel);
497 }
498 }
499 }
500
501 if (found)
502 {
503 // (*mycout) << " copy face " << facei << " from face " << fother;
504 PrintMessage (4, " copy face ", facei, " from face ", fother);
505
506 segs.SetSize(0);
507 }
508 }
509
510
511
512
513
514
515
516
Print(ostream & ost) const517 void PeriodicIdentification :: Print (ostream & ost) const
518 {
519 ost << "Periodic Identifiaction, surfaces: "
520 << s1->Name() << " - " << s2->Name() << endl;
521 s1->Print (ost);
522 ost << " - ";
523 s2->Print (ost);
524 ost << endl;
525 }
526
527
GetData(ostream & ost) const528 void PeriodicIdentification :: GetData (ostream & ost) const
529 {
530 ost << "periodic " << s1->Name() << " " << s2->Name();
531 }
532
533
534
535
536
537
538
539 CloseSurfaceIdentification ::
CloseSurfaceIdentification(int anr,const CSGeometry & ageom,const Surface * as1,const Surface * as2,const TopLevelObject * adomain,const Flags & flags)540 CloseSurfaceIdentification (int anr,
541 const CSGeometry & ageom,
542 const Surface * as1,
543 const Surface * as2,
544 const TopLevelObject * adomain,
545 const Flags & flags)
546 : Identification(anr, ageom)
547 {
548 s1 = as1;
549 s2 = as2;
550 domain = adomain;
551 ref_levels = int (flags.GetNumFlag ("reflevels", 2));
552 ref_levels_s1 = int (flags.GetNumFlag ("reflevels1", 0));
553 ref_levels_s2 = int (flags.GetNumFlag ("reflevels2", 0));
554 slices = flags.GetNumListFlag ("slices");
555 for(int i=0; i<slices.Size(); i++)
556 if((i==0 && slices[i] <= 0) ||
557 (i>0 && slices[i] <= slices[i-1]) ||
558 (slices[i] >= 1))
559 throw NgException ("slices have to be in ascending order, between 0 and 1");
560
561 // eps_n = 1e-3;
562 eps_n = 1e-6;
563
564 dom_surf_valid = 0;
565
566 if (domain)
567 for (int i = 0; i < geom.GetNTopLevelObjects(); i++)
568 if (domain == geom.GetTopLevelObject(i))
569 dom_nr = i;
570
571 usedirection = flags.NumListFlagDefined("direction");
572 if(usedirection)
573 {
574 for(int i=0; i<3; i++)
575 direction(i) = flags.GetNumListFlag("direction")[i];
576
577 direction.Normalize();
578 }
579 }
580
~CloseSurfaceIdentification()581 CloseSurfaceIdentification :: ~CloseSurfaceIdentification ()
582 {
583 ;
584 }
585
Print(ostream & ost) const586 void CloseSurfaceIdentification :: Print (ostream & ost) const
587 {
588 ost << "CloseSurface Identifiaction, surfaces: "
589 << s1->Name() << " - " << s2->Name() << endl;
590 s1->Print (ost);
591 s2->Print (ost);
592 ost << endl;
593 }
594
595
GetData(ostream & ost) const596 void CloseSurfaceIdentification :: GetData (ostream & ost) const
597 {
598 ost << "close surface " << s1->Name() << " " << s2->Name();
599 }
600
601
602 /*
603 void CloseSurfaceIdentification :: IdentifySpecialPoints
604 (ARRAY<class SpecialPoint> & points)
605 {
606 int i, j;
607 int bestj;
608 double bestval, val;
609
610 for (i = 1; i <= points.Size(); i++)
611 {
612 Point<3> p1 = points.Get(i).p;
613 Vec<3> n1;
614
615 if (!s1->PointOnSurface (p1))
616 continue;
617
618 s1->GetNormalVector (p1, n1);
619 n1 /= n1.Length();
620 if ( fabs(n1 * points.Get(i).v) > 1e-3)
621 continue;
622
623 bestval = 1e8;
624 bestj = 1;
625 for (j = 1; j <= points.Size(); j++)
626 {
627 Point<3> p2= points.Get(j).p;
628 if (!s2->PointOnSurface (p2))
629 continue;
630
631 Vec<3> n2;
632 s2->GetNormalVector (p2, n2);
633 n2 /= n2.Length();
634 if ( fabs(n2 * points.Get(j).v) > 1e-3)
635 continue;
636
637
638 Vec<3> v(p1, p2);
639 double vl = v.Length();
640 double cl = fabs (v*n1);
641
642 val = 1 - cl*cl/(vl*vl);
643
644 val += (points.Get(i).v - points.Get(j).v).Length();
645
646 if (val < bestval)
647 {
648 bestj = j;
649 bestval = val;
650 }
651 }
652
653 (*testout) << "Identify close surfaces special points: pi = "
654 << points.Get(i).p << ", vi = " << points.Get(i).v
655 << " pj = " << points.Get(bestj).p
656 << ", vj = " << points.Get(bestj).v
657 << " bestval = " << bestval << endl;
658 }
659 }
660 */
661
662 int CloseSurfaceIdentification ::
Identifyable(const SpecialPoint & sp1,const SpecialPoint & sp2,const TABLE<int> & specpoint2solid,const TABLE<int> & specpoint2surface) const663 Identifyable (const SpecialPoint & sp1, const SpecialPoint & sp2,
664 const TABLE<int> & specpoint2solid,
665 const TABLE<int> & specpoint2surface) const
666 {
667 //(*testout) << "identcheck: " << sp1.p << "; " << sp2.p << endl;
668
669 if (!dom_surf_valid)
670 {
671 const_cast<bool&> (dom_surf_valid) = 1;
672 ARRAY<int> & hsurf = const_cast<ARRAY<int>&> (domain_surfaces);
673
674 if (domain)
675 {
676 BoxSphere<3> hbox (geom.BoundingBox());
677 geom.GetIndependentSurfaceIndices (domain->GetSolid(), hbox, hsurf);
678 //(*testout) << "surfs of identification " << nr << ": " << endl << hsurf << endl;
679 }
680 else
681 {
682 hsurf.SetSize (geom.GetNSurf());
683 for (int j = 0; j < hsurf.Size(); j++)
684 hsurf[j] = j;
685 }
686 }
687
688 if (domain)
689 {
690 bool has1 = 0, has2 = 0;
691 for (int i = 0; i < specpoint2solid[sp1.nr].Size(); i++)
692 if (specpoint2solid[sp1.nr][i] == dom_nr)
693 { has1 = 1; break; }
694 for (int i = 0; i < specpoint2solid[sp2.nr].Size(); i++)
695 if (specpoint2solid[sp2.nr][i] == dom_nr)
696 { has2 = 1; break; }
697
698 if (!has1 || !has2)
699 {
700 //(*testout) << "failed at pos1" << endl;
701 return 0;
702 }
703 }
704
705 if (!s1->PointOnSurface (sp1.p))
706 {
707 //(*testout) << "failed at pos2" << endl;
708 return 0;
709 }
710
711 // (*testout) << "sp1 " << sp1.p << " sp2 " << sp2.p << endl
712 // << "specpoint2solid[sp1.nr] " << specpoint2solid[sp1.nr] << endl
713 // << "specpoint2solid[sp2.nr] " << specpoint2solid[sp2.nr] << endl;
714
715
716 Vec<3> n1 = s1->GetNormalVector (sp1.p);
717 n1.Normalize();
718 if ( fabs(n1 * sp1.v) > eps_n)
719 {
720 //(*testout) << "failed at pos3" << endl;
721 return 0;
722 }
723
724 if (!s2->PointOnSurface(sp2.p))
725 {
726 //(*testout) << "failed at pos4" << endl;
727 return 0;
728 }
729
730
731 Vec<3> n2 = s2->GetNormalVector (sp2.p);
732 n2.Normalize();
733 if ( fabs(n2 * sp2.v) > eps_n)
734 {
735 //(*testout) << "failed at pos5" << endl;
736 return 0;
737 }
738
739 // must have joint surface
740 bool joint = 0;
741
742 int j = 0, k = 0;
743 while (1)
744 {
745 int snr1 = specpoint2surface[sp1.nr][j];
746 int snr2 = specpoint2surface[sp2.nr][k];
747 if (snr1 < snr2)
748 {
749 j++;
750 if (j == specpoint2surface[sp1.nr].Size()) break;
751 }
752 else if (snr2 < snr1)
753 {
754 k++;
755 if (k == specpoint2surface[sp2.nr].Size()) break;
756 }
757 else
758 {
759 bool dom_surf = 0;
760 for (int l = 0; l < domain_surfaces.Size(); l++)
761 if (domain_surfaces[l] == snr1)
762 dom_surf = 1;
763
764 if (dom_surf)
765 {
766 Vec<3> hn1 = geom.GetSurface(snr1)->GetNormalVector (sp1.p);
767 Vec<3> hn2 = geom.GetSurface(snr1)->GetNormalVector (sp2.p);
768
769 if (hn1 * hn2 > 0)
770 {
771 joint = 1;
772 break;
773 }
774 }
775
776 j++;
777 if (j == specpoint2surface[sp1.nr].Size()) break;
778 k++;
779 if (k == specpoint2surface[sp2.nr].Size()) break;
780 }
781 }
782
783 if (!joint)
784 {
785 //(*testout) << "failed at pos6" << endl;
786 return 0;
787 }
788
789 Vec<3> v = sp2.p - sp1.p;
790 double vl = v.Length();
791 double cl = (usedirection) ? fabs(v*direction) : fabs (v*n1);
792
793
794 if(cl <= (1-eps_n*eps_n) * vl)
795 {
796 //(*testout) << "failed at pos7" << endl;
797 return 0;
798 }
799
800 double dl;
801
802 if(usedirection)
803 {
804 Vec<3> v1 = sp1.v - (sp1.v*direction)*direction; v1.Normalize();
805 Vec<3> v2 = sp2.v - (sp2.v*direction)*direction; v2.Normalize();
806
807 dl = (v1 - v2).Length();
808 }
809 else
810 dl = (sp1.v - sp2.v).Length();
811
812 if (dl < 0.1)
813 return 1;
814
815
816 //(*testout) << "failed at pos8" << endl;
817 return 0;
818 }
819
820 int CloseSurfaceIdentification ::
Identifyable(const Point<3> & p1,const Point<3> & p2) const821 Identifyable (const Point<3> & p1, const Point<3> & p2) const
822 {
823 // if (domain)
824 // if (!domain->GetSolid()->IsIn (p1) || !domain->GetSolid()->IsIn (p2))
825 // return 0;
826 return (s1->PointOnSurface (p1) && s2->PointOnSurface (p2));
827 }
828
829
830
831
832 int CloseSurfaceIdentification ::
IdentifyableCandidate(const SpecialPoint & sp1) const833 IdentifyableCandidate (const SpecialPoint & sp1) const
834 {
835 if (domain)
836 if (!domain->GetSolid()->IsIn (sp1.p))
837 return 0;
838
839 if (s1->PointOnSurface (sp1.p))
840 {
841 Vec<3> n1;
842 n1 = s1->GetNormalVector (sp1.p);
843 n1.Normalize();
844 if ( fabs(n1 * sp1.v) > eps_n)
845 return 0;
846 return 1;
847 }
848
849 if (s2->PointOnSurface (sp1.p))
850 {
851 Vec<3> n1;
852 n1 = s2->GetNormalVector (sp1.p);
853 n1.Normalize();
854 if ( fabs(n1 * sp1.v) > eps_n)
855 return 0;
856 return 1;
857 }
858 return 0;
859 }
860
861
862
863 int CloseSurfaceIdentification ::
ShortEdge(const SpecialPoint & sp1,const SpecialPoint & sp2) const864 ShortEdge (const SpecialPoint & sp1, const SpecialPoint & sp2) const
865 {
866 if ( (s1->PointOnSurface (sp1.p) && s2->PointOnSurface (sp2.p)) ||
867 (s1->PointOnSurface (sp2.p) && s2->PointOnSurface (sp1.p)) )
868 {
869 return 1;
870 }
871 return 0;
872 }
873
874
875
876 int CloseSurfaceIdentification ::
GetIdentifiedPoint(class Mesh & mesh,int pi)877 GetIdentifiedPoint (class Mesh & mesh, int pi)
878 {
879 const Surface * sold, *snew;
880 const Point<3> & p = mesh.Point (pi);
881
882 ARRAY<int,PointIndex::BASE> identmap(mesh.GetNP());
883 mesh.GetIdentifications().GetMap (nr, identmap);
884 if (identmap.Get(pi))
885 return identmap.Get(pi);
886
887
888 if (s1->PointOnSurface (p))
889 snew = s2;
890 else if (s2->PointOnSurface (p))
891 snew = s1;
892 else
893 {
894 (*testout) << "GetIdenfifiedPoint: Not possible" << endl;
895 (*testout) << "p = " << p << endl;
896 (*testout) << "surf1: " << (*s1) << endl
897 << "surf2: " << (*s2) << endl;
898
899 cerr << "GetIdenfifiedPoint: Not possible" << endl;
900 throw NgException ("GetIdenfifiedPoint: Not possible");
901 }
902
903 // project to other surface
904 Point<3> hp = p;
905 if(usedirection)
906 snew->SkewProject(hp,direction);
907 else
908 snew->Project (hp);
909
910 //(*testout) << "projecting " << p << " to " << hp << endl;
911
912 int newpi = 0;
913 for (int i = 1; i <= mesh.GetNP(); i++)
914 if (Dist2 (mesh.Point(i), hp) < 1e-12)
915 // if (Dist2 (mesh.Point(i), hp) < 1 * Dist2 (hp, p))
916 {
917 newpi = i;
918 break;
919 }
920 if (!newpi)
921 newpi = mesh.AddPoint (hp);
922
923 if (snew == s2)
924 {
925 mesh.GetIdentifications().Add (pi, newpi, nr);
926 //(*testout) << "add identification(1) " << pi << " - " << newpi << ", " << nr << endl;
927 }
928 else
929 {
930 mesh.GetIdentifications().Add (newpi, pi, nr);
931 //(*testout) << "add identification(2) " << newpi << " - " << pi << ", " << nr << endl;
932 }
933 mesh.GetIdentifications().SetType(nr,Identifications::CLOSESURFACES);
934
935
936 /*
937 (*testout) << "Identify points(closesurface), nr = " << nr << ": " << mesh.Point(pi)
938 << " and " << mesh.Point(newpi)
939 << ((snew == s2) ? "" : " inverse")
940 << endl;
941 */
942 return newpi;
943 }
944
945
946
947
948
IdentifyPoints(Mesh & mesh)949 void CloseSurfaceIdentification :: IdentifyPoints (Mesh & mesh)
950 {
951 int np = mesh.GetNP();
952
953 ARRAY<int> points_on_surf2;
954
955 for (int i2 = 1; i2 <= np; i2++)
956 if (s2->PointOnSurface (mesh.Point(i2)))
957 points_on_surf2.Append (i2);
958
959 ARRAY<int> surfs_of_p1;
960
961 for (int i1 = 1; i1 <= np; i1++)
962 {
963 Point<3> p1 = mesh.Point(i1);
964 // (*testout) << "p1 = " << i1 << " = " << p1 << endl;
965 if (domain && !domain->GetSolid()->IsIn (p1))
966 continue;
967
968 //if(domain) (*testout) << "p1 is in " << domain->GetSolid()->Name() << endl;
969
970 if (s1->PointOnSurface (p1))
971 {
972 int candi2 = 0;
973 double mindist = 1e10;
974
975 Vec<3> n1;
976 n1 = s1->GetNormalVector (p1);
977 n1.Normalize();
978
979 surfs_of_p1.SetSize(0);
980 for (int jj = 0; jj < domain_surfaces.Size(); jj++)
981 {
982 int j = domain_surfaces[jj];
983 if (geom.GetSurface(j) -> PointOnSurface(p1))
984 surfs_of_p1.Append (j);
985 }
986 //(*testout) << " surfs of p1 = " << endl << surfs_of_p1 << endl;
987
988 for (int ii2 = 0; ii2 < points_on_surf2.Size(); ii2++)
989 {
990 int i2 = points_on_surf2[ii2];
991 if (i2 == i1) continue;
992 const Point<3> p2 = mesh.Point(i2);
993
994 Vec<3> n = p2 - p1;
995 n.Normalize();
996
997 bool joint = 0;
998 for (int jj = 0; jj < surfs_of_p1.Size(); jj++)
999 {
1000 int j = surfs_of_p1[jj];
1001 if (geom.GetSurface(j) -> PointOnSurface(p2))
1002 {
1003 Vec<3> hn1 = geom.GetSurface(j)->GetNormalVector (p1);
1004 Vec<3> hn2 = geom.GetSurface(j)->GetNormalVector (p2);
1005
1006 if (hn1 * hn2 > 0)
1007 {
1008 joint = 1;
1009 break;
1010 }
1011 }
1012 }
1013
1014 if (!joint) continue;
1015
1016 if(usedirection)
1017 {
1018 if (fabs (n*direction) > 0.9)
1019 {
1020 Vec<3> p1p2 = p2-p1;
1021 double ndist = p1p2.Length2() - pow(p1p2*direction,2);
1022 if(ndist < mindist)
1023 {
1024 candi2 = i2;
1025 mindist = ndist;
1026 }
1027 }
1028
1029 }
1030 else
1031 {
1032 if (fabs (n * n1) > 0.9 &&
1033 Dist (p1, p2) < mindist)
1034 {
1035 candi2 = i2;
1036 mindist = Dist (p1, p2);
1037 }
1038 }
1039
1040 }
1041
1042 if (candi2)
1043 {
1044 //(*testout) << "identify points " << p1 << " - " << mesh.Point(candi2) << endl;
1045
1046 /*
1047 (*testout) << "Add Identification from CSI2, nr = " << nr << ", p1 = "
1048 << i1 << " = "
1049 << mesh[PointIndex(i1)] << ", p2 = " << candi2 << " = "
1050 << mesh[PointIndex(candi2)] << endl;
1051 */
1052 mesh.GetIdentifications().Add (i1, candi2, nr);
1053 mesh.GetIdentifications().SetType(nr,Identifications::CLOSESURFACES);
1054 //(*testout) << "add identification " << i1 << " - " << candi2 << ", " << nr << endl;
1055 }
1056 }
1057 }
1058 }
1059
1060
1061
IdentifyFaces(class Mesh & mesh)1062 void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh)
1063 {
1064 int fi1, fi2, side;
1065 int s1rep, s2rep;
1066
1067 for (int i = 0; i < geom.GetNSurf(); i++)
1068 {
1069 if (geom.GetSurface (i) == s1)
1070 s1rep = geom.GetSurfaceClassRepresentant(i);
1071 if (geom.GetSurface (i) == s2)
1072 s2rep = geom.GetSurfaceClassRepresentant(i);
1073 }
1074
1075 ARRAY<int> segs_on_face1, segs_on_face2;
1076
1077 identfaces.DeleteData();
1078
1079 //(*testout) << "identify faces, nr = " << nr << endl;
1080
1081 for (int i = 1; i <= mesh.GetNFD(); i++)
1082 {
1083 int surfi = mesh.GetFaceDescriptor(i).SurfNr();
1084 if (s1rep != surfi) continue;
1085
1086
1087 if (domain &&
1088 domain != geom.GetTopLevelObject (mesh.GetFaceDescriptor(i).DomainIn()-1) &&
1089 domain != geom.GetTopLevelObject (mesh.GetFaceDescriptor(i).DomainOut()-1))
1090 continue;
1091
1092 for (int j = 1; j <= mesh.GetNFD(); j++)
1093 {
1094 int surfj = mesh.GetFaceDescriptor(j).SurfNr();
1095
1096 if (surfi == surfj) continue;
1097 if (s2rep != surfj) continue;
1098
1099
1100 int idok = 1;
1101
1102 for (side = 1; side <= 2 && idok; side++)
1103 {
1104 if (side == 1)
1105 {
1106 fi1 = i;
1107 fi2 = j;
1108 }
1109 else
1110 {
1111 fi1 = j;
1112 fi2 = i;
1113 }
1114
1115
1116 segs_on_face1.SetSize(0);
1117 segs_on_face2.SetSize(0);
1118
1119 for (int k = 1; k <= mesh.GetNSeg(); k++)
1120 {
1121 if (mesh.LineSegment(k).si == fi1)
1122 segs_on_face1.Append (k);
1123 if (mesh.LineSegment(k).si == fi2)
1124 segs_on_face2.Append (k);
1125 }
1126
1127
1128 for (int k = 1; k <= mesh.GetNSeg(); k++)
1129 {
1130 const Segment & seg1 = mesh.LineSegment(k);
1131 if (seg1.si != fi1)
1132 continue;
1133
1134 int foundother = 0;
1135 /*
1136 for (int l = 1; l <= mesh.GetNSeg(); l++)
1137 {
1138 const Segment & seg2 = mesh.LineSegment(l);
1139 if (seg2.si != fi2)
1140 continue;
1141 */
1142 for (int ll = 0; ll < segs_on_face2.Size(); ll++)
1143 {
1144 int l = segs_on_face2[ll];
1145 const Segment & seg2 = mesh.LineSegment(l);
1146
1147 if (side == 1)
1148 {
1149 if (mesh.GetIdentifications().Get (seg1.p1, seg2.p1) &&
1150 mesh.GetIdentifications().Get (seg1.p2, seg2.p2))
1151 {
1152 foundother = 1;
1153 break;
1154 }
1155
1156 if (mesh.GetIdentifications().Get (seg1.p1, seg2.p2) &&
1157 mesh.GetIdentifications().Get (seg1.p2, seg2.p1))
1158 {
1159 foundother = 1;
1160 break;
1161 }
1162 }
1163 else
1164 {
1165 if (mesh.GetIdentifications().Get (seg2.p1, seg1.p1) &&
1166 mesh.GetIdentifications().Get (seg2.p2, seg1.p2))
1167 {
1168 foundother = 1;
1169 break;
1170 }
1171
1172 if (mesh.GetIdentifications().Get (seg2.p1, seg1.p2) &&
1173 mesh.GetIdentifications().Get (seg2.p2, seg1.p1))
1174 {
1175 foundother = 1;
1176 break;
1177 }
1178 }
1179 }
1180
1181 if (!foundother)
1182 {
1183 idok = 0;
1184 break;
1185 }
1186 }
1187 }
1188
1189
1190 if (idok)
1191 {
1192 //(*testout) << "Identification " << nr << ", identify faces " << i << " and " << j << endl;
1193 INDEX_2 fpair(i,j);
1194 fpair.Sort();
1195 identfaces.Set (fpair, 1);
1196 }
1197 }
1198 }
1199 }
1200
1201
1202
1203 void CloseSurfaceIdentification ::
BuildSurfaceElements(ARRAY<Segment> & segs,Mesh & mesh,const Surface * surf)1204 BuildSurfaceElements (ARRAY<Segment> & segs,
1205 Mesh & mesh, const Surface * surf)
1206 {
1207 bool found = 0;
1208 int cntquads = 0;
1209
1210 ARRAY<int,PointIndex::BASE> identmap;
1211 identmap = 0;
1212
1213 mesh.GetIdentifications().GetMap (nr, identmap);
1214
1215 for (int i = PointIndex::BASE; i < identmap.Size()+PointIndex::BASE; i++)
1216 if (identmap[i]) identmap[identmap[i]] = i;
1217
1218
1219 //(*testout) << "identification nr = " << nr << endl;
1220 //(*testout) << "surf = " << (*surf) << endl;
1221 //(*testout) << "domain = " << domain->GetSolid()->Name() << endl;
1222 //(*testout) << "segs = " << endl << segs << endl;
1223 //(*testout) << "identmap = " << endl << identmap << endl;
1224
1225 //ARRAY<bool> foundseg(segs.Size());
1226 //foundseg = false;
1227
1228 // insert quad layer:
1229 for (int i1 = 0; i1 < segs.Size(); i1++)
1230 {
1231 const Segment & s1 = segs[i1];
1232 if (identmap[s1.p1] && identmap[s1.p2])
1233 for (int i2 = 0; i2 < i1; i2++)
1234 {
1235 const Segment & s2 = segs[i2];
1236 //(*testout) << "checking " << s1 << " and " << s2 << " for ident." << endl;
1237
1238 if(domain && !((s1.domin == dom_nr ||
1239 s1.domout == dom_nr) &&
1240 (s2.domin == dom_nr ||
1241 s2.domout == dom_nr)))
1242 continue;
1243
1244 if ((mesh.GetIdentifications().Get (s1.p1, s2.p2, nr) &&
1245 mesh.GetIdentifications().Get (s1.p2, s2.p1, nr)) ||
1246 (mesh.GetIdentifications().Get (s2.p1, s1.p2, nr) &&
1247 mesh.GetIdentifications().Get (s2.p2, s1.p1, nr)))
1248 {
1249 Element2d el(s1.p1, s1.p2, s2.p1, s2.p2);
1250
1251 Vec<3> n = Cross (mesh[el[1]] - mesh[el[0]],
1252 mesh[el[3]] - mesh[el[0]]);
1253
1254 Vec<3> ns = surf->GetNormalVector (mesh[el[0]]);
1255
1256 if (n * ns < 0)
1257 {
1258 Swap (el.PNum(1), el.PNum(2));
1259 Swap (el.PNum(3), el.PNum(4));
1260 }
1261
1262 mesh.AddSurfaceElement (el);
1263 // (*testout) << "(id nr "<< nr <<") add rect element: "
1264 // << mesh.Point (el.PNum(1)) << " - "
1265 // << mesh.Point (el.PNum(2)) << " - "
1266 // << mesh.Point (el.PNum(3)) << " - "
1267 // << mesh.Point (el.PNum(4)) << endl;
1268 found = true;
1269 //foundseg[i1]=foundseg[i2] = true;
1270 cntquads++;
1271 }
1272 }
1273 }
1274 if (found)
1275 {
1276 PrintMessage(3, "insert quad layer of ", cntquads,
1277 " elements at face ", segs.Get(1).si);
1278 //ARRAY<Segment> aux;
1279 //for(int i=0; i<segs.Size();i++)
1280 // if(!foundseg[i])
1281 // aux.Append(segs[i]);
1282 segs.SetSize(0);
1283 }
1284 else
1285 {
1286 BuildSurfaceElements2 (segs, mesh, surf);
1287 }
1288 }
1289
1290
1291
1292
1293
1294
1295 void CloseSurfaceIdentification ::
BuildSurfaceElements2(ARRAY<Segment> & segs,Mesh & mesh,const Surface * surf)1296 BuildSurfaceElements2 (ARRAY<Segment> & segs,
1297 Mesh & mesh, const Surface * surf)
1298 {
1299 // copy mesh
1300
1301
1302 // (*testout) << "copy trig face, identnr = " << nr << endl;
1303 // (*testout) << "identfaces = " << endl << identfaces << endl;
1304
1305 if (!segs.Size()) return;
1306
1307 bool found = 0;
1308
1309 int fother;
1310 int facei = segs.Get(1).si;
1311 int surfnr = mesh.GetFaceDescriptor(facei).SurfNr();
1312
1313
1314 bool foundid = 0;
1315 for (INDEX_2_HASHTABLE<int>::Iterator it = identfaces.Begin();
1316 it != identfaces.End(); it++)
1317 {
1318 INDEX_2 i2;
1319 int data;
1320 identfaces.GetData (it, i2, data);
1321 if (i2.I1() == facei || i2.I2() == facei)
1322 foundid = 1;
1323 }
1324
1325 /*
1326 for (int i = 1; i <= identfaces.GetNBags(); i++)
1327 for (int j = 1; j <= identfaces.GetBagSize(i); j++)
1328 {
1329 INDEX_2 i2;
1330 int data;
1331 identfaces.GetData (i, j, i2, data);
1332 if (i2.I1() == facei || i2.I2() == facei)
1333 foundid = 1;
1334
1335 (*testout) << "identface = " << i2 << endl;
1336 (*testout) << "face " << i2.I1() << " = " << mesh.GetFaceDescriptor(i2.I1()) << endl;
1337 (*testout) << "face " << i2.I2() << " = " << mesh.GetFaceDescriptor(i2.I2()) << endl;
1338 }
1339 */
1340
1341 if (foundid)
1342 {
1343 // (*testout) << "surfaces found" << endl;
1344 // copy surface
1345 for (int i = 1; i <= mesh.GetNSE(); i++)
1346 {
1347 const Element2d & sel = mesh.SurfaceElement(i);
1348 INDEX_2 fpair (facei, sel.GetIndex());
1349 fpair.Sort();
1350 if (identfaces.Used (fpair))
1351 {
1352 found = 1;
1353 fother = sel.GetIndex();
1354
1355 // copy element
1356 Element2d newel(sel.GetType());
1357 newel.SetIndex (facei);
1358 for (int k = 1; k <= sel.GetNP(); k++)
1359 {
1360 newel.PNum(k) =
1361 GetIdentifiedPoint (mesh, sel.PNum(k));
1362 // cout << "id-point = " << sel.PNum(k) << ", np = " << newel.PNum(k) << endl;
1363 }
1364
1365 Vec<3> nt = Cross (Point<3> (mesh.Point (newel.PNum(2)))-
1366 Point<3> (mesh.Point (newel.PNum(1))),
1367 Point<3> (mesh.Point (newel.PNum(3)))-
1368 Point<3> (mesh.Point (newel.PNum(1))));
1369 Vec<3> nsurf;
1370 nsurf = geom.GetSurface (surfnr)->GetNormalVector (mesh.Point(newel.PNum(1)));
1371 if (nsurf * nt < 0)
1372 Swap (newel.PNum(2), newel.PNum(3));
1373
1374 mesh.AddSurfaceElement (newel);
1375 }
1376 }
1377 }
1378
1379 if (found)
1380 {
1381 // (*mycout) << " copy face " << facei << " from face " << fother;
1382 PrintMessage (4, " copy face ", facei, " from face ", fother);
1383 segs.SetSize(0);
1384 }
1385 }
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400 void CloseSurfaceIdentification ::
BuildVolumeElements(ARRAY<class Element2d> & surfels,class Mesh & mesh)1401 BuildVolumeElements (ARRAY<class Element2d> & surfels,
1402 class Mesh & mesh)
1403 {
1404 ;
1405 }
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419 /* ***************** Close Edges Identification ********** */
1420
1421
1422
1423 CloseEdgesIdentification ::
CloseEdgesIdentification(int anr,const CSGeometry & ageom,const Surface * afacet,const Surface * as1,const Surface * as2)1424 CloseEdgesIdentification (int anr,
1425 const CSGeometry & ageom,
1426 const Surface * afacet,
1427 const Surface * as1,
1428 const Surface * as2)
1429 : Identification(anr, ageom)
1430 {
1431 facet = afacet;
1432 s1 = as1;
1433 s2 = as2;
1434 }
1435
~CloseEdgesIdentification()1436 CloseEdgesIdentification :: ~CloseEdgesIdentification ()
1437 {
1438 ;
1439 }
1440
Print(ostream & ost) const1441 void CloseEdgesIdentification :: Print (ostream & ost) const
1442 {
1443 ost << "CloseEdges Identifiaction, facet = "
1444 << facet->Name() << ", surfaces: "
1445 << s1->Name() << " - " << s2->Name() << endl;
1446 facet->Print (ost);
1447 s1->Print (ost);
1448 s2->Print (ost);
1449 ost << endl;
1450 }
1451
1452
GetData(ostream & ost) const1453 void CloseEdgesIdentification :: GetData (ostream & ost) const
1454 {
1455 ost << "closeedges " << facet->Name() << " "
1456 << s1->Name() << " " << s2->Name();
1457 }
1458
1459
1460 /*
1461 void CloseEdgesIdentification :: IdentifySpecialPoints
1462 (ARRAY<class SpecialPoint> & points)
1463 {
1464 int i, j;
1465 int bestj;
1466 double bestval, val;
1467
1468 for (i = 1; i <= points.Size(); i++)
1469 {
1470 Point<3> p1 = points.Get(i).p;
1471 Vec<3> n1;
1472
1473 if (!s1->PointOnSurface (p1))
1474 continue;
1475
1476 s1->GetNormalVector (p1, n1);
1477 n1 /= n1.Length();
1478 if ( fabs(n1 * points.Get(i).v) > 1e-3)
1479 continue;
1480
1481 bestval = 1e8;
1482 bestj = 1;
1483 for (j = 1; j <= points.Size(); j++)
1484 {
1485 Point<3> p2= points.Get(j).p;
1486 if (!s2->PointOnSurface (p2))
1487 continue;
1488
1489 Vec<3> n2;
1490 s2->GetNormalVector (p2, n2);
1491 n2 /= n2.Length();
1492 if ( fabs(n2 * points.Get(j).v) > 1e-3)
1493 continue;
1494
1495
1496 Vec<3> v(p1, p2);
1497 double vl = v.Length();
1498 double cl = fabs (v*n1);
1499
1500 val = 1 - cl*cl/(vl*vl);
1501
1502 val += (points.Get(i).v - points.Get(j).v).Length();
1503
1504 if (val < bestval)
1505 {
1506 bestj = j;
1507 bestval = val;
1508 }
1509 }
1510
1511 (*testout) << "Identify close surfaces special points: pi = "
1512 << points.Get(i).p << ", vi = " << points.Get(i).v
1513 << " pj = " << points.Get(bestj).p
1514 << ", vj = " << points.Get(bestj).v
1515 << " bestval = " << bestval << endl;
1516 }
1517 }
1518 */
1519
1520 int CloseEdgesIdentification ::
Identifyable(const SpecialPoint & sp1,const SpecialPoint & sp2,const TABLE<int> & specpoint2solid,const TABLE<int> & specpoint2surface) const1521 Identifyable (const SpecialPoint & sp1, const SpecialPoint & sp2,
1522 const TABLE<int> & specpoint2solid,
1523 const TABLE<int> & specpoint2surface) const
1524 {
1525 int i;
1526 double val;
1527
1528 SpecialPoint hsp1 = sp1;
1529 SpecialPoint hsp2 = sp2;
1530
1531 for (i = 1; i <= 1; i++)
1532 {
1533 if (!s1->PointOnSurface (hsp1.p))
1534 continue;
1535
1536 Vec<3> n1;
1537 n1 = s1->GetNormalVector (hsp1.p);
1538 n1 /= n1.Length();
1539 if ( fabs(n1 * hsp1.v) > 1e-3)
1540 continue;
1541
1542
1543 if (!s2->PointOnSurface(hsp2.p))
1544 continue;
1545
1546 Vec<3> n2;
1547 n2 = s2->GetNormalVector (hsp2.p);
1548 n2 /= n2.Length();
1549 if ( fabs(n2 * hsp2.v) > 1e-3)
1550 continue;
1551
1552
1553 Vec<3> v = hsp2.p - hsp1.p;
1554 double vl = v.Length();
1555 double cl = fabs (v*n1);
1556
1557
1558 val = 1 - cl*cl/(vl*vl);
1559 val += (hsp1.v - hsp2.v).Length();
1560
1561 if (val < 1e-3)
1562 {
1563 return 1;
1564 }
1565 }
1566
1567 return 0;
1568 }
1569
1570
1571
1572
IdentifyPoints(Mesh & mesh)1573 void CloseEdgesIdentification :: IdentifyPoints (Mesh & mesh)
1574 {
1575 int i, j;
1576 int i1, i2;
1577
1578 int np = mesh.GetNP();
1579 for (i1 = 1; i1 <= np; i1++)
1580 for (i2 = 1; i2 <= np; i2++)
1581 {
1582 if (i2 == i1)
1583 continue;
1584
1585 const Point<3> p1 = mesh.Point(i1);
1586 const Point<3> p2 = mesh.Point(i2);
1587 Point<3> pp1 = p1;
1588 Point<3> pp2 = p2;
1589
1590 s1->Project (pp1);
1591 facet->Project (pp1);
1592 s2->Project (pp2);
1593 facet->Project (pp2);
1594
1595 if (Dist (p1, pp1) > 1e-6 || Dist (p2, pp2) > 1e-6)
1596 continue;
1597
1598 Vec<3> n1, nf, t;
1599 Vec<3> n = p2 - p1;
1600 n.Normalize();
1601
1602 n1 = s1->GetNormalVector (p1);
1603 nf = facet->GetNormalVector (p1);
1604 t = Cross (n1, nf);
1605 t /= t.Length();
1606
1607 if (fabs (n * t) < 0.5)
1608 {
1609 (*testout) << "close edges identify points " << p1 << " - " << p2 << endl;
1610 mesh.GetIdentifications().Add (i1, i2, nr);
1611 mesh.GetIdentifications().SetType(nr,Identifications::CLOSEEDGES);
1612 }
1613 }
1614 }
1615
1616 void CloseEdgesIdentification ::
BuildSurfaceElements(ARRAY<Segment> & segs,Mesh & mesh,const Surface * surf)1617 BuildSurfaceElements (ARRAY<Segment> & segs,
1618 Mesh & mesh, const Surface * surf)
1619 {
1620 int i1, i2;
1621 int found = 0;
1622 int i, j, k;
1623
1624 if (surf != facet)
1625 return;
1626
1627 for (i1 = 1; i1 <= segs.Size(); i1++)
1628 for (i2 = 1; i2 < i1; i2++)
1629 {
1630 const Segment & s1 = segs.Get(i1);
1631 const Segment & s2 = segs.Get(i2);
1632 if (mesh.GetIdentifications().Get (s1.p1, s2.p2) &&
1633 mesh.GetIdentifications().Get (s1.p2, s2.p1))
1634 {
1635 Element2d el(QUAD);
1636 el.PNum(1) = s1.p1;
1637 el.PNum(2) = s1.p2;
1638 el.PNum(3) = s2.p2;
1639 el.PNum(4) = s2.p1;
1640
1641 Vec<3> n = Cross (Point<3> (mesh.Point(el.PNum(2)))-
1642 Point<3> (mesh.Point(el.PNum(1))),
1643 Point<3> (mesh.Point(el.PNum(3)))-
1644 Point<3> (mesh.Point(el.PNum(1))));
1645 Vec<3> ns;
1646 ns = surf->GetNormalVector (mesh.Point(el.PNum(1)));
1647 //(*testout) << "n = " << n << " ns = " << ns << endl;
1648 if (n * ns < 0)
1649 {
1650 //(*testout) << "Swap the quad" << endl;
1651 Swap (el.PNum(1), el.PNum(2));
1652 Swap (el.PNum(3), el.PNum(4));
1653 }
1654
1655
1656 Swap (el.PNum(3), el.PNum(4));
1657 mesh.AddSurfaceElement (el);
1658 // (*testout) << "add rect element: "
1659 // << mesh.Point (el.PNum(1)) << " - "
1660 // << mesh.Point (el.PNum(2)) << " - "
1661 // << mesh.Point (el.PNum(3)) << " - "
1662 // << mesh.Point (el.PNum(4)) << endl;
1663 found = 1;
1664 }
1665 }
1666
1667 if (found)
1668 segs.SetSize(0);
1669 }
1670
1671 }
1672
1673