source: issm/trunk/src/c/bamg/ListofIntersectionTriangles.cpp@ 18301

Last change on this file since 18301 was 18301, checked in by Mathieu Morlighem, 11 years ago

merged trunk-jpl and trunk for revision 18299
File size: 10.9 KB
Line 
1#include <cstdio>
2#include <cstring>
3#include <cmath>
4#include <ctime>
5
6#include "./bamgobjects.h"
7#include "../shared/shared.h"
8#include "./det.h"
9
10namespace bamg {
11
12 /*Constructors Destructors*/
13 ListofIntersectionTriangles::ListofIntersectionTriangles(int n,int m)/*{{{*/
14 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
15 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
16 }
17 /*}}}*/
18 ListofIntersectionTriangles::~ListofIntersectionTriangles(){/*{{{*/
19 if (lIntTria) delete [] lIntTria,lIntTria=0;
20 if (lSegsI) delete [] lSegsI,lSegsI=0;
21 }
22 /*}}}*/
23
24 /*Methods*/
25 void ListofIntersectionTriangles::Init(void){/*{{{*/
26 state=0;
27 len=0;
28 Size=0;
29 }
30 /*}}}*/
31 double ListofIntersectionTriangles::Length(){/*{{{*/
32 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Length)*/
33
34 // computation of the length
35
36 // check Size
37 if (Size<=0){
38 _error_("Size<=0");
39 }
40
41 Metric Mx,My;
42 int ii,jj;
43 R2 x,y,xy;
44
45 SegInterpolation* SegI=lSegsI;
46 lSegsI[NbSeg].last=Size;
47 int EndSeg=Size;
48
49 y = lIntTria[0].x;
50 double sxy, s = 0;
51 lIntTria[0].s =0;
52 SegI->lBegin=s;
53
54 for (jj=0,ii=1;ii<Size;jj=ii++) {
55 // seg jj,ii
56 x = y;
57 y = lIntTria[ii].x;
58 xy = y-x;
59 Mx = lIntTria[ii].m;
60 My = lIntTria[jj].m;
61 sxy= LengthInterpole(Mx,My,xy);
62 s += sxy;
63 lIntTria[ii].s = s;
64 if (ii == EndSeg){
65 SegI->lEnd=s;
66 SegI++;
67 EndSeg=SegI->last;
68 SegI->lBegin=s;
69 }
70 }
71 len = s;
72 SegI->lEnd=s;
73
74 return s;
75 }
76 /*}}}*/
77 int ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) { /*{{{*/
78 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
79
80 int n;
81 R2 x(0,0);
82 if ( d0) x = (*tt)[0].r * d0;
83 if ( d1) x = x + (*tt)[1].r * d1;
84 if ( d2) x = x + (*tt)[2].r * d2;
85 // newer add same point
86 if(!Size || Norme2_2(lIntTria[Size-1].x-x)) {
87 if (Size==MaxSize) ReShape();
88 lIntTria[Size].t=tt;
89 lIntTria[Size].bary[0]=d0;
90 lIntTria[Size].bary[1]=d1;
91 lIntTria[Size].bary[2]=d2;
92 lIntTria[Size].x = x;
93 Metric m0,m1,m2;
94 BamgVertex * v;
95 if ((v=(*tt)(0))) m0 = v->m;
96 if ((v=(*tt)(1))) m1 = v->m;
97 if ((v=(*tt)(2))) m2 = v->m;
98 lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2);
99 n=Size++;}
100 else n=Size-1;
101 return n;
102 }
103 /*}}}*/
104 int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {/*{{{*/
105 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
106
107 int n;
108 if(!Size || Norme2_2(lIntTria[Size-1].x-A)) {
109 if (Size==MaxSize) ReShape();
110 lIntTria[Size].t=0;
111 lIntTria[Size].x=A;
112 lIntTria[Size].m=mm;
113 n=Size++;
114 }
115 else n=Size-1;
116 return n;
117 }
118 /*}}}*/
119 long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){/*{{{*/
120 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
121
122 //If length<1.5, do nothing
123 double s=Length();
124 if (s<1.5) return 0;
125
126 const long nbvold=nbv;
127 int ii = 1 ;
128 R2 y,x;
129 Metric My,Mx ;
130 double sx =0,sy;
131 int nbi=Max(2,(int) (s+0.5));
132 double sint=s/nbi;
133 double si =sint;
134
135 int EndSeg=Size;
136 SegInterpolation* SegI=NULL;
137 if (NbSeg) SegI=lSegsI,EndSeg=SegI->last;
138
139 for (int k=1;k<nbi;k++){
140 while ((ii < Size) && ( lIntTria[ii].s <= si )){
141 if (ii++ == EndSeg){
142 SegI++;
143 EndSeg=SegI->last;
144 }
145 }
146
147 int ii1=ii-1;
148 x =lIntTria[ii1].x;
149 sx =lIntTria[ii1].s;
150 Metric Mx=lIntTria[ii1].m;
151 y =lIntTria[ii].x;
152 sy =lIntTria[ii].s;
153 Metric My=lIntTria[ii].m;
154 double lxy = sy-sx;
155 double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
156
157 R2 C;
158 double cx = 1-cy;
159 C = SegI ? SegI->F(si): x * cx + y *cy;
160 //C.Echo();
161 //x.Echo();
162 //y.Echo();
163 //_printf_("cx = " << cx << ", cy=" << cy << "\n");
164
165 si += sint;
166 if ( nbv<maxnbv) {
167 vertices[nbv].r = C;
168 vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
169 }
170 else return nbv-nbvold;
171 }
172 return nbv-nbvold;
173 }
174 /*}}}*/
175 void ListofIntersectionTriangles::ReShape(){ /*{{{*/
176
177 int newsize = MaxSize*2;
178 IntersectionTriangles* nw = new IntersectionTriangles[newsize];
179 _assert_(nw);
180
181 // recopy
182 for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];
183 long int verbosity=0;
184 if(verbosity>3) _printf_(" ListofIntersectionTriangles ReShape Maxsize " << MaxSize << " -> " << MaxNbSeg << "\n");
185 MaxSize = newsize;
186 delete [] lIntTria;// remove old
187 lIntTria = nw; // copy pointer
188 }
189 /*}}}*/
190 void ListofIntersectionTriangles::SplitEdge(const Mesh & Bh, const R2 &A,const R2 &B,int nbegin) {/*{{{*/
191 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ListofIntersectionTriangles)*/
192
193 Triangle *tbegin, *t;
194
195 Icoor2 deta[3], deti,detj;
196 double ba[3];
197 int ifirst=-1,ilast;
198 int i0,i1,i2;
199 int ocut,i,j,k=-1;
200 // int OnAVertices =0;
201 Icoor2 dt[3];
202 I2 a = Bh.R2ToI2(A) ,b= Bh.R2ToI2(B);// compute the Icoor a,b
203 I2 vi,vj;
204 int iedge =-1;// not a edge
205
206 if(nbegin) {// optimisation
207 // we suppose knowing the starting triangle
208 t=tbegin=lIntTria[ilast=(Size-1)].t;
209 if (tbegin->det>=0)
210 ifirst = ilast;}
211 else {// not optimisation
212 Init();
213 t=tbegin = Bh.TriangleFindFromCoord(a,deta);
214 if( t->det>=0)
215 ilast=NewItem(t,double(deta[0])/t->det,double(deta[1])/t->det,double(deta[2])/t->det);
216 else
217 {// find the nearest boundary edge of the vertex A
218 // find a edge or such normal projection a the edge IJ is on the edge
219 // <=> IJ.IA >=0 && IJ.AJ >=0
220 ilast=ifirst;
221 double ba,bb;
222 AdjacentTriangle edge=CloseBoundaryEdge(a,t,ba,bb);
223 BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
224 NewItem(A,Metric(ba,v0,bb,v1));
225 t=edge;
226 // test if the point b is in the same side
227 if (det(v0.i,v1.i,b)>=0) {
228 AdjacentTriangle edge=CloseBoundaryEdge(a,t,ba,bb);
229 BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
230 NewItem(A,Metric(ba,v0,bb,v1));
231 return;
232 }
233 } // find the nearest boundary edge of the vertex A
234 } // end not optimisation
235 if (t->det<0) { // outside departure
236 while (t->det <0) { // intersection boundary edge and a,b,
237 k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
238 if (k<0){
239 _error_("k<0");
240 }
241 ocut = OppositeEdge[k];
242 i=VerticesOfTriangularEdge[ocut][0];
243 j=VerticesOfTriangularEdge[ocut][1];
244 vi=(*t)[i];
245 vj=(*t)[j];
246 deti = bamg::det(a,b,vi);
247 detj = bamg::det(a,b,vj);
248 if (deti>0) // go to i direction on gamma
249 ocut = PreviousEdge[ocut];
250 else if (detj<=0) // go to j direction on gamma
251 ocut = NextEdge[ocut];
252 AdjacentTriangle tadj =t->Adj(ocut);
253 t = tadj;
254 iedge= tadj;
255 if (t == tbegin) { //
256 double ba,bb;
257 AdjacentTriangle edge=CloseBoundaryEdge(a,t,ba,bb);
258 BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
259 NewItem(A,Metric(ba,v0,bb,v1));
260 return;
261 }
262 } // end while (t->det <0)
263 // theoriticaly we have: deti =<0 and detj>0
264
265 // computation of barycentric coor
266 // test if the point b is on size on t
267 // we revert vi,vj because vi,vj is def in Adj triangle
268 if ( det(vi,vj,b)>=0) {
269 t=tbegin;
270 double ba,bb;
271 AdjacentTriangle edge=CloseBoundaryEdge(b,t,ba,bb);
272 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
273 return;
274 }
275 else
276 {
277 k = OppositeVertex[iedge];
278 i=VerticesOfTriangularEdge[iedge][0];
279 j=VerticesOfTriangularEdge[iedge][1];
280 double dij = detj-deti;
281 if (i+j+k != 0 + 1 +2){
282 _error_("i+j+k != 0 + 1 +2");
283 }
284 ba[j] = detj/dij;
285 ba[i] = -deti/dij;
286 ba[k] = 0;
287 ilast=NewItem(t,ba[0],ba[1],ba[2]); }
288 } // outside departure
289
290 // recherche the intersection of [a,b] with Bh Mesh.
291 // we know a triangle ta contening the vertex a
292 // we have 2 case for intersection [a,b] with a edge [A,B] of Bh
293 // 1) the intersection point is in ]A,B[
294 // 2) is A or B
295 // first version ---
296 for (;;) {
297 // t->Draw();
298 if (iedge < 0) {
299 i0 =0;i1=1;i2=2;
300 dt[0] =bamg::det(a,b,(*t)[0]);
301 dt[1] =bamg::det(a,b,(*t)[1]);
302 dt[2] =bamg::det(a,b,(*t)[2]);}
303 else {
304 i2 = iedge;
305 i0 = NextEdge[i2];
306 i1 = NextEdge[i0];
307 dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because
308 dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side
309 dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);}
310
311 // so we have just to see the transition from - to + of the det0..2 on edge of t
312 // because we are going from a to b
313 if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
314 ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))
315 ocut =i0;
316 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
317 (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
318 ocut =i1;
319 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
320 (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
321 ocut =i2;
322 else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) &&
323 ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))
324 ocut =i0;
325 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) &&
326 (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
327 ocut =i1;
328 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) &&
329 (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
330 ocut =i2;
331 else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
332 ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0))
333 ocut =i0;
334 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
335 (dt[j=VerticesOfTriangularEdge[i1][1]] == 0))
336 ocut =i1;
337 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
338 (dt[j=VerticesOfTriangularEdge[i2][1]] == 0))
339 ocut =i2;
340 else { // On a edge (2 zero)
341 k =0;
342 if (dt[0]) ocut=0,k++;
343 if (dt[1]) ocut=1,k++;
344 if (dt[2]) ocut=2,k++;
345 if(k == 1) {
346 if (dt[ocut] >0) // triangle upper AB
347 ocut = NextEdge[ocut];
348 i= VerticesOfTriangularEdge[ocut][0];
349 j= VerticesOfTriangularEdge[ocut][1];
350 }
351 else {
352 _error_("Bug Split Edge");
353 }
354 }
355
356 k = OppositeVertex[ocut];
357
358 Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
359
360 if (detbij >= 0) { //we find the triangle contening b
361 dt[0]=bamg::det((*t)[1],(*t)[2],b);
362 dt[1]=bamg::det((*t)[2],(*t)[0],b);
363 dt[2]=bamg::det((*t)[0],(*t)[1],b);
364 double dd = t->det;
365 NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);
366 return ;}
367 else { // next triangle by adjacent by edge ocut
368 deti = dt[i];
369 detj = dt[j];
370 double dij = detj-deti;
371 ba[i] = detj/dij;
372 ba[j] = -deti/dij;
373 ba[3-i-j ] = 0;
374 ilast=NewItem(t, ba[0],ba[1],ba[2]);
375
376 AdjacentTriangle ta =t->Adj(ocut);
377 t = ta;
378 iedge= ta;
379 if (t->det <= 0) {
380 double ba,bb;
381 AdjacentTriangle edge=CloseBoundaryEdge(b,t,ba,bb);
382 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
383 return;
384 }
385 }// we go outside of omega
386 } // for(;;)
387 }
388 /*}}}*/
389
390}
Note: See TracBrowser for help on using the repository browser.