← Back to branch summary

~ubuntu-branches/ubuntu/dapper/gmsh/dapper

« back to all changes in this revision

Viewing changes to Netgen/libsrc/meshing/delaunay.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2005-09-23 23:24:14 UTC
  • mfrom: (2.1.1 sarge)
  • Revision ID: james.westby@ubuntu.com-20050923232414-iao03ey38rd8pg0m
Tags: 1.60.1-1ubuntu2
gcc4 build failure - meshtype.hpp. (Closes BTS #324310)

Show diffs side-by-side

added added

removed removed

Lines of Context:
Netgen/libsrc/meshing/delaunay.cpp
2
2
#include "meshing.hpp"
3
3
 
4
4
 
5
 
// #define TEST 
6
 
// #define TEST2
7
5
 
8
6
namespace netgen
9
7
{
 
8
 
 
9
 
 
10
  static const int deltetfaces[][3] = 
 
11
    { { 1, 2, 3 },
 
12
      { 2, 0, 3 },
 
13
      { 0, 1, 3 },
 
14
      { 1, 0, 2 } };
 
15
 
 
16
 
10
17
  class DelaunayTet
11
18
  {
12
19
    PointIndex pnums[4];
27
34
        pnums[i] = el[i];
28
35
    }
29
36
    
30
 
    PointIndex & PNum(int i) { return pnums[i-1]; }
31
 
    PointIndex PNum(int i) const { return pnums[i-1]; }
32
 
 
33
37
    PointIndex & operator[] (int i) { return pnums[i]; }
34
38
    PointIndex operator[] (int i) const { return pnums[i]; }
35
39
 
36
 
    int & NB(int i) { return nb[i-1]; }
37
 
    int NB(int i) const { return nb[i-1]; }
38
 
 
39
 
    inline void GetFace (int i, INDEX_3 & face) const;
40
 
    inline int FaceNr (INDEX_3 & face) const;  // which face nr is it ?
41
 
    inline void GetFace (int i, Element2d & face) const;
 
40
    int & NB1(int i) { return nb[i-1]; }
 
41
    int NB1(int i) const { return nb[i-1]; }
 
42
 
 
43
    int & NB(int i) { return nb[i]; }
 
44
    int NB(int i) const { return nb[i]; }
 
45
 
 
46
 
 
47
    int FaceNr (INDEX_3 & face) const  // which face nr is it ?
 
48
    {
 
49
      for (int i = 0; i < 3; i++)
 
50
        if (pnums[i] != face.I1() && 
 
51
            pnums[i] != face.I2() && 
 
52
            pnums[i] != face.I3())
 
53
          return i;
 
54
      return 3;
 
55
    }
 
56
    
 
57
    void GetFace1 (int i, INDEX_3 & face) const
 
58
    {
 
59
      face.I(1) = pnums[deltetfaces[i-1][0]];
 
60
      face.I(2) = pnums[deltetfaces[i-1][1]];
 
61
      face.I(3) = pnums[deltetfaces[i-1][2]];
 
62
    }
 
63
 
 
64
    void GetFace (int i, INDEX_3 & face) const
 
65
    {
 
66
      face.I(1) = pnums[deltetfaces[i][0]];
 
67
      face.I(2) = pnums[deltetfaces[i][1]];
 
68
      face.I(3) = pnums[deltetfaces[i][2]];
 
69
    }
 
70
      
 
71
    INDEX_3 GetFace1 (int i) const
 
72
    {
 
73
      return INDEX_3 (pnums[deltetfaces[i-1][0]],
 
74
                      pnums[deltetfaces[i-1][1]],
 
75
                      pnums[deltetfaces[i-1][2]]);
 
76
    }
 
77
 
 
78
    INDEX_3 GetFace (int i) const
 
79
    {
 
80
      return INDEX_3 (pnums[deltetfaces[i][0]],
 
81
                      pnums[deltetfaces[i][1]],
 
82
                      pnums[deltetfaces[i][2]]);
 
83
    }
 
84
     
 
85
    void GetFace1 (int i, Element2d & face) const
 
86
    {
 
87
      // face.SetType(TRIG);
 
88
      face[0] = pnums[deltetfaces[i-1][0]];
 
89
      face[1] = pnums[deltetfaces[i-1][1]];
 
90
      face[2] = pnums[deltetfaces[i-1][2]];
 
91
    }
42
92
  };
43
93
 
44
94
 
45
95
 
46
 
  inline int DelaunayTet :: FaceNr (INDEX_3 & face) const
47
 
  {
48
 
    for (int i = 0; i < 3; i++)
49
 
      if (pnums[i] != face.I1() && 
50
 
          pnums[i] != face.I2() && 
51
 
          pnums[i] != face.I3())
52
 
        return i+1;
53
 
    return 4;
54
 
  }
55
 
 
56
 
  static const int deltetfaces[][3] = 
57
 
    { { 2, 3, 4 },
58
 
      { 3, 1, 4 },
59
 
      { 1, 2, 4 },
60
 
      { 2, 1, 3 } };
61
 
 
62
 
 
63
 
  inline void DelaunayTet :: GetFace (int i, INDEX_3 & face) const
64
 
  {
65
 
    face.I(1) = PNum(deltetfaces[i-1][0]);
66
 
    face.I(2) = PNum(deltetfaces[i-1][1]);
67
 
    face.I(3) = PNum(deltetfaces[i-1][2]);
68
 
  }
69
 
 
70
 
  inline void DelaunayTet :: GetFace (int i, Element2d & face) const
71
 
  {
72
 
    face.SetType(TRIG);
73
 
    face[0] = PNum(deltetfaces[i-1][0]);
74
 
    face[1] = PNum(deltetfaces[i-1][1]);
75
 
    face[2] = PNum(deltetfaces[i-1][2]);
76
 
  }
77
 
 
78
 
 
79
 
 
80
96
 
81
97
 
82
98
 
90
106
  {
91
107
    // face nodes -> one element
92
108
    INDEX_3_CLOSED_HASHTABLE<int> faces;
93
 
    //
 
109
 
 
110
    // 
94
111
    ARRAY<DelaunayTet> & tets;
 
112
 
95
113
  public:
96
114
 
97
115
    // estimated number of points
99
117
      : faces(200), tets(atets)
100
118
    { ; }
101
119
 
102
 
 
103
120
    // add element with 4 nodes
104
121
    void Add (int elnr);
105
122
 
106
123
    // delete element with 4 nodes
107
 
    void Delete (int elnr);
 
124
    void Delete (int elnr)
 
125
    {
 
126
      DelaunayTet & el = tets.Elem(elnr);
 
127
      for (int i = 0; i < 4; i++)
 
128
        faces.Set (el.GetFace(i).Sort(), el.NB(i));
 
129
    }
108
130
 
109
131
    // get neighbour of element elnr in direction fnr 
110
132
    int GetNB (int elnr, int fnr)
111
 
    { return tets.Get(elnr).NB(fnr); }
 
133
    { 
 
134
      return tets.Get(elnr).NB1(fnr); 
 
135
    }
112
136
 
113
137
    //
114
138
    void ResetFaceHT (int size)
115
139
    {
116
140
      faces.SetSize (size);
117
141
    }
118
 
    
119
 
    void PrintMemInfo (ostream & ost) const;
120
142
  };
121
143
 
122
144
 
123
145
 
124
146
  void MeshNB :: Add (int elnr)
125
147
  {
126
 
    INDEX_3 i3, i32; 
127
 
 
128
148
    DelaunayTet & el = tets.Elem(elnr);
129
149
 
130
 
    for (int i = 1; i <= 4; i++)
 
150
    for (int i = 0; i < 4; i++)
131
151
      {
132
 
        el.GetFace (i, i3);
133
 
        i3.Sort();
134
 
      
135
 
        int bnr, posnr, othertet;
 
152
        INDEX_3 i3 = INDEX_3::Sort (el.GetFace(i));
136
153
 
 
154
        int posnr;
 
155
        
137
156
        if (!faces.PositionCreate (i3, posnr))
138
157
          {
139
158
            // face already in use
140
 
            faces.GetData (posnr, othertet);
 
159
            int othertet = faces.GetData (posnr);
141
160
 
142
161
            el.NB(i) = othertet;
143
 
          
144
162
            if (othertet)
145
163
              {
146
164
                int fnr = tets.Get(othertet).FaceNr (i3);
156
174
  }
157
175
 
158
176
 
159
 
  void MeshNB :: Delete (int elnr)
160
 
  {
161
 
    INDEX_3 i3; 
162
 
    DelaunayTet & el = tets.Elem(elnr);
163
 
 
164
 
    for (int i = 1; i <= 4; i++)
165
 
      {
166
 
        el.GetFace (i, i3);
167
 
        i3.Sort();
168
 
        faces.Set (i3, el.NB(i));
169
 
      } 
170
 
  }
171
 
 
172
 
 
173
 
  void MeshNB :: PrintMemInfo (ostream & ost) const
174
 
  {
175
 
    /*
176
 
      int uf = 0;
177
 
      for (int i = 1; i <= face2el.Size(); i++)
178
 
      if (face2el.Get(i).I1())
179
 
      uf++;
180
 
 
181
 
      ost << "MeshNB: "
182
 
      << "validels = " << validels << " totels = " << el2face.Size() << endl
183
 
      << "validfaces = " << uf << " totfaces = " << face2el.Size() << endl
184
 
      << "face2el: " << face2el.Size() * sizeof(INDEX_2) << endl
185
 
      << "el2face: " << el2face.Size() * sizeof(INDEX_4) << endl;
186
 
    */
187
 
  }
188
 
 
189
177
 
190
178
 
191
179
 
197
185
  {
198
186
    ARRAY<int> links;
199
187
  public:
200
 
    SphereList () { ; }
201
 
    inline void AddElement (int elnr);
202
 
    inline void DeleteElement (int elnr);
203
 
    inline void ConnectElement (int eli, int toi);
 
188
    SphereList () 
 
189
    { ; }
 
190
 
 
191
    void AddElement (int elnr)
 
192
    {
 
193
      if (elnr > links.Size())
 
194
        links.Append (1);
 
195
      links.Elem(elnr) = elnr;
 
196
    }
 
197
 
 
198
    void DeleteElement (int elnr)
 
199
    {
 
200
      links.Elem(elnr) = 0;
 
201
    }    
 
202
  
 
203
    void ConnectElement (int eli, int toi)
 
204
    {
 
205
      links.Elem (eli) = links.Get (toi);
 
206
      links.Elem (toi) = eli;
 
207
    }
 
208
      
204
209
    void GetList (int eli, ARRAY<int> & linked) const;
205
210
  };
206
211
 
207
 
  inline void SphereList :: AddElement (int elnr)
208
 
  {
209
 
    if (elnr > links.Size())
210
 
      links.Append (1);
211
 
    links.Elem(elnr) = elnr;
212
 
  }
213
 
 
214
 
  inline void SphereList :: DeleteElement (int elnr)
215
 
  {
216
 
    links.Elem(elnr) = 0;
217
 
  }
218
 
 
219
 
  inline void SphereList :: ConnectElement (int eli, int toi)
220
 
  {
221
 
    links.Elem (eli) = links.Get (toi);
222
 
    links.Elem (toi) = eli;
223
 
  }
224
212
 
225
213
  void SphereList :: GetList (int eli, ARRAY<int> & linked) const
226
214
  {
262
250
                         IndexSet & insphere, IndexSet & closesphere)
263
251
  {
264
252
    int i, j, k, l;
265
 
 
266
 
#ifdef TEST2
267
 
    (*testout) << endl << "add point " << newp << endl;
268
 
#endif
269
253
  
270
254
    /*
271
255
      find any sphere, such that newp is contained in
369
353
              {
370
354
                int celind = connected.Get(k);
371
355
 
372
 
                if (tempels.Get(celind).PNum(1) != -1 && 
 
356
                if (tempels.Get(celind)[0] != -1 && 
373
357
                    !insphere.IsIn (celind))
374
358
                  {
375
359
                    changed = 1;
402
386
                    }
403
387
                  else
404
388
                    {
 
389
                      /*
405
390
                      Element2d face;
406
391
                      tempels.Get(helind).GetFace (k, face);
407
392
 
408
393
                      const Point3d & p1 = mesh.Point (face.PNum(1));
409
 
                      const Point3d & p2 = mesh.Point (face.PNum(2));
410
 
                      const Point3d & p3 = mesh.Point (face.PNum(3));
 
394
                      const Point3d & p2 = mesh.Point (face[1]);
 
395
                      const Point3d & p3 = mesh.Point (face[2]);
 
396
                      */
 
397
 
 
398
                      INDEX_3 i3 = tempels.Get(helind).GetFace (k-1);
 
399
 
 
400
                      const Point3d & p1 = mesh.Point ( PointIndex (i3.I1()));
 
401
                      const Point3d & p2 = mesh.Point ( PointIndex (i3.I2()));
 
402
                      const Point3d & p3 = mesh.Point ( PointIndex (i3.I3()));
 
403
 
411
404
 
412
405
                      Vec3d v1(p1, p2);
413
406
                      Vec3d v2(p1, p3);
414
407
                      Vec3d n = Cross (v1, v2);
415
408
                      n /= n.Length();
416
409
 
417
 
                      if (n * Vec3d (p1, mesh.Point (tempels.Get(helind).PNum(k))) > 0)
 
410
                      if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k-1])) > 0)
418
411
                        n *= -1;
419
412
 
420
413
                      double dist = n * Vec3d (p1, newp);
435
428
    //      (*testout) << "newels: " << endl;
436
429
    ARRAY<Element> newels;
437
430
 
 
431
    Element2d face(TRIG);
438
432
    for (j = 1; j <= insphere.Array().Size(); j++)
439
433
      for (k = 1; k <= 4; k++)
440
434
        {
444
438
 
445
439
          if (!nbind || !insphere.IsIn (nbind))
446
440
            {
447
 
              Element2d face;
448
 
              tempels.Get (celind).GetFace (k, face);
 
441
              tempels.Get (celind).GetFace1 (k, face);
449
442
                
450
 
              Element newel(4);
 
443
              Element newel(TET);
451
444
              for (l = 1; l <= 3; l++)
452
445
                newel.PNum(l) = face.PNum(l);
453
 
              newel.PNum(4) = newpi;
 
446
              newel[3] = newpi;
454
447
 
455
448
              newels.Append (newel);
456
449
 
457
 
              Vec3d v1(mesh.Point (face.PNum(1)), mesh.Point (face.PNum(2)));
458
 
              Vec3d v2(mesh.Point (face.PNum(1)), mesh.Point (face.PNum(3)));
 
450
              Vec3d v1(mesh.Point (face[0]), mesh.Point (face[1]));
 
451
              Vec3d v2(mesh.Point (face[0]), mesh.Point (face[2]));
459
452
              Vec3d n = Cross (v1, v2);
460
453
              n /= n.Length();
461
 
              if (n * Vec3d(mesh.Point (face.PNum(1)), 
462
 
                            mesh.Point (tempels.Get(insphere.Array().Get(j)).PNum(k)))
 
454
              if (n * Vec3d(mesh.Point (face[0]), 
 
455
                            mesh.Point (tempels.Get(insphere.Array().Get(j))[k-1]))
463
456
                  > 0)
464
457
                n *= -1;
465
458
 
466
 
              double hval = n * Vec3d (mesh.Point (face.PNum(1)), newp);
 
459
              double hval = n * Vec3d (mesh.Point (face[0]), newp);
467
460
                
468
461
              if (hval > -1e-12)
469
462
                {
471
464
                  (*testout) << "vec to outer, hval = " << hval << endl;
472
465
                  (*testout) << "v1 x v2 = " << Cross (v1, v2) << endl;
473
466
                  (*testout) << "facep: "
474
 
                             << mesh.Point (face.PNum(1)) << " "
475
 
                             << mesh.Point (face.PNum(2)) << " "
476
 
                             << mesh.Point (face.PNum(3)) << endl;
 
467
                             << mesh.Point (face[0]) << " "
 
468
                             << mesh.Point (face[1]) << " "
 
469
                             << mesh.Point (face[2]) << endl;
477
470
                }
478
471
            }
479
472
        }
488
481
        meshnb.Delete (celind); 
489
482
        list.DeleteElement (celind);
490
483
          
491
 
        for (k = 1; k <= 4; k++)
492
 
          tempels.Elem(celind).PNum(k) = -1;
 
484
        for (k = 0; k < 4; k++)
 
485
          tempels.Elem(celind)[k] = -1;
493
486
 
494
487
        ((ADTree6&)tettree.Tree()).DeleteElement (celind);
495
488
        freelist.Append (celind);
600
593
 
601
594
 
602
595
    // new: local box
603
 
    mesh.GetBox (pmax, pmin);
 
596
    mesh.GetBox (pmax, pmin);   // lower bound for pmax, upper for pmin
604
597
    for (i = 1; i <= adfront->GetNF(); i++)
605
598
      {
606
599
        const Element2d & face = adfront->GetFace(i);
611
604
          }
612
605
      }
613
606
  
 
607
    for (i = 0; i < mesh.LockedPoints().Size(); i++)
 
608
      {
 
609
        pmin.SetToMin (mesh.Point (mesh.LockedPoints()[i]));
 
610
        pmax.SetToMax (mesh.Point (mesh.LockedPoints()[i]));
 
611
      }
 
612
  
 
613
 
614
614
 
615
615
    Vec3d vdiag(pmin, pmax);
616
616
    // double r1 = vdiag.Length();
631
631
 
632
632
    int np = mesh.GetNP();
633
633
 
634
 
    startel.PNum(1) = mesh.AddPoint (cp1);
635
 
    startel.PNum(2) = mesh.AddPoint (cp2);
636
 
    startel.PNum(3) = mesh.AddPoint (cp3);
637
 
    startel.PNum(4) = mesh.AddPoint (cp4);
 
634
    startel[0] = mesh.AddPoint (cp1);
 
635
    startel[1] = mesh.AddPoint (cp2);
 
636
    startel[2] = mesh.AddPoint (cp3);
 
637
    startel[3] = mesh.AddPoint (cp4);
638
638
 
639
639
    // flag points to use for Delaunay:
640
640
    BitArrayChar<PointIndex::BASE> usep(np);
674
674
    ARRAY<int> connected, treesearch;
675
675
 
676
676
 
677
 
    tpmin = tpmax = mesh.Point(startel.PNum(1));
678
 
    for (k = 2; k <= 4; k++)
 
677
    tpmin = tpmax = mesh.Point(startel[0]);
 
678
    for (k = 1; k < 4; k++)
679
679
      {
680
 
        tpmin.SetToMin (mesh.Point (startel.PNum(k)));
681
 
        tpmax.SetToMax (mesh.Point (startel.PNum(k)));
 
680
        tpmin.SetToMin (mesh.Point (startel[k]));
 
681
        tpmax.SetToMax (mesh.Point (startel[k]));
682
682
      }
683
683
    tpmax = tpmax + 0.01 * (tpmax - tpmin);
684
684
    tettree.Insert (tpmin, tpmax, 1);
686
686
 
687
687
    Point3d pc;
688
688
          
689
 
    for (k = 1; k <= 4; k++)
690
 
      pp[k-1] = &mesh.Point (startel.PNum(k));
 
689
    for (k = 0; k < 4; k++)
 
690
      pp[k] = &mesh.Point (startel[k]);
691
691
  
692
692
    CalcSphereCenter (&pp[0], pc);
693
693
 
748
748
 
749
749
 
750
750
    for (i = tempels.Size(); i >= 1; i--)
751
 
      if (tempels.Get(i).PNum(1) <= 0)
 
751
      if (tempels.Get(i)[0] <= 0)
752
752
        tempels.DeleteElement (i);
753
753
 
754
754
    PrintDot ('\n');
816
816
      for (i = 1; i <= tempels.Size(); i++)
817
817
        {   
818
818
          Element el(4);
819
 
          for (j = 1; j <= 4; j++)
820
 
            el.PNum(j) = tempels.Elem(i).PNum(j);
 
819
          for (j = 0; j < 4; j++)
 
820
            el[j] = tempels.Elem(i)[j];
821
821
 
822
822
          el.SetIndex (1);
823
 
          const Point3d & lp1 = mesh.Point (el.PNum(1));
824
 
          const Point3d & lp2 = mesh.Point (el.PNum(2));
825
 
          const Point3d & lp3 = mesh.Point (el.PNum(3));
826
 
          const Point3d & lp4 = mesh.Point (el.PNum(4));
 
823
 
 
824
          const Point3d & lp1 = mesh.Point (el[0]);
 
825
          const Point3d & lp2 = mesh.Point (el[1]);
 
826
          const Point3d & lp3 = mesh.Point (el[2]);
 
827
          const Point3d & lp4 = mesh.Point (el[3]);
827
828
          Vec3d v1(lp1, lp2);
828
829
          Vec3d v2(lp1, lp3);
829
830
          Vec3d v3(lp1, lp4);
831
832
          double vol = n * v3;
832
833
        
833
834
          if (vol > 0)
834
 
            Swap (el.PNum(3), el.PNum(4));
 
835
            swap (el[2], el[3]);
835
836
        
836
837
          tempmesh.AddVolumeElement (el);
837
838
        }
838
839
 
 
840
 
839
841
      MeshQuality3d (tempmesh);
840
842
    
841
843
      for (i = 1; i <= mesh.GetNOpenElements(); i++)
843
845
          Element2d sel = mesh.OpenElement(i);
844
846
          sel.SetIndex(1);
845
847
          tempmesh.AddSurfaceElement (sel);
846
 
          Swap (sel.PNum(2), sel.PNum(3));
 
848
          swap (sel[1], sel[2]);
847
849
          tempmesh.AddSurfaceElement (sel);
848
850
        }
