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