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