1 #include <mystdlib.h>
2 #include <myadt.hpp>
3
4 #include <linalg.hpp>
5 #include <csg.hpp>
6 #include <meshing.hpp>
7
8 namespace netgen
9 {
10
SingularEdge(double abeta,int adomnr,const CSGeometry & ageom,const Solid * asol1,const Solid * asol2,double sf,const double maxh_at_initialization)11 SingularEdge :: SingularEdge (double abeta, int adomnr,
12 const CSGeometry & ageom,
13 const Solid * asol1,
14 const Solid * asol2, double sf,
15 const double maxh_at_initialization)
16 : geom(ageom), domnr(adomnr)
17 {
18 beta = abeta;
19 maxhinit = maxh_at_initialization;
20
21 if (beta > 1)
22 {
23 beta = 1;
24 cout << "Warning: beta set to 1" << endl;
25 }
26 if (beta <= 1e-3)
27 {
28 beta = 1e-3;
29 cout << "Warning: beta set to minimal value 0.001" << endl;
30 }
31
32 sol1 = asol1;
33 sol2 = asol2;
34 factor = sf;
35 }
36
FindPointsOnEdge(class Mesh & mesh)37 void SingularEdge :: FindPointsOnEdge (class Mesh & mesh)
38 {
39 (*testout) << "find points on edge" << endl;
40 points.SetSize(0);
41 segms.SetSize(0);
42
43
44 ARRAY<int> si1, si2;
45 sol1->GetSurfaceIndices (si1);
46 sol2->GetSurfaceIndices (si2);
47
48 for (int i = 0; i < si1.Size(); i++)
49 si1[i] = geom.GetSurfaceClassRepresentant(si1[i]);
50 for (int i = 0; i < si2.Size(); i++)
51 si2[i] = geom.GetSurfaceClassRepresentant(si2[i]);
52
53
54 for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
55 {
56 INDEX_2 i2 (mesh[si].p1, mesh[si].p2);
57 /*
58
59 bool onedge = 1;
60 for (j = 1; j <= 2; j++)
61 {
62 const Point<3> p = mesh[ PointIndex (i2.I(j)) ];
63 if (sol1->IsIn (p, 1e-3) && sol2->IsIn(p, 1e-3) &&
64 !sol1->IsStrictIn (p, 1e-3) && !sol2->IsStrictIn(p, 1e-3))
65 {
66 ;
67 }
68 else
69 onedge = 0;
70 }
71 */
72
73 if (domnr != -1 && domnr != mesh[si].domin && domnr != mesh[si].domout)
74 continue;
75
76 /*
77 bool onedge = 1;
78 for (int j = 0; j < 2; j++)
79 {
80 int surfi = (j == 0) ? mesh[si].surfnr1 : mesh[si].surfnr2;
81 surfi = geom.GetSurfaceClassRepresentant(surfi);
82 if (!si1.Contains(surfi) && !si2.Contains(surfi))
83 onedge = 0;
84 }
85 */
86 int surfi1 = geom.GetSurfaceClassRepresentant(mesh[si].surfnr1);
87 int surfi2 = geom.GetSurfaceClassRepresentant(mesh[si].surfnr2);
88
89 if (si1.Contains(surfi1) && si2.Contains(surfi2) ||
90 si1.Contains(surfi2) && si2.Contains(surfi1))
91
92 // if (onedge)
93 {
94 segms.Append (i2);
95 // PrintMessage (5, "sing segment ", i2.I1(), " - ", i2.I2());
96 points.Append (mesh[ PointIndex (i2.I1())]);
97 points.Append (mesh[ PointIndex (i2.I2())]);
98 mesh[si].singedge_left = factor;
99 mesh[si].singedge_right = factor;
100 }
101 }
102
103 /*
104 (*testout) << "Singular edge points:" << endl;
105 for (int i = 0; i < points.Size(); i++)
106 (*testout) << points[i] << endl;
107 */
108
109 }
110
SetMeshSize(class Mesh & mesh,double globalh)111 void SingularEdge :: SetMeshSize (class Mesh & mesh, double globalh)
112 {
113 double hloc = pow (globalh, 1/beta);
114 if(maxhinit > 0 && maxhinit < hloc)
115 {
116 hloc = maxhinit;
117 if(points.Size() > 1)
118 {
119 for (int i = 0; i < points.Size()-1; i++)
120 mesh.RestrictLocalHLine(points[i],points[i+1],hloc);
121 }
122 else
123 {
124 for (int i = 0; i < points.Size(); i++)
125 mesh.RestrictLocalH (points[i], hloc);
126 }
127 }
128 else
129 {
130 for (int i = 0; i < points.Size(); i++)
131 mesh.RestrictLocalH (points[i], hloc);
132 }
133 }
134
135
136
SingularPoint(double abeta,const Solid * asol1,const Solid * asol2,const Solid * asol3,double sf)137 SingularPoint :: SingularPoint (double abeta,
138 const Solid * asol1,
139 const Solid * asol2,
140 const Solid * asol3, double sf)
141 {
142 beta = abeta;
143 sol1 = asol1;
144 sol2 = asol2;
145 sol3 = asol3;
146 factor = sf;
147 }
148
149
FindPoints(class Mesh & mesh)150 void SingularPoint :: FindPoints (class Mesh & mesh)
151 {
152 points.SetSize(0);
153 ARRAY<int> surfk, surf;
154
155
156 for (PointIndex pi = PointIndex::BASE;
157 pi < mesh.GetNP()+PointIndex::BASE; pi++)
158 {
159 if (mesh[pi].Type() != FIXEDPOINT) continue;
160 const Point<3> p = mesh[pi];
161
162 (*testout) << "check singular point" << p << endl;
163
164 if (sol1->IsIn (p) && sol2->IsIn(p) && sol3->IsIn(p) &&
165 !sol1->IsStrictIn (p) && !sol2->IsStrictIn(p) && !sol3->IsStrictIn(p))
166 {
167 surf.SetSize (0);
168 for (int k = 1; k <= 3; k++)
169 {
170 const Solid * solk(NULL);
171 Solid *tansol;
172 switch (k)
173 {
174 case 1: solk = sol1; break;
175 case 2: solk = sol2; break;
176 case 3: solk = sol3; break;
177 }
178
179 solk -> TangentialSolid (p, tansol, surfk, 1e-3);
180 (*testout) << "Tansol = " << *tansol << endl;
181
182 if (!tansol) continue;
183
184 ReducePrimitiveIterator rpi(Box<3> (p-Vec<3> (1e-3,1e-3,1e-3),
185 p+Vec<3> (1e-3,1e-3,1e-3)));
186 UnReducePrimitiveIterator urpi;
187
188 tansol -> IterateSolid (rpi);
189 tansol->GetSurfaceIndices (surfk);
190 tansol -> IterateSolid (urpi);
191
192 (*testout) << "surfinds = " << surfk << endl;
193
194 for (int i = 0; i < surfk.Size(); i++)
195 if (!surf.Contains (surfk[i]))
196 surf.Append (surfk[i]);
197
198 delete tansol;
199 }
200
201 if (surf.Size() < 3) continue;
202
203 points.Append (p);
204 PrintMessage (5, "Point (", p(0), ", ", p(1), ", ", p(2), ") is singular");
205 mesh[pi].Singularity(factor);
206 }
207 }
208 }
209
210
SetMeshSize(class Mesh & mesh,double globalh)211 void SingularPoint :: SetMeshSize (class Mesh & mesh, double globalh)
212 {
213 double hloc = pow (globalh, 1/beta);
214 for (int i = 1; i <= points.Size(); i++)
215 mesh.RestrictLocalH (points.Get(i), hloc);
216 }
217 }
218