849
851
 
850
852
 
851
853
      for (i = 1; i <= 4; i++)
852
854
        {
853
 
          Element2d self(3);
 
855
          Element2d self(TRIG);
854
856
          self.SetIndex (1);
855
 
          startel.GetFace (i, self);
 
857
          startel.GetFace1 (i, self);
856
858
          tempmesh.AddSurfaceElement (self);
857
859
        }
858
860
 
859
861
      
860
862
      tempmesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
861
863
      tempmesh.AddFaceDescriptor (FaceDescriptor (2, 1, 0, 0));
862
 
    
 
864
 
863
865
 
864
866
      //  for (i = mesh.GetNP() - 3; i <= mesh.GetNP(); i++)
865
867
      //    tempmesh.AddLockedPoint (i);
899
901
 
900
902
 
901
903
 
902
 
 
903
904
    // remove degenerated
904
905
 
905
906
    BitArray badnode(mesh.GetNP());
908
909
    for (i = 1; i <= tempels.Size(); i++)
909
910
      {
910
911
        Element el(4);
911
 
        for (j = 1; j <= 4; j++)
912
 
          el.PNum(j) = tempels.Elem(i).PNum(j);
 
912
        for (j = 0; j < 4; j++)
 
913
          el[j] = tempels.Elem(i)[j];
913
914
        //      Element & el = tempels.Elem(i);
914
 
        const Point3d & lp1 = mesh.Point (el.PNum(1));
915
 
        const Point3d & lp2 = mesh.Point (el.PNum(2));
916
 
        const Point3d & lp3 = mesh.Point (el.PNum(3));
917
 
        const Point3d & lp4 = mesh.Point (el.PNum(4));
 
915
        const Point3d & lp1 = mesh.Point (el[0]);
 
916
        const Point3d & lp2 = mesh.Point (el[1]);
 
917
        const Point3d & lp3 = mesh.Point (el[2]);
 
918
        const Point3d & lp4 = mesh.Point (el[3]);
918
919
        Vec3d v1(lp1, lp2);
919
920
        Vec3d v2(lp1, lp3);
920
921
        Vec3d v3(lp1, lp4);
923
924
 
924
925
        double h = v1.Length() + v2.Length() + v3.Length();
925
926
        if (fabs (vol) < 1e-8 * (h * h * h) &&
926
 
            (el.PNum(1) <= np && el.PNum(2) <= np &&
927
 
             el.PNum(3) <= np && el.PNum(4) <= np) )   // old: 1e-12
 
927
            (el[0] <= np && el[1] <= np &&
 
928
             el[2] <= np && el[3] <= np) )   // old: 1e-12
