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

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

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