![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2811
- Timestamp:
- 01/12/10 16:50:43 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 2 edited
-
Mesh2.cpp (modified) (5 diffs)
-
Triangles.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.cpp
r2810 r2811 49 49 int hinterpole=1; 50 50 51 52 long NbUnSwap =0;53 51 int ForDebugging = 0; 54 52 const Direction NoDirOfSearch = Direction(); … … 428 426 return er; 429 427 } 430 431 432 433 Metric Triangles::MetricAt (const R2 & A) const434 { //if ((vertices <= &v) && (vertices < v+nbv)) return v.m;435 I2 a = toI2(A);436 Icoor2 deta[3];437 Triangle * t =FindTriangleContening(a,deta);438 if (t->det <0) { // outside439 double ba,bb;440 TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;441 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}442 else { // inside443 Real8 aa[3];444 Real8 s = deta[0]+deta[1]+deta[2];445 aa[0]=deta[0]/s;446 aa[1]=deta[1]/s;447 aa[2]=deta[2]/s;448 return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);449 }450 }451 428 452 429 … … 988 965 989 966 990 int Triangle::swap(Int2 a,int koption){ 991 if(a/4 !=0) return 0;// arete lock or MarkUnSwap 992 993 register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent 994 register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles 995 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap 996 997 register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; 998 register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; 999 register Vertex *s1=t1->ns[OppositeVertex[a1]]; 1000 register Vertex *s2=t2->ns[OppositeVertex[a2]]; 1001 1002 Icoor2 det1=t1->det , det2=t2->det ; 1003 Icoor2 detT = det1+det2; 1004 Icoor2 detA = Abs(det1) + Abs(det2); 1005 Icoor2 detMin = Min(det1,det2); 1006 1007 int OnSwap = 0; 1008 // si 2 triangle infini (bord) => detT = -2; 1009 if (sa == 0) {// les deux triangles sont frontieres 1010 det2=bamg::det(s2->i,sb->i,s1->i); 1011 OnSwap = det2 >0;} 1012 else if (sb == 0) { // les deux triangles sont frontieres 1013 det1=bamg::det(s1->i,sa->i,s2->i); 1014 OnSwap = det1 >0;} 1015 else if(( s1 != 0) && (s2 != 0) ) { 1016 det1 = bamg::det(s1->i,sa->i,s2->i); 1017 det2 = detT - det1; 1018 OnSwap = (Abs(det1) + Abs(det2)) < detA; 1019 1020 Icoor2 detMinNew=Min(det1,det2); 1021 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test 1022 if (! OnSwap &&(detMinNew>0)) { 1023 OnSwap = detMin ==0; 1024 if (! OnSwap) { 1025 int kopt = koption; 1026 while (1) 1027 if(kopt) { 1028 // critere de Delaunay pure isotrope 1029 register Icoor2 xb1 = sb->i.x - s1->i.x, 1030 x21 = s2->i.x - s1->i.x, 1031 yb1 = sb->i.y - s1->i.y, 1032 y21 = s2->i.y - s1->i.y, 1033 xba = sb->i.x - sa->i.x, 1034 x2a = s2->i.x - sa->i.x, 1035 yba = sb->i.y - sa->i.y, 1036 y2a = s2->i.y - sa->i.y; 1037 register double 1038 cosb12 = double(xb1*x21 + yb1*y21), 1039 cosba2 = double(xba*x2a + yba*y2a) , 1040 sinb12 = double(det2), 1041 sinba2 = double(t2->det); 1042 1043 1044 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) 1045 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); 1046 // if(CurrentTh) 1047 // cout << "swap " << CurrentTh->Number(sa) << " " << CurrentTh->Number(sb) << " " ; 1048 // cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12 1049 // << " Onswap = " << OnSwap << endl; 1050 break; 1051 } 1052 else 1053 { 1054 // critere de Delaunay anisotrope 1055 Real8 som; 1056 I2 AB=(I2) *sb - (I2) *sa; 1057 I2 MAB2=((I2) *sb + (I2) *sa); 1058 R2 MAB(MAB2.x*0.5,MAB2.y*0.5); 1059 I2 A1=(I2) *s1 - (I2) *sa; 1060 I2 D = (I2) * s1 - (I2) * sb ; 1061 R2 S2(s2->i.x,s2->i.y); 1062 R2 S1(s1->i.x,s1->i.y); 1063 { 1064 Metric M=s1->m; 1065 R2 ABo = M.Orthogonal(AB); 1066 R2 A1o = M.Orthogonal(A1); 1067 // (A+B)+ x ABo = (S1+B)/2+ y A1 1068 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 1069 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); 1070 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 1071 if (Abs(d) > dd*1.e-3) { 1072 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 1073 som = M(C - S2)/M(C - S1); 1074 } else 1075 {kopt=1;continue;} 1076 1077 } 1078 { 1079 Metric M=s2->m; 1080 R2 ABo = M.Orthogonal(AB); 1081 R2 A1o = M.Orthogonal(A1); 1082 // (A+B)+ x ABo = (S1+B)/2+ y A1 1083 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 1084 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); 1085 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 1086 if(Abs(d) > dd*1.e-3) { 1087 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 1088 som += M(C - S2)/M(C - S1); 1089 } else 1090 {kopt=1;continue;} 1091 } 1092 OnSwap = som < 2; 1093 break; 1094 } 1095 1096 } // OnSwap 1097 } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) 1098 } 1099 if( OnSwap ) 1100 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); 1101 else { 1102 NbUnSwap ++; 1103 t1->SetMarkUnSwap(a1); 1104 } 1105 return OnSwap; 1106 } 967 1107 968 1108 969 Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega) … … 1213 1074 } 1214 1075 1215 1216 1217 1218 void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)1219 {1220 // -------------------------------------------1221 // s21222 // !1223 // /|\ !1224 // / | \ !1225 // / | \ !1226 // tt1 / | \ tt0 !1227 // / |s \ !1228 // / . \ !1229 // / . ` \ !1230 // / . ` \ !1231 // ---------------- !1232 // s0 tt2 s11233 //--------------------------------------------1234 1235 Triangle * tt[3]; // the 3 new Triangles1236 Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2];1237 Icoor2 det3local[3];1238 int infv = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0;1239 // infv = ordre of the infini vertex (null)1240 register int nbd0 =0; // number of zero det31241 register int izerodet=-1,iedge; // izerodet = egde contening the vertex s1242 Icoor2 detOld = t->det;1243 1244 if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) )1245 {1246 cerr << " infv " << infv << " det = " << detOld << endl;1247 cerr << Number(s) << " "<< Number(s0) << " "1248 << Number(s1) << " " << Number(s2) << endl;1249 MeshError(3);1250 }1251 1252 // if det3 do not exist then constuct det31253 if (!det3) {1254 det3 = det3local; // alloc1255 if ( infv<0 ) {1256 det3[0]=bamg::det(s ,s1,s2);1257 det3[1]=bamg::det(s0,s ,s2);1258 det3[2]=bamg::det(s0,s1,s );}1259 else {1260 // one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk1261 det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ;1262 det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ;1263 det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}}1264 1265 1266 if (!det3[0]) izerodet=0,nbd0++;1267 if (!det3[1]) izerodet=1,nbd0++;1268 if (!det3[2]) izerodet=2,nbd0++;1269 1270 if (nbd0 >0 ) // point s on a egde or on a vertex1271 if (nbd0 == 1) {1272 iedge = OppositeEdge[izerodet];1273 TriangleAdjacent ta = t->Adj(iedge);1274 1275 // the point is on the edge1276 // if the point is one the boundary1277 // add the point in outside part1278 if ( t->det >=0) { // inside triangle1279 if ((( Triangle *) ta)->det < 0 ) {1280 // add in outside triangle1281 Add(s,( Triangle *) ta);1282 return;}1283 }}1284 else {1285 cerr << " bug " << nbd0 <<endl;1286 cerr << " Bug double points in " << endl ;1287 cerr << " s = " << Number(s) << " " << s << endl;1288 cerr << " s0 = "<< Number(s0) << " " << s0 << endl;1289 cerr << " s1 = "<< Number(s1) << " " << s1 << endl;1290 cerr << " s2 = "<< Number(s2) << " " << s2 << endl;1291 MeshError(5,this);}1292 1293 // remove de MarkUnSwap edge1294 t->SetUnMarkUnSwap(0);1295 t->SetUnMarkUnSwap(1);1296 t->SetUnMarkUnSwap(2);1297 1298 tt[0]= t;1299 tt[1]= &triangles[nbt++];1300 tt[2]= &triangles[nbt++];1301 1302 if (nbt>nbtx) {1303 cerr << " No enougth triangles " << endl;1304 MeshError(999,this);1305 }1306 1307 *tt[1]= *tt[2]= *t;1308 // gestion of the link1309 tt[0]->link=tt[1];1310 tt[1]->link=tt[2];1311 1312 (* tt[0])(OppositeVertex[0])=&s;1313 (* tt[1])(OppositeVertex[1])=&s;1314 (* tt[2])(OppositeVertex[2])=&s;1315 1316 tt[0]->det=det3[0];1317 tt[1]->det=det3[1];1318 tt[2]->det=det3[2];1319 1320 // update adj des triangles externe1321 tt[0]->SetAdjAdj(0);1322 tt[1]->SetAdjAdj(1);1323 tt[2]->SetAdjAdj(2);1324 // update des adj des 3 triangle interne1325 const int i0 = 0;1326 const int i1= NextEdge[i0];1327 const int i2 = PreviousEdge[i0];1328 1329 tt[i0]->SetAdj2(i2,tt[i2],i0);1330 tt[i1]->SetAdj2(i0,tt[i0],i1);1331 tt[i2]->SetAdj2(i1,tt[i1],i2);1332 1333 tt[0]->SetTriangleContainingTheVertex();1334 tt[1]->SetTriangleContainingTheVertex();1335 tt[2]->SetTriangleContainingTheVertex();1336 1337 1338 // swap if the point s is on a edge1339 if(izerodet>=0) {1340 // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;1341 int rswap =tt[izerodet]->swap(iedge);1342 1343 if (!rswap)1344 {1345 cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;1346 }1347 assert(rswap);1348 }1349 }1350 1351 1352 Int4 Triangles::SplitInternalEdgeWithBorderVertices()1353 {1354 Int4 NbSplitEdge=0;1355 SetVertexFieldOn();1356 Int4 it;1357 Int4 nbvold=nbv;1358 long int verbosity=2;1359 for (it=0;it<nbt;it++)1360 {1361 Triangle &t=triangles[it];1362 if (t.link)1363 for (int j=0;j<3;j++)1364 if(!t.Locked(j) && !t.Hidden(j)){1365 Triangle &tt = *t.TriangleAdj(j);1366 if ( &tt && tt.link && it < Number(tt))1367 { // an internal edge1368 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];1369 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];1370 if (v0.on && v1.on)1371 {1372 R2 P= ((R2) v0 + (R2) v1)*0.5;1373 if ( nbv<nbvx) {1374 vertices[nbv].r = P;1375 vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);1376 vertices[nbv].ReferenceNumber=0;1377 vertices[nbv].DirOfSearch = NoDirOfSearch ;1378 }1379 NbSplitEdge++;1380 if (verbosity>7)1381 cout <<" Internal edge with two vertices on boundary"1382 << Number(v0) << " " << Number(v1) << " by " << endl;1383 }1384 }1385 }1386 }1387 ReMakeTriangleContainingTheVertex();1388 if (nbvold!=nbv)1389 {1390 Int4 iv = nbvold;1391 Int4 NbSwap = 0;1392 Icoor2 dete[3];1393 for (Int4 i=nbvold;i<nbv;i++)1394 {// for all the new point1395 Vertex & vi = vertices[i];1396 vi.i = toI2(vi.r);1397 vi.r = toR2(vi.i);1398 // if (!quadtree->ToClose(vi,seuil,hi,hj)) {1399 // a good new point1400 vi.ReferenceNumber=0;1401 vi.DirOfSearch =NoDirOfSearch;1402 // cout << " Add " << Number(vi) << " " << vi1403 // << " " << Number(vi) << " <--> " << Number(vi) <<endl;1404 Triangle *tcvi = FindTriangleContening(vi.i,dete);1405 if (tcvi && !tcvi->link) {1406 cout << i << " PB insert point " << Number(vi) << vi << Number(vi)1407 << " tcvi = " << tcvi << " " << tcvi->link << endl;1408 cout << (*tcvi)[1] << (*tcvi)[2] << endl;1409 tcvi = FindTriangleContening(vi.i,dete);1410 cout << (*tcvi)[1] << (*tcvi)[2] << endl;1411 MeshError(1001,this);1412 }1413 1414 1415 quadtree->Add(vi);1416 assert (tcvi && tcvi->det >= 0) ;// internal1417 Add(vi,tcvi,dete);1418 NbSwap += vi.Optim(1);1419 iv++;1420 // }1421 }1422 if (verbosity>3)1423 {1424 cout << " Nb Of New Point " << iv ;1425 cout << " Nb swap = " << NbSwap << " to split internal edges with border vertices" ;}1426 1427 nbv = iv;1428 }1429 if (NbSplitEdge > nbv-nbvold)1430 cout << " Warning not enough vertices to split all internal edges " << endl1431 << " we lost " << NbSplitEdge - ( nbv-nbvold) << " Edges Sorry " << endl;1432 if (verbosity>2)1433 cout << "SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << endl;1434 return NbSplitEdge;1435 }1436 Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)1437 {1438 long int verbosity=2;1439 Real8 seuil= 1.414/2 ;// for two close point1440 Int4 i;1441 // insertion part ---1442 1443 const Int4 nbvnew = nbv-nbvold;1444 if (verbosity>5)1445 cout << " Try to Insert the " <<nbvnew<< " new points " << endl;1446 Int4 NbSwap=0;1447 Icoor2 dete[3];1448 1449 // construction d'un ordre aleatoire1450 if (! nbvnew)1451 return 0;1452 if (nbvnew) {1453 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;1454 Int4 k3 = rand()%nbvnew ;1455 for (Int4 is3=0; is3<nbvnew; is3++) {1456 register Int4 j = nbvold +(k3 = (k3 + PrimeNumber)% nbvnew);1457 register Int4 i = nbvold+is3;1458 ordre[i]= vertices + j;1459 ordre[i]->ReferenceNumber=i;1460 }1461 // be carefull1462 Int4 iv = nbvold;1463 for (i=nbvold;i<nbv;i++)1464 {// for all the new point1465 Vertex & vi = *ordre[i];1466 vi.i = toI2(vi.r);1467 vi.r = toR2(vi.i);1468 Real4 hx,hy;1469 vi.m.Box(hx,hy);1470 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);1471 if (!quadtree->ToClose(vi,seuil,hi,hj))1472 {1473 // a good new point1474 Vertex & vj = vertices[iv];1475 Int4 j = vj.ReferenceNumber;1476 assert( &vj== ordre[j]);1477 if(i!=j)1478 { // for valgring1479 Exchange(vi,vj);1480 Exchange(ordre[j],ordre[i]);1481 }1482 vj.ReferenceNumber=0;1483 // cout << " Add " << Number(vj) << " " << vj1484 // << " " << Number(vi) << " <--> " << Number(vj) <<endl;1485 Triangle *tcvj = FindTriangleContening(vj.i,dete);1486 if (tcvj && !tcvj->link)1487 {1488 cerr << i << " PB insert point " << Number(vj) << vj << Number(vi)1489 << " tcvj = " << tcvj << " " << tcvj->link << endl;1490 cerr << (*tcvj)[1] << (*tcvj)[2] << endl;1491 tcvj = FindTriangleContening(vj.i,dete);1492 cout << (*tcvj)[1] << (*tcvj)[2] << endl;1493 MeshError(1001,this);1494 }1495 1496 1497 quadtree->Add(vj);1498 assert (tcvj && tcvj->det >= 0) ;// internal1499 Add(vj,tcvj,dete);1500 NbSwap += vj.Optim(1);1501 iv++;1502 }1503 }1504 if (verbosity>3) {1505 cout << " Nb Of New Point " << iv << " Nb Of To close Points " << nbv-iv ;1506 cout << " Nb swap = " << NbSwap << " after " ;}1507 1508 nbv = iv;1509 }1510 1511 for (i=nbvold;i<nbv;i++)1512 NbSwap += vertices[i].Optim(1);1513 if (verbosity>3)1514 cout << " NbSwap = " << NbSwap << endl;1515 1516 1517 NbTSwap += NbSwap ;1518 return nbv-nbvold;1519 1520 }1521 void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)1522 { // Triangles::NewPoints1523 long int verbosity=2;1524 Int4 nbtold(nbt),nbvold(nbv);1525 if (verbosity>2)1526 cout << " -- Triangles::NewPoints ";1527 if (verbosity>3)cout << " nbv (in) on Boundary = " << nbv <<endl;1528 Int4 i,k;1529 int j;1530 Int4 *first_np_or_next_t = new Int4[nbtx];1531 Int4 NbTSwap =0;1532 // insert old point1533 nbtold = nbt;1534 1535 if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx))1536 {1537 // Bh.SetVertexFieldOn();1538 for (i=0;i<Bh.nbv;i++)1539 {1540 Vertex & bv = Bh[i];1541 if (!bv.on) {1542 vertices[nbv].r = bv.r;1543 vertices[nbv++].m = bv.m;}1544 }1545 int nbv1=nbv;1546 Bh.ReMakeTriangleContainingTheVertex();1547 InsertNewPoints(nbvold,NbTSwap) ;1548 if (verbosity>2)1549 cout << " (Nb of Points from background mesh = "1550 << nbv-nbvold << " / " << nbv1-nbvold << ")" << endl;1551 }1552 else1553 Bh.ReMakeTriangleContainingTheVertex();1554 1555 Triangle *t;1556 // generation of the list of next Triangle1557 // at 1 time we test all the triangles1558 Int4 Headt =0,next_t;1559 for(i=0;i<nbt;i++)1560 first_np_or_next_t[i]=-(i+1);1561 // end list i >= nbt1562 // the list of test triangle is1563 // the next traingle on i is -first_np_or_next_t[i]1564 int iter=0;1565 // Big loop1566 do {1567 iter++;1568 nbtold = nbt;1569 nbvold = nbv;1570 1571 // default size of IntersectionTriangle1572 1573 i=Headt;1574 next_t=-first_np_or_next_t[i];1575 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])1576 { // for each triangle t1577 // we can change first_np_or_next_t[i]1578 // cout << " Do the triangle " << i << " Next_t=" << next_t << endl;1579 assert(i>=0 && i < nbt );1580 first_np_or_next_t[i] = iter;1581 for(j=0;j<3;j++)1582 { // for each edge1583 TriangleAdjacent tj(t,j);1584 Vertex & vA = * tj.EdgeVertex(0);1585 Vertex & vB = * tj.EdgeVertex(1);1586 1587 if (!t->link) continue;// boundary1588 if (t->det <0) continue;1589 if (t->Locked(j)) continue;1590 1591 TriangleAdjacent tadjj = t->Adj(j);1592 Triangle * ta= tadjj;1593 1594 if (ta->det <0) continue;1595 1596 R2 A = vA;1597 R2 B = vB;1598 1599 k=Number(ta);1600 1601 if(first_np_or_next_t[k]==iter) // this edge is done before1602 continue; // next edge of the triangle1603 1604 //const Int4 NbvOld = nbv;1605 lIntTria.SplitEdge(Bh,A,B);1606 lIntTria.NewPoints(vertices,nbv,nbvx);1607 } // end loop for each edge1608 1609 }// for triangle1610 1611 if (!InsertNewPoints(nbvold,NbTSwap))1612 break;1613 1614 for (i=nbtold;i<nbt;i++)1615 first_np_or_next_t[i]=iter;1616 1617 Headt = nbt; // empty list1618 for (i=nbvold;i<nbv;i++)1619 { // for all the triangle contening the vertex i1620 Vertex * s = vertices + i;1621 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);1622 Triangle * tbegin= (Triangle*) ta;1623 Int4 kt;1624 do {1625 kt = Number((Triangle*) ta);1626 if (first_np_or_next_t[kt]>0)1627 first_np_or_next_t[kt]=-Headt,Headt=kt;1628 assert( ta.EdgeVertex(0) == s);1629 ta = Next(Adj(ta));1630 } while ( (tbegin != (Triangle*) ta));1631 }1632 1633 } while (nbv!=nbvold);1634 1635 delete [] first_np_or_next_t;1636 1637 Int4 NbSwapf =0,NbSwp;1638 1639 // bofbof1640 1641 1642 NbSwp = NbSwapf;1643 for (i=0;i<nbv;i++)1644 NbSwapf += vertices[i].Optim(0);1645 /*1646 for (i=0;i<nbv;i++)1647 NbSwapf += vertices[i].Optim(0);1648 for (i=0;i<nbv;i++)1649 NbSwapf += vertices[i].Optim(0);1650 for (i=0;i<nbv;i++)1651 NbSwapf += vertices[i].Optim(0);1652 for (i=0;i<nbv;i++)1653 NbSwapf += vertices[i].Optim(0);1654 */1655 NbTSwap += NbSwapf ;1656 if (verbosity>3) cout << " " ;1657 if (verbosity>2)1658 cout << " Nb Of Vertices =" << nbv << " Nb of triangles = " << nbt-NbOutT1659 << " NbSwap final = " << NbSwapf << " Nb Total Of Swap = " << NbTSwap << endl;1660 1661 1662 }1663 1664 void Triangles::NewPointsOld(Triangles & Bh)1665 { // Triangles::NewPointsOld1666 long int verbosity=2;1667 Real8 seuil= 0.7 ;// for two neart point1668 if (verbosity>1)1669 cout << " begin : Triangles::NewPointsOld " << endl;1670 Int4 i,k;1671 int j;1672 Int4 BeginNewPoint[3];1673 Int4 EndNewPoint[3];1674 #ifdef TRACETRIANGLE1675 Int4 trace=0;1676 #endif1677 int step[3];1678 Int4 *first_np_or_next_t = new Int4[nbtx];1679 Int4 ColorEdge[3];1680 Int4 color=-1;1681 Triangle *t;1682 // generation of the list of next Triangle1683 // at 1 time we test all the triangles1684 Int4 Headt =0,next_t;1685 for(i=0;i<nbt;i++)1686 first_np_or_next_t[i]=-(i+1);1687 // end list i >= nbt1688 // the list of test triangle is1689 // the next Triangle on i is -first_np_or_next_t[i]1690 Int4 nbtold,nbvold;1691 1692 // Big loop1693 do {1694 nbtold = nbt;1695 nbvold = nbv;1696 // default size of IntersectionTriangle1697 1698 i=Headt;1699 next_t=-first_np_or_next_t[i];1700 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])1701 { // for each triangle t1702 // we can change first_np_or_next_t[i]1703 #ifdef TRACETRIANGLE1704 trace = TRACETRIANGLE <0 ? 1 : i == TRACETRIANGLE;1705 #endif1706 // cout << " Do the triangle " << i << " Next_t=" << next_t << endl;1707 assert(i>=0 && i < nbt );1708 first_np_or_next_t[i] = nbv; // to save the fist new point of triangle1709 for(j=0;j<3;j++)1710 { // for each edge1711 TriangleAdjacent tj(t,j);1712 // color++;// the color is 3i+j1713 color = 3*i + j ;;1714 ColorEdge[j]=color;1715 BeginNewPoint[j]=nbv;1716 EndNewPoint[j]=nbv-1;1717 step[j]=1;// right sens1718 Vertex & vA = * tj.EdgeVertex(0);1719 Vertex & vB = * tj.EdgeVertex(1);1720 1721 #ifdef TRACETRIANGLE1722 if(trace) {1723 cout << " i " << Number(vA) <<" j "<< Number(vB)1724 << " " << t->Locked(j) ;1725 }1726 #endif1727 if (!t->link) continue;// boundary1728 if (t->det <0) continue;1729 if (t->Locked(j)) continue;1730 1731 TriangleAdjacent tadjj = t->Adj(j);1732 Triangle * ta= tadjj;1733 1734 if (ta->det <0) continue;1735 1736 R2 A = vA;1737 R2 B = vB;1738 1739 k=Number(ta);1740 // the 2 opposite vertices1741 const Vertex & vC1 = *tj.OppositeVertex();1742 const Vertex & vC2 = *tadjj.OppositeVertex();1743 1744 #ifdef TRACETRIANGLE1745 trace = trace || k == TRACETRIANGLE;1746 if(trace) {1747 cout << "Test Arete " << i << " AB = " << A << B1748 << "i " <<Number(vA)<< "j" <<Number(vB);1749 cout << " link" <<(int)t->link << " ta=" << Number( ta)1750 << " det " <<ta->det ;1751 cout << " hA = " <<vA.m.h << " hB = " <<vB.m.h ;1752 cout << " loc " << ta->Locked(j) << endl;1753 }1754 #endif1755 1756 if(first_np_or_next_t[k]>0) { // this edge is done before1757 // find the color of the edge and begin , end of newpoint1758 register int kk = t->NuEdgeTriangleAdj(j);1759 assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1]));1760 assert ((*t)(VerticesOfTriangularEdge[j][1]) == (*ta)(VerticesOfTriangularEdge[kk][0]));1761 register Int4 kolor =3*k + kk;1762 ColorEdge[j]=kolor;1763 register Int4 kkk= 1;1764 step[j]=-1;// other sens1765 BeginNewPoint[j]=0;1766 EndNewPoint[j]=-1; // empty list1767 for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++)1768 if (vertices[iv].color > kolor)1769 break; // the color is passed1770 else if (vertices[iv].color == kolor) {1771 EndNewPoint[j]=iv;1772 if (kkk) // one time test1773 kkk=0,BeginNewPoint[j]=iv;}1774 continue; // next edge of the triangle1775 } // end if( k < i)1776 1777 1778 const Int4 NbvOld = nbv;1779 lIntTria.SplitEdge(Bh,A,B);1780 // Int4 NbvNp =1781 lIntTria.NewPoints(vertices,nbv,nbvx);1782 Int4 nbvNew=nbv;1783 nbv = NbvOld;1784 for (Int4 iv=NbvOld;iv<nbvNew;iv++) {1785 vertices[nbv].color = color;1786 vertices[nbv].ReferenceNumber=nbv;// circular link1787 R2 C = vertices[iv].r;1788 vertices[nbv].r = C;1789 vertices[nbv].m = vertices[iv].m ;1790 // test if the new point is not to close to the 2 opposite vertex1791 R2 CC1 = C-vC1 , CC2 = C-vC2;1792 if ( ( (vC1.m(CC1) + vertices[nbv].m(CC1)) > seuil)1793 && ( (vC2.m(CC2) + vertices[nbv].m(CC2)) > seuil) )1794 nbv++;1795 }1796 1797 EndNewPoint[j] = nbv-1;1798 } // end loop for each edge1799 1800 #ifdef TRACETRIANGLE1801 if(trace) {1802 // verification des point cree1803 cout << "\n ------------ " << t->link << " " << t->det1804 << " b " << BeginNewPoint[0] << " " << BeginNewPoint[1]1805 << " " << BeginNewPoint[2] << " "1806 << " e " << EndNewPoint[0] << " " << EndNewPoint[1]1807 << " " << EndNewPoint[2] << " "1808 << " s " << step[0] << " " << step[1] << " " << step[2] << " "1809 << endl;1810 }1811 #endif1812 1813 if (!t->link) continue;// boundary1814 if (t->det<=0) continue;// outside1815 // continue;1816 for(int j0=0;j0<3;j0++)1817 for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++)1818 {1819 // find the neart point one the opposite edge1820 // to compute i11821 Vertex & vi0 = vertices[i0];1822 int kstack = 0;1823 Int4 stack[10];1824 // Int4 savRef[10];1825 int j1 = j0;1826 while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge1827 // computation of the intersection of edge j1 and DOrto1828 // take the good sens1829 1830 if (BeginNewPoint[j1]> EndNewPoint[j1])1831 continue; //1832 else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) {1833 for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++)1834 stack[kstack++] = ii1;1835 continue;}1836 1837 1838 int k0,k1;1839 if (step[j1]<0) k0=1,k1=0; // reverse1840 else k0=0,k1=1;1841 R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]];1842 R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]];1843 R2 D = V11-V10;1844 Real8 c0 = vi0.m(D,(R2) vi0);1845 1846 Real8 c10 = vi0.m(D,V10);1847 Real8 c11 = vi0.m(D,V11);1848 1849 Real8 s;1850 //cout << " --i0 = " << i0 << D << V10 << V11 << endl ;1851 //cout << " c10 " << c10 << " c0 " << c0 << " c11 " << c11 << endl;1852 if (( c10 < c0 ) && (c0 < c11))1853 s = (c11-c0)/(c11-c10);1854 else if (( c11 < c0 ) && (c0 < c10))1855 s = (c11-c0) /(c11-c10);1856 else break;1857 R2 VP = V10*s + V11*(1-s);1858 int sss = (c11-c10) >0 ? 1 : -1;1859 // find the 2 point by dichotomie1860 //cout << " t =" << Number(t) << " c0 " << c0 ;1861 Int4 ii0 = BeginNewPoint[j1];1862 Int4 ii1 = EndNewPoint[j1];1863 Real8 ciii=-1,cii0=-1,cii1=-1 ;1864 if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 )1865 stack[kstack++] = ii0;//cout << " add+0 " << ii0;1866 else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 )1867 stack[kstack++] = ii1;//cout << " add+1 " << ii1;1868 else {1869 while ((ii1-ii0)> 1) {1870 Int4 iii = (ii0+ii1)/2;1871 ciii = vi0.m( D ,(R2) vertices[iii]);1872 //cout << " (iii = " << iii << " " << ciii << ") ";1873 if ( sss * (ciii - c0) <0 ) ii0 = iii;1874 else ii1 = iii;}1875 stack[kstack++] = ii0;// cout << " add0 " << ii0;1876 if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1;1877 }1878 if (kstack >5) // bug ?1879 cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl;1880 }1881 1882 stack[kstack++] = -1; // to stop1883 Int4 i1;1884 kstack =0;1885 while( (i1=stack[kstack++]) >= 0)1886 { // the two parameter is i0 and i11887 assert(i1 < nbv && i1 >= 0);1888 assert(i0 < nbv && i0 >= 0);1889 assert(i1 != i0);1890 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0];1891 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01));1892 1893 1894 #ifdef TRACETRIANGLE1895 if(trace) {1896 cout << "\n test j" << j <<" " << i01897 << " " << i1 << " d01=" << d01 <<endl;}1898 #endif1899 assert (i0 >= nbvold);1900 assert (i1 >= nbvold);1901 assert(i0 != i1);1902 if (d01 == 0)1903 break;1904 if ( d01 < seuil)1905 if (i1<nbvold) {1906 // remove all the points i0;1907 register Int4 ip,ipp;1908 for (ip=i0;i0 != (ipp = vertices[ip].ReferenceNumber);ip=ipp)1909 vertices[ip].ReferenceNumber = -1;// mark remove1910 vertices[ip].ReferenceNumber = -1;// mark remove1911 }1912 else {1913 // remove on of two points1914 register Int4 ip0, ip1, ipp0,ipp1;1915 register int kk0=1,kk1=1;1916 // count the number of common points to compute weight w0,w11917 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) kk0++;1918 for (ip1=i1;i1 != (ipp1 = vertices[ip1].ReferenceNumber);ip1=ipp1) kk1++;1919 1920 register Real8 w0 = ((Real8) kk0)/(kk0+kk1);1921 register Real8 w1 = ((Real8) kk1)/(kk0+kk1);1922 1923 // make a circular link1924 Exchange(vertices[i0].ReferenceNumber,vertices[i1].ReferenceNumber);1925 // the new coordinate1926 R2 C //= vertices[i0] ;1927 = vertices[i0].r *w0 + vertices[i1].r* w1;1928 1929 #ifdef TRACETRIANGLE1930 if(trace) {1931 cout << "\n ref = "<< vertices[i0].ref << " " <<vertices[i1].ref <<endl;1932 }1933 #endif1934 // update the new point points of the list1935 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0)1936 vertices[ip0].r = C;1937 vertices[ip0].r = C;1938 }1939 }1940 } // for (i0= ....1941 }// for triangle1942 1943 // remove of all the double points1944 1945 Int4 ip,ipp,kkk=nbvold;1946 for (i=nbvold;i<nbv;i++)1947 if (vertices[i].ReferenceNumber>=0) {// good points1948 // cout <<" i = " << i ;1949 for (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp)1950 vertices[ip].ReferenceNumber = -1;// mark remove1951 vertices[ip].ReferenceNumber = -1;// mark remove1952 // cout << i << " ---> " << kkk << endl;1953 vertices[kkk] = vertices[i];1954 vertices[kkk].i = toI2(vertices[kkk].r);1955 vertices[kkk++].ReferenceNumber = 0;1956 1957 }1958 1959 // insertion part ---1960 1961 const Int4 nbvnew = kkk-nbvold;1962 1963 cout << " Remove " << nbv - kkk << " to close vertex " ;1964 cout << " and Insert the " <<nbvnew<< " new points " << endl;1965 nbv = kkk;1966 Int4 NbSwap=0;1967 Icoor2 dete[3];1968 1969 // construction d'un ordre aleatoire1970 if (! nbvnew)1971 break;1972 if (nbvnew) {1973 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;1974 Int4 k3 = rand()%nbvnew ;1975 for (Int4 is3=0; is3<nbvnew; is3++)1976 ordre[nbvold+is3]= &vertices[nbvold +(k3 = (k3 + PrimeNumber)% nbvnew)];1977 1978 for (i=nbvold;i<nbv;i++)1979 { Vertex * vi = ordre[i];1980 Triangle *tcvi = FindTriangleContening(vi->i,dete);1981 // Vertex * nv = quadtree->NearestVertex(vi->i.x,vi->i.y);1982 // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is "1983 // << Number(nv) << nv->i << endl;1984 // Int4 kt = Number(tcvi);1985 //1986 1987 quadtree->Add(*vi); //1988 assert (tcvi->det >= 0) ;// internal1989 Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);1990 }1991 }1992 cout << " Nb swap = " << NbSwap << " after " ;1993 1994 for (i=nbvold;i<nbv;i++)1995 NbSwap += vertices[i].Optim(1);1996 cout << NbSwap << endl;1997 1998 for (i=nbtold;i<nbt;i++)1999 first_np_or_next_t[i]=1;2000 2001 Headt = nbt; // empty list2002 for (i=nbvold;i<nbv;i++)2003 { // for all the triangle contening the vertex i2004 Vertex * s = vertices + i;2005 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);2006 Triangle * tbegin= (Triangle*) ta;2007 Int4 kt;2008 do {2009 kt = Number((Triangle*) ta);2010 if (first_np_or_next_t[kt]>0)2011 first_np_or_next_t[kt]=-Headt,Headt=kt;2012 assert( ta.EdgeVertex(0) == s);2013 ta = Next(Adj(ta));2014 } while ( (tbegin != (Triangle*) ta));2015 }2016 2017 2018 } while (nbv!=nbvold);2019 delete [] first_np_or_next_t;2020 cout << " end : Triangles::NewPoints old nbv=" << nbv << endl;2021 2022 }2023 2024 void Triangles::Insert()2025 {2026 long int verbosity=2;2027 if (verbosity>2) cout << " -- Insert initial " << nbv << " vertices " << endl ;2028 Triangles * OldCurrentTh =CurrentTh;2029 2030 CurrentTh=this;2031 double time0=CPUtime(),time1,time2,time3;2032 SetIntCoor();2033 Int4 i;2034 for (i=0;i<nbv;i++)2035 ordre[i]= &vertices[i] ;2036 2037 // construction d'un ordre aleatoire2038 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;2039 Int4 k3 = rand()%nbv ;2040 for (int is3=0; is3<nbv; is3++)2041 ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv];2042 2043 2044 2045 2046 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)2047 if ( ++i >= nbv) {2048 cerr << " All the vertices are aline " << endl;2049 MeshError(998,this); }2050 2051 // echange i et 2 dans ordre afin2052 // que les 3 premiers ne soit pas aligne2053 Exchange( ordre[2], ordre[i]);2054 2055 // on ajoute un point a l'infini pour construire le maillage2056 // afin d'avoir une definition simple des aretes frontieres2057 nbt = 2;2058 2059 2060 // on construit un maillage trivale forme2061 // d'une arete et de 2 triangles2062 // construit avec le 2 aretes orientes et2063 Vertex * v0=ordre[0], *v1=ordre[1];2064 2065 triangles[0](0) = 0; // sommet pour infini2066 triangles[0](1) = v0;2067 triangles[0](2) = v1;2068 2069 triangles[1](0) = 0;// sommet pour infini2070 triangles[1](2) = v0;2071 triangles[1](1) = v1;2072 const int e0 = OppositeEdge[0];2073 const int e1 = NextEdge[e0];2074 const int e2 = PreviousEdge[e0];2075 triangles[0].SetAdj2(e0, &triangles[1] ,e0);2076 triangles[0].SetAdj2(e1, &triangles[1] ,e2);2077 triangles[0].SetAdj2(e2, &triangles[1] ,e1);2078 2079 triangles[0].det = -1; // faux triangles2080 triangles[1].det = -1; // faux triangles2081 2082 triangles[0].SetTriangleContainingTheVertex();2083 triangles[1].SetTriangleContainingTheVertex();2084 2085 triangles[0].link=&triangles[1];2086 triangles[1].link=&triangles[0];2087 2088 // nbtf = 2;2089 if ( !quadtree ) quadtree = new QuadTree(this,0);2090 quadtree->Add(*v0);2091 quadtree->Add(*v1);2092 2093 // on ajoute les sommets un Ò un2094 Int4 NbSwap=0;2095 2096 time1=CPUtime();2097 2098 if (verbosity>3) cout << " -- Begin of insertion process " << endl;2099 2100 for (Int4 icount=2; icount<nbv; icount++) {2101 Vertex *vi = ordre[icount];2102 // cout << " Insert " << Number(vi) << endl;2103 Icoor2 dete[3];2104 Triangle *tcvi = FindTriangleContening(vi->i,dete);2105 quadtree->Add(*vi);2106 Add(*vi,tcvi,dete);2107 NbSwap += vi->Optim(1,0);2108 2109 }// fin de boucle en icount2110 time2=CPUtime();2111 if (verbosity>3)2112 cout << " NbSwap of insertion " << NbSwap2113 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv2114 << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv "2115 << (float)NbUnSwap /(float) nbv2116 <<endl;2117 NbUnSwap = 0;2118 // construction d'un ordre aleatoire2119 // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ;2120 #ifdef NBLOOPOPTIM2121 2122 k3 = rand()%nbv ;2123 for (int is4=0; is4<nbv; is4++)2124 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv];2125 2126 double timeloop = time2 ;2127 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++)2128 {2129 double time000 = timeloop;2130 Int4 NbSwap = 0;2131 for (int is1=0; is1<nbv; is1++)2132 NbSwap += ordre[is1]->Optim(0,0);2133 timeloop = CPUtime();2134 if (verbosity>3)2135 cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap2136 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv2137 << " CPU=" << timeloop - time000 << " s, "2138 << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv2139 << endl;2140 NbUnSwap = 0;2141 if(!NbSwap) break;2142 }2143 ReMakeTriangleContainingTheVertex();2144 // because we break the TriangleContainingTheVertex2145 #endif2146 time3=CPUtime();2147 if (verbosity>4)2148 cout << " init " << time1 - time0 << " initialisation, "2149 << time2 - time1 << "s, insert point "2150 << time3 -time2 << "s, optim " << endl2151 << " Init Total Cpu Time = " << time3 - time0 << "s " << endl;2152 2153 CurrentTh=OldCurrentTh;2154 }2155 2156 2157 void Triangles::ForceBoundary()2158 {2159 long int verbosity=2;2160 if (verbosity > 2)2161 cout << " -- ForceBoundary nb of edge " << nbe << endl;2162 int k=0;2163 Int4 nbfe=0,nbswp=0,Nbswap=0;2164 for (Int4 t = 0; t < nbt; t++)2165 if (!triangles[t].det)2166 k++,cerr << " det T" << t << " = " << 0 << endl;2167 if (k!=0) {2168 cerr << " ther is " << k << " triangles of mes = 0 " << endl;2169 MeshError(11,this);}2170 2171 TriangleAdjacent ta(0,0);2172 for (Int4 i = 0; i < nbe; i++)2173 {2174 nbswp = ForceEdge(edges[i][0],edges[i][1],ta);2175 2176 if ( nbswp < 0) k++;2177 else Nbswap += nbswp;2178 if (nbswp) nbfe++;2179 if ( nbswp < 0 && k < 5)2180 {2181 cerr << " Missing Edge " << i << " v0 = " << Number(edges[i][0]) << edges[i][0].r2182 <<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << " " << (Triangle *) ta ;2183 if(ta.t)2184 {2185 Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1);2186 cerr << " crossing with [" << Number(aa) << ", " << Number(bb) << "]\n";2187 }2188 else cerr << endl;2189 2190 }2191 if ( nbswp >=0 && edges[i].on->Cracked())2192 ta.SetCracked();2193 }2194 2195 2196 if (k!=0) {2197 cerr << " they is " << k << " lost edges " << endl;2198 cerr << " The boundary is crossing may be!" << endl;2199 MeshError(10,this);2200 }2201 for (Int4 j=0;j<nbv;j++)2202 Nbswap += vertices[j].Optim(1,0);2203 if (verbosity > 3)2204 cout << " Nb of inforced edge = " << nbfe << " Nb of Swap " << Nbswap << endl;2205 2206 }2207 2208 void Triangles::FindSubDomain(int OutSide=0)2209 {2210 long int verbosity=0;2211 2212 if (verbosity >2)2213 {2214 if (OutSide)2215 cout << " -- Find all external sub-domain ";2216 else2217 cout << " -- Find all internal sub-domain ";2218 if(verbosity>99)2219 {2220 2221 for(int i=0;i<nbt;++i)2222 cout << i<< " " << triangles[i] << endl;2223 }2224 2225 }2226 // if (verbosity > 4) cout << " OutSide=" << OutSide << endl;2227 short * HeapArete = new short[nbt];2228 Triangle ** HeapTriangle = new Triangle* [nbt];2229 Triangle *t,*t1;2230 Int4 k,it;2231 2232 for (Int4 itt=0;itt<nbt;itt++)2233 triangles[itt].link=0; // par defaut pas de couleur2234 2235 Int4 NbSubDomTot =0;2236 for ( it=0;it<nbt;it++) {2237 if ( ! triangles[it].link ) {2238 t = triangles + it;2239 NbSubDomTot++;; // new composante connexe2240 Int4 i = 0; // niveau de la pile2241 t->link = t ; // sd forme d'un triangle cicular link2242 2243 HeapTriangle[i] =t ;2244 HeapArete[i] = 3;2245 2246 while (i >= 0) // boucle sur la pile2247 { while ( HeapArete[i]--) // boucle sur les 3 aretes2248 {2249 int na = HeapArete[i];2250 Triangle * tc = HeapTriangle[i]; // triangle courant2251 if( ! tc->Locked(na)) // arete non frontiere2252 {2253 Triangle * ta = tc->TriangleAdj(na) ; // næ triangle adjacent2254 if (ta->link == 0 ) // non deja chainer => on enpile2255 {2256 i++;2257 ta->link = t->link ; // on chaine les triangles2258 t->link = ta ; // d'un meme sous domaine2259 HeapArete[i] = 3; // pour les 3 triangles adjacents2260 HeapTriangle[i] = ta;2261 }}2262 } // deplie fin de boucle sur les 3 adjacences2263 i--;2264 }2265 }2266 }2267 2268 // supression de tous les sous domaine infini <=> contient le sommet NULL2269 it =0;2270 NbOutT = 0;2271 while (it<nbt) {2272 if (triangles[it].link)2273 {2274 if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) ))2275 {2276 // infini triangle2277 NbSubDomTot --;2278 // cout << " triangle infini " << it << triangles[it] << endl;2279 t=&triangles[it];2280 NbOutT--; // on fait un coup de trop.2281 while (t){ // cout << Number(t) << " " << endl;2282 NbOutT++;2283 t1=t;2284 t=t->link;2285 t1->link=0;}//while (t)2286 }2287 }2288 it++;} // end while (it<nbt)2289 if (nbt == NbOutT || !NbSubDomTot)2290 {2291 cout << "\n error : " << NbOutT << " " << NbSubDomTot <<" " << nbt << endl;2292 cerr << "Error: The boundary is not close => All triangles are outside " << endl;2293 MeshError(888,this);2294 }2295 2296 delete [] HeapArete;2297 delete [] HeapTriangle;2298 2299 2300 if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains )2301 { // No geom sub domain2302 Int4 i;2303 if (subdomains) delete [] subdomains;2304 subdomains = new SubDomain[ NbSubDomTot];2305 NbSubDomains= NbSubDomTot;2306 for ( i=0;i<NbSubDomains;i++) {2307 subdomains[i].head=NULL;2308 subdomains[i].ref=i+1;2309 }2310 Int4 * mark = new Int4[nbt];2311 for (it=0;it<nbt;it++)2312 mark[it]=triangles[it].link ? -1 : -2;2313 2314 it =0;2315 k = 0;2316 while (it<nbt) {2317 if (mark[it] == -1) {2318 t1 = & triangles[it];2319 t = t1->link;2320 mark[it]=k;2321 subdomains[k].head = t1;2322 // cout << " New -- " << Number(t1) << " " << it << endl;2323 do {// cout << " k " << k << " " << Number(t) << endl;2324 mark[Number(t)]=k;2325 t=t->link;2326 } while (t!=t1);2327 mark[it]=k++;}2328 // else if(mark[it] == -2 ) triangles[it].Draw(999);2329 it++;} // end white (it<nbt)2330 assert(k== NbSubDomains);2331 if(OutSide)2332 {2333 // to remove all the sub domain by parity adjacents2334 // because in this case we have only the true boundary edge2335 // so teh boundary is manifold2336 Int4 nbk = NbSubDomains;2337 while (nbk)2338 for (it=0;it<nbt && nbk ;it++)2339 for (int na=0;na<3 && nbk ;na++)2340 {2341 Triangle *ta = triangles[it].TriangleAdj(na);2342 Int4 kl = ta ? mark[Number(ta)] : -2;2343 Int4 kr = mark[it];2344 if(kr !=kl) {2345 //cout << kl << " " << kr << " rl " << subdomains[kl].ref2346 // << " rr " << subdomains[kr].ref ;2347 if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0)2348 nbk--,subdomains[kr].ref=subdomains[kl].ref-1;2349 if (kr >=0 && subdomains[kr].ref <0 && kl >=0 && subdomains[kl].ref>=0)2350 nbk--,subdomains[kl].ref=subdomains[kr].ref-1;2351 if(kr<0 && kl >=0 && subdomains[kl].ref>=0)2352 nbk--,subdomains[kl].ref=-1;2353 if(kl<0 && kr >=0 && subdomains[kr].ref>=0)2354 nbk--,subdomains[kr].ref=-1;2355 // cout << " after \t "2356 // << kl << subdomains[kl].ref << " rr " << kr2357 // << subdomains[kr].ref << endl;2358 }2359 }2360 // cout << subdomains[0].ref << subdomains[1].ref << endl;2361 Int4 j=0;2362 for ( i=0;i<NbSubDomains;i++)2363 if((-subdomains[i].ref) %2) { // good2364 //cout << " sudom ok = " << i << " " << subdomains[i].ref2365 // << " " << (-subdomains[i].ref) %2 << endl;2366 if(i != j)2367 Exchange(subdomains[i],subdomains[j]);2368 j++;}2369 else2370 { //cout << " remove sub domain " << i << endl;2371 t= subdomains[i].head;2372 while (t){// cout << Number(t) << " " << endl;2373 NbOutT++;2374 t1=t;2375 t=t->link;2376 t1->link=0;}//while (t)2377 }2378 2379 if(verbosity>4)2380 cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl;2381 NbSubDomains=j;2382 }2383 2384 delete [] mark;2385 2386 }2387 else2388 { // find the head for all sub domaine2389 if (Gh.NbSubDomains != NbSubDomains && subdomains)2390 delete [] subdomains, subdomains=0;2391 if (! subdomains )2392 subdomains = new SubDomain[ Gh.NbSubDomains];2393 NbSubDomains =Gh.NbSubDomains;2394 if(verbosity>4)2395 cout << " find the " << NbSubDomains << " sub domain " << endl;2396 Int4 err=0;2397 ReMakeTriangleContainingTheVertex();2398 Int4 * mark = new Int4[nbt];2399 Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge();2400 2401 for (it=0;it<nbt;it++)2402 mark[it]=triangles[it].link ? -1 : -2;2403 Int4 inew =0;2404 for (Int4 i=0;i<NbSubDomains;i++)2405 {2406 GeometricalEdge &eg = *Gh.subdomains[i].edge;2407 subdomains[i].ref = Gh.subdomains[i].ref;2408 int ssdlab = subdomains[i].ref;2409 // by carefull is not easy to find a edge create from a GeometricalEdge2410 // see routine MakeGeometricalEdgeToEdge2411 Edge &e = *GeometricalEdgetoEdge[Gh.Number(eg)];2412 assert(&e);2413 Vertex * v0 = e(0),*v1 = e(1);2414 Triangle *t = v0->t;2415 int sens = Gh.subdomains[i].sens;2416 // test if ge and e is in the same sens2417 // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref2418 // << " ref edge =" << eg.ref << " sens " << sens ;2419 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0)2420 sens = -sens ;2421 subdomains[i].sens = sens;2422 subdomains[i].edge = &e;2423 // cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh " << e[0].r << e[1].r << endl;2424 // cout << " v0 , v1 = " << Number(v0) << " " << Number(v1) << endl;2425 assert(t && sens);2426 2427 TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges2428 2429 while (1)2430 {2431 assert( v0 == ta.EdgeVertex(1) );2432 // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl;2433 if (ta.EdgeVertex(0) == v1) { // ok we find the edge2434 if (sens>0)2435 subdomains[i].head=t=Adj(ta);2436 else2437 subdomains[i].head=t=ta;2438 //cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl;2439 if(t<triangles || t >= triangles+nbt || t->det < 02440 || t->link == 0) // Ajoute aout 2002441 {2442 cerr << " Error in the def of sub domain "<<i2443 << " form border " << NbSubDomains - i << "/" << NbSubDomains2444 << ": Bad sens " << Gh.Number(eg) <<" "<< sens << endl;2445 err++;2446 break;}2447 Int4 it = Number(t);2448 if (mark[it] >=0) {2449 if(verbosity>10)2450 cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref2451 << " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref2452 << " skip this def " << endl;2453 break;}2454 if(i != inew)2455 Exchange(subdomains[i],subdomains[inew]);2456 inew++;2457 Triangle *tt=t;2458 Int4 kkk=0;2459 do2460 {2461 kkk++;2462 assert(mark[Number(tt)]<0);2463 mark[Number(tt)]=i;2464 tt=tt->link;2465 } while (tt!=t);2466 if(verbosity>7)2467 cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl;2468 break;}2469 ta = Previous(Adj(ta));2470 if(t == (Triangle *) ta) {2471 err++;2472 cerr << " Error in the def of sub domain " << i2473 << " edge=" << Gh.Number(eg) << " " << sens << endl;2474 break;}2475 // cout << " NB of remove subdomain " << NbSubDomTot-NbSubDomains<< endl;2476 2477 }2478 2479 }2480 if (err) MeshError(777,this);2481 2482 if (inew < NbSubDomains) {2483 if (verbosity>5)2484 cout << " Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl;2485 NbSubDomains=inew;}2486 2487 2488 for (it=0;it<nbt;it++)2489 if ( mark[it] ==-1 )2490 NbOutT++,triangles[it].link =0;2491 delete [] GeometricalEdgetoEdge;2492 delete [] mark;2493 2494 }2495 NbOutT=0;2496 for (it=0;it<nbt;it++)2497 if(!triangles[it].link) NbOutT++;2498 if (verbosity> 4)2499 cout << " " ;2500 if (verbosity> 2)2501 cout << " Nb of Sub borned Domain = " << NbSubDomTot << " NbOutTriangles = " << NbOutT <<endl;2502 2503 2504 }2505 void Triangles::ReNumberingVertex(Int4 * renu)2506 {2507 // warning be carfull because pointeur2508 // from on mesh to over mesh2509 // -- so do ReNumbering a the beginning2510 Vertex * ve = vertices+nbv;2511 Int4 it,ie,i;2512 2513 for ( it=0;it<nbt;it++)2514 triangles[it].ReNumbering(vertices,ve,renu);2515 2516 for ( ie=0;ie<nbe;ie++)2517 edges[ie].ReNumbering(vertices,ve,renu);2518 2519 for (i=0;i< NbVerticesOnGeomVertex;i++)2520 {2521 Vertex *v = VerticesOnGeomVertex[i].mv;2522 if (v>=vertices && v < ve)2523 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];2524 }2525 2526 for (i=0;i< NbVerticesOnGeomEdge;i++)2527 {2528 Vertex *v =VerticesOnGeomEdge[i].mv;2529 if (v>=vertices && v < ve)2530 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];2531 }2532 2533 for (i=0;i< NbVertexOnBThVertex;i++)2534 {2535 Vertex *v=VertexOnBThVertex[i].v;2536 if (v>=vertices && v < ve)2537 VertexOnBThVertex[i].v=vertices+renu[Number(v)];2538 }2539 2540 for (i=0;i< NbVertexOnBThEdge;i++)2541 {2542 Vertex *v=VertexOnBThEdge[i].v;2543 if (v>=vertices && v < ve)2544 VertexOnBThEdge[i].v=vertices+renu[Number(v)];2545 }2546 2547 // move the Vertices without a copy of the array2548 // be carefull not trivial code2549 Int4 j;2550 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu2551 if (renu[it] >= 0) // a new sub cycle2552 {2553 i=it;2554 Vertex ti=vertices[i],tj;2555 while ( (j=renu[i]) >= 0)2556 { // i is old, and j is new2557 renu[i] = -1-renu[i]; // mark2558 tj = vertices[j]; // save new2559 vertices[j]= ti; // new <- old2560 i=j; // next2561 ti = tj;2562 }2563 }2564 if (quadtree)2565 { delete quadtree;2566 quadtree = new QuadTree(this);2567 }2568 for ( it=0;it<nbv;it++)2569 renu[i]= -renu[i]-1;2570 2571 }2572 void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress)2573 {2574 long int verbosity=0;2575 Int4 *renu= new Int4[nbt];2576 register Triangle *t0,*t,*te=triangles+nbt;2577 register Int4 k=0,it,i,j;2578 2579 for ( it=0;it<nbt;it++)2580 renu[it]=-1; // outside triangle2581 for ( i=0;i<NbSubDomains;i++)2582 {2583 t=t0=subdomains[i].head;2584 assert(t0); // no empty sub domain2585 do {2586 Int4 kt = Number(t);2587 assert(kt>=0 && kt < nbt );2588 assert(renu[kt]==-1);2589 renu[kt]=k++;2590 }2591 while (t0 != (t=t->link));2592 }2593 if (verbosity>9)2594 cout << " number of inside triangles " << k << " nbt = " << nbt << endl;2595 // take is same numbering if possible2596 if(justcompress)2597 for ( k=0,it=0;it<nbt;it++)2598 if(renu[it] >=0 )2599 renu[it]=k++;2600 2601 // put the outside triangles at the end2602 for ( it=0;it<nbt;it++)2603 if (renu[it]==-1)2604 renu[it]=k++;2605 2606 assert(k == nbt);2607 // do the change on all the pointeur2608 for ( it=0;it<nbt;it++)2609 triangles[it].ReNumbering(triangles,te,renu);2610 2611 for ( i=0;i<NbSubDomains;i++)2612 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];2613 2614 // move the Triangles without a copy of the array2615 // be carefull not trivial code2616 for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu2617 if (renu[it] >= 0) // a new sub cycle2618 {2619 i=it;2620 Triangle ti=triangles[i],tj;2621 while ( (j=renu[i]) >= 0)2622 { // i is old, and j is new2623 renu[i] = -1; // mark2624 tj = triangles[j]; // save new2625 triangles[j]= ti; // new <- old2626 i=j; // next2627 ti = tj;2628 }2629 }2630 delete [] renu;2631 nt = nbt - NbOutT;2632 2633 }2634 Int4 Triangles::ConsRefTriangle(Int4 *reft) const2635 {2636 long int verbosity=0;2637 assert(reft);2638 register Triangle *t0,*t;2639 register Int4 k=0, num;2640 for (Int4 it=0;it<nbt;it++)2641 reft[it]=-1; // outside triangle2642 for (Int4 i=0;i<NbSubDomains;i++)2643 {2644 t=t0=subdomains[i].head;2645 assert(t0); // no empty sub domain2646 // register Int4 color=i+1;// because the color 0 is outside triangle2647 do { k++;2648 num = Number(t);2649 assert(num>=0 &&num < nbt);2650 reft[num]=i;2651 // cout << Number(t0) << " " <<Number(t)<< " " << i << endl;2652 }2653 while (t0 != (t=t->link));2654 }2655 // NbOutT = nbt - k;2656 if (verbosity>5)2657 cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k2658 << " Nbt = " << nbt << " Out Triangles = " << nbt - k << endl;2659 2660 return k;2661 2662 }2663 2664 /*2665 void Triangles::ConsLinkTriangle()2666 {2667 for (Int4 i=0;i<NbSubDomains;i++)2668 subdomains[i].head=0;2669 register Triangle * t=triangles,*tend = triangles+nbt,*hst;2670 for (;t<tend;t++)2671 {2672 if (t->link)2673 {2674 register Int4 color = t->color-1;2675 assert(color<NbSubDomains && color>=0);2676 if (hst=subdomains[color].head) {2677 t->link=hst->link;2678 hst->link=t; }2679 else {2680 subdomains[color].head = t;2681 t->link=t;}// circular link2682 }2683 }2684 {2685 for (Int4 i=0;i<NbSubDomains;i++)2686 assert(subdomains[i].head);2687 }2688 2689 }2690 */2691 /* void Triangles::RandomInit()2692 {2693 // Meshbegin = vertices;2694 // Meshend = vertices + nbvx;2695 nbv = nbvx;2696 for (int i = 0; i < nbv; i++)2697 {2698 vertices[i].r.x= rand();2699 vertices[i].r.y= rand();2700 vertices[i].ref = 0;2701 }2702 }2703 void Triangles::CubeInit(int n,int m)2704 {2705 // nbvx = n*m;2706 // Meshbegin = vertices;2707 // Meshend = vertices + nbvx;2708 nbv = n*m;2709 assert(nbv <= nbvx);2710 int k =0;2711 for (int i = 0; i < n; i++)2712 for (int j = 0; j < m; j++)2713 {2714 vertices[k].r.x= i;2715 vertices[k].r.y= j;2716 vertices[k++].ref = 0;2717 }2718 }2719 */2720 1076 Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j) 2721 1077 { return quadtree->NearestVertex(i,j); } … … 2804 1160 << " " << vertices 2805 1161 << " " << ordre << " " << triangles <<endl; 2806 2807 1162 } 2808 1163 2809 void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)2810 {2811 Gh.NbRef++;// add a ref to Gh2812 2813 2814 /*************************************************************************2815 // methode in 2 step2816 // 1 - compute the number of new edge to allocate2817 // 2 - construct the edge2818 remark:2819 in this part we suppose to have a background mesh with the same2820 geometry2821 2822 To construct the discretisation of the new mesh we have to2823 rediscretize the boundary of background Mesh2824 because we have only the pointeur from the background mesh to the geometry.2825 We need the abcisse of the background mesh vertices on geometry2826 so a vertex is2827 0 on GeometricalVertex ;2828 1 on GeometricalEdge + abcisse2829 2 internal2830 2831 *************************************************************************/2832 assert(&BTh.Gh == &Gh);2833 2834 long int verbosity=0;2835 BTh.NbRef++; // add a ref to BackGround Mesh2836 PreInit(inbvx);2837 BTh.SetVertexFieldOn();2838 int * bcurve = new int[Gh.NbOfCurves]; //2839 2840 // we have 2 ways to make the loop2841 // 1) on the geometry2842 // 2) on the background mesh2843 // if you do the loop on geometry, we don't have the pointeur on background,2844 // and if you do the loop in background we have the pointeur on geometry2845 // so do the walk on background2846 // Int4 NbVerticesOnGeomVertex;2847 // VertexOnGeom * VerticesOnGeomVertex;2848 // Int4 NbVerticesOnGeomEdge;2849 // VertexOnGeom * VerticesOnGeomEdge;2850 2851 2852 NbVerticesOnGeomVertex=0;2853 NbVerticesOnGeomEdge=0;2854 //1 copy of the Required vertex2855 int i;2856 for ( i=0;i<Gh.nbv;i++)2857 if (Gh[i].Required()) NbVerticesOnGeomVertex++;2858 2859 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];2860 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];2861 //2862 if( NbVerticesOnGeomVertex >= nbvx)2863 {2864 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;2865 MeshError(1,this);2866 }2867 assert(vertices);2868 for (i=0;i<Gh.nbv;i++)2869 if (Gh[i].Required()) {//Gh vertices Required2870 vertices[nbv] = Gh[i];2871 vertices[nbv].i = I2(0,0);2872 Gh[i].to = vertices + nbv;// save Geom -> Th2873 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);2874 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;2875 nbv++;}2876 else Gh[i].to=0;2877 //2878 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++)2879 {2880 VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i];2881 if (vog.IsRequiredVertex()) {2882 GeometricalVertex * gv = vog;2883 Vertex *bv = vog;2884 assert(gv->to);// use of Geom -> Th2885 VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv);2886 gv->to->m = bv->m; // for taking the metrix of the background mesh2887 ;}2888 }2889 assert(NbVertexOnBThVertex == NbVerticesOnGeomVertex);2890 // new stuff FH with curve2891 // find the begin of the curve in BTh2892 {2893 Gh.UnMarkEdges();2894 int bfind=0;2895 /*2896 cout << " nb curves = " << Gh.NbOfCurves << endl;2897 for(int i=0;i<Gh.NbOfCurves ;i++)2898 {2899 cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb2900 << " end e= " << Gh.Number(Gh.curves[i].ee) << " k=" << Gh.curves[i].ke << endl;2901 }*/2902 for (int i=0;i<Gh.NbOfCurves;i++)2903 {2904 bcurve[i]=-1;2905 }2906 2907 for (int iedge=0;iedge<BTh.nbe;iedge++)2908 {2909 Edge & ei = BTh.edges[iedge];2910 for(int je=0;je<2;je++) // for the 2 extremites2911 if (!ei.on->Mark() && ei[je].on->IsRequiredVertex() )2912 {2913 // a begin of curve2914 int nc = ei.on->CurveNumber;2915 2916 //cout << "curve " << nc << " v " << Gh.Number((GeometricalVertex *) *ei[je].on) << " "2917 // << " e " << " " << Gh.Number(ei.on) << " vc " << Gh.Number((*Gh.curves[nc].be)[Gh.curves[nc].kb]) << endl;2918 if(2919 ei.on==Gh.curves[nc].be &&2920 (GeometricalVertex *) *ei[je].on == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity2921 )2922 {2923 // cout << " find " << endl;2924 bcurve[nc]=iedge*2+je;2925 bfind++;2926 }2927 }2928 }2929 assert( bfind==Gh.NbOfCurves);2930 }2931 // method in 2 + 1 step2932 // 0.0) compute the length and the number of vertex to do allocation2933 // 1.0) recompute the length2934 // 1.1) compute the vertex2935 Int4 nbex=0,NbVerticesOnGeomEdgex=0;2936 for (int step=0; step <2;step++)2937 {2938 Int4 NbOfNewPoints=0;2939 Int4 NbOfNewEdge=0;2940 Int4 iedge;2941 Gh.UnMarkEdges();2942 /* add Curve loop FH2943 // find a starting back groud edges to walk2944 for (iedge=0;iedge<BTh.nbe;iedge++) {2945 Edge & ei = BTh.edges[iedge];2946 for(int jedge=0;jedge<2;jedge++) // for the 2 extremites2947 if (!ei.on->Mark() && ei[jedge].on->IsRequiredVertex() )2948 {2949 */2950 // new code FH 20042951 Real8 L=0;2952 for (int icurve=0;icurve<Gh.NbOfCurves;icurve++)2953 {2954 iedge=bcurve[icurve]/2;2955 int jedge=bcurve[icurve]%2;2956 if( ! Gh.curves[icurve].master) continue; // we skip all equi curve2957 Edge & ei = BTh.edges[iedge];2958 // warning: ei.on->Mark() can be change in2959 // loop for(jedge=0;jedge<2;jedge++)2960 // new curve2961 // good the find a starting edge2962 Real8 Lstep=0,Lcurve=0;// step between two points (phase==1)2963 Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)2964 2965 2966 // cout.precision(16);2967 for(int phase=0;phase<=step;phase++)2968 {2969 2970 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)2971 {2972 2973 int icurveequi= Gh.Number(curve);2974 2975 if( phase == 0 && icurveequi != icurve) continue;2976 2977 int k0=jedge,k1;2978 Edge * pe= BTh.edges+iedge;2979 //GeometricalEdge *ong = ei.on;2980 int iedgeequi=bcurve[icurveequi]/2;2981 int jedgeequi=bcurve[icurveequi]%2;2982 2983 int k0equi=jedgeequi,k1equi;2984 Edge * peequi= BTh.edges+iedgeequi;2985 GeometricalEdge *ongequi = peequi->on;2986 2987 Real8 sNew=Lstep;// abcisse of the new points (phase==1)2988 L=0;// length of the curve2989 Int4 i=0;// index of new points on the curve2990 register GeometricalVertex * GA0 = *(*peequi)[k0equi].on;2991 Vertex *A0;2992 A0 = GA0->to; // the vertex in new mesh2993 Vertex *A1;2994 VertexOnGeom *GA1;2995 Edge * PreviousNewEdge = 0;2996 // cout << " --------------New Curve phase " << phase2997 // << "---------- A0=" << *A0 << ei[k0] <<endl;2998 assert (A0-vertices>=0 && A0-vertices <nbv);2999 if(ongequi->Required() )3000 {3001 GeometricalVertex *GA1 = *(*peequi)[1-k0equi].on;3002 A1 = GA1->to; //3003 }3004 else3005 for(;;)3006 {3007 // assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe);3008 Edge &ee=*pe;3009 Edge &eeequi=*peequi;3010 k1 = 1-k0; // next vertex of the edge3011 k1equi= 1 - k0equi;3012 3013 assert(pe && ee.on);3014 ee.on->SetMark();3015 Vertex & v0=ee[0], & v1=ee[1];3016 R2 AB= (R2) v1 - (R2) v0;3017 Real8 L0=L,LAB;3018 LAB = LengthInterpole(v0.m,v1.m,AB);3019 L+= LAB;3020 if (phase) {// computation of the new points3021 while ((i!=NbCreatePointOnCurve) && sNew <= L) {3022 // cout << " L0= " << L0 << " L " << L << " sN="3023 // << sNew << " LAB=" << LAB << " NBPC =" <<NbCreatePointOnCurve<< " i " << i << endl;3024 assert (sNew >= L0);3025 assert(LAB);3026 3027 3028 assert(vertices && nbv<nbvx);3029 assert(edges && nbe < nbex);3030 assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge < NbVerticesOnGeomEdgex);3031 // new vertex on edge3032 A1=vertices+nbv++;3033 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;3034 Edge *e = edges + nbe++;3035 Real8 se= (sNew-L0)/LAB;3036 assert(se>=0 && se < 1.000000001);3037 se = abscisseInterpole(v0.m,v1.m,AB,se,1);3038 assert(se>=0 && se <= 1);3039 //((k1==1) != (k1==k1equi))3040 se = k1 ? se : 1. - se;3041 se = k1==k1equi ? se : 1. - se;3042 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save3043 ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);3044 A1->ReferenceNumber = eeequi.ref;3045 A1->DirOfSearch =NoDirOfSearch;3046 //cout << icurveequi << " " << i << " " << *A1 << endl;3047 e->on = ongequi;3048 e->v[0]= A0;3049 e->v[1]= A1;3050 if(verbosity>99)3051 cout << i << "+ New P "<< nbv-1 << " " <<sNew<< " L0=" << L03052 << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= "3053 << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl;3054 3055 e->ref = eeequi.ref;3056 e->adj[0]=PreviousNewEdge;3057 3058 if (PreviousNewEdge)3059 PreviousNewEdge->adj[1] = e;3060 PreviousNewEdge = e;3061 A0=A1;3062 sNew += Lstep;3063 // cout << " sNew = " << sNew << " L = " << L3064 // << " ------" <<NbCreatePointOnCurve << " " << i << endl;3065 if (++i== NbCreatePointOnCurve) break;3066 }3067 3068 }3069 assert(ee.on->CurveNumber==ei.on->CurveNumber);3070 if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl;3071 if ( ee[k1].on->IsRequiredVertex()) {3072 assert(eeequi[k1equi].on->IsRequiredVertex());3073 register GeometricalVertex * GA1 = *eeequi[k1equi].on;3074 A1=GA1->to;// the vertex in new mesh3075 assert (A1-vertices>=0 && A1-vertices <nbv);3076 break;}3077 if (!ee.adj[k1])3078 {cerr << "Error adj edge " << BTh.Number(ee) << ", nbe = " << nbe3079 << " Gh.vertices " << Gh.vertices3080 << " k1 = " << k1 << " on=" << *ee[k1].on << endl;3081 cerr << ee[k1].on->gv-Gh.vertices << endl;3082 }3083 pe = ee.adj[k1]; // next edge3084 k0 = pe->Intersection(ee);3085 peequi= eeequi.adj[k1equi]; // next edge3086 k0equi=peequi->Intersection(eeequi);3087 }// for(;;) end of the curve3088 3089 3090 if (phase) // construction of the last edge3091 {3092 Edge *e = edges + nbe++;3093 if (verbosity>10)3094 cout << " Fin curve A1" << *A1 << " " << icurve << " <=> " << icurveequi <<"-----" <<3095 NbCreatePointOnCurve << " == " <<i<<endl;3096 e->on = ongequi;3097 e->v[0]= A0;3098 e->v[1]= A1;3099 e->ref = peequi->ref;3100 e->adj[0]=PreviousNewEdge;3101 e->adj[1]=0;3102 if (PreviousNewEdge)3103 PreviousNewEdge->adj[1] = e;3104 PreviousNewEdge = e;3105 3106 assert(i==NbCreatePointOnCurve);3107 3108 }3109 } // end loop on equi curve3110 3111 if (!phase) { //3112 Int4 NbSegOnCurve = Max((Int4)(L+0.5),(Int4) 1);// nb of seg3113 Lstep = L/NbSegOnCurve;3114 Lcurve = L;3115 NbCreatePointOnCurve = NbSegOnCurve-1;3116 3117 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)3118 {3119 NbOfNewEdge += NbSegOnCurve;3120 NbOfNewPoints += NbCreatePointOnCurve;3121 }3122 if(verbosity>5)3123 cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep="3124 << Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve <<endl;3125 // do'nt3126 // if(NbCreatePointOnCurve<1) break;3127 }3128 }//for(phase;;)3129 /* modif FH add Curve class3130 }}//for (iedge=0;iedge<BTh.nbe;iedge++)3131 */3132 // new code Curve class3133 } // end of curve loop3134 // end new code3135 // do the allocation3136 if(step==0)3137 {3138 //if(!NbOfNewPoints) break;// nothing ????? bug3139 if(nbv+NbOfNewPoints > nbvx)3140 {3141 cerr << " Too much vertices on geometry " << nbv+NbOfNewPoints << " >= " << nbvx << endl;3142 MeshError(3,this);3143 }3144 //cout << " NbOfNewEdge" << NbOfNewEdge << " NbOfNewPoints " << NbOfNewPoints << endl;3145 edges = new Edge[NbOfNewEdge];3146 nbex = NbOfNewEdge;3147 if(NbOfNewPoints) { //3148 VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];3149 NbVertexOnBThEdge =NbOfNewPoints;3150 VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints];3151 NbVerticesOnGeomEdgex = NbOfNewPoints; }3152 NbOfNewPoints =0;3153 NbOfNewEdge = 0;3154 }3155 } // for(step;;)3156 assert(nbe);3157 3158 delete [] bcurve;3159 3160 3161 Insert();3162 ForceBoundary();3163 FindSubDomain();3164 3165 NewPoints(BTh,KeepBackVertices) ;3166 CurrentTh = 0;3167 }3168 3169 void Triangles::GeomToTriangles0(Int4 inbvx)3170 {3171 Gh.NbRef++;// add a ref to GH3172 3173 3174 Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve;3175 Real8 lcurve, lstep,s;3176 3177 R2 AB;3178 GeometricalVertex *a,*b;3179 Vertex *va,*vb;3180 GeometricalEdge * e;3181 PreInit(inbvx);3182 int background = &BTh != this;3183 // int SameGeom = background && (&BTh.Gh == &Gh);3184 nbv = 0;3185 NbVerticesOnGeomVertex=0;3186 NbVerticesOnGeomEdge=0;3187 for (i=0;i<Gh.nbv;i++)3188 if (Gh[i].Required() && Gh[i].IsThe() ) NbVerticesOnGeomVertex++;3189 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];3190 //3191 if( NbVerticesOnGeomVertex >= nbvx)3192 {3193 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;3194 MeshError(1,this);3195 }3196 for (i=0;i<Gh.nbv;i++)3197 if (Gh[i].Required()&& Gh[i].IsThe() ) {//Gh vertices Required3198 if (nbv < nbvx)3199 vertices[nbv] = Gh[i];3200 Gh[i].to = vertices + nbv;// save Geom -> Th3201 VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);3202 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;3203 nbv++;3204 }3205 // assert( Gh.nbv < nbvx);3206 3207 // Method in 2 step: 0 and 13208 // 1) compute de nb of edge3209 // 2) construct the edge3210 // generation of the curves3211 assert(! edges);3212 // 2 step3213 // --step=0 to compute the number of edges + alloc at end3214 // --step=1 to construct the edges3215 for (int step=0;step<2;step++)3216 {// for (int step=0;step<2;step++)3217 Int4 nbex = 0;3218 nbe = 0;3219 Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;3220 // cout << " -------------- step =" << step << endl;3221 Gh.UnMarkEdges();3222 NbOfCurves = 0;3223 for (i=0;i<Gh.nbe;i++) {3224 GeometricalEdge & ei = Gh.edges[i];3225 if (!ei.Dup()) // a good curve (not dup )3226 for(int j=0;j<2;j++)3227 if (!ei.Mark() && ei[j].Required()) {3228 // warning ei.Mark() can be change in loop for(j=0;j<2;j++)3229 // cout << " New curve = " << NbOfCurves << endl;3230 Int4 nbvend =0;3231 3232 Edge * PreviousNewEdge=0;3233 3234 lstep = -1;//to do not create points3235 if(ei.Required())3236 {3237 if (j==0)3238 if(step==0)3239 nbe++;3240 else3241 {3242 e = & ei;3243 a=ei(0)->The();3244 b=ei(1)->The();3245 assert(edges);3246 edges[nbe].v[0]=a->to;3247 edges[nbe].v[1]=b->to;;3248 edges[nbe].ref = e->ref;3249 edges[nbe].on = e;3250 edges[nbe].adj[0] = 0;3251 edges[nbe].adj[1] = 0;3252 nbe++;}3253 }3254 else3255 { // on curve ------3256 for ( int kstep=0;kstep<= step;kstep++)3257 { // begin for ( int kstep=0;kstep<= step;kstep++)3258 // if 2nd step where 2 step3259 // -- 1 compute le length of the curve3260 // -- create the points and edge3261 PreviousNewEdge=0;3262 NbNewPoints=0;3263 NbEdgeCurve=0;3264 assert(nbvend < nbvx);3265 lcurve =0;3266 s = lstep;3267 int k=j;3268 e = & ei;3269 a=ei(k)->The();3270 va = a->to;3271 e->SetMark();3272 // cout << " New curve " ;3273 3274 // if SameGeo We have go in the background geometry3275 // to find the discretisation of the curve3276 3277 for(;;)3278 {3279 k = 1-k;3280 b= (*e)(k)->The();3281 AB = b->r - a->r;3282 Metric MA = background ? BTh.MetricAt(a->r) :a->m ;3283 Metric MB = background ? BTh.MetricAt(b->r) :b->m ;3284 Real8 ledge = (MA(AB) + MB(AB))/2;3285 //3286 const int MaxSubEdge = 10;3287 int NbSubEdge = 1;3288 Real8 lSubEdge[MaxSubEdge];3289 R2 A,B;3290 if (ledge < 1.5)3291 lSubEdge[0] = ledge;3292 else {3293 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));3294 A= a->r;3295 Metric MAs =MA,MBs;3296 // cout << " lSubEdge old=" << ledge3297 // << " new " << A << MA << endl;3298 ledge = 0;3299 Real8 x =0, xstep= 1. / NbSubEdge;3300 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {3301 x += xstep;3302 B = e->F(k ? x : 1-x);3303 MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB);3304 AB = A-B;3305 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);3306 // cout << " " << lSubEdge[kk] << " x " << x3307 // << " " << A << B << MA << MB<< endl ;3308 }3309 // cout << endl;3310 }3311 3312 Real8 lcurveb = lcurve+ ledge ;3313 while (lcurve<=s && s <= lcurveb && nbv < nbvend)3314 {3315 // New points3316 3317 // Real8 aa=(lcurveb-s)/ledge;3318 // Real8 bb=(s-lcurve)/ledge;3319 3320 Real8 ss = s-lcurve;3321 // 1) find the SubEdge containing ss by dichotomie3322 int kk0=-1,kk1=NbSubEdge-1,kkk;3323 Real8 ll0=0,ll1=ledge,llk;3324 while (kk1-kk0>1)3325 {3326 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))3327 kk1=kkk,ll1=llk;3328 else3329 kk0=kkk,ll0=llk;}3330 assert(kk1 != kk0);3331 3332 Real8 sbb = (ss-ll0 )/(ll1-ll0);3333 Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb;3334 3335 // new vertex on edge3336 vb = &vertices[nbv++];3337 vb->m = Metric(aa,a->m,bb,b->m);3338 vb->ReferenceNumber = e->ref;3339 vb->DirOfSearch =NoDirOfSearch;3340 Real8 abcisse = k ? bb : aa;3341 vb->r = e->F( abcisse );3342 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);3343 3344 // to take in account the sens of the edge3345 3346 s += lstep;3347 edges[nbe].v[0]=va;3348 edges[nbe].v[1]=vb;3349 edges[nbe].ref = e->ref;3350 edges[nbe].on = e;3351 edges[nbe].adj[0] = PreviousNewEdge;3352 if(PreviousNewEdge)3353 PreviousNewEdge->adj[1] = &edges[nbe];3354 PreviousNewEdge = edges + nbe;3355 nbe++;3356 va = vb;3357 }3358 lcurve = lcurveb;3359 e->SetMark();3360 // cout << e-Gh.edges << ", " << k << " "3361 // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" "3362 // << lcurve<< ";; " ;3363 a=b;3364 if (b->Required() ) break;3365 int kprev=k;3366 k = e->SensAdj[kprev];// next vertices3367 e = e->Adj[kprev];3368 assert(e);3369 }// for(;;)3370 vb = b->to;3371 // cout << endl;3372 NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1);3373 NbNewPoints = NbEdgeCurve-1;3374 if(!kstep)3375 { NbVerticesOnGeomEdge0 += NbNewPoints;3376 NbOfCurves++;}3377 3378 nbvend=nbv+NbNewPoints;3379 3380 lstep = lcurve / NbEdgeCurve;3381 // cout <<"lstep " << lstep << " lcurve "3382 // << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " <<NbVerticesOnGeomEdge0<<" " <<NbVerticesOnGeomEdge<<" step =" <<step<< endl;3383 }3384 // end of curve --3385 if (edges) { // last edges of the curves3386 edges[nbe].v[0]=va;3387 edges[nbe].v[1]=vb;3388 edges[nbe].ref = e->ref;3389 edges[nbe].on = e;3390 edges[nbe].adj[0] = PreviousNewEdge;3391 edges[nbe].adj[1] = 0;3392 if(PreviousNewEdge)3393 PreviousNewEdge->adj[1] = & edges[nbe];3394 3395 nbe++;}3396 else3397 nbe += NbEdgeCurve;3398 } // end on curve ---3399 } // if (edges[i][j].Corner())3400 } // for (i=0;i<nbe;i++)3401 if(!step) {3402 // cout << "edges " << edges << " VerticesOnGeomEdge " <<VerticesOnGeomEdge << endl;3403 assert(!edges);3404 assert(!VerticesOnGeomEdge);3405 edges = new Edge[nbex=nbe];3406 if(NbVerticesOnGeomEdge0)3407 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];3408 assert(edges);3409 assert(VerticesOnGeomEdge || NbVerticesOnGeomEdge0 ==0);3410 // do the vertex on a geometrical vertex3411 NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;3412 }3413 else3414 assert(NbVerticesOnGeomEdge == NbVerticesOnGeomEdge0);3415 // cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe3416 // << "== " << nbex << " nbv = " << nbv << endl;3417 assert(nbex=nbe);3418 } // for (step=0;step<2;step++)3419 3420 Insert();3421 ForceBoundary();3422 FindSubDomain();3423 3424 // NewPointsOld(*this) ;3425 NewPoints(*this,0) ;3426 CurrentTh = 0;3427 }3428 3429 Edge** Triangles::MakeGeometricalEdgeToEdge()3430 {3431 assert(Gh.nbe);3432 Edge **e= new (Edge* [Gh.nbe]);3433 3434 Int4 i;3435 for ( i=0;i<Gh.nbe ; i++)3436 e[i]=NULL;3437 for ( i=0;i<nbe ; i++)3438 {3439 Edge * ei = edges+i;3440 GeometricalEdge *on = ei->on;3441 e[Gh.Number(on)] = ei;3442 }3443 for ( i=0;i<nbe ; i++)3444 for (int ii=0;ii<2;ii++) {3445 Edge * ei = edges+i;3446 GeometricalEdge *on = ei->on;3447 int j= ii;3448 while (!(*on)[j].Required()) {3449 // cout << i << " " << ii << " j= " << j << " curve = "3450 // << on->CurveNumber << " " << Gh.Number(on) << " on " << j3451 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1])3452 // << " -> " ;3453 Adj(on,j); // next geom edge3454 j=1-j;3455 // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)]3456 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl;3457 if (e[Gh.Number(on)]) break; // optimisation3458 e[Gh.Number(on)] = ei;3459 }3460 }3461 3462 int kk=0;3463 for ( i=0;i<Gh.nbe ; i++)3464 if (!e[i])3465 if(kk++<10) {3466 cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber3467 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl;3468 // assert( e[i]);3469 }3470 if(kk) MeshError(997,this);3471 3472 return e;3473 }3474 3475 Triangles::~Triangles() {3476 long int verbosity=2;3477 //assert(NbRef<=0);3478 if (CurrentTh == this) CurrentTh=0;3479 //if(vertices) delete [] vertices; //TEST crash if not commented3480 if(edges) delete [] edges;3481 if(triangles) delete [] triangles;3482 if(quadtree) delete quadtree;3483 //if(ordre) delete [] ordre; //TEST crash if not commented3484 if( subdomains) delete [] subdomains;3485 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;3486 if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex;3487 //if (name) delete [] name; //TEST crash if not commented3488 if (identity) delete [] identity;3489 if (VertexOnBThVertex) delete [] VertexOnBThVertex;3490 if (VertexOnBThEdge) delete [] VertexOnBThEdge;3491 3492 if (&Gh) {3493 if (Gh.NbRef>0) Gh.NbRef--;3494 else if (Gh.NbRef==0) delete &Gh;3495 }3496 if (&BTh && (&BTh != this)) {3497 if (BTh.NbRef>0) BTh.NbRef--;3498 else if (BTh.NbRef==0) delete &BTh;3499 }3500 PreInit(0); // set all to zero3501 }3502 3503 void Triangles::SetIntCoor(const char * strfrom)3504 {3505 pmin = vertices[0].r;3506 pmax = vertices[0].r;3507 3508 // recherche des extrema des vertices pmin,pmax3509 Int4 i;3510 for (i=0;i<nbv;i++) {3511 pmin.x = Min(pmin.x,vertices[i].r.x);3512 pmin.y = Min(pmin.y,vertices[i].r.y);3513 pmax.x = Max(pmax.x,vertices[i].r.x);3514 pmax.y = Max(pmax.y,vertices[i].r.y);3515 }3516 R2 DD = (pmax-pmin)*0.05;3517 pmin = pmin-DD;3518 pmax = pmax+DD;3519 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));3520 assert(coefIcoor >0);3521 3522 // generation of integer coord3523 for (i=0;i<nbv;i++) {3524 vertices[i].i = toI2(vertices[i].r);3525 }3526 3527 // computation of the det3528 int Nberr=0;3529 for (i=0;i<nbt;i++)3530 {3531 Vertex & v0 = triangles[i][0];3532 Vertex & v1 = triangles[i][1];3533 Vertex & v2 = triangles[i][2];3534 if ( &v0 && &v1 && &v2 ) // a good triangles;3535 {3536 triangles[i].det= det(v0,v1,v2);3537 if (triangles[i].det <=0 && Nberr++ <10)3538 {3539 if(Nberr==1)3540 if (strfrom)3541 cerr << "+++ Fatal Error " << strfrom << "(SetInCoor) Error : area of Triangle < 0 " << endl;3542 else3543 cerr << "+++ Fatal Error Triangle (in SetInCoor) area of Triangle < 0" << endl;3544 cerr << " Triangle " << i << " det (I2) = " << triangles[i].det ;3545 cerr << " (R2) " << Det(v1.r-v0.r,v2.r-v0.r);3546 cerr << "; The 3 vertices " << endl;3547 cerr << Number(v0) << " " << Number(v1) << " "3548 << Number(v2) << " : " ;3549 cerr << v0.r << v1.r << v2.r << " ; ";3550 cerr << v0.i << v1.i << v2.i << endl;3551 }3552 }3553 else3554 triangles[i].det= -1; // Null triangle;3555 }3556 if (Nberr) MeshError(899,this);3557 3558 }3559 3560 void Triangles::FillHoleInMesh() {3561 3562 Triangles* OldCurrentTh =CurrentTh;3563 CurrentTh=this;3564 3565 int verbosity=0;3566 3567 // generation of the integer coordinate3568 {3569 3570 // find extrema coordinates of vertices pmin,pmax3571 Int4 i;3572 if(verbosity>2) printf(" Filling holes in mesh of %i vertices\n",nbv);3573 3574 //initialize ordre3575 assert(ordre);3576 for (i=0;i<nbv;i++) ordre[i]=0;3577 3578 NbSubDomains =0;3579 3580 /* generation of the adjacence of the triangles*/3581 3582 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);3583 3584 //initialize st3585 Int4* st = new Int4[nbt*3];3586 for (i=0;i<nbt*3;i++) st[i]=-1;3587 3588 //check number of edges3589 Int4 kk =0;3590 for (i=0;i<nbe;i++){3591 kk=kk+(i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));3592 }3593 if (kk != nbe) {3594 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe));3595 }3596 3597 //3598 for (i=0;i<nbt;i++){3599 for (int j=0;j<3;j++) {3600 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),3601 Number(triangles[i][VerticesOfTriangularEdge[j][1]]));3602 Int4 invisible = triangles[i].Hidden(j);3603 if(st[k]==-1){3604 st[k]=3*i+j;3605 }3606 else if(st[k]>=0) {3607 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)));3608 3609 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));3610 if (invisible) triangles[i].SetHidden(j);3611 if (k<nbe) {3612 triangles[i].SetLocked(j);3613 }3614 st[k]=-2-st[k];3615 }3616 else {3617 throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",3618 Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));3619 }3620 }3621 }3622 if(verbosity>5) {3623 printf(" info on Mesh %s:\n",name);3624 printf(" - number of vertices = %i \n",nbv);3625 printf(" - number of triangles = %i \n",nbt);3626 printf(" - number of given edges = %i \n",nbe);3627 printf(" - number of all edges = %i \n" ,edge4->nb());3628 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv);3629 }3630 3631 // check the consistant of edge[].adj and the geometrical required vertex3632 Int4 k=0;3633 for (i=0;i<edge4->nb();i++){3634 if (st[i] >=0){ // edge alone3635 if (i < nbe) {3636 Int4 i0=edge4->i(i);3637 ordre[i0] = vertices+i0;3638 Int4 i1=edge4->j(i);3639 ordre[i1] = vertices+i1;3640 }3641 else {3642 k=k+1;3643 if (k <20) {3644 throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));3645 }3646 }3647 }3648 }3649 if(k != 0) {3650 throw ErrorException(__FUNCT__,exprintf("%i boundary edges are not defined as edges",k));3651 }3652 3653 /* mesh generation with boundary points*/3654 Int4 nbvb = 0;3655 for (i=0;i<nbv;i++){3656 vertices[i].t=0;3657 vertices[i].vint=0;3658 if (ordre[i]){3659 ordre[nbvb++] = ordre[i];3660 }3661 }3662 3663 Triangle* savetriangles= triangles;3664 Int4 savenbt=nbt;3665 Int4 savenbtx=nbtx;3666 SubDomain * savesubdomains = subdomains;3667 subdomains = 0;3668 3669 Int4 Nbtriafillhole = 2*nbvb;3670 Triangle* triafillhole =new Triangle[Nbtriafillhole];3671 triangles = triafillhole;3672 3673 nbt=2;3674 nbtx= Nbtriafillhole;3675 3676 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)3677 if ( ++i >= nbvb) {3678 cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;3679 MeshError(998,this); }3680 Exchange( ordre[2], ordre[i]);3681 3682 Vertex * v0=ordre[0], *v1=ordre[1];3683 3684 3685 triangles[0](0) = 0; // sommet pour infini3686 triangles[0](1) = v0;3687 triangles[0](2) = v1;3688 3689 triangles[1](0) = 0;// sommet pour infini3690 triangles[1](2) = v0;3691 triangles[1](1) = v1;3692 const int e0 = OppositeEdge[0];3693 const int e1 = NextEdge[e0];3694 const int e2 = PreviousEdge[e0];3695 triangles[0].SetAdj2(e0, &triangles[1] ,e0);3696 triangles[0].SetAdj2(e1, &triangles[1] ,e2);3697 triangles[0].SetAdj2(e2, &triangles[1] ,e1);3698 3699 triangles[0].det = -1; // faux triangles3700 triangles[1].det = -1; // faux triangles3701 3702 triangles[0].SetTriangleContainingTheVertex();3703 triangles[1].SetTriangleContainingTheVertex();3704 3705 triangles[0].link=&triangles[1];3706 triangles[1].link=&triangles[0];3707 3708 if (!quadtree )3709 delete quadtree; // ->ReInitialise();3710 3711 quadtree = new QuadTree(this,0);3712 quadtree->Add(*v0);3713 quadtree->Add(*v1);3714 3715 // We add the vertices one by one3716 Int4 NbSwap=0;3717 for (Int4 icount=2; icount<nbvb; icount++) {3718 Vertex *vi = ordre[icount];3719 // cout << " Add vertex " << Number(vi) << endl;3720 Icoor2 dete[3];3721 Triangle *tcvi = FindTriangleContening(vi->i,dete);3722 quadtree->Add(*vi);3723 Add(*vi,tcvi,dete);3724 NbSwap += vi->Optim(1,1);3725 }3726 3727 // inforce the boundary3728 TriangleAdjacent ta(0,0);3729 Int4 nbloss = 0,knbe=0;3730 for ( i = 0; i < nbe; i++){3731 if (st[i] >=0){ // edge alone => on border ... FH oct 20093732 Vertex & a=edges[i][0], & b = edges[i][1];3733 if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.3734 {3735 knbe++;3736 if (ForceEdge(a,b,ta)<0)3737 nbloss++;3738 }3739 }3740 }3741 if(nbloss) {3742 cerr << " we loss some " << nbloss << " " << " edges other " << knbe << endl;3743 MeshError(1100,this);3744 }3745 3746 FindSubDomain(1);3747 // remove all the hole3748 // remove all the good sub domain3749 Int4 krm =0;3750 for (i=0;i<nbt;i++){3751 if (triangles[i].link){ // remove triangles3752 krm++;3753 for (int j=0;j<3;j++)3754 {3755 TriangleAdjacent ta = triangles[i].Adj(j);3756 Triangle & tta = * (Triangle *) ta;3757 if(! tta.link) // edge between remove and not remove3758 { // change the link of ta;3759 int ja = ta;3760 Vertex *v0= ta.EdgeVertex(0);3761 Vertex *v1= ta.EdgeVertex(1);3762 Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);3763 assert(st[k] >=0);3764 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));3765 ta.SetLock();3766 st[k]=-2-st[k];3767 }3768 }3769 }3770 }3771 Int4 NbTfillHoll =0;3772 for (i=0;i<nbt;i++){3773 if (triangles[i].link) {3774 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);3775 triangles[i].color=-1;3776 }3777 else3778 {3779 triangles[i].color= savenbt+ NbTfillHoll++;3780 }3781 }3782 // cout << savenbt+NbTfillHoll << " " << savenbtx << endl;3783 assert(savenbt+NbTfillHoll <= savenbtx );3784 // copy of the outside triangles in saveTriangles3785 for (i=0;i<nbt;i++){3786 if(triangles[i].color>=0) {3787 savetriangles[savenbt]=triangles[i];3788 savetriangles[savenbt].link=0;3789 savenbt++;3790 }3791 }3792 // gestion of the adj3793 k =0;3794 Triangle * tmax = triangles + nbt;3795 for (i=0;i<savenbt;i++)3796 {3797 Triangle & ti = savetriangles[i];3798 for (int j=0;j<3;j++)3799 {3800 Triangle * ta = ti.TriangleAdj(j);3801 int aa = ti.NuEdgeTriangleAdj(j);3802 int lck = ti.Locked(j);3803 if (!ta) k++; // bug3804 else if ( ta >= triangles && ta < tmax)3805 {3806 ta= savetriangles + ta->color;3807 ti.SetAdj2(j,ta,aa);3808 if(lck) ti.SetLocked(j);3809 }3810 }3811 }3812 // OutSidesTriangles = triangles;3813 // Int4 NbOutSidesTriangles = nbt;3814 3815 // restore triangles;3816 nbt=savenbt;3817 nbtx=savenbtx;3818 delete [] triangles;3819 delete [] subdomains;3820 triangles = savetriangles;3821 subdomains = savesubdomains;3822 // cout << triangles << " <> " << OutSidesTriangles << endl;3823 /* k=0;3824 for (i=0;i<nbt;i++)3825 for (int j=0;j<3;j++)3826 if (!triangles[i].TriangleAdj(j))3827 k++;3828 */3829 if (k) {3830 cerr << "Error Nb of triangles edge alone = " << k << endl;3831 MeshError(9997,this);3832 }3833 FindSubDomain();3834 // cout << " NbTOld = " << NbTold << " == " << nbt - NbOutT << " " << nbt << endl;3835 3836 delete edge4;3837 delete [] st;3838 for (i=0;i<nbv;i++)3839 quadtree->Add(vertices[i]);3840 3841 SetVertexFieldOn();3842 3843 for (i=0;i<nbe;i++)3844 if(edges[i].on)3845 for(int j=0;j<2;j++)3846 if (!edges[i].adj[j])3847 if(!edges[i][j].on->IsRequiredVertex()) {3848 cerr << " Erreur adj et sommet requis edges [" << i << "][ " << j << "]= "3849 << Number(edges[i][j]) << " : " << " on = " << Gh.Number(edges[i].on) ;3850 if (edges[i][j].on->OnGeomVertex())3851 cerr << " vertex " << Gh.Number(edges[i][j].on->gv);3852 else if (edges[i][j].on->OnGeomEdge())3853 cerr << "Edges " << Gh.Number(edges[i][j].on->ge);3854 else3855 cerr << " = " << edges[i][j].on ;3856 cerr << endl;3857 }3858 }3859 CurrentTh=OldCurrentTh;3860 }3861 3862 Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) // COPY OPERATOR3863 : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this))3864 {3865 Gh.NbRef++;3866 nbvxx = Max(nbvxx,Th.nbv);3867 Int4 i;3868 // do all the allocation to be sure all the pointer existe3869 3870 char * cname = 0;3871 if (Th.name)3872 {3873 cname = new char[strlen(Th.name)+1];3874 strcpy(cname,Th.name);3875 }3876 PreInit(nbvxx,cname);// to make the allocation3877 // copy of triangles3878 nt=Th.nt;3879 nbv = Th.nbv;3880 nbt = Th.nbt;3881 nbiv = Th.nbiv;3882 nbe = Th.nbe;3883 NbSubDomains = Th.NbSubDomains;3884 NbOutT = Th.NbOutT;3885 NbOfQuad = Th.NbOfQuad ;3886 NbOfSwapTriangle =0;3887 NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex;3888 if(NbVerticesOnGeomVertex)3889 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];3890 NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge;3891 if (NbVerticesOnGeomEdge)3892 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ;3893 if (& BTh == & Th.BTh) // same back ground3894 {3895 BTh.NbRef++;3896 NbVertexOnBThVertex = Th.NbVertexOnBThVertex;3897 if(NbVertexOnBThVertex)3898 VertexOnBThVertex = new VertexOnVertex[NbVertexOnBThVertex];3899 NbVertexOnBThEdge = Th.NbVertexOnBThEdge;3900 if(NbVertexOnBThEdge)3901 VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge];3902 }3903 else3904 { // no add on back ground mesh3905 BTh.NbRef++;3906 NbVertexOnBThVertex=0;3907 VertexOnBThVertex=0;3908 NbVertexOnBThEdge=0;3909 VertexOnBThEdge=0;3910 // assert (& BTh == this); // --- a voir3911 3912 }3913 3914 3915 if(nbe)3916 edges = new Edge[nbe];3917 if(NbSubDomains)3918 subdomains = new SubDomain[NbSubDomains];3919 pmin = Th.pmin;3920 pmax = Th.pmax;3921 coefIcoor = Th.coefIcoor;3922 for(i=0;i<nbt;i++)3923 triangles[i].Set(Th.triangles[i],Th,*this);3924 for(i=0;i<nbe;i++)3925 edges[i].Set(Th,i,*this);3926 for(i=0;i<nbv;i++)3927 vertices[i].Set(Th.vertices[i],Th,*this);3928 for(i=0;i<NbSubDomains;i++)3929 subdomains[i].Set(Th,i,*this);3930 for (i=0;i<NbVerticesOnGeomVertex;i++)3931 VerticesOnGeomVertex[i].Set(Th.VerticesOnGeomVertex[i],Th,*this);3932 for (i=0;i<NbVerticesOnGeomEdge;i++)3933 VerticesOnGeomEdge[i].Set(Th.VerticesOnGeomEdge[i],Th,*this);3934 quadtree=0;3935 3936 3937 // assert(!OutSidesTriangles);3938 }3939 3940 Int4 Triangle::Optim(Int2 i,int koption)3941 {3942 // turne around in positif sens3943 Int4 NbSwap =0;3944 Triangle *t = this;3945 int k=0,j =OppositeEdge[i];3946 int jp = PreviousEdge[j];3947 // initialise tp, jp the previous triangle & edge3948 Triangle *tp= at[jp];3949 jp = aa[jp]&3;3950 do {3951 // cout << *t << " " << j << "\n\t try swap " ;3952 while (t->swap(j,koption))3953 {3954 NbSwap++;3955 assert(k++<20000);3956 t= tp->at[jp]; // set unchange t qnd j for previous triangles3957 j= NextEdge[tp->aa[jp]&3];3958 // cout << "\n\t s " << *t << " " << j << endl;3959 }3960 // end on this Triangle3961 tp = t;3962 jp = NextEdge[j];3963 3964 t= tp->at[jp]; // set unchange t qnd j for previous triangles3965 j= NextEdge[tp->aa[jp]&3];3966 3967 } while( t != this);3968 return NbSwap;3969 }3970 3971 void Triangles::SmoothingVertex(int nbiter,Real8 omega )3972 {3973 long int verbosity=0;3974 // if quatree exist remove it end reconstruct3975 if (quadtree) delete quadtree;3976 quadtree=0;3977 ReMakeTriangleContainingTheVertex();3978 Triangle vide; // a triangle to mark the boundary vertex3979 Triangle ** tstart= new Triangle* [nbv];3980 Int4 i,j,k;3981 // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide3982 if ( this == & BTh)3983 for ( i=0;i<nbv;i++)3984 tstart[i]=vertices[i].t;3985 else3986 for ( i=0;i<nbv;i++)3987 tstart[i]=0;3988 for ( j=0;j<NbVerticesOnGeomVertex;j++ )3989 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;3990 for ( k=0;k<NbVerticesOnGeomEdge;k++ )3991 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;3992 if(verbosity>2)3993 cout << " -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl;3994 for (k=0;k<nbiter;k++)3995 {3996 Int4 i,NbSwap =0;3997 Real8 delta =0;3998 for ( i=0;i<nbv;i++)3999 if (tstart[i] != &vide) // not a boundary vertex4000 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));4001 if (!NbOfQuad)4002 for ( i=0;i<nbv;i++)4003 if (tstart[i] != &vide) // not a boundary vertex4004 NbSwap += vertices[i].Optim(1);4005 if (verbosity>3)4006 cout << " Move max = " << sqrt(delta) << " iteration = "4007 << k << " Nb of Swap = " << NbSwap << endl;4008 }4009 4010 delete [] tstart;4011 if (quadtree) quadtree= new QuadTree(this);4012 }4013 void Triangles::MakeQuadTree()4014 {4015 long int verbosity=0;4016 if(verbosity>8)4017 cout << " MakeQuadTree" << endl;4018 if ( !quadtree ) quadtree = new QuadTree(this);4019 4020 }4021 void Triangles::ShowRegulaty() const// Add FH avril 20074022 {4023 const Real8 sqrt32=sqrt(3.)*0.5;4024 const Real8 aireKh=sqrt32*0.5;4025 D2 Beq(1,0),Heq(0.5,sqrt32);4026 D2xD2 Br(D2xD2(Beq,Heq).t());4027 D2xD2 B1r(Br.inv());4028 /* D2xD2 BB = Br.t()*Br;4029 cout << " BB = " << BB << " " << Br*B1r << endl;4030 MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y);4031 MatVVP2x2 VMM(MMM);4032 cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl;4033 */4034 double gammamn=1e100,hmin=1e100;4035 double gammamx=0,hmax=0;4036 double beta=1e100;4037 double beta0=0;4038 double alpha2=0;4039 double area=0,Marea=0;4040 // Real8 cf= Real8(coefIcoor);4041 // Real8 cf2= 6.*cf*cf;4042 int nt=0;4043 for (int it=0;it<nbt;it++)4044 if ( triangles[it].link)4045 {4046 nt++;4047 Triangle &K=triangles[it];4048 Real8 area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;4049 area+= area3;4050 D2xD2 B_Kt(K[0],K[1],K[2]);4051 D2xD2 B_K(B_Kt.t());4052 D2xD2 B1K = Br*B_K.inv();4053 D2xD2 BK = B_K*B1r;4054 D2xD2 B1B1 = B1K.t()*B1K;4055 MetricAnIso MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);4056 MatVVP2x2 VMK(MK);4057 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));4058 // cout << B_K << " * " << B1r << " == " << BK << " " << B_K*B_K.inv() << endl;4059 Real8 betaK=0;4060 4061 for (int j=0;j<3;j++)4062 {4063 Real8 he= Norme2(R2(K[j],K[(j+1)%3]));4064 hmin=Min(hmin,he);4065 hmax=Max(hmax,he);4066 Vertex & v=K[j];4067 D2xD2 M((MetricAnIso)v);4068 betaK += sqrt(M.det());4069 4070 D2xD2 BMB = BK.t()*M*BK;4071 MetricAnIso M1(BMB.x.x,BMB.x.y,BMB.y.y);4072 MatVVP2x2 VM1(M1);4073 //cout << B_K <<" " << M << " " << he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<< endl;4074 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);4075 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);4076 }4077 betaK *= area3;// 1/2 (somme sqrt(det))* area(K)4078 Marea+= betaK;4079 // cout << betaK << " " << area3 << " " << beta << " " << beta0 << " " << area3*3*3*3 <<endl;4080 beta=min(beta,betaK);4081 beta0=max(beta0,betaK);4082 }4083 area*=3;4084 gammamn=sqrt(gammamn);4085 gammamx=sqrt(gammamx);4086 cout << " -- adaptmesh Regulary: Nb triangles " << nt << " , h min " << hmin << " , h max " << hmax << endl;4087 cout << " area = " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl;4088 cout << " infiny-regulaty: min " << gammamn << " max " << gammamx << endl;4089 cout << " anisomax "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh)4090 << " min "<< 1./sqrt(beta0/aireKh) << endl;4091 }4092 void Triangles::ShowHistogram() const4093 {4094 4095 const Int4 kmax=10;4096 const Real8 llmin = 0.5,llmax=2;4097 const Real8 lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);4098 Int4 histo[kmax+1];4099 Int4 i,it,k, nbedges =0;4100 for (i=0;i<=kmax;i++) histo[i]=0;4101 for (it=0;it<nbt;it++)4102 if ( triangles[it].link)4103 {4104 4105 for (int j=0;j<3;j++)4106 {4107 Triangle *ta = triangles[it].TriangleAdj(j);4108 if ( !ta || !ta->link || Number(ta) >= it)4109 {4110 Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];4111 Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];4112 if ( !& vP || !&vQ) continue;4113 R2 PQ = vQ.r - vP.r;4114 Real8 l = log(LengthInterpole(vP,vQ,PQ));4115 nbedges++;4116 k = (int) ((l - lmin)*delta);4117 k = Min(Max(k,0L),kmax);4118 histo[k]++;4119 }4120 }4121 }4122 cout << " -- Histogram of the unit mesh, nb of edges" << nbedges << endl <<endl;4123 4124 cout << " length of edge in | % of edge | Nb of edges " << endl;4125 cout << " ------------------- | ---------- | ----------- " << endl;4126 for (i=0;i<=kmax;i++)4127 {4128 cout << " " ;4129 cout.width(10);4130 if (i==0) cout << " 0 " ;4131 else cout << exp(lmin+i/delta) ;4132 cout.width(); cout << "," ;4133 cout.width(10);4134 if (i==kmax) cout << " +infty " ;4135 else cout << exp(lmin+(i+1)/delta) ;4136 cout.width();cout << " | " ;4137 4138 cout.precision(4);4139 cout.width(6);4140 cout << ((long) ((10000.0 * histo[i])/ nbedges))/100.0 ;4141 cout.width();4142 cout.precision();4143 cout << " | " << histo[i] <<endl;4144 }4145 cout << " ------------------- | ---------- | ----------- " << endl <<endl;4146 4147 }4148 4149 int Triangles::Crack()4150 {4151 assert(NbCrackedEdges ==0 || NbCrackedVertices >0);4152 for (int i=0;i<NbCrackedEdges;i++)4153 CrackedEdges[i].Crack();4154 return NbCrackedEdges;4155 }4156 4157 int Triangles::UnCrack()4158 {4159 assert(NbCrackedEdges ==0 || NbCrackedVertices >0);4160 for (int i=0;i<NbCrackedEdges;i++)4161 CrackedEdges[i].UnCrack();4162 return NbCrackedEdges;4163 }4164 4165 int Triangles::CrackMesh()4166 {4167 long int verbosity=0;4168 Triangles *CurrentThOld = CurrentTh;4169 // computed the number of cracked edge4170 int i,k;4171 for (k=i=0;i<nbe;i++)4172 if(edges[i].on->Cracked()) k++;4173 if( k==0) return 0;4174 CurrentTh = this;4175 cout << " Nb of Cracked Edges = " << k << endl;4176 NbCrackedEdges =k;4177 CrackedEdges = new CrackedEdge[k];4178 // new edge4179 Edge * e = new Edge[ nbe + k];4180 4181 // copy4182 for (i=0;i<nbe;i++)4183 e[i] = edges[i];4184 delete edges;4185 edges = e;4186 4187 const int nbe0 = nbe;4188 for (k=i=0;i<nbe0;i++) // on double les arete cracked4189 if(edges[i].on->Cracked())4190 {4191 e[nbe] = e[i];4192 // return the edge4193 e[nbe].v[0] = e[i].v[1];4194 e[nbe].v[1] = e[i].v[0];4195 e[nbe].on = e[i].on->link ; // fqux4196 CrackedEdges[k++]=CrackedEdge(edges,i,nbe);4197 nbe++;4198 }4199 ReMakeTriangleContainingTheVertex() ;4200 //4201 int nbcrakev =0;4202 Vertex *vlast = vertices + nbv;4203 Vertex *vend = vertices + nbvx; // end of array4204 for (int iv=0;iv<nbv;iv++) // vertex4205 {4206 Vertex & v= vertices[iv];4207 Vertex * vv = & v;4208 int kk=0; // nb cracked4209 int kc=0;4210 int kkk =0; // nb triangle with same number4211 Triangle * tbegin = v.t;4212 int i = v.vint;4213 assert(tbegin && (i >= 0 ) && (i <3));4214 // turn around the vertex v4215 TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge4216 int k=0;4217 do {4218 int kv = VerticesOfTriangularEdge[ta][1];4219 k++;4220 Triangle * tt (ta);4221 if ( ta.Cracked() )4222 {4223 TriangleAdjacent tta=(ta.Adj());4224 assert(tta.Cracked());4225 if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart4226 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use4227 kk++;// number of cracked edge view4228 }4229 if ( tt->link ) { // if good triangles store the value4230 int it = Number(tt);4231 assert(it < nt);4232 (*tt)(kv)= vv; // Change the vertex of triangle4233 if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref4234 // tt->SetTriangleContainingTheVertex();4235 kkk++;4236 } else if (kk) { // crack + boundary4237 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use4238 }4239 4240 ta = Next(ta).Adj();4241 } while ( (tbegin != ta));4242 assert(k);4243 if (kc) nbcrakev++;4244 }4245 4246 if ( nbcrakev )4247 for (int iec =0;iec < NbCrackedEdges; iec ++)4248 CrackedEdges[iec].Set();4249 4250 // set the ref4251 cout << " set the ref " << endl ;4252 NbCrackedVertices = nbcrakev;4253 // int nbvo = nbv;4254 nbv = vlast - vertices;4255 int nbnewv = nbv - nbv; // nb of new vrtices4256 if (nbcrakev && verbosity > 1 )4257 cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nbnewv<< endl;4258 // all the new vertices are on geometry4259 // BOFBO-- A VOIR4260 if (nbnewv)4261 { //4262 Int4 n = nbnewv+NbVerticesOnGeomVertex;4263 Int4 i,j,k;4264 VertexOnGeom * vog = new VertexOnGeom[n];4265 for ( i =0; i<NbVerticesOnGeomVertex;i++)4266 vog[i]=VerticesOnGeomVertex[i];4267 delete [] VerticesOnGeomVertex;4268 VerticesOnGeomVertex = vog;4269 // loop on cracked edge4270 Vertex * LastOld = vertices + nbv - nbnewv;4271 for (int iec =0;iec < NbCrackedEdges; iec ++)4272 for (k=0;k<2;k++)4273 {4274 Edge & e = *( k ? CrackedEdges[iec].a.edge : CrackedEdges[iec].b.edge);4275 for (j=0;j<2;j++)4276 {4277 Vertex * v = e(j);4278 if ( v >= LastOld)4279 { // a new vertex4280 Int4 old = v->ReferenceNumber ; // the old same vertex4281 Int4 i = ( v - LastOld);4282 // if the old is on vertex => warning4283 // else the old is on edge => ok4284 vog[i] = vog[old];4285 // vog[i].mv = v;4286 //g[i].ge = ;4287 //og[i].abcisse = ;4288 }4289 4290 }4291 }4292 4293 NbVerticesOnGeomVertex = n;4294 }4295 SetVertexFieldOn();4296 4297 4298 if (vlast >= vend)4299 {4300 cerr << " Not enougth vertices to crack the mesh we need " << nbv << " vertices " << endl;4301 MeshError(555,this);4302 }4303 cout << " NbCrackedVertices " << NbCrackedVertices << endl;4304 CurrentTh = CurrentThOld;4305 return NbCrackedVertices;4306 }4307 4308 Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)4309 : Gh(*(new Geometry())), BTh(*this)4310 { // truncature4311 //4312 4313 char cname[] = "trunc";4314 4315 int i,k,itadj;4316 int kt=0;4317 int * kk = new int [Tho.nbv];4318 Int4 * reft = new Int4[Tho.nbt];4319 Int4 nbInT = Tho.ConsRefTriangle(reft);4320 Int4 * refv = new Int4[Tho.nbv];4321 4322 for (i=0;i<Tho.nbv;i++)4323 kk[i]=-1;4324 for (i=0;i<Tho.nbv;i++)4325 refv[i]=0;4326 int nbNewBedge =0;4327 // int nbOldBedge =0;4328 for (i=0;i<Tho.nbt;i++)4329 if( reft[i] >=0 && flag[i])4330 {4331 const Triangle & t = Tho.triangles[i];4332 kt++;4333 kk[Tho.Number(t[0])]=1;4334 kk[Tho.Number(t[1])]=1;4335 kk[Tho.Number(t[2])]=1;4336 itadj=Tho.Number(t.TriangleAdj(0));4337 if ( reft[itadj] >=0 && !flag[itadj])4338 { nbNewBedge++;4339 refv[Tho.Number(t[VerticesOfTriangularEdge[0][0]])]=bb[i];4340 refv[Tho.Number(t[VerticesOfTriangularEdge[0][1]])]=bb[i];4341 }4342 itadj=Tho.Number(t.TriangleAdj(1));4343 if ( reft[itadj] >=0 && !flag[itadj])4344 { nbNewBedge++;4345 refv[Tho.Number(t[VerticesOfTriangularEdge[1][0]])]=bb[i];4346 refv[Tho.Number(t[VerticesOfTriangularEdge[1][1]])]=bb[i];}4347 itadj=Tho.Number(t.TriangleAdj(2));4348 if ( reft[itadj] >=0 && !flag[itadj])4349 { nbNewBedge++;4350 refv[Tho.Number(t[VerticesOfTriangularEdge[2][0]])]=bb[i];4351 refv[Tho.Number(t[VerticesOfTriangularEdge[2][1]])]=bb[i];}4352 }4353 k=0;4354 for (i=0;i<Tho.nbv;i++)4355 if (kk[i]>=0)4356 kk[i]=k++;4357 cout << " number of vertices " << k << " remove = " << Tho.nbv - k << endl;4358 cout << " number of triangles " << kt << " remove = " << nbInT-kt << endl;4359 cout << " number of New boundary edge " << nbNewBedge << endl;4360 Int4 inbvx =k;4361 PreInit(inbvx,cname);4362 for (i=0;i<Tho.nbv;i++)4363 if (kk[i]>=0)4364 {4365 vertices[nbv] = Tho.vertices[i];4366 if (!vertices[nbv].ref())4367 vertices[nbv].ReferenceNumber = refv[i];4368 nbv++;4369 }4370 assert(inbvx == nbv);4371 for (i=0;i<Tho.nbt;i++)4372 if( reft[i] >=0 && flag[i])4373 {4374 const Triangle & t = Tho.triangles[i];4375 int i0 = Tho.Number(t[0]);4376 int i1 = Tho.Number(t[1]);4377 int i2 = Tho.Number(t[2]);4378 assert(i0>=0 && i1 >= 0 && i2 >= 0);4379 assert(i0<Tho.nbv && i1 <Tho.nbv && i2 <Tho.nbv);4380 // cout <<i<< " F" << flag[i] << " T " << nbt << " = " << kk[i0] << " " << kk[i1] << " " << kk[i2] ;4381 // cout << " OT " << i0 << " " << i1 << " " << i2 << " " << reft[i] << endl;4382 triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);4383 triangles[nbt].color = Tho.subdomains[reft[i]].ref;4384 nbt++;4385 }4386 assert(kt==nbt);4387 if (nbt ==0 && nbv ==0) {4388 cout << "Error all triangles was remove " << endl;4389 MeshError(999,this);4390 }4391 delete [] kk;4392 delete [] reft;4393 delete [] refv;4394 double cutoffradian = 10.0/180.0*Pi;4395 ConsGeometry(cutoffradian);4396 Gh.AfterRead();4397 SetIntCoor();4398 FillHoleInMesh();4399 4400 assert(NbSubDomains);4401 assert(subdomains[0].head && subdomains[0].head->link);4402 4403 }4404 4405 Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const { // in: B4406 // out: t4407 // out : dete[3]4408 // t the triangle and s0,s1,s2 the 3 vertices of t4409 // in dete[3] = { det(B,s1,s2) , det(s0,B,s2), det(s0,s1,B)}4410 // with det(a,b,c ) = -1 if one of 3 vertices a,b,c is NULL4411 Triangle * t=0;4412 int j,jp,jn,jj;4413 if (tstart)4414 t=tstart;4415 else4416 {4417 assert(quadtree);4418 Vertex *a = quadtree->NearestVertex(B.x,B.y) ;4419 4420 if (! a || !a->t ) {4421 if (a)4422 {cerr << " Attention PB TriangleConteningTheVertex vertex number=" << Number(a) << endl;4423 cerr << "We forget a call to ReMakeTriangleContainingTheVertex" << endl;}4424 cerr << " Pb with " << B << toR2(B) << endl;4425 MeshError(7777);4426 }4427 assert(a>= vertices && a < vertices+nbv);4428 // int k=0;4429 t = a->t;4430 assert(t>= triangles && t < triangles+nbt);4431 4432 }4433 Icoor2 detop ;4434 int kkkk =0; // number of test triangle4435 4436 while ( t->det < 0)4437 { // the initial triangles is outside4438 int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;4439 assert(k0>=0); // k0 the NULL vertex4440 int k1=NextVertex[k0],k2=PreviousVertex[k0];4441 dete[k0]=det(B,(*t)[k1],(*t)[k2]);4442 dete[k1]=dete[k2]=-1;4443 if (dete[k0] > 0) // outside B4444 return t;4445 t = t->TriangleAdj(OppositeEdge[k0]);4446 assert(kkkk++ < 2);4447 }4448 4449 jj=0;4450 detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);4451 4452 while(t->det > 0 )4453 {4454 assert( kkkk++ < 2000 );4455 j= OppositeVertex[jj];4456 dete[j] = detop; //det(*b,*s1,*s2);4457 jn = NextVertex[j];4458 jp = PreviousVertex[j];4459 dete[jp]= det(*(*t)(j),*(*t)(jn),B);4460 dete[jn] = t->det-dete[j] -dete[jp];4461 4462 // count the number k of dete <04463 int k=0,ii[3];4464 if (dete[0] < 0 ) ii[k++]=0;4465 if (dete[1] < 0 ) ii[k++]=1;4466 if (dete[2] < 0 ) ii[k++]=2;4467 // 0 => ok4468 // 1 => go in way 14469 // 2 => two way go in way 1 or 2 randomly4470 4471 if (k==0)4472 break;4473 if (k == 2 && BinaryRand())4474 Exchange(ii[0],ii[1]);4475 assert ( k < 3);4476 TriangleAdjacent t1 = t->Adj(jj=ii[0]);4477 if ((t1.det() < 0 ) && (k == 2))4478 t1 = t->Adj(jj=ii[1]);4479 t=t1;4480 j=t1;// for optimisation we now the -det[OppositeVertex[j]];4481 detop = -dete[OppositeVertex[jj]];4482 jj = j;4483 }4484 4485 if (t->det<0) // outside triangle4486 dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;4487 // NbOfTriangleSearchFind += kkkk;4488 return t;4489 }4490 1164 4491 1165 double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) { -
issm/trunk/src/c/Bamgx/Triangles.cpp
r2810 r2811 14 14 15 15 static const Direction NoDirOfSearch=Direction(); 16 long NbUnSwap =0; 16 17 17 18 /*Constructors/Destructors*/ … … 26 27 FillHoleInMesh(); 27 28 } 29 /*}}}1*/ 30 /*FUNCTION Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb){{{1*/ 31 Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb) 32 : Gh(*(new Geometry())), BTh(*this) { 33 34 char cname[] = "trunc"; 35 36 int i,k,itadj; 37 int kt=0; 38 int * kk = new int [Tho.nbv]; 39 Int4 * reft = new Int4[Tho.nbt]; 40 Int4 nbInT = Tho.ConsRefTriangle(reft); 41 Int4 * refv = new Int4[Tho.nbv]; 42 43 for (i=0;i<Tho.nbv;i++) 44 kk[i]=-1; 45 for (i=0;i<Tho.nbv;i++) 46 refv[i]=0; 47 int nbNewBedge =0; 48 // int nbOldBedge =0; 49 for (i=0;i<Tho.nbt;i++) 50 if( reft[i] >=0 && flag[i]) 51 { 52 const Triangle & t = Tho.triangles[i]; 53 kt++; 54 kk[Tho.Number(t[0])]=1; 55 kk[Tho.Number(t[1])]=1; 56 kk[Tho.Number(t[2])]=1; 57 itadj=Tho.Number(t.TriangleAdj(0)); 58 if ( reft[itadj] >=0 && !flag[itadj]) 59 { nbNewBedge++; 60 refv[Tho.Number(t[VerticesOfTriangularEdge[0][0]])]=bb[i]; 61 refv[Tho.Number(t[VerticesOfTriangularEdge[0][1]])]=bb[i]; 62 } 63 itadj=Tho.Number(t.TriangleAdj(1)); 64 if ( reft[itadj] >=0 && !flag[itadj]) 65 { nbNewBedge++; 66 refv[Tho.Number(t[VerticesOfTriangularEdge[1][0]])]=bb[i]; 67 refv[Tho.Number(t[VerticesOfTriangularEdge[1][1]])]=bb[i];} 68 itadj=Tho.Number(t.TriangleAdj(2)); 69 if ( reft[itadj] >=0 && !flag[itadj]) 70 { nbNewBedge++; 71 refv[Tho.Number(t[VerticesOfTriangularEdge[2][0]])]=bb[i]; 72 refv[Tho.Number(t[VerticesOfTriangularEdge[2][1]])]=bb[i];} 73 } 74 k=0; 75 for (i=0;i<Tho.nbv;i++) 76 if (kk[i]>=0) 77 kk[i]=k++; 78 cout << " number of vertices " << k << " remove = " << Tho.nbv - k << endl; 79 cout << " number of triangles " << kt << " remove = " << nbInT-kt << endl; 80 cout << " number of New boundary edge " << nbNewBedge << endl; 81 Int4 inbvx =k; 82 PreInit(inbvx,cname); 83 for (i=0;i<Tho.nbv;i++) 84 if (kk[i]>=0) 85 { 86 vertices[nbv] = Tho.vertices[i]; 87 if (!vertices[nbv].ref()) 88 vertices[nbv].ReferenceNumber = refv[i]; 89 nbv++; 90 } 91 assert(inbvx == nbv); 92 for (i=0;i<Tho.nbt;i++) 93 if( reft[i] >=0 && flag[i]) 94 { 95 const Triangle & t = Tho.triangles[i]; 96 int i0 = Tho.Number(t[0]); 97 int i1 = Tho.Number(t[1]); 98 int i2 = Tho.Number(t[2]); 99 assert(i0>=0 && i1 >= 0 && i2 >= 0); 100 assert(i0<Tho.nbv && i1 <Tho.nbv && i2 <Tho.nbv); 101 // cout <<i<< " F" << flag[i] << " T " << nbt << " = " << kk[i0] << " " << kk[i1] << " " << kk[i2] ; 102 // cout << " OT " << i0 << " " << i1 << " " << i2 << " " << reft[i] << endl; 103 triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]); 104 triangles[nbt].color = Tho.subdomains[reft[i]].ref; 105 nbt++; 106 } 107 assert(kt==nbt); 108 if (nbt ==0 && nbv ==0) { 109 cout << "Error all triangles was remove " << endl; 110 MeshError(999,this); 111 } 112 delete [] kk; 113 delete [] reft; 114 delete [] refv; 115 double cutoffradian = 10.0/180.0*Pi; 116 ConsGeometry(cutoffradian); 117 Gh.AfterRead(); 118 SetIntCoor(); 119 FillHoleInMesh(); 120 121 assert(NbSubDomains); 122 assert(subdomains[0].head && subdomains[0].head->link); 123 124 } 125 /*}}}1*/ 126 /*FUNCTION Triangles::~Triangles(){{{1*/ 127 Triangles::~Triangles() { 128 long int verbosity=2; 129 //assert(NbRef<=0); 130 if (CurrentTh == this) CurrentTh=0; 131 //if(vertices) delete [] vertices; //TEST crash if not commented 132 if(edges) delete [] edges; 133 if(triangles) delete [] triangles; 134 if(quadtree) delete quadtree; 135 //if(ordre) delete [] ordre; //TEST crash if not commented 136 if( subdomains) delete [] subdomains; 137 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; 138 if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex; 139 //if (name) delete [] name; //TEST crash if not commented 140 if (identity) delete [] identity; 141 if (VertexOnBThVertex) delete [] VertexOnBThVertex; 142 if (VertexOnBThEdge) delete [] VertexOnBThEdge; 143 144 if (&Gh) { 145 if (Gh.NbRef>0) Gh.NbRef--; 146 else if (Gh.NbRef==0) delete &Gh; 147 } 148 if (&BTh && (&BTh != this)) { 149 if (BTh.NbRef>0) BTh.NbRef--; 150 else if (BTh.NbRef==0) delete &BTh; 151 } 152 PreInit(0); // set all to zero 153 } 154 /*}}}1*/ 155 /*FUNCTION Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) COPY{{{1*/ 156 Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) 157 : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) { 158 Gh.NbRef++; 159 nbvxx = Max(nbvxx,Th.nbv); 160 Int4 i; 161 // do all the allocation to be sure all the pointer existe 162 163 char * cname = 0; 164 if (Th.name) 165 { 166 cname = new char[strlen(Th.name)+1]; 167 strcpy(cname,Th.name); 168 } 169 PreInit(nbvxx,cname);// to make the allocation 170 // copy of triangles 171 nt=Th.nt; 172 nbv = Th.nbv; 173 nbt = Th.nbt; 174 nbiv = Th.nbiv; 175 nbe = Th.nbe; 176 NbSubDomains = Th.NbSubDomains; 177 NbOutT = Th.NbOutT; 178 NbOfQuad = Th.NbOfQuad ; 179 NbOfSwapTriangle =0; 180 NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex; 181 if(NbVerticesOnGeomVertex) 182 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 183 NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge; 184 if (NbVerticesOnGeomEdge) 185 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ; 186 if (& BTh == & Th.BTh) // same back ground 187 { 188 BTh.NbRef++; 189 NbVertexOnBThVertex = Th.NbVertexOnBThVertex; 190 if(NbVertexOnBThVertex) 191 VertexOnBThVertex = new VertexOnVertex[NbVertexOnBThVertex]; 192 NbVertexOnBThEdge = Th.NbVertexOnBThEdge; 193 if(NbVertexOnBThEdge) 194 VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge]; 195 } 196 else 197 { // no add on back ground mesh 198 BTh.NbRef++; 199 NbVertexOnBThVertex=0; 200 VertexOnBThVertex=0; 201 NbVertexOnBThEdge=0; 202 VertexOnBThEdge=0; 203 // assert (& BTh == this); // --- a voir 204 205 } 206 207 208 if(nbe) 209 edges = new Edge[nbe]; 210 if(NbSubDomains) 211 subdomains = new SubDomain[NbSubDomains]; 212 pmin = Th.pmin; 213 pmax = Th.pmax; 214 coefIcoor = Th.coefIcoor; 215 for(i=0;i<nbt;i++) 216 triangles[i].Set(Th.triangles[i],Th,*this); 217 for(i=0;i<nbe;i++) 218 edges[i].Set(Th,i,*this); 219 for(i=0;i<nbv;i++) 220 vertices[i].Set(Th.vertices[i],Th,*this); 221 for(i=0;i<NbSubDomains;i++) 222 subdomains[i].Set(Th,i,*this); 223 for (i=0;i<NbVerticesOnGeomVertex;i++) 224 VerticesOnGeomVertex[i].Set(Th.VerticesOnGeomVertex[i],Th,*this); 225 for (i=0;i<NbVerticesOnGeomEdge;i++) 226 VerticesOnGeomEdge[i].Set(Th.VerticesOnGeomEdge[i],Th,*this); 227 quadtree=0; 228 229 230 // assert(!OutSidesTriangles); 231 } 28 232 /*}}}1*/ 29 233 … … 1648 1852 } 1649 1853 /*}}}1*/ 1854 /*FUNCTION Triangles::swap{{{1*/ 1855 int Triangle::swap(Int2 a,int koption){ 1856 if(a/4 !=0) return 0;// arete lock or MarkUnSwap 1857 1858 register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent 1859 register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles 1860 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap 1861 1862 register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; 1863 register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; 1864 register Vertex *s1=t1->ns[OppositeVertex[a1]]; 1865 register Vertex *s2=t2->ns[OppositeVertex[a2]]; 1866 1867 Icoor2 det1=t1->det , det2=t2->det ; 1868 Icoor2 detT = det1+det2; 1869 Icoor2 detA = Abs(det1) + Abs(det2); 1870 Icoor2 detMin = Min(det1,det2); 1871 1872 int OnSwap = 0; 1873 // si 2 triangle infini (bord) => detT = -2; 1874 if (sa == 0) {// les deux triangles sont frontieres 1875 det2=bamg::det(s2->i,sb->i,s1->i); 1876 OnSwap = det2 >0;} 1877 else if (sb == 0) { // les deux triangles sont frontieres 1878 det1=bamg::det(s1->i,sa->i,s2->i); 1879 OnSwap = det1 >0;} 1880 else if(( s1 != 0) && (s2 != 0) ) { 1881 det1 = bamg::det(s1->i,sa->i,s2->i); 1882 det2 = detT - det1; 1883 OnSwap = (Abs(det1) + Abs(det2)) < detA; 1884 1885 Icoor2 detMinNew=Min(det1,det2); 1886 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test 1887 if (! OnSwap &&(detMinNew>0)) { 1888 OnSwap = detMin ==0; 1889 if (! OnSwap) { 1890 int kopt = koption; 1891 while (1) 1892 if(kopt) { 1893 // critere de Delaunay pure isotrope 1894 register Icoor2 xb1 = sb->i.x - s1->i.x, 1895 x21 = s2->i.x - s1->i.x, 1896 yb1 = sb->i.y - s1->i.y, 1897 y21 = s2->i.y - s1->i.y, 1898 xba = sb->i.x - sa->i.x, 1899 x2a = s2->i.x - sa->i.x, 1900 yba = sb->i.y - sa->i.y, 1901 y2a = s2->i.y - sa->i.y; 1902 register double 1903 cosb12 = double(xb1*x21 + yb1*y21), 1904 cosba2 = double(xba*x2a + yba*y2a) , 1905 sinb12 = double(det2), 1906 sinba2 = double(t2->det); 1907 1908 1909 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) 1910 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); 1911 // if(CurrentTh) 1912 // cout << "swap " << CurrentTh->Number(sa) << " " << CurrentTh->Number(sb) << " " ; 1913 // cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12 1914 // << " Onswap = " << OnSwap << endl; 1915 break; 1916 } 1917 else 1918 { 1919 // critere de Delaunay anisotrope 1920 Real8 som; 1921 I2 AB=(I2) *sb - (I2) *sa; 1922 I2 MAB2=((I2) *sb + (I2) *sa); 1923 R2 MAB(MAB2.x*0.5,MAB2.y*0.5); 1924 I2 A1=(I2) *s1 - (I2) *sa; 1925 I2 D = (I2) * s1 - (I2) * sb ; 1926 R2 S2(s2->i.x,s2->i.y); 1927 R2 S1(s1->i.x,s1->i.y); 1928 { 1929 Metric M=s1->m; 1930 R2 ABo = M.Orthogonal(AB); 1931 R2 A1o = M.Orthogonal(A1); 1932 // (A+B)+ x ABo = (S1+B)/2+ y A1 1933 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 1934 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); 1935 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 1936 if (Abs(d) > dd*1.e-3) { 1937 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 1938 som = M(C - S2)/M(C - S1); 1939 } else 1940 {kopt=1;continue;} 1941 1942 } 1943 { 1944 Metric M=s2->m; 1945 R2 ABo = M.Orthogonal(AB); 1946 R2 A1o = M.Orthogonal(A1); 1947 // (A+B)+ x ABo = (S1+B)/2+ y A1 1948 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 1949 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); 1950 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 1951 if(Abs(d) > dd*1.e-3) { 1952 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 1953 som += M(C - S2)/M(C - S1); 1954 } else 1955 {kopt=1;continue;} 1956 } 1957 OnSwap = som < 2; 1958 break; 1959 } 1960 1961 } // OnSwap 1962 } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) 1963 } 1964 if( OnSwap ) 1965 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); 1966 else { 1967 NbUnSwap ++; 1968 t1->SetMarkUnSwap(a1); 1969 } 1970 return OnSwap; 1971 } 1972 /*}}}1*/ 1973 /*FUNCTION Triangles::MetricAt{{{1*/ 1974 Metric Triangles::MetricAt(const R2 & A) const { 1975 I2 a = toI2(A); 1976 Icoor2 deta[3]; 1977 Triangle * t =FindTriangleContening(a,deta); 1978 if (t->det <0) { // outside 1979 double ba,bb; 1980 TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ; 1981 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));} 1982 else { // inside 1983 Real8 aa[3]; 1984 Real8 s = deta[0]+deta[1]+deta[2]; 1985 aa[0]=deta[0]/s; 1986 aa[1]=deta[1]/s; 1987 aa[2]=deta[2]/s; 1988 return Metric(aa,(*t)[0],(*t)[1],(*t)[2]); 1989 } 1990 } 1991 /*}}}1*/ 1992 /*FUNCTION Triangles::Add{{{1*/ 1993 void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3) { 1994 // ------------------------------------------- 1995 // s2 1996 // ! 1997 // /|\ ! 1998 // / | \ ! 1999 // / | \ ! 2000 // tt1 / | \ tt0 ! 2001 // / |s \ ! 2002 // / . \ ! 2003 // / . ` \ ! 2004 // / . ` \ ! 2005 // ---------------- ! 2006 // s0 tt2 s1 2007 //-------------------------------------------- 2008 2009 Triangle * tt[3]; // the 3 new Triangles 2010 Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2]; 2011 Icoor2 det3local[3]; 2012 int infv = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; 2013 // infv = ordre of the infini vertex (null) 2014 register int nbd0 =0; // number of zero det3 2015 register int izerodet=-1,iedge; // izerodet = egde contening the vertex s 2016 Icoor2 detOld = t->det; 2017 2018 if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ) 2019 { 2020 cerr << " infv " << infv << " det = " << detOld << endl; 2021 cerr << Number(s) << " "<< Number(s0) << " " 2022 << Number(s1) << " " << Number(s2) << endl; 2023 MeshError(3); 2024 } 2025 2026 // if det3 do not exist then constuct det3 2027 if (!det3) { 2028 det3 = det3local; // alloc 2029 if ( infv<0 ) { 2030 det3[0]=bamg::det(s ,s1,s2); 2031 det3[1]=bamg::det(s0,s ,s2); 2032 det3[2]=bamg::det(s0,s1,s );} 2033 else { 2034 // one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk 2035 det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; 2036 det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; 2037 det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}} 2038 2039 2040 if (!det3[0]) izerodet=0,nbd0++; 2041 if (!det3[1]) izerodet=1,nbd0++; 2042 if (!det3[2]) izerodet=2,nbd0++; 2043 2044 if (nbd0 >0 ) // point s on a egde or on a vertex 2045 if (nbd0 == 1) { 2046 iedge = OppositeEdge[izerodet]; 2047 TriangleAdjacent ta = t->Adj(iedge); 2048 2049 // the point is on the edge 2050 // if the point is one the boundary 2051 // add the point in outside part 2052 if ( t->det >=0) { // inside triangle 2053 if ((( Triangle *) ta)->det < 0 ) { 2054 // add in outside triangle 2055 Add(s,( Triangle *) ta); 2056 return;} 2057 }} 2058 else { 2059 cerr << " bug " << nbd0 <<endl; 2060 cerr << " Bug double points in " << endl ; 2061 cerr << " s = " << Number(s) << " " << s << endl; 2062 cerr << " s0 = "<< Number(s0) << " " << s0 << endl; 2063 cerr << " s1 = "<< Number(s1) << " " << s1 << endl; 2064 cerr << " s2 = "<< Number(s2) << " " << s2 << endl; 2065 MeshError(5,this);} 2066 2067 // remove de MarkUnSwap edge 2068 t->SetUnMarkUnSwap(0); 2069 t->SetUnMarkUnSwap(1); 2070 t->SetUnMarkUnSwap(2); 2071 2072 tt[0]= t; 2073 tt[1]= &triangles[nbt++]; 2074 tt[2]= &triangles[nbt++]; 2075 2076 if (nbt>nbtx) { 2077 cerr << " No enougth triangles " << endl; 2078 MeshError(999,this); 2079 } 2080 2081 *tt[1]= *tt[2]= *t; 2082 // gestion of the link 2083 tt[0]->link=tt[1]; 2084 tt[1]->link=tt[2]; 2085 2086 (* tt[0])(OppositeVertex[0])=&s; 2087 (* tt[1])(OppositeVertex[1])=&s; 2088 (* tt[2])(OppositeVertex[2])=&s; 2089 2090 tt[0]->det=det3[0]; 2091 tt[1]->det=det3[1]; 2092 tt[2]->det=det3[2]; 2093 2094 // update adj des triangles externe 2095 tt[0]->SetAdjAdj(0); 2096 tt[1]->SetAdjAdj(1); 2097 tt[2]->SetAdjAdj(2); 2098 // update des adj des 3 triangle interne 2099 const int i0 = 0; 2100 const int i1= NextEdge[i0]; 2101 const int i2 = PreviousEdge[i0]; 2102 2103 tt[i0]->SetAdj2(i2,tt[i2],i0); 2104 tt[i1]->SetAdj2(i0,tt[i0],i1); 2105 tt[i2]->SetAdj2(i1,tt[i1],i2); 2106 2107 tt[0]->SetTriangleContainingTheVertex(); 2108 tt[1]->SetTriangleContainingTheVertex(); 2109 tt[2]->SetTriangleContainingTheVertex(); 2110 2111 2112 // swap if the point s is on a edge 2113 if(izerodet>=0) { 2114 // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; 2115 int rswap =tt[izerodet]->swap(iedge); 2116 2117 if (!rswap) 2118 { 2119 cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; 2120 } 2121 assert(rswap); 2122 } 2123 } 2124 /*}}}1*/ 2125 /*FUNCTION Triangles::SplitInternalEdgeWithBorderVertices{{{1*/ 2126 Int4 Triangles::SplitInternalEdgeWithBorderVertices(){ 2127 Int4 NbSplitEdge=0; 2128 SetVertexFieldOn(); 2129 Int4 it; 2130 Int4 nbvold=nbv; 2131 long int verbosity=2; 2132 for (it=0;it<nbt;it++) 2133 { 2134 Triangle &t=triangles[it]; 2135 if (t.link) 2136 for (int j=0;j<3;j++) 2137 if(!t.Locked(j) && !t.Hidden(j)){ 2138 Triangle &tt = *t.TriangleAdj(j); 2139 if ( &tt && tt.link && it < Number(tt)) 2140 { // an internal edge 2141 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 2142 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 2143 if (v0.on && v1.on) 2144 { 2145 R2 P= ((R2) v0 + (R2) v1)*0.5; 2146 if ( nbv<nbvx) { 2147 vertices[nbv].r = P; 2148 vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m); 2149 vertices[nbv].ReferenceNumber=0; 2150 vertices[nbv].DirOfSearch = NoDirOfSearch ; 2151 } 2152 NbSplitEdge++; 2153 if (verbosity>7) 2154 cout <<" Internal edge with two vertices on boundary" 2155 << Number(v0) << " " << Number(v1) << " by " << endl; 2156 } 2157 } 2158 } 2159 } 2160 ReMakeTriangleContainingTheVertex(); 2161 if (nbvold!=nbv) 2162 { 2163 Int4 iv = nbvold; 2164 Int4 NbSwap = 0; 2165 Icoor2 dete[3]; 2166 for (Int4 i=nbvold;i<nbv;i++) 2167 {// for all the new point 2168 Vertex & vi = vertices[i]; 2169 vi.i = toI2(vi.r); 2170 vi.r = toR2(vi.i); 2171 // if (!quadtree->ToClose(vi,seuil,hi,hj)) { 2172 // a good new point 2173 vi.ReferenceNumber=0; 2174 vi.DirOfSearch =NoDirOfSearch; 2175 // cout << " Add " << Number(vi) << " " << vi 2176 // << " " << Number(vi) << " <--> " << Number(vi) <<endl; 2177 Triangle *tcvi = FindTriangleContening(vi.i,dete); 2178 if (tcvi && !tcvi->link) { 2179 cout << i << " PB insert point " << Number(vi) << vi << Number(vi) 2180 << " tcvi = " << tcvi << " " << tcvi->link << endl; 2181 cout << (*tcvi)[1] << (*tcvi)[2] << endl; 2182 tcvi = FindTriangleContening(vi.i,dete); 2183 cout << (*tcvi)[1] << (*tcvi)[2] << endl; 2184 MeshError(1001,this); 2185 } 2186 2187 2188 quadtree->Add(vi); 2189 assert (tcvi && tcvi->det >= 0) ;// internal 2190 Add(vi,tcvi,dete); 2191 NbSwap += vi.Optim(1); 2192 iv++; 2193 // } 2194 } 2195 if (verbosity>3) 2196 { 2197 cout << " Nb Of New Point " << iv ; 2198 cout << " Nb swap = " << NbSwap << " to split internal edges with border vertices" ;} 2199 2200 nbv = iv; 2201 } 2202 if (NbSplitEdge > nbv-nbvold) 2203 cout << " Warning not enough vertices to split all internal edges " << endl 2204 << " we lost " << NbSplitEdge - ( nbv-nbvold) << " Edges Sorry " << endl; 2205 if (verbosity>2) 2206 cout << "SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << endl; 2207 return NbSplitEdge; 2208 } 2209 /*}}}1*/ 2210 /*FUNCTION Triangles::InsertNewPoints{{{1*/ 2211 Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) { 2212 long int verbosity=2; 2213 Real8 seuil= 1.414/2 ;// for two close point 2214 Int4 i; 2215 // insertion part --- 2216 2217 const Int4 nbvnew = nbv-nbvold; 2218 if (verbosity>5) 2219 cout << " Try to Insert the " <<nbvnew<< " new points " << endl; 2220 Int4 NbSwap=0; 2221 Icoor2 dete[3]; 2222 2223 // construction d'un ordre aleatoire 2224 if (! nbvnew) 2225 return 0; 2226 if (nbvnew) { 2227 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ; 2228 Int4 k3 = rand()%nbvnew ; 2229 for (Int4 is3=0; is3<nbvnew; is3++) { 2230 register Int4 j = nbvold +(k3 = (k3 + PrimeNumber)% nbvnew); 2231 register Int4 i = nbvold+is3; 2232 ordre[i]= vertices + j; 2233 ordre[i]->ReferenceNumber=i; 2234 } 2235 // be carefull 2236 Int4 iv = nbvold; 2237 for (i=nbvold;i<nbv;i++) 2238 {// for all the new point 2239 Vertex & vi = *ordre[i]; 2240 vi.i = toI2(vi.r); 2241 vi.r = toR2(vi.i); 2242 Real4 hx,hy; 2243 vi.m.Box(hx,hy); 2244 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor); 2245 if (!quadtree->ToClose(vi,seuil,hi,hj)) 2246 { 2247 // a good new point 2248 Vertex & vj = vertices[iv]; 2249 Int4 j = vj.ReferenceNumber; 2250 assert( &vj== ordre[j]); 2251 if(i!=j) 2252 { // for valgring 2253 Exchange(vi,vj); 2254 Exchange(ordre[j],ordre[i]); 2255 } 2256 vj.ReferenceNumber=0; 2257 // cout << " Add " << Number(vj) << " " << vj 2258 // << " " << Number(vi) << " <--> " << Number(vj) <<endl; 2259 Triangle *tcvj = FindTriangleContening(vj.i,dete); 2260 if (tcvj && !tcvj->link) 2261 { 2262 cerr << i << " PB insert point " << Number(vj) << vj << Number(vi) 2263 << " tcvj = " << tcvj << " " << tcvj->link << endl; 2264 cerr << (*tcvj)[1] << (*tcvj)[2] << endl; 2265 tcvj = FindTriangleContening(vj.i,dete); 2266 cout << (*tcvj)[1] << (*tcvj)[2] << endl; 2267 MeshError(1001,this); 2268 } 2269 2270 2271 quadtree->Add(vj); 2272 assert (tcvj && tcvj->det >= 0) ;// internal 2273 Add(vj,tcvj,dete); 2274 NbSwap += vj.Optim(1); 2275 iv++; 2276 } 2277 } 2278 if (verbosity>3) { 2279 cout << " Nb Of New Point " << iv << " Nb Of To close Points " << nbv-iv ; 2280 cout << " Nb swap = " << NbSwap << " after " ;} 2281 2282 nbv = iv; 2283 } 2284 2285 for (i=nbvold;i<nbv;i++) 2286 NbSwap += vertices[i].Optim(1); 2287 if (verbosity>3) 2288 cout << " NbSwap = " << NbSwap << endl; 2289 2290 2291 NbTSwap += NbSwap ; 2292 return nbv-nbvold; 2293 2294 } 2295 /*}}}1*/ 2296 /*FUNCTION Triangles::NewPoints{{{1*/ 2297 void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex) { 2298 long int verbosity=2; 2299 Int4 nbtold(nbt),nbvold(nbv); 2300 if (verbosity>2) 2301 cout << " -- Triangles::NewPoints "; 2302 if (verbosity>3)cout << " nbv (in) on Boundary = " << nbv <<endl; 2303 Int4 i,k; 2304 int j; 2305 Int4 *first_np_or_next_t = new Int4[nbtx]; 2306 Int4 NbTSwap =0; 2307 // insert old point 2308 nbtold = nbt; 2309 2310 if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx)) 2311 { 2312 // Bh.SetVertexFieldOn(); 2313 for (i=0;i<Bh.nbv;i++) 2314 { 2315 Vertex & bv = Bh[i]; 2316 if (!bv.on) { 2317 vertices[nbv].r = bv.r; 2318 vertices[nbv++].m = bv.m;} 2319 } 2320 int nbv1=nbv; 2321 Bh.ReMakeTriangleContainingTheVertex(); 2322 InsertNewPoints(nbvold,NbTSwap) ; 2323 if (verbosity>2) 2324 cout << " (Nb of Points from background mesh = " 2325 << nbv-nbvold << " / " << nbv1-nbvold << ")" << endl; 2326 } 2327 else 2328 Bh.ReMakeTriangleContainingTheVertex(); 2329 2330 Triangle *t; 2331 // generation of the list of next Triangle 2332 // at 1 time we test all the triangles 2333 Int4 Headt =0,next_t; 2334 for(i=0;i<nbt;i++) 2335 first_np_or_next_t[i]=-(i+1); 2336 // end list i >= nbt 2337 // the list of test triangle is 2338 // the next traingle on i is -first_np_or_next_t[i] 2339 int iter=0; 2340 // Big loop 2341 do { 2342 iter++; 2343 nbtold = nbt; 2344 nbvold = nbv; 2345 2346 // default size of IntersectionTriangle 2347 2348 i=Headt; 2349 next_t=-first_np_or_next_t[i]; 2350 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) 2351 { // for each triangle t 2352 // we can change first_np_or_next_t[i] 2353 // cout << " Do the triangle " << i << " Next_t=" << next_t << endl; 2354 assert(i>=0 && i < nbt ); 2355 first_np_or_next_t[i] = iter; 2356 for(j=0;j<3;j++) 2357 { // for each edge 2358 TriangleAdjacent tj(t,j); 2359 Vertex & vA = * tj.EdgeVertex(0); 2360 Vertex & vB = * tj.EdgeVertex(1); 2361 2362 if (!t->link) continue;// boundary 2363 if (t->det <0) continue; 2364 if (t->Locked(j)) continue; 2365 2366 TriangleAdjacent tadjj = t->Adj(j); 2367 Triangle * ta= tadjj; 2368 2369 if (ta->det <0) continue; 2370 2371 R2 A = vA; 2372 R2 B = vB; 2373 2374 k=Number(ta); 2375 2376 if(first_np_or_next_t[k]==iter) // this edge is done before 2377 continue; // next edge of the triangle 2378 2379 //const Int4 NbvOld = nbv; 2380 lIntTria.SplitEdge(Bh,A,B); 2381 lIntTria.NewPoints(vertices,nbv,nbvx); 2382 } // end loop for each edge 2383 2384 }// for triangle 2385 2386 if (!InsertNewPoints(nbvold,NbTSwap)) 2387 break; 2388 2389 for (i=nbtold;i<nbt;i++) 2390 first_np_or_next_t[i]=iter; 2391 2392 Headt = nbt; // empty list 2393 for (i=nbvold;i<nbv;i++) 2394 { // for all the triangle contening the vertex i 2395 Vertex * s = vertices + i; 2396 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]); 2397 Triangle * tbegin= (Triangle*) ta; 2398 Int4 kt; 2399 do { 2400 kt = Number((Triangle*) ta); 2401 if (first_np_or_next_t[kt]>0) 2402 first_np_or_next_t[kt]=-Headt,Headt=kt; 2403 assert( ta.EdgeVertex(0) == s); 2404 ta = Next(Adj(ta)); 2405 } while ( (tbegin != (Triangle*) ta)); 2406 } 2407 2408 } while (nbv!=nbvold); 2409 2410 delete [] first_np_or_next_t; 2411 2412 Int4 NbSwapf =0,NbSwp; 2413 2414 // bofbof 2415 2416 2417 NbSwp = NbSwapf; 2418 for (i=0;i<nbv;i++) 2419 NbSwapf += vertices[i].Optim(0); 2420 /* 2421 for (i=0;i<nbv;i++) 2422 NbSwapf += vertices[i].Optim(0); 2423 for (i=0;i<nbv;i++) 2424 NbSwapf += vertices[i].Optim(0); 2425 for (i=0;i<nbv;i++) 2426 NbSwapf += vertices[i].Optim(0); 2427 for (i=0;i<nbv;i++) 2428 NbSwapf += vertices[i].Optim(0); 2429 */ 2430 NbTSwap += NbSwapf ; 2431 if (verbosity>3) cout << " " ; 2432 if (verbosity>2) 2433 cout << " Nb Of Vertices =" << nbv << " Nb of triangles = " << nbt-NbOutT 2434 << " NbSwap final = " << NbSwapf << " Nb Total Of Swap = " << NbTSwap << endl; 2435 2436 2437 } 2438 /*}}}1*/ 2439 /*FUNCTION Triangles::NewPointsOld{{{1*/ 2440 void Triangles::NewPointsOld(Triangles & Bh) { 2441 long int verbosity=2; 2442 Real8 seuil= 0.7 ;// for two neart point 2443 if (verbosity>1) 2444 cout << " begin : Triangles::NewPointsOld " << endl; 2445 Int4 i,k; 2446 int j; 2447 Int4 BeginNewPoint[3]; 2448 Int4 EndNewPoint[3]; 2449 #ifdef TRACETRIANGLE 2450 Int4 trace=0; 2451 #endif 2452 int step[3]; 2453 Int4 *first_np_or_next_t = new Int4[nbtx]; 2454 Int4 ColorEdge[3]; 2455 Int4 color=-1; 2456 Triangle *t; 2457 // generation of the list of next Triangle 2458 // at 1 time we test all the triangles 2459 Int4 Headt =0,next_t; 2460 for(i=0;i<nbt;i++) 2461 first_np_or_next_t[i]=-(i+1); 2462 // end list i >= nbt 2463 // the list of test triangle is 2464 // the next Triangle on i is -first_np_or_next_t[i] 2465 Int4 nbtold,nbvold; 2466 2467 // Big loop 2468 do { 2469 nbtold = nbt; 2470 nbvold = nbv; 2471 // default size of IntersectionTriangle 2472 2473 i=Headt; 2474 next_t=-first_np_or_next_t[i]; 2475 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) 2476 { // for each triangle t 2477 // we can change first_np_or_next_t[i] 2478 #ifdef TRACETRIANGLE 2479 trace = TRACETRIANGLE <0 ? 1 : i == TRACETRIANGLE; 2480 #endif 2481 // cout << " Do the triangle " << i << " Next_t=" << next_t << endl; 2482 assert(i>=0 && i < nbt ); 2483 first_np_or_next_t[i] = nbv; // to save the fist new point of triangle 2484 for(j=0;j<3;j++) 2485 { // for each edge 2486 TriangleAdjacent tj(t,j); 2487 // color++;// the color is 3i+j 2488 color = 3*i + j ;; 2489 ColorEdge[j]=color; 2490 BeginNewPoint[j]=nbv; 2491 EndNewPoint[j]=nbv-1; 2492 step[j]=1;// right sens 2493 Vertex & vA = * tj.EdgeVertex(0); 2494 Vertex & vB = * tj.EdgeVertex(1); 2495 2496 #ifdef TRACETRIANGLE 2497 if(trace) { 2498 cout << " i " << Number(vA) <<" j "<< Number(vB) 2499 << " " << t->Locked(j) ; 2500 } 2501 #endif 2502 if (!t->link) continue;// boundary 2503 if (t->det <0) continue; 2504 if (t->Locked(j)) continue; 2505 2506 TriangleAdjacent tadjj = t->Adj(j); 2507 Triangle * ta= tadjj; 2508 2509 if (ta->det <0) continue; 2510 2511 R2 A = vA; 2512 R2 B = vB; 2513 2514 k=Number(ta); 2515 // the 2 opposite vertices 2516 const Vertex & vC1 = *tj.OppositeVertex(); 2517 const Vertex & vC2 = *tadjj.OppositeVertex(); 2518 2519 #ifdef TRACETRIANGLE 2520 trace = trace || k == TRACETRIANGLE; 2521 if(trace) { 2522 cout << "Test Arete " << i << " AB = " << A << B 2523 << "i " <<Number(vA)<< "j" <<Number(vB); 2524 cout << " link" <<(int)t->link << " ta=" << Number( ta) 2525 << " det " <<ta->det ; 2526 cout << " hA = " <<vA.m.h << " hB = " <<vB.m.h ; 2527 cout << " loc " << ta->Locked(j) << endl; 2528 } 2529 #endif 2530 2531 if(first_np_or_next_t[k]>0) { // this edge is done before 2532 // find the color of the edge and begin , end of newpoint 2533 register int kk = t->NuEdgeTriangleAdj(j); 2534 assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1])); 2535 assert ((*t)(VerticesOfTriangularEdge[j][1]) == (*ta)(VerticesOfTriangularEdge[kk][0])); 2536 register Int4 kolor =3*k + kk; 2537 ColorEdge[j]=kolor; 2538 register Int4 kkk= 1; 2539 step[j]=-1;// other sens 2540 BeginNewPoint[j]=0; 2541 EndNewPoint[j]=-1; // empty list 2542 for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++) 2543 if (vertices[iv].color > kolor) 2544 break; // the color is passed 2545 else if (vertices[iv].color == kolor) { 2546 EndNewPoint[j]=iv; 2547 if (kkk) // one time test 2548 kkk=0,BeginNewPoint[j]=iv;} 2549 continue; // next edge of the triangle 2550 } // end if( k < i) 2551 2552 2553 const Int4 NbvOld = nbv; 2554 lIntTria.SplitEdge(Bh,A,B); 2555 // Int4 NbvNp = 2556 lIntTria.NewPoints(vertices,nbv,nbvx); 2557 Int4 nbvNew=nbv; 2558 nbv = NbvOld; 2559 for (Int4 iv=NbvOld;iv<nbvNew;iv++) { 2560 vertices[nbv].color = color; 2561 vertices[nbv].ReferenceNumber=nbv;// circular link 2562 R2 C = vertices[iv].r; 2563 vertices[nbv].r = C; 2564 vertices[nbv].m = vertices[iv].m ; 2565 // test if the new point is not to close to the 2 opposite vertex 2566 R2 CC1 = C-vC1 , CC2 = C-vC2; 2567 if ( ( (vC1.m(CC1) + vertices[nbv].m(CC1)) > seuil) 2568 && ( (vC2.m(CC2) + vertices[nbv].m(CC2)) > seuil) ) 2569 nbv++; 2570 } 2571 2572 EndNewPoint[j] = nbv-1; 2573 } // end loop for each edge 2574 2575 #ifdef TRACETRIANGLE 2576 if(trace) { 2577 // verification des point cree 2578 cout << "\n ------------ " << t->link << " " << t->det 2579 << " b " << BeginNewPoint[0] << " " << BeginNewPoint[1] 2580 << " " << BeginNewPoint[2] << " " 2581 << " e " << EndNewPoint[0] << " " << EndNewPoint[1] 2582 << " " << EndNewPoint[2] << " " 2583 << " s " << step[0] << " " << step[1] << " " << step[2] << " " 2584 << endl; 2585 } 2586 #endif 2587 2588 if (!t->link) continue;// boundary 2589 if (t->det<=0) continue;// outside 2590 // continue; 2591 for(int j0=0;j0<3;j0++) 2592 for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++) 2593 { 2594 // find the neart point one the opposite edge 2595 // to compute i1 2596 Vertex & vi0 = vertices[i0]; 2597 int kstack = 0; 2598 Int4 stack[10]; 2599 // Int4 savRef[10]; 2600 int j1 = j0; 2601 while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge 2602 // computation of the intersection of edge j1 and DOrto 2603 // take the good sens 2604 2605 if (BeginNewPoint[j1]> EndNewPoint[j1]) 2606 continue; // 2607 else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) { 2608 for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++) 2609 stack[kstack++] = ii1; 2610 continue;} 2611 2612 2613 int k0,k1; 2614 if (step[j1]<0) k0=1,k1=0; // reverse 2615 else k0=0,k1=1; 2616 R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]]; 2617 R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]]; 2618 R2 D = V11-V10; 2619 Real8 c0 = vi0.m(D,(R2) vi0); 2620 2621 Real8 c10 = vi0.m(D,V10); 2622 Real8 c11 = vi0.m(D,V11); 2623 2624 Real8 s; 2625 //cout << " --i0 = " << i0 << D << V10 << V11 << endl ; 2626 //cout << " c10 " << c10 << " c0 " << c0 << " c11 " << c11 << endl; 2627 if (( c10 < c0 ) && (c0 < c11)) 2628 s = (c11-c0)/(c11-c10); 2629 else if (( c11 < c0 ) && (c0 < c10)) 2630 s = (c11-c0) /(c11-c10); 2631 else break; 2632 R2 VP = V10*s + V11*(1-s); 2633 int sss = (c11-c10) >0 ? 1 : -1; 2634 // find the 2 point by dichotomie 2635 //cout << " t =" << Number(t) << " c0 " << c0 ; 2636 Int4 ii0 = BeginNewPoint[j1]; 2637 Int4 ii1 = EndNewPoint[j1]; 2638 Real8 ciii=-1,cii0=-1,cii1=-1 ; 2639 if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 ) 2640 stack[kstack++] = ii0;//cout << " add+0 " << ii0; 2641 else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 ) 2642 stack[kstack++] = ii1;//cout << " add+1 " << ii1; 2643 else { 2644 while ((ii1-ii0)> 1) { 2645 Int4 iii = (ii0+ii1)/2; 2646 ciii = vi0.m( D ,(R2) vertices[iii]); 2647 //cout << " (iii = " << iii << " " << ciii << ") "; 2648 if ( sss * (ciii - c0) <0 ) ii0 = iii; 2649 else ii1 = iii;} 2650 stack[kstack++] = ii0;// cout << " add0 " << ii0; 2651 if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1; 2652 } 2653 if (kstack >5) // bug ? 2654 cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl; 2655 } 2656 2657 stack[kstack++] = -1; // to stop 2658 Int4 i1; 2659 kstack =0; 2660 while( (i1=stack[kstack++]) >= 0) 2661 { // the two parameter is i0 and i1 2662 assert(i1 < nbv && i1 >= 0); 2663 assert(i0 < nbv && i0 >= 0); 2664 assert(i1 != i0); 2665 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0]; 2666 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01)); 2667 2668 2669 #ifdef TRACETRIANGLE 2670 if(trace) { 2671 cout << "\n test j" << j <<" " << i0 2672 << " " << i1 << " d01=" << d01 <<endl;} 2673 #endif 2674 assert (i0 >= nbvold); 2675 assert (i1 >= nbvold); 2676 assert(i0 != i1); 2677 if (d01 == 0) 2678 break; 2679 if ( d01 < seuil) 2680 if (i1<nbvold) { 2681 // remove all the points i0; 2682 register Int4 ip,ipp; 2683 for (ip=i0;i0 != (ipp = vertices[ip].ReferenceNumber);ip=ipp) 2684 vertices[ip].ReferenceNumber = -1;// mark remove 2685 vertices[ip].ReferenceNumber = -1;// mark remove 2686 } 2687 else { 2688 // remove on of two points 2689 register Int4 ip0, ip1, ipp0,ipp1; 2690 register int kk0=1,kk1=1; 2691 // count the number of common points to compute weight w0,w1 2692 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) kk0++; 2693 for (ip1=i1;i1 != (ipp1 = vertices[ip1].ReferenceNumber);ip1=ipp1) kk1++; 2694 2695 register Real8 w0 = ((Real8) kk0)/(kk0+kk1); 2696 register Real8 w1 = ((Real8) kk1)/(kk0+kk1); 2697 2698 // make a circular link 2699 Exchange(vertices[i0].ReferenceNumber,vertices[i1].ReferenceNumber); 2700 // the new coordinate 2701 R2 C //= vertices[i0] ; 2702 = vertices[i0].r *w0 + vertices[i1].r* w1; 2703 2704 #ifdef TRACETRIANGLE 2705 if(trace) { 2706 cout << "\n ref = "<< vertices[i0].ref << " " <<vertices[i1].ref <<endl; 2707 } 2708 #endif 2709 // update the new point points of the list 2710 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) 2711 vertices[ip0].r = C; 2712 vertices[ip0].r = C; 2713 } 2714 } 2715 } // for (i0= .... 2716 }// for triangle 2717 2718 // remove of all the double points 2719 2720 Int4 ip,ipp,kkk=nbvold; 2721 for (i=nbvold;i<nbv;i++) 2722 if (vertices[i].ReferenceNumber>=0) {// good points 2723 // cout <<" i = " << i ; 2724 for (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp) 2725 vertices[ip].ReferenceNumber = -1;// mark remove 2726 vertices[ip].ReferenceNumber = -1;// mark remove 2727 // cout << i << " ---> " << kkk << endl; 2728 vertices[kkk] = vertices[i]; 2729 vertices[kkk].i = toI2(vertices[kkk].r); 2730 vertices[kkk++].ReferenceNumber = 0; 2731 2732 } 2733 2734 // insertion part --- 2735 2736 const Int4 nbvnew = kkk-nbvold; 2737 2738 cout << " Remove " << nbv - kkk << " to close vertex " ; 2739 cout << " and Insert the " <<nbvnew<< " new points " << endl; 2740 nbv = kkk; 2741 Int4 NbSwap=0; 2742 Icoor2 dete[3]; 2743 2744 // construction d'un ordre aleatoire 2745 if (! nbvnew) 2746 break; 2747 if (nbvnew) { 2748 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ; 2749 Int4 k3 = rand()%nbvnew ; 2750 for (Int4 is3=0; is3<nbvnew; is3++) 2751 ordre[nbvold+is3]= &vertices[nbvold +(k3 = (k3 + PrimeNumber)% nbvnew)]; 2752 2753 for (i=nbvold;i<nbv;i++) 2754 { Vertex * vi = ordre[i]; 2755 Triangle *tcvi = FindTriangleContening(vi->i,dete); 2756 // Vertex * nv = quadtree->NearestVertex(vi->i.x,vi->i.y); 2757 // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is " 2758 // << Number(nv) << nv->i << endl; 2759 // Int4 kt = Number(tcvi); 2760 // 2761 2762 quadtree->Add(*vi); // 2763 assert (tcvi->det >= 0) ;// internal 2764 Add(*vi,tcvi,dete),NbSwap += vi->Optim(1); 2765 } 2766 } 2767 cout << " Nb swap = " << NbSwap << " after " ; 2768 2769 for (i=nbvold;i<nbv;i++) 2770 NbSwap += vertices[i].Optim(1); 2771 cout << NbSwap << endl; 2772 2773 for (i=nbtold;i<nbt;i++) 2774 first_np_or_next_t[i]=1; 2775 2776 Headt = nbt; // empty list 2777 for (i=nbvold;i<nbv;i++) 2778 { // for all the triangle contening the vertex i 2779 Vertex * s = vertices + i; 2780 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]); 2781 Triangle * tbegin= (Triangle*) ta; 2782 Int4 kt; 2783 do { 2784 kt = Number((Triangle*) ta); 2785 if (first_np_or_next_t[kt]>0) 2786 first_np_or_next_t[kt]=-Headt,Headt=kt; 2787 assert( ta.EdgeVertex(0) == s); 2788 ta = Next(Adj(ta)); 2789 } while ( (tbegin != (Triangle*) ta)); 2790 } 2791 2792 2793 } while (nbv!=nbvold); 2794 delete [] first_np_or_next_t; 2795 cout << " end : Triangles::NewPoints old nbv=" << nbv << endl; 2796 2797 } 2798 /*}}}1*/ 2799 /*FUNCTION Triangles::Insert{{{1*/ 2800 void Triangles::Insert() { 2801 long int verbosity=2; 2802 if (verbosity>2) cout << " -- Insert initial " << nbv << " vertices " << endl ; 2803 Triangles * OldCurrentTh =CurrentTh; 2804 2805 CurrentTh=this; 2806 double time0=CPUtime(),time1,time2,time3; 2807 SetIntCoor(); 2808 Int4 i; 2809 for (i=0;i<nbv;i++) 2810 ordre[i]= &vertices[i] ; 2811 2812 // construction d'un ordre aleatoire 2813 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ; 2814 Int4 k3 = rand()%nbv ; 2815 for (int is3=0; is3<nbv; is3++) 2816 ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 2817 2818 2819 2820 2821 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 2822 if ( ++i >= nbv) { 2823 cerr << " All the vertices are aline " << endl; 2824 MeshError(998,this); } 2825 2826 // echange i et 2 dans ordre afin 2827 // que les 3 premiers ne soit pas aligne 2828 Exchange( ordre[2], ordre[i]); 2829 2830 // on ajoute un point a l'infini pour construire le maillage 2831 // afin d'avoir une definition simple des aretes frontieres 2832 nbt = 2; 2833 2834 2835 // on construit un maillage trivale forme 2836 // d'une arete et de 2 triangles 2837 // construit avec le 2 aretes orientes et 2838 Vertex * v0=ordre[0], *v1=ordre[1]; 2839 2840 triangles[0](0) = 0; // sommet pour infini 2841 triangles[0](1) = v0; 2842 triangles[0](2) = v1; 2843 2844 triangles[1](0) = 0;// sommet pour infini 2845 triangles[1](2) = v0; 2846 triangles[1](1) = v1; 2847 const int e0 = OppositeEdge[0]; 2848 const int e1 = NextEdge[e0]; 2849 const int e2 = PreviousEdge[e0]; 2850 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 2851 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 2852 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 2853 2854 triangles[0].det = -1; // faux triangles 2855 triangles[1].det = -1; // faux triangles 2856 2857 triangles[0].SetTriangleContainingTheVertex(); 2858 triangles[1].SetTriangleContainingTheVertex(); 2859 2860 triangles[0].link=&triangles[1]; 2861 triangles[1].link=&triangles[0]; 2862 2863 // nbtf = 2; 2864 if ( !quadtree ) quadtree = new QuadTree(this,0); 2865 quadtree->Add(*v0); 2866 quadtree->Add(*v1); 2867 2868 // on ajoute les sommets un Ò un 2869 Int4 NbSwap=0; 2870 2871 time1=CPUtime(); 2872 2873 if (verbosity>3) cout << " -- Begin of insertion process " << endl; 2874 2875 for (Int4 icount=2; icount<nbv; icount++) { 2876 Vertex *vi = ordre[icount]; 2877 // cout << " Insert " << Number(vi) << endl; 2878 Icoor2 dete[3]; 2879 Triangle *tcvi = FindTriangleContening(vi->i,dete); 2880 quadtree->Add(*vi); 2881 Add(*vi,tcvi,dete); 2882 NbSwap += vi->Optim(1,0); 2883 2884 }// fin de boucle en icount 2885 time2=CPUtime(); 2886 if (verbosity>3) 2887 cout << " NbSwap of insertion " << NbSwap 2888 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2889 << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv " 2890 << (float)NbUnSwap /(float) nbv 2891 <<endl; 2892 NbUnSwap = 0; 2893 // construction d'un ordre aleatoire 2894 // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ; 2895 #ifdef NBLOOPOPTIM 2896 2897 k3 = rand()%nbv ; 2898 for (int is4=0; is4<nbv; is4++) 2899 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 2900 2901 double timeloop = time2 ; 2902 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++) 2903 { 2904 double time000 = timeloop; 2905 Int4 NbSwap = 0; 2906 for (int is1=0; is1<nbv; is1++) 2907 NbSwap += ordre[is1]->Optim(0,0); 2908 timeloop = CPUtime(); 2909 if (verbosity>3) 2910 cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap 2911 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2912 << " CPU=" << timeloop - time000 << " s, " 2913 << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv 2914 << endl; 2915 NbUnSwap = 0; 2916 if(!NbSwap) break; 2917 } 2918 ReMakeTriangleContainingTheVertex(); 2919 // because we break the TriangleContainingTheVertex 2920 #endif 2921 time3=CPUtime(); 2922 if (verbosity>4) 2923 cout << " init " << time1 - time0 << " initialisation, " 2924 << time2 - time1 << "s, insert point " 2925 << time3 -time2 << "s, optim " << endl 2926 << " Init Total Cpu Time = " << time3 - time0 << "s " << endl; 2927 2928 CurrentTh=OldCurrentTh; 2929 } 2930 /*}}}1*/ 2931 /*FUNCTION Triangles::ForceBoundary{{{1*/ 2932 void Triangles::ForceBoundary() { 2933 long int verbosity=2; 2934 if (verbosity > 2) 2935 cout << " -- ForceBoundary nb of edge " << nbe << endl; 2936 int k=0; 2937 Int4 nbfe=0,nbswp=0,Nbswap=0; 2938 for (Int4 t = 0; t < nbt; t++) 2939 if (!triangles[t].det) 2940 k++,cerr << " det T" << t << " = " << 0 << endl; 2941 if (k!=0) { 2942 cerr << " ther is " << k << " triangles of mes = 0 " << endl; 2943 MeshError(11,this);} 2944 2945 TriangleAdjacent ta(0,0); 2946 for (Int4 i = 0; i < nbe; i++) 2947 { 2948 nbswp = ForceEdge(edges[i][0],edges[i][1],ta); 2949 2950 if ( nbswp < 0) k++; 2951 else Nbswap += nbswp; 2952 if (nbswp) nbfe++; 2953 if ( nbswp < 0 && k < 5) 2954 { 2955 cerr << " Missing Edge " << i << " v0 = " << Number(edges[i][0]) << edges[i][0].r 2956 <<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << " " << (Triangle *) ta ; 2957 if(ta.t) 2958 { 2959 Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1); 2960 cerr << " crossing with [" << Number(aa) << ", " << Number(bb) << "]\n"; 2961 } 2962 else cerr << endl; 2963 2964 } 2965 if ( nbswp >=0 && edges[i].on->Cracked()) 2966 ta.SetCracked(); 2967 } 2968 2969 2970 if (k!=0) { 2971 cerr << " they is " << k << " lost edges " << endl; 2972 cerr << " The boundary is crossing may be!" << endl; 2973 MeshError(10,this); 2974 } 2975 for (Int4 j=0;j<nbv;j++) 2976 Nbswap += vertices[j].Optim(1,0); 2977 if (verbosity > 3) 2978 cout << " Nb of inforced edge = " << nbfe << " Nb of Swap " << Nbswap << endl; 2979 2980 } 2981 /*}}}1*/ 2982 /*FUNCTION Triangles::FindSubDomain{{{1*/ 2983 void Triangles::FindSubDomain(int OutSide=0) { 2984 long int verbosity=0; 2985 2986 if (verbosity >2) 2987 { 2988 if (OutSide) 2989 cout << " -- Find all external sub-domain "; 2990 else 2991 cout << " -- Find all internal sub-domain "; 2992 if(verbosity>99) 2993 { 2994 2995 for(int i=0;i<nbt;++i) 2996 cout << i<< " " << triangles[i] << endl; 2997 } 2998 2999 } 3000 // if (verbosity > 4) cout << " OutSide=" << OutSide << endl; 3001 short * HeapArete = new short[nbt]; 3002 Triangle ** HeapTriangle = new Triangle* [nbt]; 3003 Triangle *t,*t1; 3004 Int4 k,it; 3005 3006 for (Int4 itt=0;itt<nbt;itt++) 3007 triangles[itt].link=0; // par defaut pas de couleur 3008 3009 Int4 NbSubDomTot =0; 3010 for ( it=0;it<nbt;it++) { 3011 if ( ! triangles[it].link ) { 3012 t = triangles + it; 3013 NbSubDomTot++;; // new composante connexe 3014 Int4 i = 0; // niveau de la pile 3015 t->link = t ; // sd forme d'un triangle cicular link 3016 3017 HeapTriangle[i] =t ; 3018 HeapArete[i] = 3; 3019 3020 while (i >= 0) // boucle sur la pile 3021 { while ( HeapArete[i]--) // boucle sur les 3 aretes 3022 { 3023 int na = HeapArete[i]; 3024 Triangle * tc = HeapTriangle[i]; // triangle courant 3025 if( ! tc->Locked(na)) // arete non frontiere 3026 { 3027 Triangle * ta = tc->TriangleAdj(na) ; // næ triangle adjacent 3028 if (ta->link == 0 ) // non deja chainer => on enpile 3029 { 3030 i++; 3031 ta->link = t->link ; // on chaine les triangles 3032 t->link = ta ; // d'un meme sous domaine 3033 HeapArete[i] = 3; // pour les 3 triangles adjacents 3034 HeapTriangle[i] = ta; 3035 }} 3036 } // deplie fin de boucle sur les 3 adjacences 3037 i--; 3038 } 3039 } 3040 } 3041 3042 // supression de tous les sous domaine infini <=> contient le sommet NULL 3043 it =0; 3044 NbOutT = 0; 3045 while (it<nbt) { 3046 if (triangles[it].link) 3047 { 3048 if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) )) 3049 { 3050 // infini triangle 3051 NbSubDomTot --; 3052 // cout << " triangle infini " << it << triangles[it] << endl; 3053 t=&triangles[it]; 3054 NbOutT--; // on fait un coup de trop. 3055 while (t){ // cout << Number(t) << " " << endl; 3056 NbOutT++; 3057 t1=t; 3058 t=t->link; 3059 t1->link=0;}//while (t) 3060 } 3061 } 3062 it++;} // end while (it<nbt) 3063 if (nbt == NbOutT || !NbSubDomTot) 3064 { 3065 cout << "\n error : " << NbOutT << " " << NbSubDomTot <<" " << nbt << endl; 3066 cerr << "Error: The boundary is not close => All triangles are outside " << endl; 3067 MeshError(888,this); 3068 } 3069 3070 delete [] HeapArete; 3071 delete [] HeapTriangle; 3072 3073 3074 if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains ) 3075 { // No geom sub domain 3076 Int4 i; 3077 if (subdomains) delete [] subdomains; 3078 subdomains = new SubDomain[ NbSubDomTot]; 3079 NbSubDomains= NbSubDomTot; 3080 for ( i=0;i<NbSubDomains;i++) { 3081 subdomains[i].head=NULL; 3082 subdomains[i].ref=i+1; 3083 } 3084 Int4 * mark = new Int4[nbt]; 3085 for (it=0;it<nbt;it++) 3086 mark[it]=triangles[it].link ? -1 : -2; 3087 3088 it =0; 3089 k = 0; 3090 while (it<nbt) { 3091 if (mark[it] == -1) { 3092 t1 = & triangles[it]; 3093 t = t1->link; 3094 mark[it]=k; 3095 subdomains[k].head = t1; 3096 // cout << " New -- " << Number(t1) << " " << it << endl; 3097 do {// cout << " k " << k << " " << Number(t) << endl; 3098 mark[Number(t)]=k; 3099 t=t->link; 3100 } while (t!=t1); 3101 mark[it]=k++;} 3102 // else if(mark[it] == -2 ) triangles[it].Draw(999); 3103 it++;} // end white (it<nbt) 3104 assert(k== NbSubDomains); 3105 if(OutSide) 3106 { 3107 // to remove all the sub domain by parity adjacents 3108 // because in this case we have only the true boundary edge 3109 // so teh boundary is manifold 3110 Int4 nbk = NbSubDomains; 3111 while (nbk) 3112 for (it=0;it<nbt && nbk ;it++) 3113 for (int na=0;na<3 && nbk ;na++) 3114 { 3115 Triangle *ta = triangles[it].TriangleAdj(na); 3116 Int4 kl = ta ? mark[Number(ta)] : -2; 3117 Int4 kr = mark[it]; 3118 if(kr !=kl) { 3119 //cout << kl << " " << kr << " rl " << subdomains[kl].ref 3120 // << " rr " << subdomains[kr].ref ; 3121 if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0) 3122 nbk--,subdomains[kr].ref=subdomains[kl].ref-1; 3123 if (kr >=0 && subdomains[kr].ref <0 && kl >=0 && subdomains[kl].ref>=0) 3124 nbk--,subdomains[kl].ref=subdomains[kr].ref-1; 3125 if(kr<0 && kl >=0 && subdomains[kl].ref>=0) 3126 nbk--,subdomains[kl].ref=-1; 3127 if(kl<0 && kr >=0 && subdomains[kr].ref>=0) 3128 nbk--,subdomains[kr].ref=-1; 3129 // cout << " after \t " 3130 // << kl << subdomains[kl].ref << " rr " << kr 3131 // << subdomains[kr].ref << endl; 3132 } 3133 } 3134 // cout << subdomains[0].ref << subdomains[1].ref << endl; 3135 Int4 j=0; 3136 for ( i=0;i<NbSubDomains;i++) 3137 if((-subdomains[i].ref) %2) { // good 3138 //cout << " sudom ok = " << i << " " << subdomains[i].ref 3139 // << " " << (-subdomains[i].ref) %2 << endl; 3140 if(i != j) 3141 Exchange(subdomains[i],subdomains[j]); 3142 j++;} 3143 else 3144 { //cout << " remove sub domain " << i << endl; 3145 t= subdomains[i].head; 3146 while (t){// cout << Number(t) << " " << endl; 3147 NbOutT++; 3148 t1=t; 3149 t=t->link; 3150 t1->link=0;}//while (t) 3151 } 3152 3153 if(verbosity>4) 3154 cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl; 3155 NbSubDomains=j; 3156 } 3157 3158 delete [] mark; 3159 3160 } 3161 else 3162 { // find the head for all sub domaine 3163 if (Gh.NbSubDomains != NbSubDomains && subdomains) 3164 delete [] subdomains, subdomains=0; 3165 if (! subdomains ) 3166 subdomains = new SubDomain[ Gh.NbSubDomains]; 3167 NbSubDomains =Gh.NbSubDomains; 3168 if(verbosity>4) 3169 cout << " find the " << NbSubDomains << " sub domain " << endl; 3170 Int4 err=0; 3171 ReMakeTriangleContainingTheVertex(); 3172 Int4 * mark = new Int4[nbt]; 3173 Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge(); 3174 3175 for (it=0;it<nbt;it++) 3176 mark[it]=triangles[it].link ? -1 : -2; 3177 Int4 inew =0; 3178 for (Int4 i=0;i<NbSubDomains;i++) 3179 { 3180 GeometricalEdge &eg = *Gh.subdomains[i].edge; 3181 subdomains[i].ref = Gh.subdomains[i].ref; 3182 int ssdlab = subdomains[i].ref; 3183 // by carefull is not easy to find a edge create from a GeometricalEdge 3184 // see routine MakeGeometricalEdgeToEdge 3185 Edge &e = *GeometricalEdgetoEdge[Gh.Number(eg)]; 3186 assert(&e); 3187 Vertex * v0 = e(0),*v1 = e(1); 3188 Triangle *t = v0->t; 3189 int sens = Gh.subdomains[i].sens; 3190 // test if ge and e is in the same sens 3191 // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref 3192 // << " ref edge =" << eg.ref << " sens " << sens ; 3193 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) 3194 sens = -sens ; 3195 subdomains[i].sens = sens; 3196 subdomains[i].edge = &e; 3197 // cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh " << e[0].r << e[1].r << endl; 3198 // cout << " v0 , v1 = " << Number(v0) << " " << Number(v1) << endl; 3199 assert(t && sens); 3200 3201 TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges 3202 3203 while (1) 3204 { 3205 assert( v0 == ta.EdgeVertex(1) ); 3206 // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl; 3207 if (ta.EdgeVertex(0) == v1) { // ok we find the edge 3208 if (sens>0) 3209 subdomains[i].head=t=Adj(ta); 3210 else 3211 subdomains[i].head=t=ta; 3212 //cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl; 3213 if(t<triangles || t >= triangles+nbt || t->det < 0 3214 || t->link == 0) // Ajoute aout 200 3215 { 3216 cerr << " Error in the def of sub domain "<<i 3217 << " form border " << NbSubDomains - i << "/" << NbSubDomains 3218 << ": Bad sens " << Gh.Number(eg) <<" "<< sens << endl; 3219 err++; 3220 break;} 3221 Int4 it = Number(t); 3222 if (mark[it] >=0) { 3223 if(verbosity>10) 3224 cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref 3225 << " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref 3226 << " skip this def " << endl; 3227 break;} 3228 if(i != inew) 3229 Exchange(subdomains[i],subdomains[inew]); 3230 inew++; 3231 Triangle *tt=t; 3232 Int4 kkk=0; 3233 do 3234 { 3235 kkk++; 3236 assert(mark[Number(tt)]<0); 3237 mark[Number(tt)]=i; 3238 tt=tt->link; 3239 } while (tt!=t); 3240 if(verbosity>7) 3241 cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl; 3242 break;} 3243 ta = Previous(Adj(ta)); 3244 if(t == (Triangle *) ta) { 3245 err++; 3246 cerr << " Error in the def of sub domain " << i 3247 << " edge=" << Gh.Number(eg) << " " << sens << endl; 3248 break;} 3249 // cout << " NB of remove subdomain " << NbSubDomTot-NbSubDomains<< endl; 3250 3251 } 3252 3253 } 3254 if (err) MeshError(777,this); 3255 3256 if (inew < NbSubDomains) { 3257 if (verbosity>5) 3258 cout << " Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl; 3259 NbSubDomains=inew;} 3260 3261 3262 for (it=0;it<nbt;it++) 3263 if ( mark[it] ==-1 ) 3264 NbOutT++,triangles[it].link =0; 3265 delete [] GeometricalEdgetoEdge; 3266 delete [] mark; 3267 3268 } 3269 NbOutT=0; 3270 for (it=0;it<nbt;it++) 3271 if(!triangles[it].link) NbOutT++; 3272 if (verbosity> 4) 3273 cout << " " ; 3274 if (verbosity> 2) 3275 cout << " Nb of Sub borned Domain = " << NbSubDomTot << " NbOutTriangles = " << NbOutT <<endl; 3276 3277 3278 } 3279 /*}}}1*/ 3280 /*FUNCTION Triangles::ReNumberingVertex{{{1*/ 3281 void Triangles::ReNumberingVertex(Int4 * renu) { 3282 // warning be carfull because pointeur 3283 // from on mesh to over mesh 3284 // -- so do ReNumbering a the beginning 3285 Vertex * ve = vertices+nbv; 3286 Int4 it,ie,i; 3287 3288 for ( it=0;it<nbt;it++) 3289 triangles[it].ReNumbering(vertices,ve,renu); 3290 3291 for ( ie=0;ie<nbe;ie++) 3292 edges[ie].ReNumbering(vertices,ve,renu); 3293 3294 for (i=0;i< NbVerticesOnGeomVertex;i++) 3295 { 3296 Vertex *v = VerticesOnGeomVertex[i].mv; 3297 if (v>=vertices && v < ve) 3298 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)]; 3299 } 3300 3301 for (i=0;i< NbVerticesOnGeomEdge;i++) 3302 { 3303 Vertex *v =VerticesOnGeomEdge[i].mv; 3304 if (v>=vertices && v < ve) 3305 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)]; 3306 } 3307 3308 for (i=0;i< NbVertexOnBThVertex;i++) 3309 { 3310 Vertex *v=VertexOnBThVertex[i].v; 3311 if (v>=vertices && v < ve) 3312 VertexOnBThVertex[i].v=vertices+renu[Number(v)]; 3313 } 3314 3315 for (i=0;i< NbVertexOnBThEdge;i++) 3316 { 3317 Vertex *v=VertexOnBThEdge[i].v; 3318 if (v>=vertices && v < ve) 3319 VertexOnBThEdge[i].v=vertices+renu[Number(v)]; 3320 } 3321 3322 // move the Vertices without a copy of the array 3323 // be carefull not trivial code 3324 Int4 j; 3325 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu 3326 if (renu[it] >= 0) // a new sub cycle 3327 { 3328 i=it; 3329 Vertex ti=vertices[i],tj; 3330 while ( (j=renu[i]) >= 0) 3331 { // i is old, and j is new 3332 renu[i] = -1-renu[i]; // mark 3333 tj = vertices[j]; // save new 3334 vertices[j]= ti; // new <- old 3335 i=j; // next 3336 ti = tj; 3337 } 3338 } 3339 if (quadtree) 3340 { delete quadtree; 3341 quadtree = new QuadTree(this); 3342 } 3343 for ( it=0;it<nbv;it++) 3344 renu[i]= -renu[i]-1; 3345 3346 } 3347 /*}}}1*/ 3348 /*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/ 3349 void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) { 3350 long int verbosity=0; 3351 Int4 *renu= new Int4[nbt]; 3352 register Triangle *t0,*t,*te=triangles+nbt; 3353 register Int4 k=0,it,i,j; 3354 3355 for ( it=0;it<nbt;it++) 3356 renu[it]=-1; // outside triangle 3357 for ( i=0;i<NbSubDomains;i++) 3358 { 3359 t=t0=subdomains[i].head; 3360 assert(t0); // no empty sub domain 3361 do { 3362 Int4 kt = Number(t); 3363 assert(kt>=0 && kt < nbt ); 3364 assert(renu[kt]==-1); 3365 renu[kt]=k++; 3366 } 3367 while (t0 != (t=t->link)); 3368 } 3369 if (verbosity>9) 3370 cout << " number of inside triangles " << k << " nbt = " << nbt << endl; 3371 // take is same numbering if possible 3372 if(justcompress) 3373 for ( k=0,it=0;it<nbt;it++) 3374 if(renu[it] >=0 ) 3375 renu[it]=k++; 3376 3377 // put the outside triangles at the end 3378 for ( it=0;it<nbt;it++) 3379 if (renu[it]==-1) 3380 renu[it]=k++; 3381 3382 assert(k == nbt); 3383 // do the change on all the pointeur 3384 for ( it=0;it<nbt;it++) 3385 triangles[it].ReNumbering(triangles,te,renu); 3386 3387 for ( i=0;i<NbSubDomains;i++) 3388 subdomains[i].head=triangles+renu[Number(subdomains[i].head)]; 3389 3390 // move the Triangles without a copy of the array 3391 // be carefull not trivial code 3392 for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu 3393 if (renu[it] >= 0) // a new sub cycle 3394 { 3395 i=it; 3396 Triangle ti=triangles[i],tj; 3397 while ( (j=renu[i]) >= 0) 3398 { // i is old, and j is new 3399 renu[i] = -1; // mark 3400 tj = triangles[j]; // save new 3401 triangles[j]= ti; // new <- old 3402 i=j; // next 3403 ti = tj; 3404 } 3405 } 3406 delete [] renu; 3407 nt = nbt - NbOutT; 3408 3409 } 3410 /*}}}1*/ 3411 /*FUNCTION Triangles::ConsRefTriangle{{{1*/ 3412 Int4 Triangles::ConsRefTriangle(Int4 *reft) const { 3413 long int verbosity=0; 3414 assert(reft); 3415 register Triangle *t0,*t; 3416 register Int4 k=0, num; 3417 for (Int4 it=0;it<nbt;it++) 3418 reft[it]=-1; // outside triangle 3419 for (Int4 i=0;i<NbSubDomains;i++) 3420 { 3421 t=t0=subdomains[i].head; 3422 assert(t0); // no empty sub domain 3423 // register Int4 color=i+1;// because the color 0 is outside triangle 3424 do { k++; 3425 num = Number(t); 3426 assert(num>=0 &&num < nbt); 3427 reft[num]=i; 3428 // cout << Number(t0) << " " <<Number(t)<< " " << i << endl; 3429 } 3430 while (t0 != (t=t->link)); 3431 } 3432 // NbOutT = nbt - k; 3433 if (verbosity>5) 3434 cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k 3435 << " Nbt = " << nbt << " Out Triangles = " << nbt - k << endl; 3436 3437 return k; 3438 3439 } 3440 /*}}}1*/ 3441 /*FUNCTION Triangles::GeomToTriangles1{{{1*/ 3442 void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices) { 3443 Gh.NbRef++;// add a ref to Gh 3444 3445 3446 /************************************************************************* 3447 // methode in 2 step 3448 // 1 - compute the number of new edge to allocate 3449 // 2 - construct the edge 3450 remark: 3451 in this part we suppose to have a background mesh with the same 3452 geometry 3453 3454 To construct the discretisation of the new mesh we have to 3455 rediscretize the boundary of background Mesh 3456 because we have only the pointeur from the background mesh to the geometry. 3457 We need the abcisse of the background mesh vertices on geometry 3458 so a vertex is 3459 0 on GeometricalVertex ; 3460 1 on GeometricalEdge + abcisse 3461 2 internal 3462 3463 *************************************************************************/ 3464 assert(&BTh.Gh == &Gh); 3465 3466 long int verbosity=0; 3467 BTh.NbRef++; // add a ref to BackGround Mesh 3468 PreInit(inbvx); 3469 BTh.SetVertexFieldOn(); 3470 int * bcurve = new int[Gh.NbOfCurves]; // 3471 3472 // we have 2 ways to make the loop 3473 // 1) on the geometry 3474 // 2) on the background mesh 3475 // if you do the loop on geometry, we don't have the pointeur on background, 3476 // and if you do the loop in background we have the pointeur on geometry 3477 // so do the walk on background 3478 // Int4 NbVerticesOnGeomVertex; 3479 // VertexOnGeom * VerticesOnGeomVertex; 3480 // Int4 NbVerticesOnGeomEdge; 3481 // VertexOnGeom * VerticesOnGeomEdge; 3482 3483 3484 NbVerticesOnGeomVertex=0; 3485 NbVerticesOnGeomEdge=0; 3486 //1 copy of the Required vertex 3487 int i; 3488 for ( i=0;i<Gh.nbv;i++) 3489 if (Gh[i].Required()) NbVerticesOnGeomVertex++; 3490 3491 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 3492 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex]; 3493 // 3494 if( NbVerticesOnGeomVertex >= nbvx) 3495 { 3496 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl; 3497 MeshError(1,this); 3498 } 3499 assert(vertices); 3500 for (i=0;i<Gh.nbv;i++) 3501 if (Gh[i].Required()) {//Gh vertices Required 3502 vertices[nbv] = Gh[i]; 3503 vertices[nbv].i = I2(0,0); 3504 Gh[i].to = vertices + nbv;// save Geom -> Th 3505 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); 3506 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl; 3507 nbv++;} 3508 else Gh[i].to=0; 3509 // 3510 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++) 3511 { 3512 VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i]; 3513 if (vog.IsRequiredVertex()) { 3514 GeometricalVertex * gv = vog; 3515 Vertex *bv = vog; 3516 assert(gv->to);// use of Geom -> Th 3517 VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv); 3518 gv->to->m = bv->m; // for taking the metrix of the background mesh 3519 ;} 3520 } 3521 assert(NbVertexOnBThVertex == NbVerticesOnGeomVertex); 3522 // new stuff FH with curve 3523 // find the begin of the curve in BTh 3524 { 3525 Gh.UnMarkEdges(); 3526 int bfind=0; 3527 /* 3528 cout << " nb curves = " << Gh.NbOfCurves << endl; 3529 for(int i=0;i<Gh.NbOfCurves ;i++) 3530 { 3531 cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb 3532 << " end e= " << Gh.Number(Gh.curves[i].ee) << " k=" << Gh.curves[i].ke << endl; 3533 }*/ 3534 for (int i=0;i<Gh.NbOfCurves;i++) 3535 { 3536 bcurve[i]=-1; 3537 } 3538 3539 for (int iedge=0;iedge<BTh.nbe;iedge++) 3540 { 3541 Edge & ei = BTh.edges[iedge]; 3542 for(int je=0;je<2;je++) // for the 2 extremites 3543 if (!ei.on->Mark() && ei[je].on->IsRequiredVertex() ) 3544 { 3545 // a begin of curve 3546 int nc = ei.on->CurveNumber; 3547 3548 //cout << "curve " << nc << " v " << Gh.Number((GeometricalVertex *) *ei[je].on) << " " 3549 // << " e " << " " << Gh.Number(ei.on) << " vc " << Gh.Number((*Gh.curves[nc].be)[Gh.curves[nc].kb]) << endl; 3550 if( 3551 ei.on==Gh.curves[nc].be && 3552 (GeometricalVertex *) *ei[je].on == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity 3553 ) 3554 { 3555 // cout << " find " << endl; 3556 bcurve[nc]=iedge*2+je; 3557 bfind++; 3558 } 3559 } 3560 } 3561 assert( bfind==Gh.NbOfCurves); 3562 } 3563 // method in 2 + 1 step 3564 // 0.0) compute the length and the number of vertex to do allocation 3565 // 1.0) recompute the length 3566 // 1.1) compute the vertex 3567 Int4 nbex=0,NbVerticesOnGeomEdgex=0; 3568 for (int step=0; step <2;step++) 3569 { 3570 Int4 NbOfNewPoints=0; 3571 Int4 NbOfNewEdge=0; 3572 Int4 iedge; 3573 Gh.UnMarkEdges(); 3574 /* add Curve loop FH 3575 // find a starting back groud edges to walk 3576 for (iedge=0;iedge<BTh.nbe;iedge++) { 3577 Edge & ei = BTh.edges[iedge]; 3578 for(int jedge=0;jedge<2;jedge++) // for the 2 extremites 3579 if (!ei.on->Mark() && ei[jedge].on->IsRequiredVertex() ) 3580 { 3581 */ 3582 // new code FH 2004 3583 Real8 L=0; 3584 for (int icurve=0;icurve<Gh.NbOfCurves;icurve++) 3585 { 3586 iedge=bcurve[icurve]/2; 3587 int jedge=bcurve[icurve]%2; 3588 if( ! Gh.curves[icurve].master) continue; // we skip all equi curve 3589 Edge & ei = BTh.edges[iedge]; 3590 // warning: ei.on->Mark() can be change in 3591 // loop for(jedge=0;jedge<2;jedge++) 3592 // new curve 3593 // good the find a starting edge 3594 Real8 Lstep=0,Lcurve=0;// step between two points (phase==1) 3595 Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 3596 3597 3598 // cout.precision(16); 3599 for(int phase=0;phase<=step;phase++) 3600 { 3601 3602 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next) 3603 { 3604 3605 int icurveequi= Gh.Number(curve); 3606 3607 if( phase == 0 && icurveequi != icurve) continue; 3608 3609 int k0=jedge,k1; 3610 Edge * pe= BTh.edges+iedge; 3611 //GeometricalEdge *ong = ei.on; 3612 int iedgeequi=bcurve[icurveequi]/2; 3613 int jedgeequi=bcurve[icurveequi]%2; 3614 3615 int k0equi=jedgeequi,k1equi; 3616 Edge * peequi= BTh.edges+iedgeequi; 3617 GeometricalEdge *ongequi = peequi->on; 3618 3619 Real8 sNew=Lstep;// abcisse of the new points (phase==1) 3620 L=0;// length of the curve 3621 Int4 i=0;// index of new points on the curve 3622 register GeometricalVertex * GA0 = *(*peequi)[k0equi].on; 3623 Vertex *A0; 3624 A0 = GA0->to; // the vertex in new mesh 3625 Vertex *A1; 3626 VertexOnGeom *GA1; 3627 Edge * PreviousNewEdge = 0; 3628 // cout << " --------------New Curve phase " << phase 3629 // << "---------- A0=" << *A0 << ei[k0] <<endl; 3630 assert (A0-vertices>=0 && A0-vertices <nbv); 3631 if(ongequi->Required() ) 3632 { 3633 GeometricalVertex *GA1 = *(*peequi)[1-k0equi].on; 3634 A1 = GA1->to; // 3635 } 3636 else 3637 for(;;) 3638 { 3639 // assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe); 3640 Edge &ee=*pe; 3641 Edge &eeequi=*peequi; 3642 k1 = 1-k0; // next vertex of the edge 3643 k1equi= 1 - k0equi; 3644 3645 assert(pe && ee.on); 3646 ee.on->SetMark(); 3647 Vertex & v0=ee[0], & v1=ee[1]; 3648 R2 AB= (R2) v1 - (R2) v0; 3649 Real8 L0=L,LAB; 3650 LAB = LengthInterpole(v0.m,v1.m,AB); 3651 L+= LAB; 3652 if (phase) {// computation of the new points 3653 while ((i!=NbCreatePointOnCurve) && sNew <= L) { 3654 // cout << " L0= " << L0 << " L " << L << " sN=" 3655 // << sNew << " LAB=" << LAB << " NBPC =" <<NbCreatePointOnCurve<< " i " << i << endl; 3656 assert (sNew >= L0); 3657 assert(LAB); 3658 3659 3660 assert(vertices && nbv<nbvx); 3661 assert(edges && nbe < nbex); 3662 assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge < NbVerticesOnGeomEdgex); 3663 // new vertex on edge 3664 A1=vertices+nbv++; 3665 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge; 3666 Edge *e = edges + nbe++; 3667 Real8 se= (sNew-L0)/LAB; 3668 assert(se>=0 && se < 1.000000001); 3669 se = abscisseInterpole(v0.m,v1.m,AB,se,1); 3670 assert(se>=0 && se <= 1); 3671 //((k1==1) != (k1==k1equi)) 3672 se = k1 ? se : 1. - se; 3673 se = k1==k1equi ? se : 1. - se; 3674 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 3675 ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 3676 A1->ReferenceNumber = eeequi.ref; 3677 A1->DirOfSearch =NoDirOfSearch; 3678 //cout << icurveequi << " " << i << " " << *A1 << endl; 3679 e->on = ongequi; 3680 e->v[0]= A0; 3681 e->v[1]= A1; 3682 if(verbosity>99) 3683 cout << i << "+ New P "<< nbv-1 << " " <<sNew<< " L0=" << L0 3684 << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= " 3685 << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl; 3686 3687 e->ref = eeequi.ref; 3688 e->adj[0]=PreviousNewEdge; 3689 3690 if (PreviousNewEdge) 3691 PreviousNewEdge->adj[1] = e; 3692 PreviousNewEdge = e; 3693 A0=A1; 3694 sNew += Lstep; 3695 // cout << " sNew = " << sNew << " L = " << L 3696 // << " ------" <<NbCreatePointOnCurve << " " << i << endl; 3697 if (++i== NbCreatePointOnCurve) break; 3698 } 3699 3700 } 3701 assert(ee.on->CurveNumber==ei.on->CurveNumber); 3702 if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl; 3703 if ( ee[k1].on->IsRequiredVertex()) { 3704 assert(eeequi[k1equi].on->IsRequiredVertex()); 3705 register GeometricalVertex * GA1 = *eeequi[k1equi].on; 3706 A1=GA1->to;// the vertex in new mesh 3707 assert (A1-vertices>=0 && A1-vertices <nbv); 3708 break;} 3709 if (!ee.adj[k1]) 3710 {cerr << "Error adj edge " << BTh.Number(ee) << ", nbe = " << nbe 3711 << " Gh.vertices " << Gh.vertices 3712 << " k1 = " << k1 << " on=" << *ee[k1].on << endl; 3713 cerr << ee[k1].on->gv-Gh.vertices << endl; 3714 } 3715 pe = ee.adj[k1]; // next edge 3716 k0 = pe->Intersection(ee); 3717 peequi= eeequi.adj[k1equi]; // next edge 3718 k0equi=peequi->Intersection(eeequi); 3719 }// for(;;) end of the curve 3720 3721 3722 if (phase) // construction of the last edge 3723 { 3724 Edge *e = edges + nbe++; 3725 if (verbosity>10) 3726 cout << " Fin curve A1" << *A1 << " " << icurve << " <=> " << icurveequi <<"-----" << 3727 NbCreatePointOnCurve << " == " <<i<<endl; 3728 e->on = ongequi; 3729 e->v[0]= A0; 3730 e->v[1]= A1; 3731 e->ref = peequi->ref; 3732 e->adj[0]=PreviousNewEdge; 3733 e->adj[1]=0; 3734 if (PreviousNewEdge) 3735 PreviousNewEdge->adj[1] = e; 3736 PreviousNewEdge = e; 3737 3738 assert(i==NbCreatePointOnCurve); 3739 3740 } 3741 } // end loop on equi curve 3742 3743 if (!phase) { // 3744 Int4 NbSegOnCurve = Max((Int4)(L+0.5),(Int4) 1);// nb of seg 3745 Lstep = L/NbSegOnCurve; 3746 Lcurve = L; 3747 NbCreatePointOnCurve = NbSegOnCurve-1; 3748 3749 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next) 3750 { 3751 NbOfNewEdge += NbSegOnCurve; 3752 NbOfNewPoints += NbCreatePointOnCurve; 3753 } 3754 if(verbosity>5) 3755 cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep=" 3756 << Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve <<endl; 3757 // do'nt 3758 // if(NbCreatePointOnCurve<1) break; 3759 } 3760 }//for(phase;;) 3761 /* modif FH add Curve class 3762 }}//for (iedge=0;iedge<BTh.nbe;iedge++) 3763 */ 3764 // new code Curve class 3765 } // end of curve loop 3766 // end new code 3767 // do the allocation 3768 if(step==0) 3769 { 3770 //if(!NbOfNewPoints) break;// nothing ????? bug 3771 if(nbv+NbOfNewPoints > nbvx) 3772 { 3773 cerr << " Too much vertices on geometry " << nbv+NbOfNewPoints << " >= " << nbvx << endl; 3774 MeshError(3,this); 3775 } 3776 //cout << " NbOfNewEdge" << NbOfNewEdge << " NbOfNewPoints " << NbOfNewPoints << endl; 3777 edges = new Edge[NbOfNewEdge]; 3778 nbex = NbOfNewEdge; 3779 if(NbOfNewPoints) { // 3780 VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints]; 3781 NbVertexOnBThEdge =NbOfNewPoints; 3782 VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints]; 3783 NbVerticesOnGeomEdgex = NbOfNewPoints; } 3784 NbOfNewPoints =0; 3785 NbOfNewEdge = 0; 3786 } 3787 } // for(step;;) 3788 assert(nbe); 3789 3790 delete [] bcurve; 3791 3792 3793 Insert(); 3794 ForceBoundary(); 3795 FindSubDomain(); 3796 3797 NewPoints(BTh,KeepBackVertices) ; 3798 CurrentTh = 0; 3799 } 3800 /*}}}1*/ 3801 /*FUNCTION Triangles::GeomToTriangles0{{{1*/ 3802 void Triangles::GeomToTriangles0(Int4 inbvx) { 3803 Gh.NbRef++;// add a ref to GH 3804 3805 3806 Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve; 3807 Real8 lcurve, lstep,s; 3808 3809 R2 AB; 3810 GeometricalVertex *a,*b; 3811 Vertex *va,*vb; 3812 GeometricalEdge * e; 3813 PreInit(inbvx); 3814 int background = &BTh != this; 3815 // int SameGeom = background && (&BTh.Gh == &Gh); 3816 nbv = 0; 3817 NbVerticesOnGeomVertex=0; 3818 NbVerticesOnGeomEdge=0; 3819 for (i=0;i<Gh.nbv;i++) 3820 if (Gh[i].Required() && Gh[i].IsThe() ) NbVerticesOnGeomVertex++; 3821 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 3822 // 3823 if( NbVerticesOnGeomVertex >= nbvx) 3824 { 3825 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl; 3826 MeshError(1,this); 3827 } 3828 for (i=0;i<Gh.nbv;i++) 3829 if (Gh[i].Required()&& Gh[i].IsThe() ) {//Gh vertices Required 3830 if (nbv < nbvx) 3831 vertices[nbv] = Gh[i]; 3832 Gh[i].to = vertices + nbv;// save Geom -> Th 3833 VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]); 3834 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl; 3835 nbv++; 3836 } 3837 // assert( Gh.nbv < nbvx); 3838 3839 // Method in 2 step: 0 and 1 3840 // 1) compute de nb of edge 3841 // 2) construct the edge 3842 // generation of the curves 3843 assert(! edges); 3844 // 2 step 3845 // --step=0 to compute the number of edges + alloc at end 3846 // --step=1 to construct the edges 3847 for (int step=0;step<2;step++) 3848 {// for (int step=0;step<2;step++) 3849 Int4 nbex = 0; 3850 nbe = 0; 3851 Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; 3852 // cout << " -------------- step =" << step << endl; 3853 Gh.UnMarkEdges(); 3854 NbOfCurves = 0; 3855 for (i=0;i<Gh.nbe;i++) { 3856 GeometricalEdge & ei = Gh.edges[i]; 3857 if (!ei.Dup()) // a good curve (not dup ) 3858 for(int j=0;j<2;j++) 3859 if (!ei.Mark() && ei[j].Required()) { 3860 // warning ei.Mark() can be change in loop for(j=0;j<2;j++) 3861 // cout << " New curve = " << NbOfCurves << endl; 3862 Int4 nbvend =0; 3863 3864 Edge * PreviousNewEdge=0; 3865 3866 lstep = -1;//to do not create points 3867 if(ei.Required()) 3868 { 3869 if (j==0) 3870 if(step==0) 3871 nbe++; 3872 else 3873 { 3874 e = & ei; 3875 a=ei(0)->The(); 3876 b=ei(1)->The(); 3877 assert(edges); 3878 edges[nbe].v[0]=a->to; 3879 edges[nbe].v[1]=b->to;; 3880 edges[nbe].ref = e->ref; 3881 edges[nbe].on = e; 3882 edges[nbe].adj[0] = 0; 3883 edges[nbe].adj[1] = 0; 3884 nbe++;} 3885 } 3886 else 3887 { // on curve ------ 3888 for ( int kstep=0;kstep<= step;kstep++) 3889 { // begin for ( int kstep=0;kstep<= step;kstep++) 3890 // if 2nd step where 2 step 3891 // -- 1 compute le length of the curve 3892 // -- create the points and edge 3893 PreviousNewEdge=0; 3894 NbNewPoints=0; 3895 NbEdgeCurve=0; 3896 assert(nbvend < nbvx); 3897 lcurve =0; 3898 s = lstep; 3899 int k=j; 3900 e = & ei; 3901 a=ei(k)->The(); 3902 va = a->to; 3903 e->SetMark(); 3904 // cout << " New curve " ; 3905 3906 // if SameGeo We have go in the background geometry 3907 // to find the discretisation of the curve 3908 3909 for(;;) 3910 { 3911 k = 1-k; 3912 b= (*e)(k)->The(); 3913 AB = b->r - a->r; 3914 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; 3915 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; 3916 Real8 ledge = (MA(AB) + MB(AB))/2; 3917 // 3918 const int MaxSubEdge = 10; 3919 int NbSubEdge = 1; 3920 Real8 lSubEdge[MaxSubEdge]; 3921 R2 A,B; 3922 if (ledge < 1.5) 3923 lSubEdge[0] = ledge; 3924 else { 3925 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); 3926 A= a->r; 3927 Metric MAs =MA,MBs; 3928 // cout << " lSubEdge old=" << ledge 3929 // << " new " << A << MA << endl; 3930 ledge = 0; 3931 Real8 x =0, xstep= 1. / NbSubEdge; 3932 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { 3933 x += xstep; 3934 B = e->F(k ? x : 1-x); 3935 MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB); 3936 AB = A-B; 3937 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); 3938 // cout << " " << lSubEdge[kk] << " x " << x 3939 // << " " << A << B << MA << MB<< endl ; 3940 } 3941 // cout << endl; 3942 } 3943 3944 Real8 lcurveb = lcurve+ ledge ; 3945 while (lcurve<=s && s <= lcurveb && nbv < nbvend) 3946 { 3947 // New points 3948 3949 // Real8 aa=(lcurveb-s)/ledge; 3950 // Real8 bb=(s-lcurve)/ledge; 3951 3952 Real8 ss = s-lcurve; 3953 // 1) find the SubEdge containing ss by dichotomie 3954 int kk0=-1,kk1=NbSubEdge-1,kkk; 3955 Real8 ll0=0,ll1=ledge,llk; 3956 while (kk1-kk0>1) 3957 { 3958 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) 3959 kk1=kkk,ll1=llk; 3960 else 3961 kk0=kkk,ll0=llk;} 3962 assert(kk1 != kk0); 3963 3964 Real8 sbb = (ss-ll0 )/(ll1-ll0); 3965 Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb; 3966 3967 // new vertex on edge 3968 vb = &vertices[nbv++]; 3969 vb->m = Metric(aa,a->m,bb,b->m); 3970 vb->ReferenceNumber = e->ref; 3971 vb->DirOfSearch =NoDirOfSearch; 3972 Real8 abcisse = k ? bb : aa; 3973 vb->r = e->F( abcisse ); 3974 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); 3975 3976 // to take in account the sens of the edge 3977 3978 s += lstep; 3979 edges[nbe].v[0]=va; 3980 edges[nbe].v[1]=vb; 3981 edges[nbe].ref = e->ref; 3982 edges[nbe].on = e; 3983 edges[nbe].adj[0] = PreviousNewEdge; 3984 if(PreviousNewEdge) 3985 PreviousNewEdge->adj[1] = &edges[nbe]; 3986 PreviousNewEdge = edges + nbe; 3987 nbe++; 3988 va = vb; 3989 } 3990 lcurve = lcurveb; 3991 e->SetMark(); 3992 // cout << e-Gh.edges << ", " << k << " " 3993 // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" " 3994 // << lcurve<< ";; " ; 3995 a=b; 3996 if (b->Required() ) break; 3997 int kprev=k; 3998 k = e->SensAdj[kprev];// next vertices 3999 e = e->Adj[kprev]; 4000 assert(e); 4001 }// for(;;) 4002 vb = b->to; 4003 // cout << endl; 4004 NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1); 4005 NbNewPoints = NbEdgeCurve-1; 4006 if(!kstep) 4007 { NbVerticesOnGeomEdge0 += NbNewPoints; 4008 NbOfCurves++;} 4009 4010 nbvend=nbv+NbNewPoints; 4011 4012 lstep = lcurve / NbEdgeCurve; 4013 // cout <<"lstep " << lstep << " lcurve " 4014 // << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " <<NbVerticesOnGeomEdge0<<" " <<NbVerticesOnGeomEdge<<" step =" <<step<< endl; 4015 } 4016 // end of curve -- 4017 if (edges) { // last edges of the curves 4018 edges[nbe].v[0]=va; 4019 edges[nbe].v[1]=vb; 4020 edges[nbe].ref = e->ref; 4021 edges[nbe].on = e; 4022 edges[nbe].adj[0] = PreviousNewEdge; 4023 edges[nbe].adj[1] = 0; 4024 if(PreviousNewEdge) 4025 PreviousNewEdge->adj[1] = & edges[nbe]; 4026 4027 nbe++;} 4028 else 4029 nbe += NbEdgeCurve; 4030 } // end on curve --- 4031 } // if (edges[i][j].Corner()) 4032 } // for (i=0;i<nbe;i++) 4033 if(!step) { 4034 // cout << "edges " << edges << " VerticesOnGeomEdge " <<VerticesOnGeomEdge << endl; 4035 assert(!edges); 4036 assert(!VerticesOnGeomEdge); 4037 edges = new Edge[nbex=nbe]; 4038 if(NbVerticesOnGeomEdge0) 4039 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0]; 4040 assert(edges); 4041 assert(VerticesOnGeomEdge || NbVerticesOnGeomEdge0 ==0); 4042 // do the vertex on a geometrical vertex 4043 NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge; 4044 } 4045 else 4046 assert(NbVerticesOnGeomEdge == NbVerticesOnGeomEdge0); 4047 // cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe 4048 // << "== " << nbex << " nbv = " << nbv << endl; 4049 assert(nbex=nbe); 4050 } // for (step=0;step<2;step++) 4051 4052 Insert(); 4053 ForceBoundary(); 4054 FindSubDomain(); 4055 4056 // NewPointsOld(*this) ; 4057 NewPoints(*this,0) ; 4058 CurrentTh = 0; 4059 } 4060 /*}}}1*/ 4061 /*FUNCTION Triangles::MakeGeometricalEdgeToEdge{{{1*/ 4062 Edge** Triangles::MakeGeometricalEdgeToEdge() { 4063 assert(Gh.nbe); 4064 Edge **e= new (Edge* [Gh.nbe]); 4065 4066 Int4 i; 4067 for ( i=0;i<Gh.nbe ; i++) 4068 e[i]=NULL; 4069 for ( i=0;i<nbe ; i++) 4070 { 4071 Edge * ei = edges+i; 4072 GeometricalEdge *on = ei->on; 4073 e[Gh.Number(on)] = ei; 4074 } 4075 for ( i=0;i<nbe ; i++) 4076 for (int ii=0;ii<2;ii++) { 4077 Edge * ei = edges+i; 4078 GeometricalEdge *on = ei->on; 4079 int j= ii; 4080 while (!(*on)[j].Required()) { 4081 // cout << i << " " << ii << " j= " << j << " curve = " 4082 // << on->CurveNumber << " " << Gh.Number(on) << " on " << j 4083 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) 4084 // << " -> " ; 4085 Adj(on,j); // next geom edge 4086 j=1-j; 4087 // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)] 4088 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl; 4089 if (e[Gh.Number(on)]) break; // optimisation 4090 e[Gh.Number(on)] = ei; 4091 } 4092 } 4093 4094 int kk=0; 4095 for ( i=0;i<Gh.nbe ; i++) 4096 if (!e[i]) 4097 if(kk++<10) { 4098 cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber 4099 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl; 4100 // assert( e[i]); 4101 } 4102 if(kk) MeshError(997,this); 4103 4104 return e; 4105 } 4106 /*}}}1*/ 4107 /*FUNCTION Triangles::SetIntCoor{{{1*/ 4108 void Triangles::SetIntCoor(const char * strfrom) { 4109 pmin = vertices[0].r; 4110 pmax = vertices[0].r; 4111 4112 // recherche des extrema des vertices pmin,pmax 4113 Int4 i; 4114 for (i=0;i<nbv;i++) { 4115 pmin.x = Min(pmin.x,vertices[i].r.x); 4116 pmin.y = Min(pmin.y,vertices[i].r.y); 4117 pmax.x = Max(pmax.x,vertices[i].r.x); 4118 pmax.y = Max(pmax.y,vertices[i].r.y); 4119 } 4120 R2 DD = (pmax-pmin)*0.05; 4121 pmin = pmin-DD; 4122 pmax = pmax+DD; 4123 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 4124 assert(coefIcoor >0); 4125 4126 // generation of integer coord 4127 for (i=0;i<nbv;i++) { 4128 vertices[i].i = toI2(vertices[i].r); 4129 } 4130 4131 // computation of the det 4132 int Nberr=0; 4133 for (i=0;i<nbt;i++) 4134 { 4135 Vertex & v0 = triangles[i][0]; 4136 Vertex & v1 = triangles[i][1]; 4137 Vertex & v2 = triangles[i][2]; 4138 if ( &v0 && &v1 && &v2 ) // a good triangles; 4139 { 4140 triangles[i].det= det(v0,v1,v2); 4141 if (triangles[i].det <=0 && Nberr++ <10) 4142 { 4143 if(Nberr==1) 4144 if (strfrom) 4145 cerr << "+++ Fatal Error " << strfrom << "(SetInCoor) Error : area of Triangle < 0 " << endl; 4146 else 4147 cerr << "+++ Fatal Error Triangle (in SetInCoor) area of Triangle < 0" << endl; 4148 cerr << " Triangle " << i << " det (I2) = " << triangles[i].det ; 4149 cerr << " (R2) " << Det(v1.r-v0.r,v2.r-v0.r); 4150 cerr << "; The 3 vertices " << endl; 4151 cerr << Number(v0) << " " << Number(v1) << " " 4152 << Number(v2) << " : " ; 4153 cerr << v0.r << v1.r << v2.r << " ; "; 4154 cerr << v0.i << v1.i << v2.i << endl; 4155 } 4156 } 4157 else 4158 triangles[i].det= -1; // Null triangle; 4159 } 4160 if (Nberr) MeshError(899,this); 4161 4162 } 4163 /*}}}1*/ 4164 /*FUNCTION Triangles::FillHoleInMesh{{{1*/ 4165 void Triangles::FillHoleInMesh() { 4166 4167 Triangles* OldCurrentTh =CurrentTh; 4168 CurrentTh=this; 4169 4170 int verbosity=0; 4171 4172 // generation of the integer coordinate 4173 { 4174 4175 // find extrema coordinates of vertices pmin,pmax 4176 Int4 i; 4177 if(verbosity>2) printf(" Filling holes in mesh of %i vertices\n",nbv); 4178 4179 //initialize ordre 4180 assert(ordre); 4181 for (i=0;i<nbv;i++) ordre[i]=0; 4182 4183 NbSubDomains =0; 4184 4185 /* generation of the adjacence of the triangles*/ 4186 4187 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 4188 4189 //initialize st 4190 Int4* st = new Int4[nbt*3]; 4191 for (i=0;i<nbt*3;i++) st[i]=-1; 4192 4193 //check number of edges 4194 Int4 kk =0; 4195 for (i=0;i<nbe;i++){ 4196 kk=kk+(i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1]))); 4197 } 4198 if (kk != nbe) { 4199 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe)); 4200 } 4201 4202 // 4203 for (i=0;i<nbt;i++){ 4204 for (int j=0;j<3;j++) { 4205 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), 4206 Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 4207 Int4 invisible = triangles[i].Hidden(j); 4208 if(st[k]==-1){ 4209 st[k]=3*i+j; 4210 } 4211 else if(st[k]>=0) { 4212 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))); 4213 4214 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); 4215 if (invisible) triangles[i].SetHidden(j); 4216 if (k<nbe) { 4217 triangles[i].SetLocked(j); 4218 } 4219 st[k]=-2-st[k]; 4220 } 4221 else { 4222 throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles", 4223 Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]))); 4224 } 4225 } 4226 } 4227 if(verbosity>5) { 4228 printf(" info on Mesh %s:\n",name); 4229 printf(" - number of vertices = %i \n",nbv); 4230 printf(" - number of triangles = %i \n",nbt); 4231 printf(" - number of given edges = %i \n",nbe); 4232 printf(" - number of all edges = %i \n" ,edge4->nb()); 4233 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); 4234 } 4235 4236 // check the consistant of edge[].adj and the geometrical required vertex 4237 Int4 k=0; 4238 for (i=0;i<edge4->nb();i++){ 4239 if (st[i] >=0){ // edge alone 4240 if (i < nbe) { 4241 Int4 i0=edge4->i(i); 4242 ordre[i0] = vertices+i0; 4243 Int4 i1=edge4->j(i); 4244 ordre[i1] = vertices+i1; 4245 } 4246 else { 4247 k=k+1; 4248 if (k <20) { 4249 throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i))); 4250 } 4251 } 4252 } 4253 } 4254 if(k != 0) { 4255 throw ErrorException(__FUNCT__,exprintf("%i boundary edges are not defined as edges",k)); 4256 } 4257 4258 /* mesh generation with boundary points*/ 4259 Int4 nbvb = 0; 4260 for (i=0;i<nbv;i++){ 4261 vertices[i].t=0; 4262 vertices[i].vint=0; 4263 if (ordre[i]){ 4264 ordre[nbvb++] = ordre[i]; 4265 } 4266 } 4267 4268 Triangle* savetriangles= triangles; 4269 Int4 savenbt=nbt; 4270 Int4 savenbtx=nbtx; 4271 SubDomain * savesubdomains = subdomains; 4272 subdomains = 0; 4273 4274 Int4 Nbtriafillhole = 2*nbvb; 4275 Triangle* triafillhole =new Triangle[Nbtriafillhole]; 4276 triangles = triafillhole; 4277 4278 nbt=2; 4279 nbtx= Nbtriafillhole; 4280 4281 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 4282 if ( ++i >= nbvb) { 4283 cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl; 4284 MeshError(998,this); } 4285 Exchange( ordre[2], ordre[i]); 4286 4287 Vertex * v0=ordre[0], *v1=ordre[1]; 4288 4289 4290 triangles[0](0) = 0; // sommet pour infini 4291 triangles[0](1) = v0; 4292 triangles[0](2) = v1; 4293 4294 triangles[1](0) = 0;// sommet pour infini 4295 triangles[1](2) = v0; 4296 triangles[1](1) = v1; 4297 const int e0 = OppositeEdge[0]; 4298 const int e1 = NextEdge[e0]; 4299 const int e2 = PreviousEdge[e0]; 4300 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 4301 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 4302 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 4303 4304 triangles[0].det = -1; // faux triangles 4305 triangles[1].det = -1; // faux triangles 4306 4307 triangles[0].SetTriangleContainingTheVertex(); 4308 triangles[1].SetTriangleContainingTheVertex(); 4309 4310 triangles[0].link=&triangles[1]; 4311 triangles[1].link=&triangles[0]; 4312 4313 if (!quadtree ) 4314 delete quadtree; // ->ReInitialise(); 4315 4316 quadtree = new QuadTree(this,0); 4317 quadtree->Add(*v0); 4318 quadtree->Add(*v1); 4319 4320 // We add the vertices one by one 4321 Int4 NbSwap=0; 4322 for (Int4 icount=2; icount<nbvb; icount++) { 4323 Vertex *vi = ordre[icount]; 4324 // cout << " Add vertex " << Number(vi) << endl; 4325 Icoor2 dete[3]; 4326 Triangle *tcvi = FindTriangleContening(vi->i,dete); 4327 quadtree->Add(*vi); 4328 Add(*vi,tcvi,dete); 4329 NbSwap += vi->Optim(1,1); 4330 } 4331 4332 // inforce the boundary 4333 TriangleAdjacent ta(0,0); 4334 Int4 nbloss = 0,knbe=0; 4335 for ( i = 0; i < nbe; i++){ 4336 if (st[i] >=0){ // edge alone => on border ... FH oct 2009 4337 Vertex & a=edges[i][0], & b = edges[i][1]; 4338 if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1. 4339 { 4340 knbe++; 4341 if (ForceEdge(a,b,ta)<0) 4342 nbloss++; 4343 } 4344 } 4345 } 4346 if(nbloss) { 4347 cerr << " we loss some " << nbloss << " " << " edges other " << knbe << endl; 4348 MeshError(1100,this); 4349 } 4350 4351 FindSubDomain(1); 4352 // remove all the hole 4353 // remove all the good sub domain 4354 Int4 krm =0; 4355 for (i=0;i<nbt;i++){ 4356 if (triangles[i].link){ // remove triangles 4357 krm++; 4358 for (int j=0;j<3;j++) 4359 { 4360 TriangleAdjacent ta = triangles[i].Adj(j); 4361 Triangle & tta = * (Triangle *) ta; 4362 if(! tta.link) // edge between remove and not remove 4363 { // change the link of ta; 4364 int ja = ta; 4365 Vertex *v0= ta.EdgeVertex(0); 4366 Vertex *v1= ta.EdgeVertex(1); 4367 Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv); 4368 assert(st[k] >=0); 4369 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); 4370 ta.SetLock(); 4371 st[k]=-2-st[k]; 4372 } 4373 } 4374 } 4375 } 4376 Int4 NbTfillHoll =0; 4377 for (i=0;i<nbt;i++){ 4378 if (triangles[i].link) { 4379 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL); 4380 triangles[i].color=-1; 4381 } 4382 else 4383 { 4384 triangles[i].color= savenbt+ NbTfillHoll++; 4385 } 4386 } 4387 // cout << savenbt+NbTfillHoll << " " << savenbtx << endl; 4388 assert(savenbt+NbTfillHoll <= savenbtx ); 4389 // copy of the outside triangles in saveTriangles 4390 for (i=0;i<nbt;i++){ 4391 if(triangles[i].color>=0) { 4392 savetriangles[savenbt]=triangles[i]; 4393 savetriangles[savenbt].link=0; 4394 savenbt++; 4395 } 4396 } 4397 // gestion of the adj 4398 k =0; 4399 Triangle * tmax = triangles + nbt; 4400 for (i=0;i<savenbt;i++) 4401 { 4402 Triangle & ti = savetriangles[i]; 4403 for (int j=0;j<3;j++) 4404 { 4405 Triangle * ta = ti.TriangleAdj(j); 4406 int aa = ti.NuEdgeTriangleAdj(j); 4407 int lck = ti.Locked(j); 4408 if (!ta) k++; // bug 4409 else if ( ta >= triangles && ta < tmax) 4410 { 4411 ta= savetriangles + ta->color; 4412 ti.SetAdj2(j,ta,aa); 4413 if(lck) ti.SetLocked(j); 4414 } 4415 } 4416 } 4417 // OutSidesTriangles = triangles; 4418 // Int4 NbOutSidesTriangles = nbt; 4419 4420 // restore triangles; 4421 nbt=savenbt; 4422 nbtx=savenbtx; 4423 delete [] triangles; 4424 delete [] subdomains; 4425 triangles = savetriangles; 4426 subdomains = savesubdomains; 4427 // cout << triangles << " <> " << OutSidesTriangles << endl; 4428 /* k=0; 4429 for (i=0;i<nbt;i++) 4430 for (int j=0;j<3;j++) 4431 if (!triangles[i].TriangleAdj(j)) 4432 k++; 4433 */ 4434 if (k) { 4435 cerr << "Error Nb of triangles edge alone = " << k << endl; 4436 MeshError(9997,this); 4437 } 4438 FindSubDomain(); 4439 // cout << " NbTOld = " << NbTold << " == " << nbt - NbOutT << " " << nbt << endl; 4440 4441 delete edge4; 4442 delete [] st; 4443 for (i=0;i<nbv;i++) 4444 quadtree->Add(vertices[i]); 4445 4446 SetVertexFieldOn(); 4447 4448 for (i=0;i<nbe;i++) 4449 if(edges[i].on) 4450 for(int j=0;j<2;j++) 4451 if (!edges[i].adj[j]) 4452 if(!edges[i][j].on->IsRequiredVertex()) { 4453 cerr << " Erreur adj et sommet requis edges [" << i << "][ " << j << "]= " 4454 << Number(edges[i][j]) << " : " << " on = " << Gh.Number(edges[i].on) ; 4455 if (edges[i][j].on->OnGeomVertex()) 4456 cerr << " vertex " << Gh.Number(edges[i][j].on->gv); 4457 else if (edges[i][j].on->OnGeomEdge()) 4458 cerr << "Edges " << Gh.Number(edges[i][j].on->ge); 4459 else 4460 cerr << " = " << edges[i][j].on ; 4461 cerr << endl; 4462 } 4463 } 4464 CurrentTh=OldCurrentTh; 4465 } 4466 /*}}}1*/ 4467 /*FUNCTION Triangles::Optim{{{1*/ 4468 Int4 Triangle::Optim(Int2 i,int koption) { 4469 // turne around in positif sens 4470 Int4 NbSwap =0; 4471 Triangle *t = this; 4472 int k=0,j =OppositeEdge[i]; 4473 int jp = PreviousEdge[j]; 4474 // initialise tp, jp the previous triangle & edge 4475 Triangle *tp= at[jp]; 4476 jp = aa[jp]&3; 4477 do { 4478 // cout << *t << " " << j << "\n\t try swap " ; 4479 while (t->swap(j,koption)) 4480 { 4481 NbSwap++; 4482 assert(k++<20000); 4483 t= tp->at[jp]; // set unchange t qnd j for previous triangles 4484 j= NextEdge[tp->aa[jp]&3]; 4485 // cout << "\n\t s " << *t << " " << j << endl; 4486 } 4487 // end on this Triangle 4488 tp = t; 4489 jp = NextEdge[j]; 4490 4491 t= tp->at[jp]; // set unchange t qnd j for previous triangles 4492 j= NextEdge[tp->aa[jp]&3]; 4493 4494 } while( t != this); 4495 return NbSwap; 4496 } 4497 /*}}}1*/ 4498 /*FUNCTION Triangles::SmoothingVertex{{{1*/ 4499 void Triangles::SmoothingVertex(int nbiter,Real8 omega ) { 4500 long int verbosity=0; 4501 // if quatree exist remove it end reconstruct 4502 if (quadtree) delete quadtree; 4503 quadtree=0; 4504 ReMakeTriangleContainingTheVertex(); 4505 Triangle vide; // a triangle to mark the boundary vertex 4506 Triangle ** tstart= new Triangle* [nbv]; 4507 Int4 i,j,k; 4508 // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide 4509 if ( this == & BTh) 4510 for ( i=0;i<nbv;i++) 4511 tstart[i]=vertices[i].t; 4512 else 4513 for ( i=0;i<nbv;i++) 4514 tstart[i]=0; 4515 for ( j=0;j<NbVerticesOnGeomVertex;j++ ) 4516 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide; 4517 for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 4518 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide; 4519 if(verbosity>2) 4520 cout << " -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl; 4521 for (k=0;k<nbiter;k++) 4522 { 4523 Int4 i,NbSwap =0; 4524 Real8 delta =0; 4525 for ( i=0;i<nbv;i++) 4526 if (tstart[i] != &vide) // not a boundary vertex 4527 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega)); 4528 if (!NbOfQuad) 4529 for ( i=0;i<nbv;i++) 4530 if (tstart[i] != &vide) // not a boundary vertex 4531 NbSwap += vertices[i].Optim(1); 4532 if (verbosity>3) 4533 cout << " Move max = " << sqrt(delta) << " iteration = " 4534 << k << " Nb of Swap = " << NbSwap << endl; 4535 } 4536 4537 delete [] tstart; 4538 if (quadtree) quadtree= new QuadTree(this); 4539 } 4540 /*}}}1*/ 4541 /*FUNCTION Triangles::MakeQuadTree{{{1*/ 4542 void Triangles::MakeQuadTree() { 4543 long int verbosity=0; 4544 if(verbosity>8) 4545 cout << " MakeQuadTree" << endl; 4546 if ( !quadtree ) quadtree = new QuadTree(this); 4547 4548 } 4549 /*}}}1*/ 4550 /*FUNCTION Triangles::ShowRegulaty{{{1*/ 4551 void Triangles::ShowRegulaty() const { 4552 const Real8 sqrt32=sqrt(3.)*0.5; 4553 const Real8 aireKh=sqrt32*0.5; 4554 D2 Beq(1,0),Heq(0.5,sqrt32); 4555 D2xD2 Br(D2xD2(Beq,Heq).t()); 4556 D2xD2 B1r(Br.inv()); 4557 /* D2xD2 BB = Br.t()*Br; 4558 cout << " BB = " << BB << " " << Br*B1r << endl; 4559 MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y); 4560 MatVVP2x2 VMM(MMM); 4561 cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl; 4562 */ 4563 double gammamn=1e100,hmin=1e100; 4564 double gammamx=0,hmax=0; 4565 double beta=1e100; 4566 double beta0=0; 4567 double alpha2=0; 4568 double area=0,Marea=0; 4569 // Real8 cf= Real8(coefIcoor); 4570 // Real8 cf2= 6.*cf*cf; 4571 int nt=0; 4572 for (int it=0;it<nbt;it++) 4573 if ( triangles[it].link) 4574 { 4575 nt++; 4576 Triangle &K=triangles[it]; 4577 Real8 area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.; 4578 area+= area3; 4579 D2xD2 B_Kt(K[0],K[1],K[2]); 4580 D2xD2 B_K(B_Kt.t()); 4581 D2xD2 B1K = Br*B_K.inv(); 4582 D2xD2 BK = B_K*B1r; 4583 D2xD2 B1B1 = B1K.t()*B1K; 4584 MetricAnIso MK(B1B1.x.x,B1B1.x.y,B1B1.y.y); 4585 MatVVP2x2 VMK(MK); 4586 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1)); 4587 // cout << B_K << " * " << B1r << " == " << BK << " " << B_K*B_K.inv() << endl; 4588 Real8 betaK=0; 4589 4590 for (int j=0;j<3;j++) 4591 { 4592 Real8 he= Norme2(R2(K[j],K[(j+1)%3])); 4593 hmin=Min(hmin,he); 4594 hmax=Max(hmax,he); 4595 Vertex & v=K[j]; 4596 D2xD2 M((MetricAnIso)v); 4597 betaK += sqrt(M.det()); 4598 4599 D2xD2 BMB = BK.t()*M*BK; 4600 MetricAnIso M1(BMB.x.x,BMB.x.y,BMB.y.y); 4601 MatVVP2x2 VM1(M1); 4602 //cout << B_K <<" " << M << " " << he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<< endl; 4603 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2); 4604 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2); 4605 } 4606 betaK *= area3;// 1/2 (somme sqrt(det))* area(K) 4607 Marea+= betaK; 4608 // cout << betaK << " " << area3 << " " << beta << " " << beta0 << " " << area3*3*3*3 <<endl; 4609 beta=min(beta,betaK); 4610 beta0=max(beta0,betaK); 4611 } 4612 area*=3; 4613 gammamn=sqrt(gammamn); 4614 gammamx=sqrt(gammamx); 4615 cout << " -- adaptmesh Regulary: Nb triangles " << nt << " , h min " << hmin << " , h max " << hmax << endl; 4616 cout << " area = " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl; 4617 cout << " infiny-regulaty: min " << gammamn << " max " << gammamx << endl; 4618 cout << " anisomax "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh) 4619 << " min "<< 1./sqrt(beta0/aireKh) << endl; 4620 } 4621 /*}}}1*/ 4622 /*FUNCTION Triangles::ShowHistogram{{{1*/ 4623 void Triangles::ShowHistogram() const { 4624 4625 const Int4 kmax=10; 4626 const Real8 llmin = 0.5,llmax=2; 4627 const Real8 lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin); 4628 Int4 histo[kmax+1]; 4629 Int4 i,it,k, nbedges =0; 4630 for (i=0;i<=kmax;i++) histo[i]=0; 4631 for (it=0;it<nbt;it++) 4632 if ( triangles[it].link) 4633 { 4634 4635 for (int j=0;j<3;j++) 4636 { 4637 Triangle *ta = triangles[it].TriangleAdj(j); 4638 if ( !ta || !ta->link || Number(ta) >= it) 4639 { 4640 Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]]; 4641 Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]]; 4642 if ( !& vP || !&vQ) continue; 4643 R2 PQ = vQ.r - vP.r; 4644 Real8 l = log(LengthInterpole(vP,vQ,PQ)); 4645 nbedges++; 4646 k = (int) ((l - lmin)*delta); 4647 k = Min(Max(k,0L),kmax); 4648 histo[k]++; 4649 } 4650 } 4651 } 4652 cout << " -- Histogram of the unit mesh, nb of edges" << nbedges << endl <<endl; 4653 4654 cout << " length of edge in | % of edge | Nb of edges " << endl; 4655 cout << " ------------------- | ---------- | ----------- " << endl; 4656 for (i=0;i<=kmax;i++) 4657 { 4658 cout << " " ; 4659 cout.width(10); 4660 if (i==0) cout << " 0 " ; 4661 else cout << exp(lmin+i/delta) ; 4662 cout.width(); cout << "," ; 4663 cout.width(10); 4664 if (i==kmax) cout << " +infty " ; 4665 else cout << exp(lmin+(i+1)/delta) ; 4666 cout.width();cout << " | " ; 4667 4668 cout.precision(4); 4669 cout.width(6); 4670 cout << ((long) ((10000.0 * histo[i])/ nbedges))/100.0 ; 4671 cout.width(); 4672 cout.precision(); 4673 cout << " | " << histo[i] <<endl; 4674 } 4675 cout << " ------------------- | ---------- | ----------- " << endl <<endl; 4676 4677 } 4678 /*}}}1*/ 4679 /*FUNCTION Triangles::Crack{{{1*/ 4680 int Triangles::Crack() { 4681 assert(NbCrackedEdges ==0 || NbCrackedVertices >0); 4682 for (int i=0;i<NbCrackedEdges;i++) 4683 CrackedEdges[i].Crack(); 4684 return NbCrackedEdges; 4685 } 4686 /*}}}1*/ 4687 /*FUNCTION Triangles::UnCrack{{{1*/ 4688 int Triangles::UnCrack() { 4689 assert(NbCrackedEdges ==0 || NbCrackedVertices >0); 4690 for (int i=0;i<NbCrackedEdges;i++) 4691 CrackedEdges[i].UnCrack(); 4692 return NbCrackedEdges; 4693 } 4694 /*}}}1*/ 4695 /*FUNCTION Triangles::CrackMesh{{{1*/ 4696 int Triangles::CrackMesh() { 4697 long int verbosity=0; 4698 Triangles *CurrentThOld = CurrentTh; 4699 // computed the number of cracked edge 4700 int i,k; 4701 for (k=i=0;i<nbe;i++) 4702 if(edges[i].on->Cracked()) k++; 4703 if( k==0) return 0; 4704 CurrentTh = this; 4705 cout << " Nb of Cracked Edges = " << k << endl; 4706 NbCrackedEdges =k; 4707 CrackedEdges = new CrackedEdge[k]; 4708 // new edge 4709 Edge * e = new Edge[ nbe + k]; 4710 4711 // copy 4712 for (i=0;i<nbe;i++) 4713 e[i] = edges[i]; 4714 delete edges; 4715 edges = e; 4716 4717 const int nbe0 = nbe; 4718 for (k=i=0;i<nbe0;i++) // on double les arete cracked 4719 if(edges[i].on->Cracked()) 4720 { 4721 e[nbe] = e[i]; 4722 // return the edge 4723 e[nbe].v[0] = e[i].v[1]; 4724 e[nbe].v[1] = e[i].v[0]; 4725 e[nbe].on = e[i].on->link ; // fqux 4726 CrackedEdges[k++]=CrackedEdge(edges,i,nbe); 4727 nbe++; 4728 } 4729 ReMakeTriangleContainingTheVertex() ; 4730 // 4731 int nbcrakev =0; 4732 Vertex *vlast = vertices + nbv; 4733 Vertex *vend = vertices + nbvx; // end of array 4734 for (int iv=0;iv<nbv;iv++) // vertex 4735 { 4736 Vertex & v= vertices[iv]; 4737 Vertex * vv = & v; 4738 int kk=0; // nb cracked 4739 int kc=0; 4740 int kkk =0; // nb triangle with same number 4741 Triangle * tbegin = v.t; 4742 int i = v.vint; 4743 assert(tbegin && (i >= 0 ) && (i <3)); 4744 // turn around the vertex v 4745 TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge 4746 int k=0; 4747 do { 4748 int kv = VerticesOfTriangularEdge[ta][1]; 4749 k++; 4750 Triangle * tt (ta); 4751 if ( ta.Cracked() ) 4752 { 4753 TriangleAdjacent tta=(ta.Adj()); 4754 assert(tta.Cracked()); 4755 if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart 4756 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use 4757 kk++;// number of cracked edge view 4758 } 4759 if ( tt->link ) { // if good triangles store the value 4760 int it = Number(tt); 4761 assert(it < nt); 4762 (*tt)(kv)= vv; // Change the vertex of triangle 4763 if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref 4764 // tt->SetTriangleContainingTheVertex(); 4765 kkk++; 4766 } else if (kk) { // crack + boundary 4767 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use 4768 } 4769 4770 ta = Next(ta).Adj(); 4771 } while ( (tbegin != ta)); 4772 assert(k); 4773 if (kc) nbcrakev++; 4774 } 4775 4776 if ( nbcrakev ) 4777 for (int iec =0;iec < NbCrackedEdges; iec ++) 4778 CrackedEdges[iec].Set(); 4779 4780 // set the ref 4781 cout << " set the ref " << endl ; 4782 NbCrackedVertices = nbcrakev; 4783 // int nbvo = nbv; 4784 nbv = vlast - vertices; 4785 int nbnewv = nbv - nbv; // nb of new vrtices 4786 if (nbcrakev && verbosity > 1 ) 4787 cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nbnewv<< endl; 4788 // all the new vertices are on geometry 4789 // BOFBO-- A VOIR 4790 if (nbnewv) 4791 { // 4792 Int4 n = nbnewv+NbVerticesOnGeomVertex; 4793 Int4 i,j,k; 4794 VertexOnGeom * vog = new VertexOnGeom[n]; 4795 for ( i =0; i<NbVerticesOnGeomVertex;i++) 4796 vog[i]=VerticesOnGeomVertex[i]; 4797 delete [] VerticesOnGeomVertex; 4798 VerticesOnGeomVertex = vog; 4799 // loop on cracked edge 4800 Vertex * LastOld = vertices + nbv - nbnewv; 4801 for (int iec =0;iec < NbCrackedEdges; iec ++) 4802 for (k=0;k<2;k++) 4803 { 4804 Edge & e = *( k ? CrackedEdges[iec].a.edge : CrackedEdges[iec].b.edge); 4805 for (j=0;j<2;j++) 4806 { 4807 Vertex * v = e(j); 4808 if ( v >= LastOld) 4809 { // a new vertex 4810 Int4 old = v->ReferenceNumber ; // the old same vertex 4811 Int4 i = ( v - LastOld); 4812 // if the old is on vertex => warning 4813 // else the old is on edge => ok 4814 vog[i] = vog[old]; 4815 // vog[i].mv = v; 4816 //g[i].ge = ; 4817 //og[i].abcisse = ; 4818 } 4819 4820 } 4821 } 4822 4823 NbVerticesOnGeomVertex = n; 4824 } 4825 SetVertexFieldOn(); 4826 4827 4828 if (vlast >= vend) 4829 { 4830 cerr << " Not enougth vertices to crack the mesh we need " << nbv << " vertices " << endl; 4831 MeshError(555,this); 4832 } 4833 cout << " NbCrackedVertices " << NbCrackedVertices << endl; 4834 CurrentTh = CurrentThOld; 4835 return NbCrackedVertices; 4836 } 4837 /*}}}1*/ 4838 /*FUNCTION Triangles::FindTriangleContening{{{1*/ 4839 Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const { 4840 Triangle * t=0; 4841 int j,jp,jn,jj; 4842 if (tstart) 4843 t=tstart; 4844 else 4845 { 4846 assert(quadtree); 4847 Vertex *a = quadtree->NearestVertex(B.x,B.y) ; 4848 4849 if (! a || !a->t ) { 4850 if (a) 4851 {cerr << " Attention PB TriangleConteningTheVertex vertex number=" << Number(a) << endl; 4852 cerr << "We forget a call to ReMakeTriangleContainingTheVertex" << endl;} 4853 cerr << " Pb with " << B << toR2(B) << endl; 4854 MeshError(7777); 4855 } 4856 assert(a>= vertices && a < vertices+nbv); 4857 // int k=0; 4858 t = a->t; 4859 assert(t>= triangles && t < triangles+nbt); 4860 4861 } 4862 Icoor2 detop ; 4863 int kkkk =0; // number of test triangle 4864 4865 while ( t->det < 0) 4866 { // the initial triangles is outside 4867 int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; 4868 assert(k0>=0); // k0 the NULL vertex 4869 int k1=NextVertex[k0],k2=PreviousVertex[k0]; 4870 dete[k0]=det(B,(*t)[k1],(*t)[k2]); 4871 dete[k1]=dete[k2]=-1; 4872 if (dete[k0] > 0) // outside B 4873 return t; 4874 t = t->TriangleAdj(OppositeEdge[k0]); 4875 assert(kkkk++ < 2); 4876 } 4877 4878 jj=0; 4879 detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B); 4880 4881 while(t->det > 0 ) 4882 { 4883 assert( kkkk++ < 2000 ); 4884 j= OppositeVertex[jj]; 4885 dete[j] = detop; //det(*b,*s1,*s2); 4886 jn = NextVertex[j]; 4887 jp = PreviousVertex[j]; 4888 dete[jp]= det(*(*t)(j),*(*t)(jn),B); 4889 dete[jn] = t->det-dete[j] -dete[jp]; 4890 4891 // count the number k of dete <0 4892 int k=0,ii[3]; 4893 if (dete[0] < 0 ) ii[k++]=0; 4894 if (dete[1] < 0 ) ii[k++]=1; 4895 if (dete[2] < 0 ) ii[k++]=2; 4896 // 0 => ok 4897 // 1 => go in way 1 4898 // 2 => two way go in way 1 or 2 randomly 4899 4900 if (k==0) 4901 break; 4902 if (k == 2 && BinaryRand()) 4903 Exchange(ii[0],ii[1]); 4904 assert ( k < 3); 4905 TriangleAdjacent t1 = t->Adj(jj=ii[0]); 4906 if ((t1.det() < 0 ) && (k == 2)) 4907 t1 = t->Adj(jj=ii[1]); 4908 t=t1; 4909 j=t1;// for optimisation we now the -det[OppositeVertex[j]]; 4910 detop = -dete[OppositeVertex[jj]]; 4911 jj = j; 4912 } 4913 4914 if (t->det<0) // outside triangle 4915 dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop; 4916 // NbOfTriangleSearchFind += kkkk; 4917 return t; 4918 } 4919 /*}}}1*/ 1650 4920 1651 4921 } // end of namespace bamg
Note:
See TracChangeset
for help on using the changeset viewer.