928
929
          {
929
 
            badnode.Set(el.PNum(1));
930
 
            badnode.Set(el.PNum(2));
931
 
            badnode.Set(el.PNum(3));
932
 
            badnode.Set(el.PNum(4));
 
930
            badnode.Set(el[0]);
 
931
            badnode.Set(el[1]);
 
932
            badnode.Set(el[2]);
 
933
            badnode.Set(el[3]);
933
934
            ndeg++;
934
935
            (*testout) << "vol = " << vol << " h = " << h << endl;
935
936
          }
936
937
 
937
938
        if (vol > 0)
938
 
          Swap (el.PNum(3), el.PNum(4));
 
939
          Swap (el[2], el[3]);
939
940
      }
940
941
 
941
942
    ne = tempels.Size();
942
943
    for (i = ne; i >= 1; i--)
943
944
      {
944
945
        const DelaunayTet & el = tempels.Get(i);
945
 
        if (badnode.Test(el.PNum(1)) ||
946
 
            badnode.Test(el.PNum(2)) ||
947
 
            badnode.Test(el.PNum(3)) ||
948
 
            badnode.Test(el.PNum(4)) )
 
946
        if (badnode.Test(el[0]) ||
 
947
            badnode.Test(el[1]) ||
 
948
            badnode.Test(el[2]) ||
 
949
            badnode.Test(el[3]) )
949
950
          tempels.DeleteElement(i);
950
951
      }
