1 // Author : Gaetan Bricteux
2
3 #ifndef INTEGRATION3D_CC
4 #define INTEGRATION3D_CC
5
6 #include "Integration3D.h"
7 #include "recurCut.h"
8 #include "GaussIntegration.h"
9 #include "polynomialBasis.h"
10 #include "GmshDefines.h"
11
12 #define ZERO_LS_TOL 1.e-9
13 #define EQUALITY_TOL 1.e-15
14
15 // cross product
cross(const double * v1,const double * v2,double * v)16 inline void cross(const double *v1, const double *v2, double *v) {
17 v[0] = v1[1] * v2[2] - v2[1] * v1[2];
18 v[1] = v2[0] * v1[2] - v1[0] * v2[2];
19 v[2] = v1[0] * v2[1] - v2[0] * v1[1];
20 }
21 // dot product
dot(const double * v1,const double * v2)22 inline double dot(const double *v1, const double *v2) {
23 return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
24 }
25 // norm
norm(const double * v)26 inline double norm(const double *v) {
27 return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
28 }
29 // substract
vec(const double * x1,const double * x2,double * v)30 inline void vec (const double *x1, const double *x2, double *v) {
31 v[0] = x2[0] - x1[0]; v[1] = x2[1] - x1[1]; v[2] = x2[2] - x1[2];
32 }
vec(const DI_Point * p1,const DI_Point * p2,double * v)33 inline void vec (const DI_Point *p1, const DI_Point *p2, double *v) {
34 v[0] = p2->x() - p1->x(); v[1] = p2->y() - p1->y(); v[2] = p2->z() - p1->z();
35 }
sq2(const double a)36 inline double sq2 (const double a) {return a * a;}
37
38 // determinant
det3(double d11,double d12,double d13,double d21,double d22,double d23,double d31,double d32,double d33)39 inline double det3(double d11, double d12, double d13,
40 double d21, double d22, double d23,
41 double d31, double d32, double d33) {
42 return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22);
43 }
det4(double d11,double d12,double d13,double d14,double d21,double d22,double d23,double d24,double d31,double d32,double d33,double d34,double d41,double d42,double d43,double d44)44 inline double det4 (double d11, double d12, double d13, double d14,
45 double d21, double d22, double d23, double d24,
46 double d31, double d32, double d33, double d34,
47 double d41, double d42, double d43, double d44) {
48 return d11 * det3(d22, d23, d24, d32, d33, d34, d42, d43, d44)
49 - d21 * det3(d12, d13, d14, d32, d33, d34, d42, d43, d44)
50 + d31 * det3(d12, d13, d14, d22, d23, d24, d42, d43, d44)
51 - d41 * det3(d12, d13, d14, d22, d23, d24, d32, d33, d34);
52 }
53
54 // distance
distance(const DI_Point * p1,const DI_Point * p2)55 inline double distance(const DI_Point *p1, const DI_Point *p2) {
56 return sqrt((p1->x() - p2->x()) * (p1->x() - p2->x())
57 + (p1->y() - p2->y()) * (p1->y() - p2->y())
58 + (p1->z() - p2->z()) * (p1->z() - p2->z()));
59 }
distance(const DI_CuttingPoint * p1,const DI_CuttingPoint * p2)60 inline double distance(const DI_CuttingPoint *p1, const DI_CuttingPoint *p2) {
61 return sqrt((p1->x() - p2->x()) * (p1->x() - p2->x())
62 + (p1->y() - p2->y()) * (p1->y() - p2->y())
63 + (p1->z() - p2->z()) * (p1->z() - p2->z()));
64 }
65
66 // middle of 2 points
middle(const DI_Point * p1,const DI_Point * p2)67 inline DI_Point* middle (const DI_Point *p1, const DI_Point *p2) {
68 return new DI_Point ((p1->x() + p2->x()) / 2, (p1->y() + p2->y()) / 2, (p1->z() + p2->z()) / 2);
69 }
middle(const DI_Point * p1,const DI_Point * p2,const DI_Element * e,const std::vector<gLevelset * > & RPNi)70 inline DI_Point* middle (const DI_Point *p1, const DI_Point *p2,
71 const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
72 return new DI_Point ((p1->x() + p2->x()) / 2, (p1->y() + p2->y()) / 2, (p1->z() + p2->z()) / 2, e, RPNi);
73 }
74
75 // barycentre
barycenter(const DI_Element * el,const DI_Element * e,const std::vector<gLevelset * > & RPN)76 inline DI_Point* barycenter (const DI_Element *el,
77 const DI_Element *e, const std::vector<gLevelset *> &RPN) {
78 double x[3] = {0., 0., 0.};
79 int n;
80 for (n = 0; n < el->nbVert(); n++) {
81 x[0] += el->x(n);
82 x[1] += el->y(n);
83 x[2] += el->z(n);
84 }
85 return new DI_Point(x[0] / n, x[1] / n, x[2] / n, e, RPN);
86 }
87
88 // swap two integers
swap(int & a,int & b)89 inline void swap(int &a, int &b) {
90 int t = a;
91 a = b; b = t;
92 }
93 // swap two doubles
swap(double & a,double & b)94 inline void swap(double &a, double &b) {
95 double t = a;
96 a = b; b = t;
97 }
98 // swap two points
swap(DI_Point ** p1,DI_Point ** p2)99 inline void swap (DI_Point **p1, DI_Point **p2) {
100 DI_Point *pt = *p1;
101 *p1 = *p2; *p2 = pt;
102 }
103
104 // return true if the 4 points are planar
isPlanar(const DI_Point * p1,const DI_Point * p2,const DI_Point * p3,const DI_Point * p4)105 inline bool isPlanar (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3, const DI_Point *p4) {
106 double v12[3]; vec(p1, p2, v12);
107 double v13[3]; vec(p1, p3, v13);
108 double v14[3]; vec(p1, p4, v14);
109 double c1[3]; cross(v12, v13, c1);
110 double c2[3]; cross(v12, v14, c2);
111 double c[3]; cross(c1, c2, c);
112 return (c[0] == 0. && c[1] == 0. && c[2] == 0.);
113 }
114 // return true if the 3 points are linear
isLinear(const DI_Point * p1,const DI_Point * p2,const DI_Point * p3)115 inline bool isLinear (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3) {
116 double v12[3]; vec(p1, p2, v12);
117 double v13[3]; vec(p1, p3, v13);
118 double c[3]; cross(v12, v13, c);
119 return (c[0] == 0. && c[1] == 0. && c[2] == 0.);
120 }
121 // return true if the 4 points form a quad in good ordering
ordered4Nodes(const DI_Point * p1,const DI_Point * p2,const DI_Point * p3,const DI_Point * p4)122 inline bool ordered4Nodes (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3, const DI_Point *p4) {
123 double v21[3]; vec(p2, p1, v21);
124 double v23[3]; vec(p2, p3, v23);
125 double v24[3]; vec(p2, p4, v24);
126 double norm21 = norm(v21);
127 double norm23 = norm(v23);
128 double norm24 = norm(v24);
129 double dot213 = dot(v21, v23);
130 double dot214 = dot(v21, v24);
131 double theta213 = acos(dot213 / (norm21 * norm23));
132 double theta214 = acos(dot214 / (norm21 * norm24));
133 if(theta214 > theta213) return false;
134 return true;
135 }
136
isInQE(const DI_Triangle * t,const DI_QualError * QE)137 bool isInQE (const DI_Triangle *t, const DI_QualError *QE) {
138 int b = 0;
139 for(int i = 0; i < 3; i++) {
140 for(int j = 0; j < 4; j++)
141 if (t->pt(i)->equal(QE->pt(j))) {b++; break;}
142 }
143 return (b == 3);
144 }
belongsTo(const DI_Element * e,const DI_Element * E)145 bool belongsTo (const DI_Element *e, const DI_Element *E) {
146 int b = 0;
147 for(int k = 0; k < E->nbVert(); k++) {
148 for(int l = 0; l < e->nbVert(); l++)
149 if((e->pt(l))->equal(E->pt(k))) {b++; break;}
150 if(b == e->nbVert()) return true;
151 }
152 return false;
153 }
isLastLnInV(const std::vector<DI_Line * > & v,const int i1=0)154 bool isLastLnInV (const std::vector<DI_Line *> &v, const int i1 = 0) {
155 int b;
156 for (int i = i1; i < (int)v.size() - 1; i++) {
157 b = 0;
158 for (int k = 0; k < 2; k++)
159 for(int l = 0; l < 2; l++)
160 if(v[i]->pt(k)->equal(v[v.size() - 1]->pt(l))) {b++; break;}
161 if(b == 2) return true;
162 }
163 return false;
164 }
isLastTrInV(const std::vector<DI_Triangle * > & v,const int i1=0)165 bool isLastTrInV (const std::vector<DI_Triangle *> &v, const int i1 = 0) {
166 int b;
167 for (int i = i1; i < (int)v.size() - 1; i++) {
168 b = 0;
169 for (int k = 0; k < 3; k++)
170 for(int l = 0; l < 3; l++)
171 if(v[i]->pt(k)->equal(v[v.size() - 1]->pt(l))) {b++; break;}
172 if(b == 3) return true;
173 }
174 return false;
175 }
isLastQInV(const std::vector<DI_Quad * > & v,const int i1=0)176 bool isLastQInV (const std::vector<DI_Quad *> &v, const int i1 = 0) {
177 int b;
178 for (int i = i1; i < (int)v.size() - 1; i++) {
179 b = 0;
180 for(int k = 0; k < 4; k++)
181 for(int l = 0; l < 4; l++)
182 if(v[i]->pt(k)->equal(v[v.size() - 1]->pt(l))) {b++; break;}
183 if(b == 4) return true;
184 }
185 return false;
186 }
187
188 // return common vertex of 2 edges (s11,s12) and (s21,s22) // s11 != s12 and s21 != s22
commonV(int & s11,int & s12,int & s21,int & s22)189 int commonV (int &s11, int &s12, int &s21, int &s22) {
190 if(s11 == s21 || s11 == s22) return s11;
191 if(s12 == s21 || s12 == s22) return s12;
192 printf("no common summit, %d,%d,%d,%d\n", s11, s12, s21, s22);
193 return 0;
194 }
195
adjustLs(double ls)196 inline double adjustLs (double ls) {
197 return (fabs(ls) < ZERO_LS_TOL) ? 0. : ls;
198 }
isCrossed(const DI_Point * p1,const DI_Point * p2)199 inline bool isCrossed (const DI_Point *p1, const DI_Point *p2) {
200 double ls1 = p1->ls(), ls2 = p2->ls();
201 return (ls1 * ls2 < 0.);
202 }
203
204 // return the index of the point with minimum x,y and z
minimum(double * x,double * y,double * z,const int num)205 int minimum(double *x, double *y, double *z, const int num) {
206 double xm = x[0];
207 for(int i = 1; i < num; i++) if(x[i] < xm) xm = x[i];
208 std::vector<int> INDx(num); int countx = 0;
209 for(int i = 0; i < num; i++) if(x[i] == xm) INDx[countx++] = i;
210 if(countx == 1) return INDx[0];
211
212 double ym = y[INDx[0]];
213 for(int i = 1; i < countx; i++) if(y[INDx[i]] < ym) ym = y[INDx[i]];
214 std::vector<int> INDy(countx); int county = 0;
215 for(int i = 0; i < countx; i++) if(y[INDx[i]] == ym) INDy[county++] = INDx[i];
216 if(county == 1) return INDy[0];
217
218 double zm = z[INDy[0]];
219 for(int i = 1; i < county; i++) if(z[INDy[i]] < zm) zm = z[INDy[i]];
220 std::vector<int> INDz(county); int countz = 0;
221 for(int i = 0; i < county; i++) if(z[INDy[i]] == zm) INDz[countz++] = INDy[i];
222 return INDz[0];
223 }
224
225 // Return the quality of a triangle
qualityTri(const DI_Point * p0,const DI_Point * p1,const DI_Point * p2)226 double qualityTri(const DI_Point *p0, const DI_Point *p1, const DI_Point *p2) {
227 double a = distance(p0, p1);
228 double b = distance(p0, p2);
229 double c = distance(p1, p2);
230
231 double rhoO = a * b * c / sqrt((a + b + c) * (b + c - a) * (c + a - b) * (a + b - c));
232 double rhoI = a * b * c / (2 * (a + b + c) * rhoO);
233 return 2 * rhoI / rhoO;
234 }
235
236 // Return the best quality for a quad cut into 2 triangles
bestQuality(const DI_Point * p1,const DI_Point * p2,const DI_Point * p3,const DI_Point * p4,DI_Triangle ** t1,DI_Triangle ** t2)237 int bestQuality(const DI_Point *p1, const DI_Point *p2, const DI_Point *p3, const DI_Point *p4,
238 DI_Triangle **t1, DI_Triangle **t2) {
239 int cut = 0;
240 double qual11 = qualityTri(p1, p2, p3);
241 double qual12 = qualityTri(p1, p3, p4);
242 double qual21 = qualityTri(p1, p2, p4);
243 double qual22 = qualityTri(p2, p3, p4);
244 double worst1 = std::min(qual11, qual12);
245 double worst2 = std::min(qual21, qual22);
246 if(worst1 - worst2 > EQUALITY_TOL)
247 cut = 1;
248 else if(worst2 - worst1 > EQUALITY_TOL)
249 cut = 2;
250 else {
251 worst1 = std::max(qual11, qual12);
252 worst2 = std::max(qual21, qual22);
253 if(worst1 - worst2 > EQUALITY_TOL)
254 cut = 1;
255 else if(worst2 - worst1 > EQUALITY_TOL)
256 cut = 2;
257 else { // same min and max qualities
258 double x[4] = {p1->x(), p2->x(), p3->x(), p4->x()};
259 double y[4] = {p1->y(), p2->y(), p3->y(), p4->y()};
260 double z[4] = {p1->z(), p2->z(), p3->z(), p4->z()};
261 int IND = minimum(x, y, z, 4);
262 if(IND == 0 || IND == 2) cut = 1;
263 else cut = 2;
264 }
265 }
266 if(cut == 1) {
267 *t1 = new DI_Triangle(p1, p2, p3);
268 *t2 = new DI_Triangle(p1, p3, p4);
269 return 1;
270 }
271 else {
272 *t1 = new DI_Triangle(p1, p2, p4);
273 *t2 = new DI_Triangle(p2, p3, p4);
274 return 2;
275 }
276 }
277
278 // Return the quality of a tetrahedron
qualityTet(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3,double x4,double y4,double z4)279 inline double qualityTet (double x1, double y1, double z1, double x2, double y2, double z2,
280 double x3, double y3, double z3, double x4, double y4, double z4){
281 double Dx = det4(x1 * x1 + y1 * y1 + z1 * z1, y1, z1, 1., x2 * x2 + y2 * y2 + z2 * z2, y2, z2, 1.,
282 x3 * x3 + y3 * y3 + z3 * z3, y3, z3, 1., x4 * x4 + y4 * y4 + z4 * z4, y4, z4, 1.);
283 double Dy = -det4(x1 * x1 + y1 * y1 + z1 * z1, x1, z1, 1., x2 * x2 + y2 * y2 + z2 * z2, x2, z2, 1.,
284 x3 * x3 + y3 * y3 + z3 * z3, x3, z3, 1., x4 * x4 + y4 * y4 + z4 * z4, x4, z4, 1.);
285 double Dz = det4(x1 * x1 + y1 * y1 + z1 * z1, x1, y1, 1., x2 * x2 + y2 * y2 + z2 * z2, x2, y2, 1.,
286 x3 * x3 + y3 * y3 + z3 * z3, x3, y3, 1., x4 * x4 + y4 * y4 + z4 * z4, x4, y4, 1.);
287 double a = det4(x1, y1, z1, 1., x2, y2, z2, 1., x3, y3, z3, 1., x4, y4, z4, 1.);
288 double c = det4(x1 * x1 + y1 * y1 + z1 * z1, x1, y1, z1, x2 * x2 + y2 * y2 + z2 * z2, x2, y2, z2,
289 x3 * x3 + y3 * y3 + z3 * z3, x3, y3, z3, x4 * x4 + y4 * y4 + z4 * z4, x4, y4, z4);
290 double rhoO = sqrt(Dx * Dx + Dy * Dy + Dz * Dz - 4. * a * c) / (2. * fabs(a));
291
292 double res[4][3];
293 double v12[3] = {x2 - x1, y2 - y1, z2 - z1};
294 double v13[3] = {x3 - x1, y3 - y1, z3 - z1};
295 double v14[3] = {x4 - x1, y4 - y1, z4 - z1};
296 double v32[3] = {x2 - x3, y2 - y3, z2 - z3};
297 double v43[3] = {x3 - x4, y3 - y4, z3 - z4};
298 double v34[3] = {x4 - x3, y4 - y3, z4 - z3};
299 cross(v12, v13, res[0]);
300 cross(v13, v14, res[1]);
301 cross(v14, v12, res[2]);
302 cross(v32, v43, res[3]);
303 double V = dot(res[0], v34) / 6.;
304 double SA = norm(res[0]) / 2. + norm(res[1]) / 2.
305 + norm(res[2]) / 2. + norm(res[3]) / 2.;
306 double rhoI = 3 * V / SA;
307 return 3 * rhoI / rhoO;
308 }
309
310 // Return the cutting of a pyramid cut into 2 tetras with the best quality faces
311 // (x1,x2,x3,x4 form the base in ccw and x5 is the summit)
bestQuality(const DI_Point * p1,const DI_Point * p2,const DI_Point * p3,const DI_Point * p4,const DI_Point * p5,DI_Tetra ** t1,DI_Tetra ** t2)312 int bestQuality(const DI_Point *p1, const DI_Point *p2, const DI_Point *p3,
313 const DI_Point *p4, const DI_Point *p5, DI_Tetra **t1, DI_Tetra **t2) {
314 int cut = 0;
315 double qual11 = qualityTri(p1, p2, p3);
316 double qual12 = qualityTri(p1, p3, p4);
317 double qual21 = qualityTri(p1, p2, p4);
318 double qual22 = qualityTri(p2, p3, p4);
319
320 double worst1 = std::min(qual11, qual12);
321 double worst2 = std::min(qual21, qual22);
322 if(worst1 - worst2 > EQUALITY_TOL)
323 cut = 1;
324 else if(worst2 - worst1 > EQUALITY_TOL)
325 cut = 2;
326 else {
327 worst1 = std::max(qual11, qual12);
328 worst2 = std::max(qual21, qual22);
329 if(worst1 - worst2 > EQUALITY_TOL)
330 cut = 1;
331 else if(worst2 - worst1 > EQUALITY_TOL)
332 cut = 2;
333 else { // same min and max qualities
334 double x[4] = {p1->x(), p2->x(), p3->x(), p4->x()};
335 double y[4] = {p1->y(), p2->y(), p3->y(), p4->y()};
336 double z[4] = {p1->z(), p2->z(), p3->z(), p4->z()};
337 int IND = minimum(x, y, z, 4);
338 if(IND == 0 || IND == 2) cut = 1;
339 else cut = 2;
340 }
341 }
342 if(cut == 1) {
343 *t1 = new DI_Tetra(p1, p2, p3, p5);
344 *t2 = new DI_Tetra(p1, p3, p4, p5);
345 return 1;
346 }
347 else {
348 *t1 = new DI_Tetra(p1, p2, p4, p5);
349 *t2 = new DI_Tetra(p2, p3, p4, p5);
350 return 2;
351 }
352 }
353
354 // Return the cutting of a prism into 3 tetras with the best quality faces
bestQuality(DI_Point * p0,DI_Point * p1,DI_Point * p2,DI_Point * p3,DI_Point * p4,DI_Point * p5,DI_Tetra ** t1,DI_Tetra ** t2,DI_Tetra ** t3,std::vector<DI_QualError * > & QError)355 int bestQuality (DI_Point *p0, DI_Point *p1, DI_Point *p2,
356 DI_Point *p3, DI_Point *p4, DI_Point *p5,
357 DI_Tetra **t1, DI_Tetra **t2, DI_Tetra **t3, std::vector<DI_QualError *> &QError) {
358 int cut3[3] = {0, 0, 0}; int cut = -1;
359 int fa[3][4] = {{0, 3, 4, 1}, {0, 2, 5, 3}, {1, 4, 5, 2}}; //faces
360 DI_Point* pt[6] = {p0, p1, p2, p3, p4, p5};
361 for(int f = 0; f < 3; f++){
362 int i1 = fa[f][0], i2 = fa[f][1], i3 = fa[f][2], i4 = fa[f][3];
363 double qual11 = qualityTri(pt[i1], pt[i2], pt[i3]);
364 double qual12 = qualityTri(pt[i1], pt[i3], pt[i4]);
365 double qual21 = qualityTri(pt[i1], pt[i2], pt[i4]);
366 double qual22 = qualityTri(pt[i2], pt[i3], pt[i4]);
367 //printf("face%d, q11=%.18f,q12=%.18f,q21=%.18f,q22=%.18f\n",f,qual11,qual12,qual21,qual22);
368 double worst1 = std::min(qual11, qual12);
369 double worst2 = std::min(qual21, qual22);
370 if(worst1 - worst2 > EQUALITY_TOL)
371 cut3[f] = 1;
372 else if(worst2 - worst1 > EQUALITY_TOL)
373 cut3[f] = 2;
374 else {
375 worst1 = std::max(qual11, qual12);
376 worst2 = std::max(qual21, qual22);
377 if(worst1 - worst2 > EQUALITY_TOL)
378 cut3[f] = 1;
379 else if(worst2 - worst1 > EQUALITY_TOL)
380 cut3[f] = 2;
381 else { // same min and max qualities
382 double xc[4] = {pt[i1]->x(), pt[i2]->x(), pt[i3]->x(), pt[i4]->x()};
383 double yc[4] = {pt[i1]->y(), pt[i2]->y(), pt[i3]->y(), pt[i4]->y()};
384 double zc[4] = {pt[i1]->z(), pt[i2]->z(), pt[i3]->z(), pt[i4]->z()};
385 int IND = minimum(xc, yc, zc, 4);
386 if(IND == 0 || IND == 2) cut3[f] = 1;
387 else cut3[f] = 2;
388 }
389 }
390 }
391
392 if(cut3[0] == 1){
393 if(cut3[1] == 1){
394 if(cut3[2] == 1) cut = 1;
395 else cut = 2;
396 }
397 else{
398 if(cut3[2] == 1) QError.push_back(new DI_QualError(p1, p5, p2, p4));
399 cut = 3;
400 }
401 }
402 else {
403 if(cut3[1] == 1){
404 if(cut3[2] == 2) QError.push_back(new DI_QualError(p1, p5, p2, p4));
405 cut = 0;
406 }
407 else{
408 if(cut3[2] == 1) cut = 5;
409 else cut = 4;
410 }
411 }
412
413 if(cut == 0) {
414 *t1 = new DI_Tetra(p0, p1, p2, p5);
415 *t2 = new DI_Tetra(p0, p1, p5, p3);
416 *t3 = new DI_Tetra(p1, p5, p3, p4);
417 return 1;
418 }
419 else if(cut == 1) {
420 *t1 = new DI_Tetra(p0, p1, p2, p5);
421 *t2 = new DI_Tetra(p0, p1, p5, p4);
422 *t3 = new DI_Tetra(p0, p4, p5, p3);
423 return 2;
424 }
425 else if(cut == 2) {
426 *t1 = new DI_Tetra(p0, p1, p2, p4);
427 *t2 = new DI_Tetra(p0, p4, p2, p5);
428 *t3 = new DI_Tetra(p0, p4, p5, p3);
429 return 3;
430 }
431 else if(cut == 3) {
432 *t1 = new DI_Tetra(p0, p1, p2, p4);
433 *t2 = new DI_Tetra(p0, p4, p2, p3);
434 *t3 = new DI_Tetra(p2, p3, p4, p5);
435 return 4;
436 }
437 else if(cut == 4) {
438 *t1 = new DI_Tetra(p0, p1, p2, p3);
439 *t2 = new DI_Tetra(p1, p2, p3, p4);
440 *t3 = new DI_Tetra(p2, p3, p4, p5);
441 return 5;
442 }
443 else if(cut == 5) {
444 *t1 = new DI_Tetra(p0, p1, p2, p3);
445 *t2 = new DI_Tetra(p1, p2, p3, p5);
446 *t3 = new DI_Tetra(p1, p5, p3, p4);
447 return 6;
448 }
449 return 0;
450 }
451
452 // computes the intersection between a level set and a linear edge
Newton(const DI_Point * p1,const DI_Point * p2,const DI_Element * e,const std::vector<gLevelset * > & RPNi)453 DI_Point* Newton(const DI_Point *p1, const DI_Point *p2,
454 const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
455 double eps = -p1->ls() / (p2->ls() - p1->ls());
456 DI_Point* p = new DI_Point(p1->x() + eps * (p2->x() - p1->x()), p1->y() + eps * (p2->y() - p1->y()),
457 p1->z() + eps * (p2->z() - p1->z()));
458 double pls = e->evalLs(p->x(), p->y(), p->z());
459 // Newton iteration to find the intersection between the level set and the edge
460 while(fabs(pls) > EQUALITY_TOL){
461 if(pls * p1->ls() < 0.) {
462 eps = -pls / (p1->ls() - pls);
463 p->move(p->x() + eps * (p1->x() - p->x()), p->y() + eps * (p1->y() - p->y()), p->z() + eps * (p1->z() - p->z()));
464 }
465 else {
466 eps = -pls / (p2->ls() - pls);
467 p->move(p->x() + eps * (p2->x() - p->x()), p->y() + eps * (p2->y() - p->y()), p->z() + eps * (p2->z() - p->z()));
468 }
469 pls = e->evalLs(p->x(), p->y(), p->z());
470 }
471 p->computeLs(e, RPNi);
472 return p;
473 }
474
475 // compute the position of the middle of a quadratic edge
476 //(intersection between the bisector of the linear edge and the levelset)
quadMidNode(const DI_Point * p1,const DI_Point * p2,const DI_Point * pf,const DI_Element * e,const std::vector<gLevelset * > & RPNi)477 DI_Point* quadMidNode(const DI_Point *p1, const DI_Point *p2, const DI_Point *pf,
478 const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
479 // middle of the edge between the 2 cutting points
480 DI_Point midEN((p1->x() + p2->x()) / 2., (p1->y() + p2->y()) / 2., (p1->z() + p2->z()) / 2.);
481 midEN.addLs(e);
482 // the bisector is computed in the plane of the face of the tetra [(x1,x2)x(x1,xf)]x(x1,x2)
483 double v12[3]; vec(p1, p2, v12);
484 double v1f[3]; vec(p1, pf, v1f);
485 double nf[3]; cross(v12, v1f, nf);
486 double bisector[3]; cross(nf, v12, bisector);
487 double normB = norm(bisector);
488 if(normB) {
489 for (int i = 0; i < 3; i++)
490 bisector[i] = bisector[i] / normB;
491 }
492 // raise the length of bisector if needed .........
493 DI_Point pt(midEN.x() + bisector[0], midEN.y() + bisector[1], midEN.z() + bisector[2]);
494 pt.addLs(e);
495 if (pt.ls() * midEN.ls() > 0.) {
496 pt.move(midEN.x() - bisector[0], midEN.y() - bisector[1], midEN.z() - bisector[2]);
497 pt.changeLs(e->evalLs(pt.x(), pt.y(), pt.z()));
498 }
499 DI_Point* qmn = Newton(&midEN, &pt, e, RPNi);
500 return qmn;
501 }
502
503 // ------------------------------------------------------------------------------------------------
504 // ------------------------------------------------------------------------------------------------
505
506 // DI_Point methods -------------------------------------------------------------------------------
DI_Point(double x,double y,double z,gLevelset * ls)507 DI_Point::DI_Point (double x, double y, double z, gLevelset *ls) : x_(x), y_(y), z_(z) {
508 addLs((*ls)(x, y, z));
509 }
operator =(const DI_Point & rhs)510 DI_Point & DI_Point::operator= (const DI_Point &rhs) {
511 if(this != &rhs){
512 x_ = rhs.x(); y_ = rhs.y(); z_ = rhs.z(); Ls = rhs.Ls;
513 }
514 return *this;
515 }
addLs(const double ls)516 void DI_Point::addLs (const double ls) {
517 Ls.push_back(adjustLs(ls));
518 }
addLs(const DI_Element * e)519 void DI_Point::addLs (const DI_Element *e) {
520 addLs(e->evalLs(x_, y_, z_));
521 }
chooseLs(const gLevelset * Lsi)522 void DI_Point::chooseLs (const gLevelset *Lsi) {
523 if(Ls.size() < 2) return;
524 double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
525 double ls = Lsi->choose(ls1, ls2);
526 Ls.pop_back(); Ls.pop_back();
527 addLs(ls);
528 }
computeLs(const DI_Element * e,const std::vector<gLevelset * > & RPNi)529 void DI_Point::computeLs (const DI_Element *e, const std::vector<gLevelset*> &RPNi) {
530 int prim = 0; Ls.clear();
531 if(e->sizeLs() == 0) return;
532 for(int l = 0; l < (int)RPNi.size(); l++) {
533 if (RPNi[l]->isPrimitive())
534 addLs(e->evalLs(x_, y_, z_, prim++));
535 else
536 chooseLs(RPNi[l]);
537 }
538 }
computeLs(const gLevelset * ls)539 void DI_Point::computeLs (const gLevelset *ls) {
540 Ls.clear();
541 addLs((*ls)(x_, y_, z_));
542 }
equal(const DI_Point * p) const543 bool DI_Point::equal(const DI_Point *p) const {
544 return (fabs(x() - p->x()) < EQUALITY_TOL && fabs(y() - p->y()) < EQUALITY_TOL && fabs(z() - p->z()) < EQUALITY_TOL);
545 }
546
547 // DI_IntegrationPoint methods --------------------------------------------------------------------
computeLs(const DI_Element * e,const std::vector<gLevelset * > & RPN)548 void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<gLevelset*> &RPN) {
549 int prim = 0; std::vector<double> Ls;
550 for(int l = 0; l < (int)RPN.size(); l++) {
551 if(RPN[l]->isPrimitive())
552 Ls.push_back(adjustLs(e->evalLs(x_, y_, z_, prim++)));
553 else {
554 double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
555 double ls = RPN[l]->choose(ls1, ls2);
556 Ls.pop_back(); Ls.pop_back();
557 Ls.push_back(adjustLs(ls));
558 }
559 }
560 setLs(Ls.back());
561 }
562
563 // DI_CuttingPoint methods ------------------------------------------------------------------------
DI_CuttingPoint(const DI_Point * pt)564 DI_CuttingPoint::DI_CuttingPoint(const DI_Point *pt) : DI_Point(pt->x(), pt->y(), pt->z()),
565 xl_(pt->x()), yl_(pt->y()), zl_(pt->z()) {
566 for(int i = 0; i < pt->sizeLs(); i++)
567 addLs(pt->ls(i));
568 }
569
570 // DI_PointLessThan methods -----------------------------------------------------------------------
571 double DI_PointLessThan::tolerance = 1.e-6;
operator ()(const DI_Point * p1,const DI_Point * p2) const572 bool DI_PointLessThan::operator()(const DI_Point *p1, const DI_Point *p2) const
573 {
574 if(p1->x() - p2->x() > tolerance) return true;
575 if(p1->x() - p2->x() < -tolerance) return false;
576 if(p1->y() - p2->y() > tolerance) return true;
577 if(p1->y() - p2->y() < -tolerance) return false;
578 if(p1->z() - p2->z() > tolerance) return true;
579 return false;
580 }
581
582 // DI_Element methods -----------------------------------------------------------------------------
DI_Element(const DI_Element & cp)583 DI_Element::DI_Element(const DI_Element &cp) : lsTag_(cp.lsTag()), polOrder_(cp.getPolynomialOrder()),
584 integral_(cp.integral()) {
585 pts_ = new DI_Point [cp.nbVert()];
586 for(int i = 0; i < cp.nbVert(); i++)
587 pts_[i] = DI_Point(*cp.pt(i));
588 if(cp.nbMid()) {
589 mid_ = new DI_Point [cp.nbMid()];
590 for(int i = 0; i < cp.nbMid(); i++)
591 mid_[i] = DI_Point(*cp.mid(i));
592 }
593 else
594 mid_ = NULL;
595 }
operator =(const DI_Element & rhs)596 DI_Element & DI_Element::operator= (const DI_Element &rhs){
597 if(type() != rhs.type()) {
598 printf("Error : try to assign element of different type!\n");
599 return *this;
600 }
601 if(this != &rhs) {
602 delete [] pts_;
603 pts_ = new DI_Point[rhs.nbVert()];
604 for(int i = 0; i < nbVert(); i++) {
605 pts_[i] = DI_Point(*rhs.pt(i));
606 }
607 if(rhs.nbMid()) {
608 if(mid_) delete [] mid_;
609 mid_ = new DI_Point[rhs.nbMid()];
610 for(int i = 0; i < rhs.nbMid(); i++)
611 mid_[i] = DI_Point(*(rhs.mid(i)));
612 }
613 else mid_ = NULL;
614 polOrder_ = rhs.getPolynomialOrder();
615 integral_ = rhs.integral();
616 lsTag_ = rhs.lsTag();
617 }
618 return *this;
619 }
getShapeFunctions(double u,double v,double w,double s[],int o) const620 void DI_Element::getShapeFunctions(double u, double v, double w, double s[], int o) const
621 {
622 //printf("type elem =%d order =%d\n", type(), o);
623 const nodalBasis* fs = getFunctionSpace(o);
624 if(fs) fs->f(u, v, w, s);
625 else Msg::Error("Function space not implemented for this type of element");
626 }
627
getGradShapeFunctions(double u,double v,double w,double s[][3],int o) const628 void DI_Element::getGradShapeFunctions(double u, double v, double w, double s[][3],
629 int o) const
630 {
631 const nodalBasis* fs = getFunctionSpace(o);
632 if(fs) fs->df(u, v, w, s);
633 else Msg::Error("Function space not implemented for this type of element");
634 }
setPolynomialOrder(int o)635 void DI_Element::setPolynomialOrder (int o) {
636 if(polOrder_ == o) return;
637 if (mid_){
638 delete [] mid_;
639 mid_ = NULL;
640 }
641 polOrder_ = o;
642
643 if(o == 1) return;
644
645 const nodalBasis *fs = getFunctionSpace(o);
646 if(!fs) Msg::Error("Function space not implemented for this type of element");
647
648 mid_ = new DI_Point[nbMid()];
649 int j = nbVert();
650 int dim = getDim();
651 double u, v, w;
652 double xyz[3];
653 for(int i = 0; i < nbMid(); i++, j++) {
654 u = fs->points(j,0);
655 v = (dim > 1) ? fs->points(j,1) : 0.;
656 w = (dim > 2) ? fs->points(j,2) : 0.;
657 evalC(u, v, w, xyz, 1);
658 mid_[i] = DI_Point(xyz[0], xyz[1], xyz[2]);
659 }
660 }
setPolynomialOrder(int o,const DI_Element * e,const std::vector<gLevelset * > & RPNi)661 void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
662 if(polOrder_ == o) return;
663 if (mid_){
664 delete [] mid_;
665 mid_ = NULL;
666 }
667 polOrder_ = o;
668
669 if(o == 1) return;
670
671 const nodalBasis *fs = getFunctionSpace(o);
672 if(!fs) Msg::Error("Function space not implemented for this type of element");
673
674 mid_ = new DI_Point[nbMid()];
675 int j = nbVert();
676 int dim = getDim();
677 double u, v, w;
678 double xyz[3];
679 for(int i = 0; i < nbMid(); i++, j++) {
680 u = fs->points(j,0);
681 v = (dim > 1) ? fs->points(j,1) : 0.;
682 w = (dim > 2) ? fs->points(j,2) : 0.;
683 evalC(u, v, w, xyz, 1);
684 mid_[i] = DI_Point(xyz[0], xyz[1], xyz[2], e, RPNi);
685 }
686 }
addLs(const double * ls)687 void DI_Element::addLs (const double *ls) {
688 for(int i = 0; i < nbVert() + nbMid(); i++)
689 pt(i)->addLs(ls[i]);
690 }
addLs(const DI_Element * e)691 void DI_Element::addLs (const DI_Element *e) {
692 if(e->sizeLs() < 1) return;
693 for(int i = 0; i < nbVert() + nbMid(); i++)
694 pt(i)->addLs(e);
695 }
addLs(const DI_Element * e,const gLevelset * Ls)696 void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls) {
697 if(type() != e->type()) {
698 printf("Error : addLs with element of different type\n");
699 }
700 for(int i = 0; i < nbVert() + nbMid(); ++i) {
701 double ls = (*Ls)(e->x(i), e->y(i), e->z(i));
702 pt(i)->addLs(ls);
703 }
704 }
addLs(const DI_Element * e,const gLevelset * Ls,int iLs,double ** nodeLs)705 void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls, int iLs, double **nodeLs) {
706 if(!nodeLs || !nodeLs[0][iLs]) {
707 addLs(e, Ls);
708 return;
709 }
710 for(int i = 0; i < nbVert() + nbMid(); i++)
711 pt(i)->addLs(nodeLs[i][iLs]);
712 }
chooseLs(const gLevelset * Lsi)713 void DI_Element::chooseLs (const gLevelset *Lsi) {
714 if(sizeLs() < 2)
715 printf("chooseLs with element ls size < 2 : typeEl=%d\n", type());
716 for(int i = 0; i < nbVert() + nbMid(); i++)
717 pt(i)->chooseLs(Lsi);
718 }
clearLs()719 void DI_Element::clearLs() {
720 for(int i = 0; i < nbVert() + nbMid(); i++)
721 pt(i)->clearLs();
722 }
addQuadEdge(int edge,DI_Point * xm,const DI_Element * e,const std::vector<gLevelset * > & RPNi)723 bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
724 const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
725 if(edge >= nbEdg()) {printf("wrong number (%d) for quadratic edge for a ", edge); print(); return false;}
726 int s1, s2; vert(edge, s1, s2);
727 int order0 = getPolynomialOrder();
728 if(order0 == 1) setPolynomialOrder(2, e, RPNi);
729 double dist = distance(mid(edge), xm);
730 double sideLength = distance(pt(s1), pt(s2));
731 if (dist / sideLength < 1.e-5) {
732 if(order0 == 1) setPolynomialOrder(1);
733 printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quadratic edge\n", dist, sideLength, dist/sideLength);
734 return true; // do not add the quadratic edge if xm is very close from the middle of the edge
735 }
736 DI_Point *p0 = &mid_[edge];
737 mid_[edge].move(xm->x(), xm->y(), xm->z());//mid_[edge] = xm;
738 // check if the quadratic edge will produce a twist in the element (detJ<0)
739 if(!testDetJ()){
740 // reinitialize quad edges if the element was not quad at the begining
741 if(order0 == 1) setPolynomialOrder(1);
742 else mid_[edge].move(p0->x(), p0->y(), p0->z());//mid_[edge] = p0;
743 printf("detJ<0 when trying to add a quadratic edge in "); print();
744 return false;
745 }
746 printf("in add quad edge \n");
747 computeIntegral();
748 return true;
749 }
contain(const DI_Point * p) const750 bool DI_Element::contain (const DI_Point *p) const {
751 for(int i = 0; i < nbVert(); i++){
752 if(p->equal(pt(i)))
753 return true;
754 }
755 switch(getDim()){
756 case 1 :
757 if(!isLinear(pt(0), pt(1), p)) return false;
758 if(distance(p, pt(0)) > integral() || distance(p, pt(1)) > integral()) return false;
759 return true;
760 case 2 :
761 if(!isPlanar(pt(0), pt(1), pt(2), p)) return false;
762 for(int i = 0; i < nbVert(); i++){
763 double v1[3]; vec(pt(i), pt((i + 1) % nbVert()), v1);
764 double v2[3]; vec(pt(i), pt((i + 2) % nbVert()), v2);
765 double vp[3]; vec(pt(i), p, vp);
766 double c1[3]; cross(v1, v2, c1);
767 double c2[3]; cross(v1, vp, c2);
768 if(dot(c1, c2) < 0.) return false; // assert pt is on the same side of (pt(i),pt(i+1)) than pt(i+2)
769 }
770 return true;
771 case 3 :
772 for(int i = 0; i < nbVert(); i++){
773 double v1[3]; vec(pt(i), pt((i + 1) % nbVert()), v1);
774 double v2[3]; vec(pt(i), pt((i + 2) % nbVert()), v2);
775 double v3[3]; vec(pt(i), pt((i + 3) % nbVert()), v3);
776 double vp[3]; vec(pt(i), p, vp);
777 double v1v2[3]; cross(v1, v2, v1v2);
778 double c1 = dot(v1v2, v3);
779 double c2 = dot(v1v2, vp);
780 if(c1 * c2 < 0.) return false; // assert pt is on the same side of (pt(i),pt(i+1),pt(i+2)) than pt(i+3)
781 }
782 return true;
783 default :
784 return false;
785 }
786 }
contain(const DI_Element * e) const787 bool DI_Element::contain (const DI_Element *e) const {
788 for(int p = 0; p < e->nbVert(); p++)
789 if(!contain(e->pt(p))) return false;
790 return true;
791 }
mappingP(DI_Point * pt) const792 void DI_Element::mappingP (DI_Point *pt) const {
793 double e[3]; evalC(pt->x(), pt->y(), pt->z(), e);
794 pt->move(e[0], e[1], e[2]);
795 }
mappingIP(DI_IntegrationPoint * in) const796 void DI_Element::mappingIP (DI_IntegrationPoint *in) const {
797 double e[3]; evalC(in->x(), in->y(), in->z(), e);
798 double w = in->weight() * integral() / refIntegral();
799 in->move(e[0], e[1], e[2]);
800 in->setWeight(w);
801 }
mappingCP(DI_CuttingPoint * cp) const802 void DI_Element::mappingCP (DI_CuttingPoint *cp) const {
803 double e[3]; evalC(cp->x(), cp->y(), cp->z(), e);
804 cp->move(e[0], e[1], e[2]);
805 }
mappingEl(DI_Element * el) const806 void DI_Element::mappingEl (DI_Element *el) const {
807 double s[3];
808 for (int i = 0; i < el->nbVert() + el->nbMid(); i++) {
809 evalC(el->x(i), el->y(i), el->z(i), s);
810 el->pt(i)->move(s[0], s[1], s[2]);
811 }
812 el->computeIntegral();
813 }
integrationPoints(const int polynomialOrder,const DI_Element * loc,const DI_Element * e,const std::vector<gLevelset * > & RPN,std::vector<DI_IntegrationPoint * > & ip) const814 void DI_Element::integrationPoints (const int polynomialOrder, const DI_Element *loc,
815 const DI_Element *e, const std::vector<gLevelset *> &RPN,
816 std::vector<DI_IntegrationPoint *> &ip) const {
817 std::vector<DI_IntegrationPoint *> ip_ref;
818 getRefIntegrationPoints(polynomialOrder, ip_ref);
819 for (int j = 0; j < (int)ip_ref.size(); j++) {
820 DI_IntegrationPoint IPloc (*ip_ref[j]);
821 loc->mappingIP(&IPloc);
822 mappingIP(ip_ref[j]);
823 ip_ref[j]->addLocC(IPloc.x(), IPloc.y(), IPloc.z());
824 DI_IntegrationPoint pp(IPloc); pp.computeLs(e, RPN);
825 //ip_ref[j]->setLs(loc->evalLs(IPloc.x(), IPloc.y(), IPloc.z()));//levelset computed in the subelement FIXME check sign!
826 //ip_ref[j]->setLs(evalLs(ip_ref[j]->x(), ip_ref[j]->y(), ip_ref[j]->z()));
827 ip_ref[j]->setLs(pp.ls());
828 /*printf("pt loc : "); IPloc.print(); printf("pt reel : "); ip_ref[j]->print();
829 printf("el loc : "); loc->printls(); printf("el reel : "); printls(); printf("el e : "); e->printls();
830 printf("ls_loc=%g ls_reel=%g ls_e=%g\n\n", loc->evalLs(IPloc.x(), IPloc.y(), IPloc.z()),
831 evalLs(ip_ref[j]->x(), ip_ref[j]->y(), ip_ref[j]->z()), pp.ls());*/
832 ip.push_back(ip_ref[j]);
833 }
834 }
getCuttingPoints(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_CuttingPoint * > & cp) const835 void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
836 std::vector<DI_CuttingPoint *> &cp) const {
837 int s1, s2;
838 for(int i = 0; i < nbEdg(); i++){
839 vert(i, s1, s2);
840 if (isCrossed(pt(s1), pt(s2))) {
841 DI_Point* p = Newton(pt(s1), pt(s2), e, RPNi);
842 cp.push_back(new DI_CuttingPoint(p));
843 delete p;
844 }
845 }
846 for(int i = 0; i < nbVert(); i++){
847 if(ls(i) == 0.)
848 cp.push_back(new DI_CuttingPoint(pt(i)));
849 }
850 }
evalC(const double u,const double v,const double w,double * ev,int order) const851 void DI_Element::evalC (const double u, const double v, const double w, double *ev, int order) const {
852 int nbV = nbVert() + nbMid();
853 std::vector<double> s(nbV);
854 ev[0] = 0; ev[1] = 0; ev[2] = 0;
855 getShapeFunctions (u, v, w, &s[0], order);
856 for(int i = 0; i < nbV; i++){
857 ev[0] += x(i) * s[i];
858 ev[1] += y(i) * s[i];
859 ev[2] += z(i) * s[i];
860 }
861 }
evalLs(const double u,const double v,const double w,int iLs,int order) const862 double DI_Element::evalLs (const double u, const double v, const double w, int iLs, int order) const{
863 int nbV = nbVert() + nbMid();
864 if(iLs == -1) iLs = sizeLs() - 1; //last ls value
865 double vls = 0;
866 std::vector<double> s(nbV);
867 getShapeFunctions (u, v, w, &s[0], order);
868 for(int i = 0; i < nbV; i++)
869 vls += ls(i, iLs) * s[i];
870 return adjustLs(vls);
871 }
testDetJ() const872 bool DI_Element::testDetJ() const {
873 double detJ0 = detJ(x(0), y(0), z(0));
874 for(int i = 1; i < nbVert() + nbMid(); i++)
875 if(detJ0 * detJ(x(i), y(i), z(i)) < 0.) return false;
876 return true;
877 }
detJ(const double u,const double v,const double w) const878 double DI_Element::detJ (const double u, const double v, const double w) const {
879 double J[3][3];
880 J[0][0] = J[0][1] = J[0][2] = 0.;
881 J[1][0] = J[1][1] = J[1][2] = 0.;
882 J[2][0] = J[2][1] = J[2][2] = 0.;
883 /*double **s=new double*[nbVert()+ nbMid()];
884 for(int i = 0; i < nbVert()+ nbMid(); ++i){
885 s[i] = new double[3];//s[i] = new double[getDim()];
886 }*/
887 typedef double d3[3];
888 int nbV = nbVert() + nbMid();
889 d3 *s = new d3[nbV];
890 getGradShapeFunctions(u, v, w, s);
891 switch(getDim()){
892 case 3 :
893 for(int i = 0; i < nbV; i++) {
894 J[0][0] += x(i) * s[i][0]; J[0][1] += y(i) * s[i][0]; J[0][2] += z(i) * s[i][0];
895 J[1][0] += x(i) * s[i][1]; J[1][1] += y(i) * s[i][1]; J[1][2] += z(i) * s[i][1];
896 J[2][0] += x(i) * s[i][2]; J[2][1] += y(i) * s[i][2]; J[2][2] += z(i) * s[i][2];
897 // delete [] s[i];
898 }
899 delete [] s;
900 return det3(J[0][0], J[0][1], J[0][2], J[1][0], J[1][1], J[1][2], J[2][0], J[2][1], J[2][2]);
901 case 2 :
902 for(int i = 0; i < nbV; i++) {
903 J[0][0] += x(i) * s[i][0]; J[0][1] += y(i) * s[i][0]; J[0][2] += z(i) * s[i][0];
904 J[1][0] += x(i) * s[i][1]; J[1][1] += y(i) * s[i][1]; J[1][2] += z(i) * s[i][1];
905 // delete [] s[i];
906 }
907 delete [] s;
908 return sqrt(sq2(J[0][0] * J[1][1] - J[0][1] * J[1][0]) +
909 sq2(J[0][2] * J[1][0] - J[0][0] * J[1][2]) +
910 sq2(J[0][1] * J[1][2] - J[0][2] * J[1][1])); // allways > 0 => does'nt allow to check if twist !!!!
911 case 1 :
912 for(int i = 0; i < nbV; i++) {
913 J[0][0] += x(i) * s[i][0]; J[0][1] += y(i) * s[i][0]; J[0][2] += z(i) * s[i][0];
914 // delete [] s[i];
915 }
916 delete [] s;
917 return sqrt(sq2(J[0][0]) + sq2(J[0][1]) + sq2(J[0][2])); // allways > 0 !!!!
918 default :
919 delete [] s;
920 return 1.;
921 }
922 }
923 // set the lsTag to +1 if the element is inside the domain (lsTag is -1 by default)
computeLsTagDom(const DI_Element * e,const std::vector<gLevelset * > & RPN)924 void DI_Element::computeLsTagDom(const DI_Element *e, const std::vector<gLevelset *> &RPN) {
925 for(int i = 0; i < nbVert(); i++) {
926 if(pt(i)->isOutsideDomain())
927 return;
928 if(pt(i)->isInsideDomain())
929 {setLsTag(1); return;}
930 }
931 DI_Point* bar = barycenter(this, e, RPN);
932 if(bar->isOutsideDomain())
933 {delete bar; return;}
934 if(bar->isInsideDomain())
935 {setLsTag(1); delete bar; return;}
936 for(int i = 0; i < nbVert(); i++) {
937 DI_Point* mid = middle (pt(i), bar, e, RPN);
938 if(mid->isOutsideDomain())
939 {delete mid; delete bar; return;}
940 if(mid->isInsideDomain())
941 {setLsTag(1); delete mid; delete bar; return;}
942 delete mid;
943 }
944 delete bar;
945 printf("Error : Unable to determine the sign of the element : \n");
946 printf(" - Parent element : "); e->printls();
947 printf(" - Element : "); printls();
948 }
949 // set the lsTag to -1 if the element is not on the border of the domain
computeLsTagBound(std::vector<DI_Hexa * > & hexas,std::vector<DI_Tetra * > & tetras)950 void DI_Element::computeLsTagBound(std::vector<DI_Hexa *> &hexas, std::vector<DI_Tetra *> &tetras){
951 for(int i = 0; i < nbVert(); i++) {
952 if(!pt(i)->isOnBorder()) {
953 setLsTag(-1);
954 return;
955 }
956 }
957
958 /*DI_Element *e1 = NULL, *e2 = NULL;
959 for(int i = 0; i < (int)tetras.size(); i++){
960 if(belongsTo(this, tetras[i])) {
961 if(!e1) e1 = tetras[i];
962 else {e2 = tetras[i]; break;}
963 }
964 }
965 if(e1 && e2) {
966 if(e1->lsTag() * e2->lsTag() > 0.) setLsTag(-1);
967 return;
968 }
969 for(int i = 0; i < (int)hexas.size(); i++){
970 if(belongsTo(this, hexas[i])) {
971 if(!e1) e1 = hexas[i];
972 else {e2 = hexas[i]; break;}
973 }
974 }
975 if(e1 && e2 && e1->lsTag() * e2->lsTag() > 0.) setLsTag(-1);*/
976 }
computeLsTagBound(std::vector<DI_Quad * > & quads,std::vector<DI_Triangle * > & triangles)977 void DI_Element::computeLsTagBound(std::vector<DI_Quad *> &quads, std::vector<DI_Triangle *> &triangles){
978 for(int i = 0; i < nbVert(); i++) {
979 if(ls(i) != 0.) {
980 setLsTag(-1);
981 return;
982 }
983 }
984
985 /*DI_Element *e1 = NULL, *e2 = NULL, *et = NULL;
986 int NT = triangles.size(), NQ = quads.size();
987 for(int i = 0; i < NT + NQ; i++){
988 if((i < NT && belongsTo(this, triangles[i])) ||
989 (i >= NT && belongsTo(this, quads[i-NT]))) {
990 if(i < NT) et = triangles[i];
991 else et = quads[i-NT];
992 if(!e1) e1 = et;
993 else {e2 = et; break;}
994 }
995 }
996 if(e1 && e2 && e1->lsTag() * e2->lsTag() > 0.) setLsTag(-1);*/
997 }
998
print() const999 void DI_Element::print() const {
1000 switch(type()) {
1001 case DI_LIN : printf("Line"); break;
1002 case DI_TRI : printf("Triangle"); break;
1003 case DI_QUA : printf("Quad"); break;
1004 case DI_TET : printf("Tetra"); break;
1005 case DI_HEX : printf("Hexa"); break;
1006 default : printf("Element");
1007 }
1008 printf("%d ", getPolynomialOrder());
1009 for(int i = 0; i < nbVert() + nbMid(); i++)
1010 printf("(%g,%g,%g) ", x(i), y(i), z(i));
1011 printf("tag=%d\n", lsTag());
1012 }
printls() const1013 void DI_Element::printls() const {
1014 switch(type()) {
1015 case DI_LIN : printf("Line"); break;
1016 case DI_TRI : printf("Triangle"); break;
1017 case DI_QUA : printf("Quad"); break;
1018 case DI_TET : printf("Tetra"); break;
1019 case DI_HEX : printf("Hexa"); break;
1020 default : printf("Element");
1021 }
1022 printf("%d ", getPolynomialOrder());
1023 for(int i = 0; i < nbVert() + nbMid(); i++) {
1024 printf("(%g,%g,%g) ls=(", x(i), y(i), z(i));
1025 for(int j = 0; j < sizeLs(); j++) printf("%g,", ls(i,j));
1026 printf("); ");
1027 }
1028 printf("tag=%d\n", lsTag());
1029 }
1030
1031 double DI_ElementLessThan::tolerance = 1.e-6;
operator ()(const DI_Element * e1,const DI_Element * e2) const1032 bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2) const
1033 {
1034 for(int i = 0; i < e1->nbVert(); i++) {
1035 if(e1->x(i) - e2->x(i) > tolerance) return true;
1036 if(e1->x(i) - e2->x(i) < -tolerance) return false;
1037 if(e1->y(i) - e2->y(i) > tolerance) return true;
1038 if(e1->y(i) - e2->y(i) < -tolerance) return false;
1039 if(e1->z(i) - e2->z(i) > tolerance) return true;
1040 }
1041 return false;
1042 }
1043
1044 // DI_Line methods --------------------------------------------------------------------------------
getFunctionSpace(int o) const1045 const nodalBasis* DI_Line::getFunctionSpace(int o) const{
1046 int order = (o == -1) ? getPolynomialOrder() : o;
1047 int tag = ElementType::getType(TYPE_LIN, order);
1048 return BasisFactory::getNodalBasis(tag);
1049 }
1050
computeIntegral()1051 void DI_Line::computeIntegral() {
1052 integral_ = LineLength(pt(0), pt(1));
1053 // switch (getPolynomialOrder()) {
1054 // case 1 :
1055 // integral_ = LineLength(pt(0), pt(1));
1056 // break;
1057 // case 2 :
1058 // integral_ = LineLength(pt(0), mid(0)) + LineLength(mid(0), pt(1));
1059 // break;
1060 // default :
1061 // printf("Order %d line function space not implemented ", getPolynomialOrder()); print();
1062 // }
1063 }
1064 // DI_Triangle methods ----------------------------------------------------------------------------
getFunctionSpace(int o) const1065 const nodalBasis* DI_Triangle::getFunctionSpace(int o) const
1066 {
1067 int order = (o == -1) ? getPolynomialOrder() : o;
1068 int tag = ElementType::getType(TYPE_TRI, order);
1069 return BasisFactory::getNodalBasis(tag);
1070 }
computeIntegral()1071 void DI_Triangle::computeIntegral() {
1072 integral_ = TriSurf(pt(0), pt(1), pt(2));
1073 // switch (getPolynomialOrder()) {
1074 // case 1 :
1075 // integral_ = TriSurf(pt(0), pt(1), pt(2));
1076 // break;
1077 // case 2 :
1078 // integral_ = TriSurf(pt(0), mid(0), mid(2)) + TriSurf(pt(1), mid(0), mid(1))
1079 // + TriSurf(pt(2), mid(1), mid(2)) + TriSurf(mid(0), mid(1), mid(2));
1080 // break;
1081 // default :
1082 // printf("Order %d triangle function space not implemented ", getPolynomialOrder()); print();
1083 // }
1084 }
quality() const1085 double DI_Triangle::quality() const {
1086 return qualityTri(pt(0), pt(1), pt(2));
1087 }
1088
1089 // DI_Quad methods --------------------------------------------------------------------------------
getFunctionSpace(int o) const1090 const nodalBasis* DI_Quad::getFunctionSpace(int o) const{
1091 int order = (o == -1) ? getPolynomialOrder() : o;
1092 int tag = ElementType::getType(TYPE_QUA, order);
1093 return BasisFactory::getNodalBasis(tag);
1094 }
1095
computeIntegral()1096 void DI_Quad::computeIntegral() {
1097 integral_ = TriSurf(pt(0), pt(1), pt(2)) + TriSurf(pt(0), pt(2), pt(3));
1098 // switch (getPolynomialOrder()) {
1099 // case 1 :
1100 // integral_ = TriSurf(pt(0), pt(1), pt(2)) + TriSurf(pt(0), pt(2), pt(3));
1101 // break;
1102 // case 2 :
1103 // integral_ = TriSurf(pt(0), mid(0), mid(4)) + TriSurf(pt(0), mid(4), mid(3))
1104 // + TriSurf(pt(1), mid(1), mid(4)) + TriSurf(pt(1), mid(4), mid(0))
1105 // + TriSurf(pt(2), mid(2), mid(4)) + TriSurf(pt(2), mid(4), mid(1))
1106 // + TriSurf(pt(3), mid(3), mid(4)) + TriSurf(pt(3), mid(4), mid(2));
1107 // break;
1108 // default :
1109 // printf("Order %d quadrangle function space not implemented ", getPolynomialOrder()); print();
1110 // }
1111 }
1112
1113 // DI_Tetra methods -------------------------------------------------------------------------------
getFunctionSpace(int o) const1114 const nodalBasis* DI_Tetra::getFunctionSpace(int o) const{
1115 int order = (o == -1) ? getPolynomialOrder() : o;
1116 int tag = ElementType::getType(TYPE_TET, order);
1117 return BasisFactory::getNodalBasis(tag);
1118 }
1119
computeIntegral()1120 void DI_Tetra::computeIntegral() {
1121 integral_ = TetraVol(pt(0), pt(1), pt(2), pt(3));
1122 }
quality() const1123 double DI_Tetra::quality() const {
1124 return qualityTet(x(0), y(0), z(0), x(1), y(1), z(1), x(2), y(2), z(2), x(3), y(3), z(3));
1125 }
1126
1127
1128 // Hexahedron methods -----------------------------------------------------------------------------
getFunctionSpace(int o) const1129 const nodalBasis* DI_Hexa::getFunctionSpace(int o) const{
1130 int order = (o == -1) ? getPolynomialOrder() : o;
1131 int tag = ElementType::getType(TYPE_HEX, order);
1132 return BasisFactory::getNodalBasis(tag);
1133 }
computeIntegral()1134 void DI_Hexa::computeIntegral() {
1135 integral_ = TetraVol(pt(0), pt(1), pt(3), pt(4)) + TetraVol(pt(1), pt(4), pt(5), pt(7))
1136 + TetraVol(pt(1), pt(3), pt(4), pt(7)) + TetraVol(pt(2), pt(5), pt(6), pt(7))
1137 + TetraVol(pt(1), pt(2), pt(3), pt(7)) + TetraVol(pt(1), pt(5), pt(2), pt(7));
1138 }
1139
1140 // ----------------------------------------------------------------------------
1141 // -------------------------- SPLITTING ---------------------------------------
1142 // ----------------------------------------------------------------------------
1143
1144 // Split a reference line cut by a level set into sublines
selfSplit(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Line * > & subLines,std::vector<DI_CuttingPoint * > & cuttingPoints) const1145 void DI_Line::selfSplit (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
1146 std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
1147 {
1148 if (!isCrossed(pt(0), pt(1))) {
1149 subLines.push_back(new DI_Line(*this));
1150 return;
1151 }
1152
1153 // compute the intersection between the line and the level set
1154 DI_Point *U = Newton(pt(0), pt(1), e, RPNi);
1155 cuttingPoints.push_back(new DI_CuttingPoint(U));
1156
1157 // line cut => split into 2 sublines
1158 subLines.push_back(new DI_Line(pt(0), U, lsTag()));
1159 subLines.push_back(new DI_Line(U, pt(1), lsTag()));
1160 delete U;
1161 }
1162
1163 // Split a triangle into 4 sub-triangles
splitIntoSubTriangles(std::vector<DI_Triangle * > & subTriangles) const1164 void DI_Triangle::splitIntoSubTriangles (std::vector<DI_Triangle *> &subTriangles) const {
1165 DI_Point *p01 = middle(pt(0), pt(1));
1166 DI_Point *p02 = middle(pt(0), pt(2));
1167 DI_Point *p12 = middle(pt(1), pt(2));
1168
1169 subTriangles.push_back(new DI_Triangle(p01, p02, p12)); // 01-02-12
1170 subTriangles.push_back(new DI_Triangle(pt(0), p01, p02)); // 0-01-02
1171 subTriangles.push_back(new DI_Triangle(pt(1), p01, p12)); // 1-01-12
1172 subTriangles.push_back(new DI_Triangle(pt(2), p02, p12)); // 2-02-12
1173 delete p01; delete p02; delete p12;
1174 }
1175 // Split a reference triangle cut by a level set into subtriangles
selfSplit(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Quad * > & subQuads,std::vector<DI_Triangle * > & subTriangles,std::vector<DI_Line * > & surfLines,std::vector<DI_CuttingPoint * > & cuttingPoints) const1176 void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
1177 std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
1178 std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
1179 {
1180 //bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // if true, use quadratic sub triangles
1181 bool quadT = false;
1182 int LStag = RPNi.back()->getTag();
1183
1184 int nbZe = 0, ze[3];
1185 for(int i = 0; i < 3; i++)
1186 if(pt(i)->isOnBorder())
1187 ze[nbZe++] = i;
1188 for(int i = 0; i < nbZe; i++)
1189 cuttingPoints.push_back(new DI_CuttingPoint(pt(ze[i])));
1190 // !! we add these points several times => remove later
1191
1192 if (!isCrossed(pt(0), pt(1)) && !isCrossed(pt(0), pt(2)) && !isCrossed(pt(1), pt(2))) {
1193 subTriangles.push_back(new DI_Triangle(*this));
1194 if(nbZe == 2)
1195 surfLines.push_back(new DI_Line(pt(ze[0]), pt(ze[1]), LStag));
1196 // !! we add these lines twice => remove the second later
1197 return;
1198 }
1199
1200 DI_Point *U[2]; // cutting points (max2)
1201 int IND[2]; // ids of edges cut
1202 int COUNT = 0; // number of edges cut
1203
1204 // compute the intersections between the edges of the triangle and the level set
1205 for(int i = 0; i < 3; i++){
1206 int n1 = i, n2 = (i + 1) % 3;
1207 if (isCrossed(pt(n1), pt(n2))) {
1208 U[COUNT] = Newton(pt(n1), pt(n2), e, RPNi);
1209 IND[COUNT++] = i;
1210 }
1211 }
1212
1213 for(int i = 0; i < COUNT; i++)
1214 cuttingPoints.push_back(new DI_CuttingPoint(U[i]));
1215
1216 // Adjust the indices of the nodes in order to have the same pattern for each case with the same number of edges cut;
1217 // compute the position of the middle nodes on the quadratic edges of the sub triangles;
1218 // and create the sub triangles
1219 switch (COUNT) {
1220 case 1 : {// 1 edge cut => split into 2 triangles
1221 int s1 = IND[0];
1222 int s2 = (s1 + 1) % 3;
1223 int s3 = (s2 + 1) % 3;
1224
1225 DI_Triangle *t1 = new DI_Triangle(pt(s3), pt(s1), U[0], lsTag());
1226 DI_Triangle *t2 = new DI_Triangle(pt(s2), pt(s3), U[0], lsTag());
1227 DI_Line *ln = new DI_Line(U[0], pt(s3), LStag);
1228
1229 if(quadT){
1230 DI_Point *midEN2 = quadMidNode(U[0], pt(s3), pt(s1), e, RPNi); // intersection between ls and the bissector between the cutting points
1231 cuttingPoints.push_back(new DI_CuttingPoint(midEN2));
1232 t1->addQuadEdge (2, midEN2, e, RPNi);
1233 t2->addQuadEdge (1, midEN2, e, RPNi);
1234 ln->addQuadEdge (0, midEN2, e, RPNi);
1235 delete midEN2;
1236 }
1237 subTriangles.push_back(t1);
1238 subTriangles.push_back(t2);
1239 surfLines.push_back(ln);
1240 delete U[0];
1241 break;
1242 }
1243 case 2: {// 2 edges cut => 1 triangle + 1 quad => split into 3 triangles
1244 int s1 = 2 * IND[0] - IND[1] + 2; // arbitrary formulation to find the good vertex !!!
1245 int s2 = (s1 + 1) % 3;
1246 int s3 = (s1 + 2) % 3;
1247 if(s1 == 0) { // (0,2) => swap U[0] and U[1]
1248 swap(&U[0], &U[1]);
1249 }
1250
1251 bool useQuads = false; // if true, split the triangle into 1 triangle and 1 quad
1252
1253 DI_Triangle *t1 = new DI_Triangle(U[0], pt(s1), U[1], lsTag());
1254 DI_Line *ln = new DI_Line(U[0], U[1], LStag);
1255 DI_Quad *q1; DI_Triangle *t2, *t3;
1256 if(useQuads) q1 = new DI_Quad(U[0], U[1], pt(s2), pt(s3), lsTag());
1257 else bestQuality(U[0], U[1], pt(s2), pt(s3), &t2, &t3);
1258 t2->setLsTag(lsTag()); t3->setLsTag(lsTag());
1259
1260 if(quadT){
1261 DI_Point *midEN2 = quadMidNode(U[0], U[1], pt(s2), e, RPNi); // intersection between ls and the bissector between the cutting points
1262 cuttingPoints.push_back(new DI_CuttingPoint(midEN2));
1263 //midEN2.printls();
1264 t1->addQuadEdge(2, midEN2, e, RPNi); //t1.printls();
1265 ln->addQuadEdge(0, midEN2, e, RPNi); //ln.printls();
1266 if(useQuads) q1->addQuadEdge(0, midEN2, e, RPNi);
1267 else t2->addQuadEdge(0, midEN2, e, RPNi); //t2.printls();
1268 delete midEN2;
1269 }
1270 subTriangles.push_back(t1);
1271 surfLines.push_back(ln);
1272 if(useQuads) subQuads.push_back(q1);
1273 else {
1274 subTriangles.push_back(t2);
1275 subTriangles.push_back(t3);
1276 }
1277 delete U[0]; delete U[1];
1278 break;
1279 }
1280 default:
1281 printf("Error : %d edge(s) cut in the triangle (ls : %g %g %g)\n",
1282 COUNT, ls(0), ls(1), ls(2));
1283 }
1284 }
1285
1286 // Split a reference Quadrangle into 2 triangles
splitIntoTriangles(std::vector<DI_Triangle * > & triangles) const1287 void DI_Quad::splitIntoTriangles(std::vector<DI_Triangle *> &triangles) const
1288 {
1289 triangles.push_back(new DI_Triangle(pt(0), pt(1), pt(3), lsTag())); // 013
1290 triangles.push_back(new DI_Triangle(pt(1), pt(2), pt(3), lsTag())); // 123
1291 }
1292
1293 // Split a reference tetrahedron cut by a level set (defined in a hex) into sub tetrahedra
1294 // and create triangles on the surface of the level set
selfSplit(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Tetra * > & subTetras,std::vector<DI_Triangle * > & surfTriangles,std::vector<DI_CuttingPoint * > & cuttingPoints,std::vector<DI_QualError * > & QError) const1295 void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
1296 std::vector<DI_Tetra *> &subTetras, std::vector<DI_Triangle *> &surfTriangles,
1297 std::vector<DI_CuttingPoint *> &cuttingPoints, std::vector<DI_QualError *> &QError) const
1298 {
1299 //bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // use of quadratic surf triangles and quadratic sub tetrahedra
1300 bool quadT = false;
1301 int tag = RPNi.back()->getTag();
1302
1303 int nbZe = 0, ze[4];
1304 for(int i = 0; i < 4; i++)
1305 if(pt(i)->isOnBorder())
1306 ze[nbZe++] = i;
1307 for(int i = 0; i < nbZe; i++)
1308 cuttingPoints.push_back(new DI_CuttingPoint(pt(ze[i]))); // !! we add these points several times => remove later
1309
1310 // case 0 : the tetrahedron is not cut by the level set
1311 if (!isCrossed(pt(0), pt(1)) && !isCrossed(pt(0), pt(2)) && !isCrossed(pt(1), pt(2)) &&
1312 !isCrossed(pt(0), pt(3)) && !isCrossed(pt(1), pt(3)) && !isCrossed(pt(2), pt(3))) {
1313 subTetras.push_back(new DI_Tetra(*this));
1314 if(nbZe == 3)
1315 surfTriangles.push_back(new DI_Triangle(pt(ze[0]), pt(ze[1]), pt(ze[2]), tag));
1316 // !! we add these triangles twice => remove the second later
1317 return;
1318 }
1319
1320 DI_Point *U[4]; // cutting points
1321 int COUNT = 0; // number of edges cut
1322 int IND[4]; // ids of edges cut
1323
1324 // edges nodes and opposite edges nodes
1325 int tab [6][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 2, 0, 3}, {2, 3, 0, 1}, {3, 1, 0, 2}};
1326
1327 // compute the intersections between the edges of the tetra and the level set.
1328 for(int i = 0; i < 6; i++){
1329 int n1 = tab[i][0], n2 = tab[i][1];
1330 if (isCrossed(pt(n1), pt(n2))) {
1331 U[COUNT] = Newton(pt(n1), pt(n2), e, RPNi);
1332 IND[COUNT++] = i;
1333 }
1334 }
1335 for(int i = 0; i < COUNT; i++)
1336 cuttingPoints.push_back(new DI_CuttingPoint(U[i]));
1337
1338 // Adjust the indices of the nodes in order to have the same pattern for each case with the same number of edges cut;
1339 // compute the position of the middle nodes on the quadratic edges of the sub tetras;
1340 // and create the sub tetras and sub triangles on the level set surface
1341 switch (COUNT) {
1342 case 1 : {// 1 edge cut => split into 2 tetras
1343 int i1 = IND[0]; //printf("case 1 : ind : %d\n",i1);
1344
1345 int s1 = tab[i1][0];
1346 int s2 = tab[i1][1];
1347 int s3 = tab[i1][2];
1348 int s4 = tab[i1][3];
1349
1350 DI_Tetra *t1 = new DI_Tetra(U[0], pt(s4), pt(s3), pt(s1));
1351 DI_Tetra *t2 = new DI_Tetra(pt(s4), U[0], pt(s3), pt(s2));
1352 DI_Triangle *tr = new DI_Triangle(pt(s4), U[0], pt(s3), tag);
1353
1354 if(quadT){
1355 DI_Point *midEN2[2]; // intersection between ls and the bissector between the cutting points
1356 midEN2[0] = quadMidNode(U[0], pt(s4), pt(s1), e, RPNi);
1357 midEN2[1] = quadMidNode(U[0], pt(s3), pt(s1), e, RPNi);
1358 tr->addQuadEdge (0, midEN2[0], e, RPNi);
1359 t1->addQuadEdge (0, midEN2[0], e, RPNi);
1360 t2->addQuadEdge (0, midEN2[0], e, RPNi);
1361 tr->addQuadEdge (1, midEN2[1], e, RPNi);
1362 t1->addQuadEdge (1, midEN2[1], e, RPNi);
1363 t2->addQuadEdge (3, midEN2[1], e, RPNi);
1364 delete midEN2[0]; delete midEN2[1];
1365 }
1366 subTetras.push_back(t1);
1367 subTetras.push_back(t2);
1368 surfTriangles.push_back(tr);
1369 delete U[0];
1370 break;
1371 }
1372 case 2 : {// 2 edges cut => 1 tetra + 1 pyramid => split into 3 tetras
1373 int i1 = IND[0], i2 = IND[1];
1374
1375 if((i1 == 0 && (i2 == 2 ||i2 == 3)) || (i1 == 1 && i2 == 4) || ((i1 == 2 || i1 == 3) && i2 == 5)) {
1376 swap(&U[0], &U[1]); swap(i1, i2);
1377 }
1378 //printf("case 2 : ind : %d,%d\n",i1,i2);
1379
1380 int s1 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]);
1381 int s2 = commonV(tab[i1][0], tab[i1][1], tab[i2][0], tab[i2][1]);
1382 int s3 = commonV(tab[i1][2], tab[i1][3], tab[i2][0], tab[i2][1]);
1383 int s4 = commonV(tab[i1][0], tab[i1][1], tab[i2][2], tab[i2][3]);
1384
1385 DI_Tetra *t1 = new DI_Tetra(U[0], U[1], pt(s2), pt(s1));
1386 DI_Triangle *tr = new DI_Triangle(U[0], U[1], pt(s1), tag);
1387 DI_Tetra *t2, *t3;
1388 int qual = bestQuality(U[1], U[0], pt(s4), pt(s3), pt(s1), &t2, &t3);
1389
1390 if(quadT){
1391 DI_Point *midEN2[3]; // intersection between ls and the bissector between the cutting points
1392 midEN2[0] = quadMidNode(U[1], pt(s1), pt(s2), e, RPNi);
1393 midEN2[1] = quadMidNode(U[0], pt(s1), pt(s2), e, RPNi);
1394 midEN2[2] = quadMidNode(U[1], U[0], pt(s4), e, RPNi);
1395 tr->addQuadEdge(1, midEN2[0], e, RPNi);
1396 t1->addQuadEdge(5, midEN2[0], e, RPNi);
1397 t2->addQuadEdge(2, midEN2[0], e, RPNi);
1398 t3->addQuadEdge(2, midEN2[0], e, RPNi);
1399
1400 tr->addQuadEdge(2, midEN2[1], e, RPNi);
1401 t1->addQuadEdge(2, midEN2[1], e, RPNi);
1402 t2->addQuadEdge(5, midEN2[1], e, RPNi);
1403
1404 tr->addQuadEdge(0, midEN2[2], e, RPNi);
1405 if(t2->addQuadEdge(0, midEN2[2], e, RPNi)) {
1406 t1->addQuadEdge(0, midEN2[2], e, RPNi);
1407 }
1408 else {
1409 DI_Point* mid = middle(pt(s3), pt(s4), e, RPNi);
1410 DI_Point* quad = middle(mid, midEN2[2], e, RPNi);
1411 if(qual == 1){
1412 t2->addQuadEdge(1, quad, e, RPNi);
1413 t3->addQuadEdge(0, quad, e, RPNi);
1414 }
1415 else {
1416 t2->addQuadEdge(3, quad, e, RPNi);
1417 t3->addQuadEdge(1, quad, e, RPNi);
1418 }
1419 if(t2->addQuadEdge(0, midEN2[2], e, RPNi)) {
1420 t1->addQuadEdge(0, midEN2[2], e, RPNi);
1421 }
1422 else printf("unable to add quadratic edge U0U1 to the subtetrahedra cas 2.\n");
1423 delete mid; delete quad;
1424 }
1425 delete midEN2[0]; delete midEN2[1]; delete midEN2[2];
1426 }
1427
1428 subTetras.push_back(t1);
1429 subTetras.push_back(t2);
1430 subTetras.push_back(t3);
1431 surfTriangles.push_back(tr);
1432 delete U[0]; delete U[1];
1433 break;
1434 }
1435 case 3 : {// 3 edges cut => 1 tetra + 1 prism => split into 4 tetras
1436 int i1 = IND[0], i2 = IND[1], i3 = IND[2];
1437
1438 if(i1 == 0 && i2 == 3) {
1439 swap(&U[1], &U[2]); swap(i2, i3);
1440 }
1441 //printf("case 3 : ind : %d,%d,%d\n",i1,i2,i3);
1442
1443 int s1 = commonV(tab[i1][0], tab[i1][1], tab[i2][0], tab[i2][1]);
1444 int s2 = commonV(tab[i2][2], tab[i2][3], tab[i3][2], tab[i3][3]);
1445 int s3 = commonV(tab[i1][2], tab[i1][3], tab[i3][2], tab[i3][3]);
1446 int s4 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]);
1447
1448 DI_Tetra *t1 = new DI_Tetra(pt(s1), U[0], U[1], U[2]);
1449 DI_Triangle *tr = new DI_Triangle(U[0], U[1], U[2], tag);
1450 DI_Tetra *t2, *t3, *t4;
1451 bestQuality(U[0], U[1], U[2], pt(s2), pt(s3), pt(s4), &t2, &t3, &t4, QError);
1452 /*t2 = DI_Tetra(U[0],V[0],W[0], x(s2),y(s2),z(s2), U[1],V[1],W[1], U[2],V[2],W[2]);
1453 t3 = DI_Tetra(U[2],V[2],W[2], U[1],V[1],W[1], x(s3),y(s3),z(s3), x(s2),y(s2),z(s2));
1454 t4 = DI_Tetra(U[2],V[2],W[2], x(s3),y(s3),z(s3), x(s4),y(s4),z(s4), x(s2),y(s2),z(s2));*/
1455 if(quadT){
1456 DI_Point *midEN2[3]; // intersection between ls and the bissector between the cutting points
1457 midEN2[0] = quadMidNode(U[2], U[1], pt(s3), e, RPNi);
1458 midEN2[1] = quadMidNode(U[1], U[0], pt(s2), e, RPNi);
1459 midEN2[2] = quadMidNode(U[2], U[0], pt(s2), e, RPNi);
1460 tr->addQuadEdge(1, midEN2[0], e, RPNi);
1461 if(t3->addQuadEdge(0, midEN2[0], e, RPNi)) {
1462 t1->addQuadEdge(5, midEN2[0], e, RPNi);
1463 t2->addQuadEdge(4, midEN2[0], e, RPNi);
1464 }
1465 else {
1466 DI_Point *mid = middle(pt(s3), pt(s4), e, RPNi);
1467 DI_Point *quad = middle(mid, midEN2[0], e, RPNi);
1468 t3->addQuadEdge(1, quad, e, RPNi);
1469 t4->addQuadEdge(0, quad, e, RPNi);
1470 if(t3->addQuadEdge(0, midEN2[0], e, RPNi)) {
1471 t1->addQuadEdge(5, midEN2[0], e, RPNi);
1472 t2->addQuadEdge(4, midEN2[0], e, RPNi);
1473 }
1474 else printf("unable to add quadratic edge U1U2 to the subtetrahedra cas 3.\n");
1475 delete mid; delete quad;
1476 }
1477
1478 tr->addQuadEdge(0, midEN2[1], e, RPNi);
1479 if(t2->addQuadEdge(1, midEN2[1], e, RPNi)) {
1480 t1->addQuadEdge(0, midEN2[1], e, RPNi);
1481 }
1482 else {
1483 DI_Point *mid = middle(pt(s2), pt(s3), e, RPNi);
1484 DI_Point *quad = middle(mid, midEN2[1], e, RPNi);
1485 t2->addQuadEdge(3, quad, e, RPNi);
1486 t3->addQuadEdge(5, quad, e, RPNi);
1487 if(t2->addQuadEdge(1, midEN2[1], e, RPNi)) {
1488 t1->addQuadEdge(0, midEN2[1], e, RPNi);
1489 }
1490 else printf("unable to add quadratic edge U1U2 to the subtetrahedra cas 3.\n");
1491 delete mid; delete quad;
1492 }
1493
1494 tr->addQuadEdge(2, midEN2[2], e, RPNi);
1495 if(t2->addQuadEdge(2, midEN2[2], e, RPNi)) {
1496 t1->addQuadEdge(2, midEN2[2], e, RPNi);
1497 }
1498 else {
1499 DI_Point *mid = middle(pt(s2), pt(s4), e, RPNi);
1500 DI_Point *quad = middle(mid, midEN2[2], e, RPNi);
1501 t2->addQuadEdge(5, quad, e, RPNi);
1502 t3->addQuadEdge(2, quad, e, RPNi);
1503 t4->addQuadEdge(2, quad, e, RPNi);
1504 if(t2->addQuadEdge(2, midEN2[2], e, RPNi)) {
1505 t1->addQuadEdge(2, midEN2[2], e, RPNi);
1506 }
1507 else printf("unable to add quadratic edge U0U2 to the subtetrahedra cas 3.\n");
1508 delete mid; delete quad;
1509 }
1510 delete midEN2[0]; delete midEN2[1]; delete midEN2[2];
1511 }
1512
1513 subTetras.push_back(t1);
1514 subTetras.push_back(t2);
1515 subTetras.push_back(t3);
1516 subTetras.push_back(t4);
1517 surfTriangles.push_back(tr);
1518 delete U[0]; delete U[1]; delete U[2];
1519 break;
1520 }
1521 case 4 : {// 4 edges cut => 2 prisms => split into 6 tetras
1522 int i1 = IND[0], i2 = IND[1], i3 = IND[2], i4 = IND[3];
1523
1524 if(i1 == 0 && i2 == 2) {
1525 swap(&U[0], &U[1]); swap(i1, i2);
1526 }
1527 else if(i1 == 1 && i2 == 2) {
1528 swap(&U[2], &U[3]); swap(i3, i4);
1529 }
1530
1531 int s1 = commonV(tab[i1][0], tab[i1][1], tab[i2][0], tab[i2][1]);
1532 int s2 = commonV(tab[i1][0], tab[i1][1], tab[i2][2], tab[i2][3]);
1533 int s3 = commonV(tab[i1][2], tab[i1][3], tab[i2][0], tab[i2][1]);
1534 int s4 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]);
1535
1536 DI_Tetra *t1, *t2, *t3, *t4, *t5, *t6;
1537 DI_Triangle *tr1, *tr2;
1538 bestQuality(U[0], U[1], U[2], U[3], &tr1, &tr2);
1539 tr1->setLsTag(tag); tr2->setLsTag(tag);
1540 bestQuality(pt(s1), U[0], U[1], pt(s4), U[3], U[2], &t1, &t2, &t3, QError);
1541 bestQuality(pt(s2), U[0], U[3], pt(s3), U[1], U[2], &t4, &t5, &t6, QError);
1542 /*t1 = DI_Tetra(U[0],V[0],W[0], x(s2),y(s2),z(s2), U[1],V[1],W[1], U[3],V[3],W[3]);
1543 t2 = DI_Tetra(U[1],V[1],W[1], U[3],V[3],W[3], x(s2),y(s2),z(s2), U[2],V[2],W[2]);
1544 t3 = DI_Tetra(U[1],V[1],W[1], x(s2),y(s2),z(s2), x(s3),y(s3),z(s3), U[2],V[2],W[2]);
1545 t4 = DI_Tetra(U[0],V[0],W[0], U[1],V[1],W[1], x(s1),y(s1),z(s1), x(s4),y(s4),z(s4));
1546 t5 = DI_Tetra(U[0],V[0],W[0], U[1],V[1],W[1], x(s4),y(s4),z(s4), U[3],V[3],W[3]);
1547 t6 = DI_Tetra(U[3],V[3],W[3], U[1],V[1],W[1], x(s4),y(s4),z(s4), U[2],V[2],W[2]);
1548 tr1 = DI_Triangle(U[0],V[0],W[0], U[1],V[1],W[1], U[3],V[3],W[3]);
1549 tr2 = DI_Triangle(U[1],V[1],W[1], U[2],V[2],W[2], U[3],V[3],W[3]);*/
1550
1551 if(quadT){
1552 DI_Point *midEN2[5]; // intersection between ls and the bissector between the cutting points
1553 midEN2[0] = quadMidNode(U[1], U[2], pt(s4), e, RPNi);
1554 midEN2[1] = quadMidNode(U[1], U[0], pt(s2), e, RPNi);
1555 midEN2[2] = quadMidNode(U[0], U[3], pt(s4), e, RPNi);
1556 midEN2[3] = quadMidNode(U[2], U[3], pt(s2), e, RPNi);
1557 midEN2[4] = quadMidNode(U[1], U[3], pt(s2), e, RPNi);
1558 //quadMidNode(U[1], U[3], (pt(s2)+pt(s3))/2, (pt(s4)+pt(s1))/2, e, RPNi);
1559 tr2->addQuadEdge(0, midEN2[0], e, RPNi);
1560 if(t6->addQuadEdge(5, midEN2[0], e, RPNi)) {
1561 t2->addQuadEdge(2, midEN2[0], e, RPNi);
1562 t3->addQuadEdge(2, midEN2[0], e, RPNi);
1563 }
1564 else {
1565 DI_Point *mid = middle(pt(s1), pt(s4), e, RPNi);
1566 DI_Point *quad = middle(mid, midEN2[0], e, RPNi);
1567 t4->addQuadEdge(5, quad, e, RPNi);
1568 t5->addQuadEdge(3, quad, e, RPNi);
1569 t6->addQuadEdge(3, quad, e, RPNi);
1570 if(t6->addQuadEdge(5, midEN2[0], e, RPNi)) {
1571 t2->addQuadEdge(2, midEN2[0], e, RPNi);
1572 t3->addQuadEdge(2, midEN2[0], e, RPNi);
1573 }
1574 else printf("unable to add quadratic edge U1U2 to the subtetrahedra cas 4.\n");
1575 delete mid; delete quad;
1576 }
1577
1578 tr1->addQuadEdge(0, midEN2[1], e, RPNi);
1579 if(t1->addQuadEdge(1, midEN2[1], e, RPNi)) {
1580 t4->addQuadEdge(0, midEN2[1], e, RPNi);
1581 t5->addQuadEdge(0, midEN2[1], e, RPNi);
1582 }
1583 else {
1584 DI_Point *mid = middle(pt(s2), pt(s3), e, RPNi);
1585 DI_Point *quad = middle(mid, midEN2[1], e, RPNi);
1586 t1->addQuadEdge(3, quad, e, RPNi);
1587 t2->addQuadEdge(1, quad, e, RPNi);
1588 t3->addQuadEdge(0, quad, e, RPNi);
1589 if(t1->addQuadEdge(1, midEN2[1], e, RPNi)) {
1590 t4->addQuadEdge(0, midEN2[1], e, RPNi);
1591 t5->addQuadEdge(0, midEN2[1], e, RPNi);
1592 }
1593 else printf("unable to add quadratic edge U0U1 to the subtetrahedra cas 4.\n");
1594 delete mid; delete quad;
1595 }
1596
1597 tr1->addQuadEdge(2, midEN2[2], e, RPNi);
1598 if(t5->addQuadEdge(2, midEN2[2], e, RPNi)) {
1599 t1->addQuadEdge(2, midEN2[2], e, RPNi);
1600 }
1601 else {
1602 DI_Point *mid = middle(pt(s1), pt(s4), e, RPNi);
1603 DI_Point *quad = middle(mid, midEN2[2], e, RPNi);
1604 t4->addQuadEdge(2, quad, e, RPNi);
1605 t5->addQuadEdge(1, quad, e, RPNi);
1606 if(t5->addQuadEdge(2, midEN2[2], e, RPNi)) {
1607 t1->addQuadEdge(2, midEN2[2], e, RPNi);
1608 }
1609 else printf("unable to add quadratic edge U0U3 to the subtetrahedra cas 4.\n");
1610 delete mid; delete quad;
1611 }
1612
1613 tr2->addQuadEdge(1, midEN2[3], e, RPNi);
1614 if(t2->addQuadEdge(5, midEN2[3], e, RPNi)) {
1615 t6->addQuadEdge(2, midEN2[3], e, RPNi);
1616 }
1617 else {
1618 DI_Point *mid = middle(pt(s2), pt(s3), e, RPNi);
1619 DI_Point *quad = middle(mid, midEN2[3], e, RPNi);
1620 t2->addQuadEdge(4, quad, e, RPNi);
1621 t3->addQuadEdge(5, quad, e, RPNi);
1622 if(t2->addQuadEdge(5, midEN2[3], e, RPNi)) {
1623 t6->addQuadEdge(2, midEN2[3], e, RPNi);
1624 }
1625 else printf("unable to add quadratic edge U2U3 to the subtetrahedra cas 4.\n");
1626 delete mid; delete quad;
1627 }
1628
1629 tr1->addQuadEdge(1, midEN2[4], e, RPNi);
1630 tr2->addQuadEdge(2, midEN2[4], e, RPNi);
1631 t1->addQuadEdge(4, midEN2[4], e, RPNi);
1632 t2->addQuadEdge(0, midEN2[4], e, RPNi);
1633 t5->addQuadEdge(5, midEN2[4], e, RPNi);
1634 t6->addQuadEdge(0, midEN2[4], e, RPNi);
1635 for(int i = 0; i < 5; i++) delete midEN2[i];
1636 }
1637
1638 subTetras.push_back(t1);
1639 subTetras.push_back(t2);
1640 subTetras.push_back(t3);
1641 subTetras.push_back(t4);
1642 subTetras.push_back(t5);
1643 subTetras.push_back(t6);
1644 surfTriangles.push_back(tr1);
1645 surfTriangles.push_back(tr2);
1646 for(int i = 0; i < 4; i++) delete U[i];
1647 break;
1648 }
1649 default:
1650 printf("Error : %d edge(s) cut in the tetrahedron (ls : %g %g %g %g)\n",
1651 COUNT, ls(0), ls(1), ls(2), ls(3));
1652 }
1653 }
1654
1655 // Split a reference Hexahedron into 6 tetrahedra
splitIntoTetras(std::vector<DI_Tetra * > & tetras) const1656 void DI_Hexa::splitIntoTetras(std::vector<DI_Tetra *> &tetras) const
1657 {
1658 tetras.push_back(new DI_Tetra(pt(0), pt(1), pt(3), pt(4))); // 0134
1659 tetras.push_back(new DI_Tetra(pt(1), pt(4), pt(5), pt(7))); // 1457
1660 tetras.push_back(new DI_Tetra(pt(1), pt(3), pt(4), pt(7))); // 1347
1661 tetras.push_back(new DI_Tetra(pt(2), pt(5), pt(6), pt(7))); // 2567
1662 tetras.push_back(new DI_Tetra(pt(1), pt(2), pt(3), pt(7))); // 1237
1663 tetras.push_back(new DI_Tetra(pt(1), pt(5), pt(2), pt(7))); // 1527
1664 }
1665
1666 // -----------------------------------------------------------------------------
1667 // -------------------INTEGRATION POINTS --------------------------------------
1668 // -----------------------------------------------------------------------------
1669
1670 // return the integration points in the reference line
getRefIntegrationPoints(const int polynomialOrder,std::vector<DI_IntegrationPoint * > & ip) const1671 void DI_Line::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
1672 int N = getNGQLPts(polynomialOrder);
1673 IntPt* intpt (getGQLPts(polynomialOrder));
1674 for (int i = 0; i < N; ++i){
1675 ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
1676 }
1677 }
1678
1679 // return the integration points in the reference triangle
getRefIntegrationPoints(const int polynomialOrder,std::vector<DI_IntegrationPoint * > & ip) const1680 void DI_Triangle::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
1681 int pO = polynomialOrder;
1682 if(pO == 11 || pO == 16 || pO == 18 || pO == 20) pO++;
1683 if(pO == 15) pO = 17;
1684 int N = getNGQTPts(pO);
1685 IntPt* intpt (getGQTPts(pO));
1686 for (int i = 0; i < N; ++i){
1687 ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
1688 }
1689 }
1690
1691 // return the integration points in the reference triangle
getRefIntegrationPoints(const int polynomialOrder,std::vector<DI_IntegrationPoint * > & ip) const1692 void DI_Quad::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
1693 int N = getNGQQPts(polynomialOrder);
1694 IntPt* intpt (getGQQPts(polynomialOrder));
1695 for (int i = 0; i < N; ++i){
1696 ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
1697 }
1698 }
1699
1700 // return the integration points in the reference tetra
getRefIntegrationPoints(const int polynomialOrder,std::vector<DI_IntegrationPoint * > & ip) const1701 void DI_Tetra::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
1702 int pO = polynomialOrder;
1703 if(pO == 9) pO++;
1704 int N = getNGQTetPts(pO);
1705 IntPt* intpt (getGQTetPts(pO));
1706 for (int i = 0; i < N; ++i){
1707 ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
1708 }
1709 }
1710
1711 // return the integration points in the reference Cube
getRefIntegrationPoints(const int polynomialOrder,std::vector<DI_IntegrationPoint * > & ip) const1712 void DI_Hexa::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
1713 int N = getNGQHPts(polynomialOrder);
1714 IntPt* intpt (getGQHPts(polynomialOrder));
1715 for (int i = 0; i < N; ++i){
1716 ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
1717 }
1718 }
1719
1720 // -----------------------------------------------------------------------------
1721 // ----------------------------- CUTTING --------------------------------------
1722 // -----------------------------------------------------------------------------
1723
1724 // cut a line into sublines along the levelset curve
cut(std::vector<gLevelset * > & RPN,std::vector<DI_IntegrationPoint * > & ip,DI_Point::Container & cp,const int polynomialOrder,std::vector<DI_Line * > & lines,int recurLevel,double ** nodeLs) const1725 bool DI_Line::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
1726 DI_Point::Container &cp, const int polynomialOrder,
1727 std::vector<DI_Line *> &lines, int recurLevel, double **nodeLs) const
1728 {
1729 std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
1730
1731 DI_Line *ll = new DI_Line(-1, 0, 0, 1, 0, 0, 2.); //reference line
1732 ll->setPolynomialOrder(getPolynomialOrder());
1733 std::vector<DI_Line *> ll_subLines;
1734 std::vector<DI_CuttingPoint *> ll_cp;
1735
1736 RecurElement *re = new RecurElement(ll);
1737 bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
1738 pushSubElements(re, ll_subLines);
1739 delete re;
1740
1741 int iPrim = 0;
1742 if(signChange){
1743 for(int l = 0; l < (int)RPN.size(); l++) {
1744 gLevelset *Lsi = RPN[l];
1745 RPNi.push_back(Lsi);
1746 if(Lsi->isPrimitive()) {
1747 ll->addLs(this, Lsi, iPrim++, nodeLs);
1748 int nbSubLn = ll_subLines.size();
1749 int nbCP = ll_cp.size();
1750 for(int i = 0; i < nbSubLn; i++) ll_subLines[i]->addLs(ll);
1751 for(int i = 0; i < nbCP; i++) ll_cp[i]->addLs(ll);
1752
1753 for(int ln = 0; ln < nbSubLn; ln++){
1754 DI_Line *lnt = ll_subLines[0];
1755 ll_subLines.erase(ll_subLines.begin());
1756 bool c = lnt->cut(ll, RPNi, ll_subLines, ll_cp);
1757 if(c) delete lnt;
1758 }
1759 }
1760 else {
1761 for(int l = 0; l < (int)ll_subLines.size(); l++)
1762 ll_subLines[l]->chooseLs(Lsi);
1763 for(int p = 0; p < (int)ll_cp.size(); p++)
1764 ll_cp[p]->chooseLs(Lsi);
1765 }
1766 }
1767 }
1768 else{
1769 for(int l = 0; l < (int)RPN.size(); l++) {
1770 gLevelset *Lsi = RPN[l];
1771 if(Lsi->isPrimitive()) {
1772 ll->addLs(this, Lsi, iPrim++, nodeLs);
1773 }
1774 }
1775 if(iPrim == 1) iPrim--;
1776 ll_subLines[0]->addLs(this, RPN.back(), iPrim, nodeLs);
1777 }
1778
1779 for(int l = 0; l < (int)ll_subLines.size(); l++) {
1780 ll_subLines[l]->computeLsTagDom(ll, RPN);
1781 DI_Line *ll_subLn = new DI_Line(*ll_subLines[l]);
1782 mappingEl(ll_subLines[l]);
1783 ll_subLines[l]->integrationPoints(polynomialOrder, ll_subLn, ll, RPN, ip);
1784 lines.push_back(ll_subLines[l]);
1785 delete ll_subLn;
1786 }
1787 for(int p = 0; p < (int)ll_cp.size(); p++) {
1788 if(ll_cp[p]->ls() != 0) {delete ll_cp[p]; continue;}
1789 mappingCP(ll_cp[p]);
1790 std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(ll_cp[p]);
1791 if(!isIn.second) delete ll_cp[p];
1792 }
1793 delete ll;
1794 return signChange;
1795 }
1796
1797 // cut a line into sublines along one levelset primitive and return true if it is cut
cut(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Line * > & subLines,std::vector<DI_CuttingPoint * > & cp)1798 bool DI_Line::cut(const DI_Element *e, const std::vector<gLevelset *> &RPNi,
1799 std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint *> &cp)
1800 {
1801 // check if the line is cut by the level set
1802 int on = 0, pos = 0, neg = 0, ze[2];
1803 for (int i = 0; i < 2; i++){
1804 if(pt(i)->isOnBorder()) ze[on++] = i;
1805 else if (ls(i) > 0.) pos++;
1806 else neg++;
1807 }
1808 if(pos && neg)
1809 selfSplit(e, RPNi, subLines, cp);
1810 else {
1811 for(int i = 0; i < on; i++)
1812 cp.push_back(new DI_CuttingPoint(pt(ze[i])));
1813 subLines.push_back(this);
1814 }
1815 return (pos && neg);
1816 }
1817
1818 // cut a triangle into subtriangles along the levelset curve
cut(std::vector<gLevelset * > & RPN,std::vector<DI_IntegrationPoint * > & ip,std::vector<DI_IntegrationPoint * > & ipS,DI_Point::Container & cp,const int polynomialOrderQ,const int polynomialOrderTr,const int polynomialOrderL,std::vector<DI_Quad * > & subQuads,std::vector<DI_Triangle * > & subTriangles,std::vector<DI_Line * > & surfLines,int recurLevel,double ** nodeLs) const1819 bool DI_Triangle::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
1820 std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
1821 const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
1822 std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
1823 std::vector<DI_Line *> &surfLines, int recurLevel, double **nodeLs) const
1824 {
1825 std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
1826
1827 DI_Triangle *tt = new DI_Triangle(0, 0, 0, 1, 0, 0, 0, 1, 0, 0.5); //reference triangle
1828 tt->setPolynomialOrder(getPolynomialOrder());
1829 std::vector<DI_Quad *> tt_subQuads;
1830 std::vector<DI_Triangle *> tt_subTriangles;
1831 std::vector<DI_Line *> tt_surfLines;
1832 std::vector<DI_CuttingPoint *> tt_cp;
1833
1834 RecurElement *re = new RecurElement(tt);
1835 bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
1836 pushSubElements(re, tt_subTriangles);
1837 delete re;
1838
1839 int iPrim = 0;
1840 if(signChange){
1841 for(int l = 0; l < (int)RPN.size(); l++) {
1842 gLevelset *Lsi = RPN[l];
1843 RPNi.push_back(Lsi);
1844 if(Lsi->isPrimitive()) {
1845 tt->addLs(this, Lsi, iPrim++, nodeLs);
1846 int nbSubQ = tt_subQuads.size(), nbSubQ1 = nbSubQ;
1847 int nbSubTr = tt_subTriangles.size(), nbSubTr1 = nbSubTr;
1848 int nbSurfLn = tt_surfLines.size(), nbSurfLn1 = nbSurfLn;
1849 int nbCP = tt_cp.size();
1850 for(int i = 0; i < nbSubQ; i++) tt_subQuads[i]->addLs(tt);
1851 for(int i = 0; i < nbSubTr; i++) tt_subTriangles[i]->addLs(tt);
1852 for(int i = 0; i < nbSurfLn; i++) tt_surfLines[i]->addLs(tt);
1853 for(int i = 0; i < nbCP; i++) tt_cp[i]->addLs(tt);
1854 int ze[3], cz = 0;
1855 for (int i = 0; i < 3; i++) if(tt->pt(i)->isOnBorder()) ze[cz++] = i;
1856
1857 for(int q = 0; q < nbSubQ; q++) {
1858 nbSubQ1 = tt_subQuads.size();
1859 nbSurfLn1 = tt_surfLines.size();
1860 DI_Quad *qt = tt_subQuads[0];
1861 tt_subQuads.erase(tt_subQuads.begin());
1862 bool c = qt->cut(tt, RPNi, tt_subQuads, tt_subTriangles, tt_surfLines, tt_cp);
1863 if(c) delete qt;
1864 if((int)tt_surfLines.size() - nbSurfLn1 == 1 && (int)tt_subQuads.size() == nbSubQ1){ // 1 line created on boundary of the quad
1865 if (cz == 2) { // 1 line created on boundary of the big triangle => check among surfLines
1866 DI_Line lf (pt(ze[0]), pt(ze[1]));
1867 for(int i = (int)surfLines.size() - 1; i >= 0; i--){
1868 if (lf.contain(surfLines[i])){ // line allready created on another surface => pop the new one
1869 delete tt_surfLines.back(); tt_surfLines.pop_back(); break;
1870 }
1871 }
1872 }
1873 else { // 1 line created inside the big triangle => check among tt_surfLines
1874 if(isLastLnInV(tt_surfLines)) {
1875 delete tt_surfLines.back(); tt_surfLines.pop_back();
1876 }
1877 }
1878 }
1879 }
1880 for(int t = 0; t < nbSubTr; t++){
1881 nbSubTr1 = tt_subTriangles.size();
1882 nbSurfLn1 = tt_surfLines.size();
1883 DI_Triangle *trt = tt_subTriangles[0];
1884 tt_subTriangles.erase(tt_subTriangles.begin());
1885 bool c = trt->cut(tt, RPNi, tt_subQuads, tt_subTriangles, tt_surfLines, tt_cp);
1886 if(c) delete trt;
1887 if((int)tt_surfLines.size() - nbSurfLn1 == 1 && (int)tt_subTriangles.size() == nbSubTr1) { // 1 line created on boundary of the triangle
1888 if(cz == 2) { // 1 line created on boundary of the big triangle => check among surfLines
1889 DI_Line lf(pt(ze[0]), pt(ze[1]));
1890 for(int i = (int)surfLines.size() - 1; i >= 0; i--){
1891 if(lf.contain(surfLines[i])) {
1892 delete tt_surfLines.back(); tt_surfLines.pop_back(); break;
1893 }
1894 }
1895 }
1896 else{ // 1 line created inside the big triangle => check among tt_surfLines
1897 if(isLastLnInV(tt_surfLines)) {
1898 delete tt_surfLines.back(); tt_surfLines.pop_back();
1899 }
1900 }
1901 }
1902 }
1903 for(int ln = 0; ln < nbSurfLn; ln++){
1904 DI_Line* lnt = tt_surfLines[0];
1905 tt_surfLines.erase(tt_surfLines.begin());
1906 bool c = lnt->cut(tt, RPNi, tt_surfLines, tt_cp);
1907 if(c) delete lnt;
1908 }
1909 }
1910 else {
1911 for(int q = 0; q < (int)tt_subQuads.size(); q++)
1912 tt_subQuads[q]->chooseLs(Lsi);
1913 for(int t = 0; t < (int)tt_subTriangles.size(); t++)
1914 tt_subTriangles[t]->chooseLs(Lsi);
1915 for(int l = 0; l < (int)tt_surfLines.size(); l++)
1916 tt_surfLines[l]->chooseLs(Lsi);
1917 for(int p = 0; p < (int)tt_cp.size(); p++)
1918 tt_cp[p]->chooseLs(Lsi);
1919 }
1920 }
1921 }
1922 else{
1923 for(int l = 0; l < (int)RPN.size(); l++) {
1924 gLevelset *Lsi = RPN[l];
1925 if(Lsi->isPrimitive()) {
1926 tt->addLs(this, Lsi, iPrim++, nodeLs);
1927 }
1928 }
1929 if(iPrim == 1) iPrim--;
1930 tt_subTriangles[0]->addLs(this, RPN.back(), iPrim, nodeLs);
1931 }
1932
1933 for(int q = 0; q < (int)tt_subQuads.size(); q++) {
1934 tt_subQuads[q]->computeLsTagDom(tt, RPN);
1935 DI_Quad *tt_subQ = new DI_Quad(*tt_subQuads[q]);
1936 mappingEl(tt_subQuads[q]);
1937 tt_subQuads[q]->integrationPoints(polynomialOrderQ, tt_subQ, tt, RPN, ip);
1938 subQuads.push_back(tt_subQuads[q]);
1939 delete tt_subQ;
1940 }
1941 for(int t = 0; t < (int)tt_subTriangles.size(); t++) {
1942 tt_subTriangles[t]->computeLsTagDom(tt, RPN);
1943 DI_Triangle *tt_subTr = new DI_Triangle(*tt_subTriangles[t]);
1944 mappingEl(tt_subTriangles[t]);
1945 tt_subTriangles[t]->integrationPoints(polynomialOrderTr, tt_subTr, tt, RPN, ip);
1946 subTriangles.push_back(tt_subTriangles[t]);
1947 delete tt_subTr;
1948 }
1949 for(int l = 0; l < (int)tt_surfLines.size(); l++) {
1950 tt_surfLines[l]->computeLsTagBound(tt_subQuads, tt_subTriangles);
1951 if(tt_surfLines[l]->lsTag() == -1) {delete tt_surfLines[l]; continue;}
1952 DI_Line *tt_surfLn = new DI_Line(*tt_surfLines[l]);
1953 mappingEl(tt_surfLines[l]);
1954 tt_surfLines[l]->integrationPoints(polynomialOrderL, tt_surfLn, tt, RPN, ipS);
1955 surfLines.push_back(tt_surfLines[l]);
1956 delete tt_surfLn;
1957 }
1958 for(int p = 0; p < (int)tt_cp.size(); p++) {
1959 if(tt_cp[p]->ls() != 0) {delete tt_cp[p]; continue;}
1960 mappingCP(tt_cp[p]);
1961 std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(tt_cp[p]);
1962 if(!isIn.second) delete tt_cp[p];
1963 }
1964 delete tt;
1965 return signChange;
1966 }
1967
1968 // cut a triangle into subtriangles along one levelset primitive
cut(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Quad * > & subQuads,std::vector<DI_Triangle * > & subTriangles,std::vector<DI_Line * > & surfLines,std::vector<DI_CuttingPoint * > & cp)1969 bool DI_Triangle::cut (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
1970 std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
1971 std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cp)
1972 {
1973 // check if the triangle is cut by the level set
1974 int on = 0, pos = 0, neg = 0, ze[3];
1975 for (int i = 0; i < 3; i++){
1976 if(pt(i)->isOnBorder()) ze[on++] = i;
1977 else if (pt(i)->ls() > 0.) pos++;
1978 else neg++;
1979 }
1980 if(pos && neg)
1981 selfSplit(e, RPNi, subQuads, subTriangles, surfLines, cp);
1982 else {
1983 if(on == 2){ // the level set is zero on an edge of the triangle
1984 surfLines.push_back(new DI_Line(pt(ze[0]), pt(ze[1]), RPNi.back()->getTag()));
1985 // add the line twice if the edge belongs to 2 triangles => remove it later!
1986 }
1987 if(on == 3)
1988 printf("Warning : triangle with zero levelset on every vertex.\n");
1989 for(int i = 0; i < on; i++)
1990 cp.push_back(new DI_CuttingPoint(pt(ze[i])));
1991 subTriangles.push_back(this);
1992 }
1993 return (pos && neg);
1994 }
1995
1996 // cut a quadrangle into subtriangles along the levelset curve
cut(std::vector<gLevelset * > & RPN,std::vector<DI_IntegrationPoint * > & ip,std::vector<DI_IntegrationPoint * > & ipS,DI_Point::Container & cp,const int polynomialOrderQ,const int polynomialOrderTr,const int polynomialOrderL,std::vector<DI_Quad * > & subQuads,std::vector<DI_Triangle * > & subTriangles,std::vector<DI_Line * > & surfLines,int recurLevel,double ** nodeLs) const1997 bool DI_Quad::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
1998 std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
1999 const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
2000 std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
2001 std::vector<DI_Line *> &surfLines, int recurLevel, double **nodeLs) const
2002 {
2003 std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
2004
2005 DI_Quad *qq = new DI_Quad(-1, -1, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, 4.); // reference quad
2006 qq->setPolynomialOrder(getPolynomialOrder());
2007 std::vector<DI_Quad *> qq_subQuads;
2008 std::vector<DI_Triangle *> qq_subTriangles;
2009 std::vector<DI_Line *> qq_surfLines;
2010 std::vector<DI_CuttingPoint *> qq_cp;
2011
2012 RecurElement *re = new RecurElement(qq);
2013 bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
2014 pushSubElements(re, qq_subQuads);
2015 delete re;
2016
2017 int iPrim = 0;
2018 if(signChange) {
2019 for(int l = 0; l < (int)RPN.size(); l++) {
2020 gLevelset *Lsi = RPN[l];
2021 RPNi.push_back(Lsi);
2022 if(Lsi->isPrimitive()) {
2023 qq->addLs(this, Lsi, iPrim++, nodeLs);
2024 int nbSubQ = qq_subQuads.size(), nbSubQ1 = nbSubQ;
2025 int nbSubTr = qq_subTriangles.size(), nbSubTr1 = nbSubTr;
2026 int nbSurfLn = qq_surfLines.size(), nbSurfLn1 = nbSurfLn;
2027 int nbCP = qq_cp.size();
2028 for(int i = 0; i < nbSubQ; i++) qq_subQuads[i]->addLs(qq);
2029 for(int i = 0; i < nbSubTr; i++) qq_subTriangles[i]->addLs(qq);
2030 for(int i = 0; i < nbSurfLn; i++) qq_surfLines[i]->addLs(qq);
2031 for(int i = 0; i < nbCP; i++) qq_cp[i]->addLs(qq);
2032 int pos = 0, neg = 0, ze[4], cz = 0;
2033 for (int i = 0; i < 4; i++){
2034 if(qq->pt(i)->isOnBorder()) ze[cz++] = i;
2035 else if (qq->ls(i) > 0.) pos++;
2036 else neg++;
2037 }
2038
2039 for(int q = 0; q < nbSubQ; q++) {
2040 nbSubQ1 = qq_subQuads.size();
2041 nbSurfLn1 = qq_surfLines.size();
2042 DI_Quad *qt = qq_subQuads[0];
2043 qq_subQuads.erase(qq_subQuads.begin());
2044 bool c = qt->cut(qq, RPNi, qq_subQuads, qq_subTriangles, qq_surfLines, qq_cp);
2045 if(c) delete qt;
2046 if((int)qq_surfLines.size() - nbSurfLn1 == 1 && (int)qq_subQuads.size() == nbSubQ1){ // 1 line created on boundary of the quad
2047 if (cz == 2 && !(pos && neg)) { // 1 line created on boundary of the big quad => check among surfLines
2048 DI_Line lf (pt(ze[0]), pt(ze[1]));
2049 for(int i = (int)surfLines.size() - 1; i >= 0; i--){
2050 if (lf.contain(surfLines[i])){ // line allready created on another surface => pop the new one
2051 delete qq_surfLines.back(); qq_surfLines.pop_back(); break;
2052 }
2053 }
2054 }
2055 else { // 1 line created inside the big quad => check among qq_surfLines
2056 if(isLastLnInV(qq_surfLines)) {
2057 delete qq_surfLines.back(); qq_surfLines.pop_back();
2058 }
2059 }
2060 }
2061 }
2062 for(int t = 0; t < nbSubTr; t++){
2063 nbSubTr1 = qq_subTriangles.size();
2064 nbSurfLn1 = qq_surfLines.size();
2065 DI_Triangle *trt = qq_subTriangles[0];
2066 qq_subTriangles.erase(qq_subTriangles.begin());
2067 bool c = trt->cut(qq, RPNi, qq_subQuads, qq_subTriangles, qq_surfLines, qq_cp);
2068 if(c) delete trt;
2069 if((int)qq_surfLines.size() - nbSurfLn1 == 1 && (int)qq_subTriangles.size() == nbSubTr1) { // 1 line created on boundary of the triangle
2070 if(cz == 2 && !(pos && neg)) { // 1 line created on boundary of the big quad => check among surfLines
2071 DI_Line lf(pt(ze[0]), pt(ze[1]));
2072 for(int i = (int)surfLines.size() - 1; i >= 0; i--){
2073 if(lf.contain(surfLines[i])) {
2074 delete qq_surfLines.back(); qq_surfLines.pop_back(); break;
2075 }
2076 }
2077 }
2078 else{ // 1 line created inside the big quad => check among qq_surfLines
2079 if(isLastLnInV(qq_surfLines)) {
2080 delete qq_surfLines.back(); qq_surfLines.pop_back();
2081 }
2082 }
2083 }
2084 }
2085 for(int ln = 0; ln < nbSurfLn; ln++){
2086 DI_Line *lnt = qq_surfLines[0];
2087 qq_surfLines.erase(qq_surfLines.begin());
2088 bool c = lnt->cut(qq, RPNi, qq_surfLines, qq_cp);
2089 if(c) delete lnt;
2090 }
2091 }
2092 else {
2093 for(int q = 0; q < (int)qq_subQuads.size(); q++)
2094 qq_subQuads[q]->chooseLs(Lsi);
2095 for(int t = 0; t < (int)qq_subTriangles.size(); t++)
2096 qq_subTriangles[t]->chooseLs(Lsi);
2097 for(int l = 0; l < (int)qq_surfLines.size(); l++)
2098 qq_surfLines[l]->chooseLs(Lsi);
2099 for(int p = 0; p < (int)qq_cp.size(); p++)
2100 qq_cp[p]->chooseLs(Lsi);
2101 }
2102 }
2103 }
2104 else {
2105 for(int l = 0; l < (int)RPN.size(); l++) {
2106 gLevelset *Lsi = RPN[l];
2107 if(Lsi->isPrimitive()) {
2108 qq->addLs(this, Lsi, iPrim++, nodeLs);
2109 }
2110 }
2111 if(iPrim == 1) iPrim--;
2112 qq_subQuads[0]->addLs(this, RPN.back(), iPrim, nodeLs);
2113 }
2114
2115 for(int q = 0; q < (int)qq_subQuads.size(); q++) {
2116 qq_subQuads[q]->computeLsTagDom(qq, RPN);
2117 DI_Quad *qq_subQ = new DI_Quad(*qq_subQuads[q]);
2118 mappingEl(qq_subQuads[q]);
2119 qq_subQuads[q]->integrationPoints(polynomialOrderQ, qq_subQ, qq, RPN, ip);
2120 subQuads.push_back(qq_subQuads[q]);
2121 delete qq_subQ;
2122 }
2123 for(int t = 0; t < (int)qq_subTriangles.size(); t++) {
2124 qq_subTriangles[t]->computeLsTagDom(qq, RPN);
2125 DI_Triangle *qq_subTr = new DI_Triangle(*qq_subTriangles[t]);
2126 mappingEl(qq_subTriangles[t]);
2127 qq_subTriangles[t]->integrationPoints(polynomialOrderTr, qq_subTr, qq, RPN, ip);
2128 subTriangles.push_back(qq_subTriangles[t]);
2129 delete qq_subTr;
2130 }
2131 for(int l = 0; l < (int)qq_surfLines.size(); l++) {
2132 qq_surfLines[l]->computeLsTagBound(qq_subQuads, qq_subTriangles);
2133 if(qq_surfLines[l]->lsTag() == -1) {delete qq_surfLines[l]; continue;} //FIXME
2134 DI_Line *qq_surfLn = new DI_Line(*qq_surfLines[l]);
2135 mappingEl(qq_surfLines[l]);
2136 qq_surfLines[l]->integrationPoints(polynomialOrderL, qq_surfLn, qq, RPN, ipS);
2137 surfLines.push_back(qq_surfLines[l]);
2138 delete qq_surfLn;
2139 }
2140 for(int p = 0; p < (int)qq_cp.size(); p++) {
2141 if(qq_cp[p]->ls() != 0) {delete qq_cp[p]; continue;}
2142 mappingCP(qq_cp[p]);
2143 std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(qq_cp[p]);
2144 if(!isIn.second) delete qq_cp[p];
2145 }
2146 delete qq;
2147 return signChange;
2148 }
2149
2150 // cut a quadrangle into subtriangles along the levelset curve
cut(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Quad * > & subQuads,std::vector<DI_Triangle * > & subTriangles,std::vector<DI_Line * > & surfLines,std::vector<DI_CuttingPoint * > & cp)2151 bool DI_Quad::cut (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
2152 std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
2153 std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cp)
2154 {
2155 // check if the quadrangle is cut by the level set
2156 int on = 0, pos = 0, neg = 0, ze[4];
2157 for (int i = 0; i < 4; i++){
2158 if(pt(i)->isOnBorder()) ze[on++] = i;
2159 else if (pt(i)->ls() > 0.) pos++;
2160 else neg++;
2161 }
2162 if (pos && neg) {
2163 std::vector<DI_Triangle *> triangles;
2164 splitIntoTriangles (triangles); // Split the quad into 2 sub triangles
2165 int nl0 = surfLines.size(), nt1, nl1, nt2, nl2;
2166 for (int t = 0; t < (int)triangles.size(); t++) {
2167 nt1 = subTriangles.size(); nl1 = surfLines.size();
2168 triangles[t]->selfSplit(e, RPNi, subQuads, subTriangles, surfLines, cp);
2169 nt2 = subTriangles.size(); nl2 = surfLines.size();
2170 if((nt2 - nt1) == 1 && (nl2 - nl1) == 1){ // only 1 line created on an edge of a triangle => check if it was not yet created
2171 if(isLastLnInV(surfLines, nl0)) {delete surfLines.back(); surfLines.pop_back(); nl2--;}
2172 }
2173 }
2174 }
2175 else {
2176 if(on == 2){ // the level set is zero on an edge of the quad
2177 surfLines.push_back(new DI_Line(pt(ze[0]), pt(ze[1]), RPNi.back()->getTag()));
2178 // we add the line twice if the edge belongs to 2 quads => remove it later!
2179 }
2180 if(on == 4)
2181 printf("Warning : quadrangle with zero levelset on every vertex.\n");
2182 for(int i = 0; i < on; i++)
2183 cp.push_back(new DI_CuttingPoint(pt(ze[i])));
2184 subQuads.push_back(this);
2185 }
2186 return (pos && neg);
2187 }
2188
2189 // cut a tetrahedron into subtetrahedra along the levelset surface
cut(std::vector<gLevelset * > & RPN,std::vector<DI_IntegrationPoint * > & ip,std::vector<DI_IntegrationPoint * > & ipS,DI_Point::Container & cp,const int polynomialOrderT,const int polynomialOrderQ,const int polynomialOrderTr,std::vector<DI_Tetra * > & subTetras,std::vector<DI_Quad * > & surfQuads,std::vector<DI_Triangle * > & surfTriangles,int recurLevel,double ** nodeLs) const2190 bool DI_Tetra::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
2191 std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
2192 const int polynomialOrderT, const int polynomialOrderQ, const int polynomialOrderTr,
2193 std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
2194 std::vector<DI_Triangle *> &surfTriangles, int recurLevel, double **nodeLs) const
2195 {
2196 std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
2197
2198 DI_Tetra *tt = new DI_Tetra(0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1./6.); // reference tetra
2199 tt->setPolynomialOrder(getPolynomialOrder());
2200 std::vector<DI_Hexa *> tt_subHexas;
2201 std::vector<DI_Tetra *> tt_subTetras;
2202 std::vector<DI_Quad *> tt_surfQuads;
2203 std::vector<DI_Triangle *> tt_surfTriangles;
2204 std::vector<DI_Line *> tt_surfLines; // not used ...
2205 std::vector<DI_CuttingPoint *> tt_cp;
2206 std::vector<DI_QualError *> QError;
2207
2208 RecurElement *re = new RecurElement(tt);
2209 bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
2210 pushSubElements(re, tt_subTetras);
2211 delete re;
2212
2213 int iPrim = 0;
2214 if(signChange) {
2215 for(int l = 0; l < (int)RPN.size(); l++) {
2216 gLevelset *Lsi = RPN[l];
2217 RPNi.push_back(Lsi);
2218 if(Lsi->isPrimitive()) {
2219 tt->addLs(this, Lsi, iPrim++, nodeLs);
2220 int nbSubT = tt_subTetras.size(), nbSubT1 = nbSubT;
2221 int nbSurfQ = tt_surfQuads.size();
2222 int nbSurfTr = tt_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
2223 int nbCP = tt_cp.size();
2224 for(int i = 0; i < nbSubT; i++) tt_subTetras[i]->addLs(tt);
2225 for(int i = 0; i < nbSurfQ; i++) tt_surfQuads[i]->addLs(tt);
2226 for(int i = 0; i < nbSurfTr; i++) tt_surfTriangles[i]->addLs(tt);
2227 for(int i = 0; i < nbCP; i++) tt_cp[i]->addLs(tt);
2228 int ze[4], cz = 0;
2229 for (int i = 0; i < 4; i++) if(tt->pt(i)->isOnBorder()) ze[cz++] = i;
2230
2231 for(int t = 0; t < nbSubT; t++) {
2232 nbSurfTr1 = tt_surfTriangles.size();
2233 nbSubT1 = tt_subTetras.size();
2234 DI_Tetra *tet = tt_subTetras[0];
2235 tt_subTetras.erase(tt_subTetras.begin());
2236 bool c = tet->cut(tt, RPNi, tt_subTetras, tt_surfQuads, tt_surfTriangles, tt_cp, QError);
2237 if(c) delete tet;
2238 if((int)tt_surfTriangles.size() - nbSurfTr1 == 1 && (int)tt_subTetras.size() == nbSubT1) { // 1 triangle created on surface of the subtetra
2239 // check among the tt_surfTriangles
2240 if(isLastTrInV(tt_surfTriangles)) {
2241 delete tt_surfTriangles.back(); tt_surfTriangles.pop_back();
2242 }
2243 else if(cz == 3) { // 1 triangle created on surface of the reference tet => check among surfTriangles
2244 DI_Triangle tf(pt(ze[0]), pt(ze[1]), pt(ze[2]));
2245 for(int i = (int)surfTriangles.size() - 1; i >= 0; i--){
2246 if(tf.contain(surfTriangles[i])) {
2247 delete tt_surfTriangles.back(); tt_surfTriangles.pop_back(); break;
2248 }
2249 }
2250 }
2251 }
2252 }
2253 for(int q = 0; q < nbSurfQ; q++) {
2254 DI_Quad *qt = tt_surfQuads[0];
2255 tt_surfQuads.erase(tt_surfQuads.begin());
2256 bool c = qt->cut(tt, RPNi, tt_surfQuads, tt_surfTriangles, tt_surfLines, tt_cp);
2257 if(c) delete qt;
2258 }
2259 for(int t = 0; t < nbSurfTr; t++){
2260 DI_Triangle *trt = tt_surfTriangles[0];
2261 tt_surfTriangles.erase(tt_surfTriangles.begin());
2262 bool c = trt->cut(tt, RPNi, tt_surfQuads, tt_surfTriangles, tt_surfLines, tt_cp);
2263 if(c) delete trt;
2264 }
2265 }
2266 else {
2267 for(int t = 0; t < (int)tt_subTetras.size(); t++)
2268 tt_subTetras[t]->chooseLs(Lsi);
2269 for(int q = 0; q < (int)tt_surfQuads.size(); q++)
2270 tt_surfQuads[q]->chooseLs(Lsi);
2271 for(int t = 0; t < (int)tt_surfTriangles.size(); t++)
2272 tt_surfTriangles[t]->chooseLs(Lsi);
2273 for(int p = 0; p < (int)tt_cp.size(); p++)
2274 tt_cp[p]->chooseLs(Lsi);
2275 }
2276 }
2277 }
2278 else {
2279 for(int l = 0; l < (int)RPN.size(); l++) {
2280 gLevelset *Lsi = RPN[l];
2281 if(Lsi->isPrimitive()) {
2282 tt->addLs(this, Lsi, iPrim++, nodeLs);
2283 }
2284 }
2285 if(iPrim == 1) iPrim--;
2286 tt_subTetras[0]->addLs(this, RPN.back(), iPrim, nodeLs);
2287 }
2288
2289 for(int i = 0; i < (int)QError.size(); i++) {
2290 QError[i]->print(this);
2291 delete QError[i];
2292 }
2293
2294 for(int t = 0; t < (int)tt_subTetras.size(); t++) {
2295 tt_subTetras[t]->computeLsTagDom(tt, RPN);
2296 DI_Tetra *tt_subT = new DI_Tetra(*tt_subTetras[t]);
2297 mappingEl(tt_subTetras[t]);
2298 tt_subTetras[t]->integrationPoints(polynomialOrderT, tt_subT, tt, RPN, ip);
2299 subTetras.push_back(tt_subTetras[t]);
2300 delete tt_subT;
2301 }
2302 for(int q = 0; q < (int)tt_surfQuads.size(); q++) {
2303 tt_surfQuads[q]->computeLsTagBound(tt_subHexas, tt_subTetras);
2304 if(tt_surfQuads[q]->lsTag() == -1) {delete tt_surfQuads[q]; continue;}
2305 DI_Quad *tt_surfQ = new DI_Quad(*tt_surfQuads[q]);
2306 mappingEl(tt_surfQuads[q]);
2307 tt_surfQuads[q]->integrationPoints(polynomialOrderQ, tt_surfQ, tt, RPN, ipS);
2308 surfQuads.push_back(tt_surfQuads[q]);
2309 delete tt_surfQ;
2310 }
2311 for(int t = 0; t < (int)tt_surfTriangles.size(); t++) {
2312 tt_surfTriangles[t]->computeLsTagBound(tt_subHexas, tt_subTetras);
2313 if(tt_surfTriangles[t]->lsTag() == -1) {delete tt_surfTriangles[t]; continue;}
2314 DI_Triangle *tt_surfTr = new DI_Triangle(*tt_surfTriangles[t]);
2315 mappingEl(tt_surfTriangles[t]);
2316 tt_surfTriangles[t]->integrationPoints(polynomialOrderTr, tt_surfTr, tt, RPN, ipS);
2317 surfTriangles.push_back(tt_surfTriangles[t]);
2318 delete tt_surfTr;
2319 }
2320 for(int p = 0; p < (int)tt_cp.size(); p++) {
2321 if(tt_cp[p]->ls() != 0) {delete tt_cp[p]; continue;}
2322 mappingCP(tt_cp[p]);
2323 std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(tt_cp[p]);
2324 if(!isIn.second) delete tt_cp[p];
2325 }
2326 delete tt;
2327 return signChange;
2328 }
2329
2330 // cut a tetrahedron into subtetrahedra along the levelset surface
cut(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Tetra * > & subTetras,std::vector<DI_Quad * > & surfQuads,std::vector<DI_Triangle * > & surfTriangles,std::vector<DI_CuttingPoint * > & cp,std::vector<DI_QualError * > & QError)2331 bool DI_Tetra::cut (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
2332 std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
2333 std::vector<DI_Triangle *> &surfTriangles, std::vector<DI_CuttingPoint *> &cp,
2334 std::vector<DI_QualError *> &QError)
2335 {
2336 // check if the tetra is cut by the level set
2337 int on = 0, pos = 0, neg = 0, ze[4];
2338 for (int i = 0; i < 4; i++){
2339 if(pt(i)->isOnBorder()) ze[on++] = i;
2340 else if (pt(i)->ls() > 0.) pos++;
2341 else neg++;
2342 }
2343 if (pos && neg)
2344 selfSplit(e, RPNi, subTetras, surfTriangles, cp, QError);
2345 else{
2346 if(on == 3){ // the level set is zero on a face of the tetra
2347 surfTriangles.push_back(new DI_Triangle(pt(ze[0]), pt(ze[1]), pt(ze[2]), RPNi.back()->getTag()));
2348 // add the triangle twice if the face belongs to 2 tetras => remove it later!
2349 }
2350 if(on == 4)
2351 printf("Warning : tetrahedron with zero levelset on every vertex.\n");
2352 for(int i = 0; i < on; i++)
2353 cp.push_back(new DI_CuttingPoint(pt(ze[i])));
2354 subTetras.push_back(this); // add the tetra to subTetras if it is not cut
2355 }
2356 return (pos && neg);
2357 }
2358
2359 // cut a hex into subtetras along the levelset surface
cut(std::vector<gLevelset * > & RPN,std::vector<DI_IntegrationPoint * > & ip,std::vector<DI_IntegrationPoint * > & ipS,DI_Point::Container & cp,const int polynomialOrderH,const int polynomialOrderT,const int polynomialOrderQ,const int polynomialOrderTr,std::vector<DI_Hexa * > & subHexas,std::vector<DI_Tetra * > & subTetras,std::vector<DI_Quad * > & surfQuads,std::vector<DI_Triangle * > & surfTriangles,std::vector<DI_Line * > & frontLines,int recurLevel,double ** nodeLs) const2360 bool DI_Hexa::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
2361 std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
2362 const int polynomialOrderH, const int polynomialOrderT,
2363 const int polynomialOrderQ, const int polynomialOrderTr,
2364 std::vector<DI_Hexa *> &subHexas, std::vector<DI_Tetra *> &subTetras,
2365 std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
2366 std::vector<DI_Line *> &frontLines, int recurLevel, double **nodeLs) const
2367 {
2368 std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
2369
2370 // reference hexa
2371 DI_Hexa *hh = new DI_Hexa(-1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 8.);
2372 hh->setPolynomialOrder(getPolynomialOrder());
2373 std::vector<DI_Hexa *> hh_subHexas;
2374 std::vector<DI_Tetra *> hh_subTetras;
2375 std::vector<DI_Quad *> hh_surfQuads;
2376 std::vector<DI_Triangle *> hh_surfTriangles;
2377 std::vector<DI_Line *> hh_frontLines;
2378 std::vector<DI_CuttingPoint *> hh_cp;
2379 std::vector<DI_QualError *> QError;
2380
2381 RecurElement *re = new RecurElement(hh);
2382 bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
2383 pushSubElements(re, hh_subHexas);
2384 delete re;
2385
2386 int iPrim = 0;
2387 if(signChange){
2388 for(int l = 0; l < (int)RPN.size(); l++) {
2389 gLevelset *Lsi = RPN[l];
2390 RPNi.push_back(Lsi);
2391 if(Lsi->isPrimitive()) {
2392 hh->addLs(this, Lsi, iPrim++, nodeLs);
2393 int nbSubH = hh_subHexas.size();
2394 int nbSubT = hh_subTetras.size(), nbSubT1 = nbSubT;
2395 int nbSurfQ = hh_surfQuads.size(), nbSurfQ1 = nbSurfQ;
2396 int nbSurfTr = hh_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
2397 int nbFrontLn = hh_frontLines.size();
2398 int nbCP = hh_cp.size();
2399 for(int i = 0; i < nbSubH; i++) hh_subHexas[i]->addLs(hh);
2400 for(int i = 0; i < nbSubT; i++) hh_subTetras[i]->addLs(hh);
2401 for(int i = 0; i < nbSurfQ; i++) hh_surfQuads[i]->addLs(hh);
2402 for(int i = 0; i < nbSurfTr; i++) hh_surfTriangles[i]->addLs(hh);
2403 for(int i = 0; i < nbFrontLn; i++) hh_frontLines[i]->addLs(hh);
2404 for(int i = 0; i < nbCP; i++) hh_cp[i]->addLs(hh);
2405 int pos = 0, neg = 0, ze[8], cz = 0;
2406 for (int i = 0; i < 8; i++){
2407 if(hh->pt(i)->isOnBorder()) ze[cz++] = i;
2408 else if (hh->ls(i) > 0.) pos++;
2409 else neg++;
2410 }
2411
2412 for(int h = 0; h < nbSubH; h++) {
2413 nbSurfQ1 = hh_surfQuads.size();
2414 DI_Hexa *ht = hh_subHexas[0];
2415 hh_subHexas.erase(hh_subHexas.begin());
2416 bool c = ht->cut(hh, RPNi, hh_subHexas, hh_subTetras, hh_surfQuads, hh_surfTriangles, hh_cp, QError);
2417 if(c) delete ht;
2418 if((int)hh_surfQuads.size() > nbSurfQ1){ // 1 quad created on surface of the subHexa
2419 if(isLastQInV(hh_surfQuads)) {
2420 delete hh_surfQuads.back(); hh_surfQuads.pop_back();
2421 }
2422 else if (cz == 4) {
2423 if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3])))
2424 swap(ze[2], ze[3]);
2425 DI_Quad qf(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]));
2426 for(int i = (int)surfQuads.size() - 1; i >= 0; i--){
2427 if (qf.contain(surfQuads[i])){
2428 delete hh_surfQuads.back(); hh_surfQuads.pop_back(); break;
2429 }
2430 }
2431 for(int i = (int)surfTriangles.size() - 1; i >= 0; i--){
2432 if (qf.contain(surfTriangles[i])) {
2433 delete hh_surfQuads.back(); hh_surfQuads.pop_back(); break;
2434 }
2435 }
2436 }
2437 }
2438 }
2439 for(int t = 0; t < nbSubT; t++) {
2440 nbSurfTr1 = hh_surfTriangles.size();
2441 nbSubT1 = hh_subTetras.size();
2442 DI_Tetra *tet = hh_subTetras[0];
2443 hh_subTetras.erase(hh_subTetras.begin());
2444 bool c = tet->cut(hh, RPNi, hh_subTetras, hh_surfQuads, hh_surfTriangles, hh_cp, QError);
2445 if(c) delete tet;
2446 if((int)hh_surfTriangles.size() - nbSurfTr1 == 1 && (int)hh_subTetras.size() == nbSubT1) { // 1 triangle created on surface of the subtetra
2447 // check among hh_surfTriangles
2448 if(isLastTrInV(hh_surfTriangles)) {
2449 delete hh_surfTriangles.back(); hh_surfTriangles.pop_back();
2450 }
2451 else if (cz == 4 && !(pos && neg)){ // 1 triangle created on surface of the hex => check among surfTriangles
2452 if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3])))
2453 swap(ze[2], ze[3]);
2454 DI_Quad qt(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]));
2455 for(int i = (int)surfTriangles.size() - 1; i >= 0; i--){
2456 if(qt.contain(surfTriangles[i])) {
2457 delete hh_surfTriangles.back(); hh_surfTriangles.pop_back(); break;
2458 }
2459 }
2460 }
2461 else { // 1 triangle created inside the hex and not in hh_surfTriangles
2462 // check among the quads with quality error
2463 for(int i = 0; i < (int)QError.size(); i++){
2464 if(isInQE(hh_surfTriangles[nbSurfTr1], QError[i])) {
2465 delete hh_surfTriangles.back(); hh_surfTriangles.pop_back();
2466 break;
2467 }
2468 }
2469 }
2470 }
2471 }
2472 for(int q = 0; q < nbSurfQ; q++) {
2473 DI_Quad *qt = hh_surfQuads[0];
2474 hh_surfQuads.erase(hh_surfQuads.begin());
2475 bool c = qt->cut(hh, RPNi, hh_surfQuads, hh_surfTriangles, hh_frontLines, hh_cp);
2476 if(c) delete qt;
2477 }
2478 for(int t = 0; t < nbSurfTr; t++){
2479 DI_Triangle *trt = hh_surfTriangles[0];
2480 hh_surfTriangles.erase(hh_surfTriangles.begin());
2481 bool c = trt->cut(hh, RPNi, hh_surfQuads, hh_surfTriangles, hh_frontLines, hh_cp);
2482 if(c) delete trt;
2483 }
2484 for(int l = 0; l < nbFrontLn; l++){
2485 DI_Line *lnt = hh_frontLines[0];
2486 hh_frontLines.erase(hh_frontLines.begin());
2487 bool c = lnt->cut(hh, RPNi, hh_frontLines, hh_cp);
2488 if(c) delete lnt;
2489 }
2490 }
2491 else {
2492 for(int h = 0; h < (int)hh_subHexas.size(); h++)
2493 hh_subHexas[h]->chooseLs(Lsi);
2494 for(int t = 0; t < (int)hh_subTetras.size(); t++)
2495 hh_subTetras[t]->chooseLs(Lsi);
2496 for(int q = 0; q < (int)hh_surfQuads.size(); q++)
2497 hh_surfQuads[q]->chooseLs(Lsi);
2498 for(int t = 0; t < (int)hh_surfTriangles.size(); t++)
2499 hh_surfTriangles[t]->chooseLs(Lsi);
2500 for(int l = 0; l < (int)hh_frontLines.size(); l++)
2501 hh_frontLines[l]->chooseLs(Lsi);
2502 for(int p = 0; p < (int)hh_cp.size(); p++)
2503 hh_cp[p]->chooseLs(Lsi);
2504 }
2505 }
2506 }
2507 else {
2508 for(int l = 0; l < (int)RPN.size(); l++) {
2509 gLevelset *Lsi = RPN[l];
2510 if(Lsi->isPrimitive()) {
2511 hh->addLs(this, Lsi, iPrim++, nodeLs);
2512 }
2513 }
2514 if(iPrim == 1) iPrim--;
2515 hh_subHexas[0]->addLs(this, RPN.back(), iPrim, nodeLs);
2516 }
2517
2518 for(int i = 0; i < (int)QError.size(); i++) {
2519 QError[i]->print(this);
2520 delete QError[i];
2521 }
2522
2523 for(int h = 0; h < (int)hh_subHexas.size(); h++) {
2524 hh_subHexas[h]->computeLsTagDom(hh, RPN);
2525 DI_Hexa *hh_subH = new DI_Hexa(*hh_subHexas[h]);
2526 mappingEl(hh_subHexas[h]);
2527 hh_subHexas[h]->integrationPoints(polynomialOrderH, hh_subH, hh, RPN, ip);
2528 subHexas.push_back(hh_subHexas[h]);
2529 delete hh_subH;
2530 }
2531 for(int t = 0; t < (int)hh_subTetras.size(); t++) {
2532 hh_subTetras[t]->computeLsTagDom(hh, RPN);
2533 DI_Tetra *hh_subT = new DI_Tetra(*hh_subTetras[t]);
2534 mappingEl(hh_subTetras[t]);
2535 hh_subTetras[t]->integrationPoints(polynomialOrderT, hh_subT, hh, RPN, ip);
2536 subTetras.push_back(hh_subTetras[t]);
2537 delete hh_subT;
2538 }
2539 for(int q = 0; q < (int)hh_surfQuads.size(); q++) {
2540 hh_surfQuads[q]->computeLsTagBound(hh_subHexas, hh_subTetras);
2541 if(hh_surfQuads[q]->lsTag() == -1) {delete hh_surfQuads[q]; continue;}
2542 DI_Quad *hh_surfQ = new DI_Quad(*hh_surfQuads[q]);
2543 mappingEl(hh_surfQuads[q]);
2544 hh_surfQuads[q]->integrationPoints(polynomialOrderQ, hh_surfQ, hh, RPN, ipS);
2545 surfQuads.push_back(hh_surfQuads[q]);
2546 delete hh_surfQ;
2547 }
2548 for(int t = 0; t < (int)hh_surfTriangles.size(); t++) {
2549 hh_surfTriangles[t]->computeLsTagBound(hh_subHexas, hh_subTetras);
2550 if(hh_surfTriangles[t]->lsTag() == -1) {delete hh_surfTriangles[t]; continue;}
2551 DI_Triangle *hh_surfTr = new DI_Triangle(*hh_surfTriangles[t]);
2552 mappingEl(hh_surfTriangles[t]);
2553 hh_surfTriangles[t]->integrationPoints(polynomialOrderTr, hh_surfTr, hh, RPN, ipS);
2554 surfTriangles.push_back(hh_surfTriangles[t]);
2555 delete hh_surfTr;
2556 }
2557 for(int l = 0; l < (int)hh_frontLines.size(); l++) {
2558 hh_frontLines[l]->computeLsTagBound(hh_surfQuads, hh_surfTriangles);
2559 if(hh_frontLines[l]->lsTag() == -1) {delete hh_frontLines[l]; continue;}
2560 DI_Line *hh_frontLn = new DI_Line(*hh_frontLines[l]);
2561 mappingEl(hh_frontLines[l]);
2562 //hh_frontLines[l]->integrationPoints(polynomialOrderL, tt_surfLn, tt, RPN, ipS);
2563 frontLines.push_back(hh_frontLines[l]);
2564 delete hh_frontLn;
2565 }
2566 for(int p = 0; p < (int)hh_cp.size(); p++) {
2567 if(hh_cp[p]->ls() != 0) {delete hh_cp[p]; continue;} // returns only the cutting points with ls==0
2568 mappingCP(hh_cp[p]);
2569 std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(hh_cp[p]);
2570 if(!isIn.second) delete hh_cp[p];
2571 }
2572 delete hh;
2573 return signChange;
2574 }
2575
2576 // cut a hex into subtetras along the levelset surface
cut(const DI_Element * e,const std::vector<gLevelset * > & RPNi,std::vector<DI_Hexa * > & Hexas,std::vector<DI_Tetra * > & subTetras,std::vector<DI_Quad * > & surfQuads,std::vector<DI_Triangle * > & surfTriangles,std::vector<DI_CuttingPoint * > & cp,std::vector<DI_QualError * > & QError)2577 bool DI_Hexa::cut (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
2578 std::vector<DI_Hexa *> &Hexas, std::vector<DI_Tetra *> &subTetras,
2579 std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
2580 std::vector<DI_CuttingPoint *> &cp, std::vector<DI_QualError *> &QError)
2581 {
2582 // check if the hex is cut by the level set
2583 int on = 0, pos = 0, neg = 0, ze[8];
2584 for (int i = 0; i < 8; i++){
2585 if(pt(i)->isOnBorder()) ze[on++] = i;
2586 else if (pt(i)->ls() > 0.) pos++;
2587 else neg++;
2588 }
2589 if (pos && neg) {
2590 std::vector<DI_Tetra *> tetras; tetras.reserve(6);
2591 // Split the quad into sub tetras
2592 splitIntoTetras (tetras);
2593 // each of the sub tetras is split again into sub tetras with respect to the levelset
2594 int nt0 = surfTriangles.size(), nT1, nt1, nT2, nt2;
2595 for (int t = 0; t < (int)tetras.size(); t++) {
2596 nT1 = subTetras.size(); nt1 = surfTriangles.size();
2597 tetras[t]->selfSplit(e, RPNi, subTetras, surfTriangles, cp, QError);
2598 nT2 = subTetras.size(); nt2 = surfTriangles.size();
2599 if((nT2 - nT1) == 1 && (nt2 - nt1) == 1){ // only 1 triangle created on a face of a tetra => check if it was not yet created
2600 if(isLastTrInV(surfTriangles, nt0)) {delete surfTriangles.back(); surfTriangles.pop_back(); nt2--;}
2601 }
2602 }
2603 }
2604
2605 else{
2606 if(on == 4){ // the level set is zero on a face of the hex
2607 // assert the nodes are in the same plane
2608 if(!isPlanar(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) {
2609 printf("Error : The 4 nodes with zero levelset are not planar!\n");}
2610 else {
2611 // order the 4 nodes
2612 if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) swap(ze[2], ze[3]);
2613 // add the quad twice if the face belongs to 2 hexas => remove it later!
2614 if(ze[0] == 2) surfQuads.push_back(new DI_Quad(pt(ze[1]), pt(ze[2]), pt(ze[3]), pt(ze[0]), RPNi.back()->getTag()));
2615 // to assert that the quad will be split into triangles in the same manner as the hexa into tetras.
2616 else surfQuads.push_back(new DI_Quad(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]), RPNi.back()->getTag()));
2617 }
2618 }
2619 for(int i = 0; i < on; i++)
2620 cp.push_back(new DI_CuttingPoint(pt(ze[i])));
2621 Hexas.push_back(this);
2622 }
2623 return (pos && neg);
2624 }
2625
2626 #endif
2627