951
952
 
959
960
    for (i = 1; i <= mesh.GetNOpenElements(); i++)
960
961
      {
961
962
        const Element2d & tri = mesh.OpenElement(i);
962
 
        INDEX_3 i3(tri.PNum(1), tri.PNum(2), tri.PNum(3));
 
963
        INDEX_3 i3(tri[0], tri[1], tri[2]);
963
964
        i3.Sort();
964
965
        openeltab.Set (i3, i);
965
966
      }
966
967
 
967
968
    for (i = 1; i <= tempels.Size(); i++)
968
969
      {
969
 
        for (j = 1; j <= 4; j++)
 
970
        for (j = 0; j < 4; j++)
970
971
          {
971
 
            Element2d face;
972
 
            tempels.Get(i).GetFace (j, face);
973
 
            INDEX_3 i3(face.PNum(1), face.PNum(2), face.PNum(3));
 
972
            INDEX_3 i3 = tempels.Get(i).GetFace (j);
974
973
            i3.Sort();
975
974
            if (openeltab.Used(i3))
976
975
              openeltab.Set (i3, 0);
1028
1027
          {
1029
1028
            switch (j)
1030
1029
              {
1031
 
              case 1: i2 = INDEX_2(el.PNum(1), el.PNum(2)); break;
1032
 
              case 2: i2 = INDEX_2(el.PNum(1), el.PNum(3)); break;
1033
 
              case 3: i2 = INDEX_2(el.PNum(1), el.PNum(4)); break;
1034
 
              case 4: i2 = INDEX_2(el.PNum(2), el.PNum(3)); break;
1035
 
              case 5: i2 = INDEX_2(el.PNum(2), el.PNum(4)); break;
1036
 
              case 6: i2 = INDEX_2(el.PNum(3), el.PNum(4)); break;
 
1030
              case 1: i2 = INDEX_2(el[0], el[1]); break;
 
1031
              case 2: i2 = INDEX_2(el[0], el[2]); break;
 
1032
              case 3: i2 = INDEX_2(el[0], el[3]); break;
 
1033
              case 4: i2 = INDEX_2(el[1], el[2]); break;
 
1034
              case 5: i2 = INDEX_2(el[1], el[3]); break;
 
1035
              case 6: i2 = INDEX_2(el[2], el[3]); break;
1037
1036
              }
1038
1037
            i2.Sort();
1039
1038
            tetedges.Set (i2, 1);
1042
1041
    //  cout << "tetedges:";
1043
1042
    //  tetedges.PrintMemInfo (cout);
1044
1043
 
 
1044
 
 
1045
    for (INDEX_2_HASHTABLE<INDEX_2>::Iterator it = twotrias.Begin();
 
1046
         it != twotrias.End(); it++)
 
1047
      {
 
1048
        INDEX_2 hi2, hi3;
 
1049
        twotrias.GetData (it, hi2, hi3);
 
1050
        hi3.Sort();
 
1051
        if (tetedges.Used (hi3))
 
1052
          {
 
1053
            const Point3d & p1 = mesh.Point ( PointIndex (hi2.I1()));
 
1054
            const Point3d & p2 = mesh.Point ( PointIndex (hi2.I2()));
 
1055
            const Point3d & p3 = mesh.Point ( PointIndex (hi3.I1()));
 
1056
            const Point3d & p4 = mesh.Point ( PointIndex (hi3.I2()));
 
1057
            Vec3d v1(p1, p2);
 
1058
            Vec3d v2(p1, p3);
 
1059
            Vec3d v3(p1, p4);
 
1060
            Vec3d n = Cross (v1, v2);
 
1061
            double vol = n * v3;
 
1062
            
 
1063
            double h = v1.Length() + v2.Length() + v3.Length();
 
1064
            if (fabs (vol) < 1e-4 * (h * h * h))   // old: 1e-12
 
1065
              {
 
1066
                badnode.Set(hi3.I1());  
 
1067
                badnode.Set(hi3.I2());  
 
1068
              }
 
1069
          }
 
1070
      }
 
1071
 
 
1072
    /*
1045
1073
    for (i = 1; i <= twotrias.GetNBags(); i++)
1046
1074
      for (j = 1; j <= twotrias.GetBagSize (i); j++)
1047
1075
        {
1068
1096
                }
1069
1097
            }
1070
1098
        }
1071
 
 
 
1099
    */
1072
1100
 
1073
1101
    ne = tempels.Size();
1074
1102
    for (i = ne; i >= 1; i--)
1075
1103
      {
1076
1104
        const DelaunayTet & el = tempels.Get(i);
1077
 
        if (badnode.Test(el.PNum(1)) ||
1078
 
            badnode.Test(el.PNum(2)) ||
1079
 
            badnode.Test(el.PNum(3)) ||
1080
 
            badnode.Test(el.PNum(4)) )
 
1105
        if (badnode.Test(el[0]) ||
 
1106
            badnode.Test(el[1]) ||
 
1107
            badnode.Test(el[2]) ||
 
1108
            badnode.Test(el[3]) )
1081
1109
          tempels.DeleteElement(i);
1082
1110
      }
1083
1111
 
1103
1131
              {
1104
1132
                const Element2d & tri = mesh.OpenElement(fnr);
1105
1133
              
1106
 
                Point3d ltpmin (mesh.Point(tri.PNum(1)));
 
1134
                Point3d ltpmin (mesh.Point(tri[0]));
1107
1135
                Point3d ltpmax (tpmin);
1108
1136
              
1109
1137
                for (k = 2; k <= 3; k++)
1124
1152
          
1125
1153
            int intersect = 0;
1126
1154
          
1127
 
            for (j = 1; j <= 4; j++)
 
1155
            for (j = 0; j < 4; j++)
1128
1156
              {
1129
 
                pp[j-1] = &mesh.Point(el.PNum(j));
1130
 
                tetpi[j-1] = el.PNum(j);
 
1157
                pp[j] = &mesh.Point(el[j]);
 
1158
                tetpi[j] = el[j];
1131
1159
              }
1132
1160
          
1133
1161
            Point3d tetpmin(*pp[0]);
1161
1189
              
1162
1190
                if (IntersectTetTriangle (&pp[0], &tripp[0], tetpi, tripi))
1163
1191
                  {
 
1192
                    /*
1164
1193
                    int il1, il2;
1165
1194
                    (*testout) << "intersect !" << endl;
1166
1195
                    (*testout) << "triind: ";
1180
1209
                    for (il2 = 0; il2 < 4; il2++)
1181
1210
                      (*testout) << " " << *pp[il2];
1182
1211
                    (*testout) << endl;
1183
 
              
 
1212
                    */
1184
1213
                  
1185
1214
                  
1186
1215
                    intersect = 1;
1220
1249
    for (i = 1; i <= mesh.GetNOpenElements(); i++)
1221
1250
      {
1222
1251
        const Element2d & tri = mesh.OpenElement(i);
1223
 
        INDEX_3 i3 (tri.PNum(1), tri.PNum(2), tri.PNum(3));
 
1252
        INDEX_3 i3 (tri[0], tri[1], tri[2]);
1224
1253
        i3.Sort();
1225
1254
        boundaryfaces.PrepareSet (i3);
1226
1255
      }
1228
1257
    for (i = 1; i <= mesh.GetNOpenElements(); i++)
1229
1258
      {
1230
1259
        const Element2d & tri = mesh.OpenElement(i);
1231
 
        INDEX_3 i3 (tri.PNum(1), tri.PNum(2), tri.PNum(3));
 
1260
        INDEX_3 i3 (tri[0], tri[1], tri[2]);
1232
1261
        i3.Sort();
1233
1262
        boundaryfaces.Set (i3, 1);
1234
1263
      }
1235
1264
 
1236
 
    for (i = 1; i <= tempels.Size(); i++)
1237
 
      for (j = 1; j <= 4; j++)
1238
 
        tempels.Elem(i).NB(j) = 0;
 
1265
    for (i = 0; i < tempels.Size(); i++)
 
1266
      for (j = 0; j < 4; j++)
 
1267
        tempels[i].NB(j) = 0;
1239
1268
  
1240
1269
    TABLE<int,PointIndex::BASE> elsonpoint(mesh.GetNP());
1241
1270
    for (i = 0; i < tempels.Size(); i++)
1263
1292
 
1264
1293
    INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100);   
1265
1294
  
1266
 
    Element2d hel(3);
 
1295
    Element2d hel(TRIG);
1267
1296
    for (pi = PointIndex::BASE; 
1268
1297
         pi < mesh.GetNP()+PointIndex::BASE; pi++)
1269
1298
      {
1275
1304
 
1276
1305
            for (j = 1; j <= 4; j++)
1277
1306
              {
1278
 
                el.GetFace (j, hel);
 
1307
                el.GetFace1 (j, hel);
1279
1308
                hel.Invert();
1280
1309
                hel.NormalizeNumbering();
1281
1310
              
1282
 
                if (hel.PNum(1) == pi)
 
1311
                if (hel[0] == pi)
1283
1312
                  {
1284
 
                    INDEX_3 i3(hel.PNum(1), hel.PNum(2), hel.PNum(3));
 
1313
                    INDEX_3 i3(hel[0], hel[1], hel[2]);
1285
1314
                  
1286
1315
                    if (!boundaryfaces.Used (i3))
1287
1316
                      {
1289
1318
                          {
1290
1319
                            INDEX_2 i2 = faceht.Get(i3);
1291
1320
                          
1292
 
                            tempels.Elem(i).NB(j) = i2.I1();
1293
 
                            tempels.Elem(i2.I1()).NB(i2.I2()) = i;
 
1321
                            tempels.Elem(i).NB1(j) = i2.I1();
 
1322
                            tempels.Elem(i2.I1()).NB1(i2.I2()) = i;
1294
1323
                          }
1295
1324
                        else
1296
1325
                          {
1297
1326
                            hel.Invert();
1298
1327
                            hel.NormalizeNumbering();
1299
 
                            INDEX_3 i3(hel.PNum(1), hel.PNum(2), hel.PNum(3));
 
1328
                            INDEX_3 i3(hel[0], hel[1], hel[2]);
1300
1329
                            INDEX_2 i2(i, j);
1301
1330
                            faceht.Set (i3, i2);
1302
1331
                          }
1314
1343
      {
1315
1344
      INDEX_3 i3;
1316
1345
      Element2d face;
1317
 
      el.GetFace (j, face);
 
1346
      el.GetFace1 (j, face);
1318
1347
      for (int kk = 1; kk <= 3; kk++)
1319
1348
      i3.I(kk) = face.PNum(kk);
1320
1349
 
1346
1375
      {
1347
1376
      (*testout) << i << " ";
1348
1377
      for (j = 1; j <= 4; j++)
1349
 
      (*testout) << tempels.Get(i).NB(j) << " ";
 
1378
      (*testout) << tempels.Get(i).NB1(j) << " ";
1350
1379
      (*testout) << endl;
1351
1380
      }
1352
1381
  
1415
1444
        if (done) break;
1416
1445
      
1417
1446
        const DelaunayTet & el = tempels.Get(i);
1418
 
        const Point3d & p1 = mesh.Point (el.PNum(1));
1419
 
        const Point3d & p2 = mesh.Point (el.PNum(2));
1420
 
        const Point3d & p3 = mesh.Point (el.PNum(3));
1421
 
        const Point3d & p4 = mesh.Point (el.PNum(4));
 
1447
        const Point3d & p1 = mesh.Point (el[0]);
 
1448
        const Point3d & p2 = mesh.Point (el[1]);
 
1449
        const Point3d & p3 = mesh.Point (el[2]);
 
1450
        const Point3d & p4 = mesh.Point (el[3]);
1422
1451
      
1423
1452
        Point3d ci = Center (p1, p2, p3, p4);
1424
1453
 
1449
1478
 
1450
1479
                for (j = 1; j <= 4; j++)
1451
1480
                  {
1452
 
                    INDEX_3 i3;
 
1481
                    INDEX_3 i3 = tempels.Get(ei).GetFace1(j);
 
1482
                    /*
1453
1483
                    Element2d face;
1454
1484
                    tempels.Get(ei).GetFace(j, face);
1455
1485
                    for (int kk = 1; kk <= 3; kk++)
1456
1486
                      i3.I(kk) = face.PNum(kk);
 
1487
                    */
1457
1488
                    i3.Sort();
1458
1489
                  
1459
1490
 
1460
 
                    if (tempels.Get(ei).NB(j))
1461
 
                      elstack.Append (tempels.Get(ei).NB(j));
 
1491
                    if (tempels.Get(ei).NB1(j))
 
1492
                      elstack.Append (tempels.Get(ei).NB1(j));
1462
1493
 
1463
1494
                    /*
1464
1495
                      if (innerfaces.Used(i3))
1466
1497
                      INDEX_2 i2 = innerfaces.Get(i3);
1467
1498
                      int other = i2.I1() + i2.I2() - ei;
1468
1499
 
1469
 
                      if (other != tempels.Get(ei).NB(j))
 
1500
                      if (other != tempels.Get(ei).NB1(j))
1470
1501
                      cerr << "different1 !!" << endl;
1471
1502
 
1472
1503
                      if (other)
1475
1506
                      }
1476
1507
                      }
1477
1508
                      else
1478
 
                      if (tempels.Get(ei).NB(j))
 
1509
                      if (tempels.Get(ei).NB1(j))
1479
1510
                      cerr << "different2 !!" << endl;
1480
1511
                    */
1481
1512
 
1492
1523
        for (i = 1; i <= ne; i++)
1493
1524
          {
1494
1525
            const DelaunayTet & el = tempels.Get(i);
1495
 
            const Point3d & p1 = mesh.Point (el.PNum(1));
1496
 
            const Point3d & p2 = mesh.Point (el.PNum(2));
1497
 
            const Point3d & p3 = mesh.Point (el.PNum(3));
1498
 
            const Point3d & p4 = mesh.Point (el.PNum(4));
 
1526
            const Point3d & p1 = mesh.Point (el[0]);
 
1527
            const Point3d & p2 = mesh.Point (el[1]);
 
1528
            const Point3d & p3 = mesh.Point (el[2]);
 
1529
            const Point3d & p4 = mesh.Point (el[3]);
1499
1530
          
1500
1531
            Point3d ci = Center (p1, p2, p3, p4);
1501
1532
          
1636
1667
    for (i = 1; i <= tempels.Size(); i++)
1637
1668
      {
1638
1669
        Element el(4);
1639
 
        for (j = 1; j <= 4; j++)
1640
 
          el.PNum(j) = tempels.Elem(i).PNum(j);
 
1670
        for (j = 0; j < 4; j++)
 
1671
          el[j] = tempels.Elem(i)[j];
1641
1672
        mesh.AddVolumeElement (el);
1642
1673
      }
1643
1674