![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3231
- Timestamp:
- 03/09/10 09:52:02 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Bamgx.cpp
r3217 r3231 143 143 if (bamgopts->field){ 144 144 if (verbosity>1) printf(" Generating Metric from solution field...\n"); 145 BTh. IntersectConsMetric(bamgopts);145 BTh.AddMetric(bamgopts); 146 146 } 147 147 } 148 148 149 BTh. IntersectGeomMetric(bamgopts);149 BTh.AddGeometryMetric(bamgopts); 150 150 if (verbosity>1) printf(" Smoothing Metric..."); 151 151 if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation); -
issm/trunk/src/c/Bamgx/Mesh2.h
r3230 r3231 705 705 return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y); 706 706 } 707 void Add (Vertex & s,Triangle * t,Icoor2 * =0) ;707 void AddVertex(Vertex & s,Triangle * t,Icoor2 * =0) ; 708 708 void Insert(); 709 709 void ForceBoundary(); 710 710 void Renumber(BamgOpts* bamgopts); 711 711 void FindSubDomain(int OutSide=0); 712 Int4 ConsRefTriangle(Int4 *) const;712 Int4 TriangleReferenceList(Int4 *) const; 713 713 void ShowHistogram() const; 714 714 void ShowRegulaty() const; … … 748 748 void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef); 749 749 void WriteMetric(BamgOpts* bamgopts); 750 void IntersectConsMetric(BamgOpts* bamgopts);750 void AddMetric(BamgOpts* bamgopts); 751 751 void BuildMetric0(BamgOpts* bamgopts); 752 752 void BuildMetric1(BamgOpts* bamgopts); 753 void IntersectGeomMetric(BamgOpts* bamgopts);753 void AddGeometryMetric(BamgOpts* bamgopts); 754 754 int isCracked() const {return NbCrackedVertices != 0;} 755 755 int Crack(); 756 756 int UnCrack(); 757 757 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); 758 void FillHoleInMesh() ;758 void GenerateMeshProperties() ; 759 759 int CrackMesh(); 760 760 -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3230 r3231 27 27 ReadMesh(bamgmesh,bamgopts); 28 28 SetIntCoor(); 29 FillHoleInMesh();29 GenerateMeshProperties(); 30 30 } 31 31 /*}}}1*/ … … 36 36 ReadMesh(index,x,y,nods,nels); 37 37 SetIntCoor(); 38 FillHoleInMesh();38 GenerateMeshProperties(); 39 39 } 40 40 /*}}}1*/ … … 46 46 int * kk = new int [Tho.nbv]; 47 47 Int4 * reft = new Int4[Tho.nbt]; 48 Int4 nbInT = Tho. ConsRefTriangle(reft);48 Int4 nbInT = Tho.TriangleReferenceList(reft); 49 49 Int4 * refv = new Int4[Tho.nbv]; 50 50 … … 130 130 Gh.AfterRead(); 131 131 SetIntCoor(); 132 FillHoleInMesh();132 GenerateMeshProperties(); 133 133 134 134 if (!NbSubDomains){ … … 552 552 //Build reft that holds the number the subdomain number of each triangle 553 553 Int4* reft = new Int4[nbt]; 554 Int4 nbInT = ConsRefTriangle(reft);554 Int4 nbInT = TriangleReferenceList(reft); 555 555 556 556 //Vertices … … 862 862 863 863 /*Methods*/ 864 /*FUNCTION Triangles::Add{{{1*/ 865 void Triangles::Add( Vertex &s,Triangle* t, Icoor2* det3) { 864 /*FUNCTION Triangles::AddGeometryMetric{{{1*/ 865 void Triangles::AddGeometryMetric(BamgOpts* bamgopts){ 866 867 /*Get options*/ 868 int verbosity=bamgopts->verbose; 869 double anisomax =bamgopts->anisomax; 870 double errg =bamgopts->errg; 871 872 Real8 ss[2]={0.00001,0.99999}; 873 Real8 errC = 2*sqrt(2*errg); 874 Real8 hmax = Gh.MaximalHmax(); 875 Real8 hmin = Gh.MinimalHmin(); 876 877 //check that hmax is positive 878 if (hmax<=0){ 879 throw ErrorException(__FUNCT__,exprintf("hmax<=0")); 880 } 881 882 //errC cannot be higher than 1 883 if (errC > 1) errC = 1; 884 885 //Set all vertices to "on" 886 SetVertexFieldOn(); 887 888 //loop over all the vertices on edges 889 for (Int4 i=0;i<nbe;i++){ 890 for (int j=0;j<2;j++){ 891 892 Vertex V; 893 VertexOnGeom GV; 894 Gh.ProjectOnCurve(edges[i],ss[j],V,GV); 895 896 GeometricalEdge * eg = GV; 897 Real8 s = GV; 898 R2 tg; 899 Real8 R1= eg->R1tg(s,tg); 900 Real8 ht=hmax; 901 // err relative to the length of the edge 902 if (R1>1.0e-20) { 903 ht = Min(Max(errC/R1,hmin),hmax); 904 } 905 Real8 hn=Min(hmax,ht*anisomax); 906 if (ht<=0 || hn<=0){ 907 throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0")); 908 } 909 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg); 910 Metric MVp(Vp); 911 edges[i][j].m.IntersectWith(MVp); 912 } 913 } 914 // the problem is for the vertex on vertex 915 } 916 /*}}}1*/ 917 /*FUNCTION Triangles::AddMetric{{{1*/ 918 void Triangles::AddMetric(BamgOpts* bamgopts){ 919 // Hessiantype = 0 => H is computed using double P2 projection 920 // Hessiantype = 1 => H is computed with green formula 921 922 /*Options*/ 923 int Hessiantype=bamgopts->Hessiantype; 924 925 if (Hessiantype==0){ 926 BuildMetric0(bamgopts); 927 } 928 else if (Hessiantype==1){ 929 BuildMetric1(bamgopts); 930 } 931 else{ 932 throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (1->use Green formula, 0-> double P2 projection)",Hessiantype)); 933 } 934 } 935 /*}}}1*/ 936 /*FUNCTION Triangles::AddVertex{{{1*/ 937 void Triangles::AddVertex( Vertex &s,Triangle* t, Icoor2* det3) { 866 938 // ------------------------------------------- 867 939 // s2 … … 935 1007 if ((( Triangle *) ta)->det < 0 ) { 936 1008 // add in outside triangle 937 Add (s,( Triangle *) ta);1009 AddVertex(s,( Triangle *) ta); 938 1010 return; 939 1011 } … … 1978 2050 } 1979 2051 /*}}}1*/ 1980 /*FUNCTION Triangles::ConsRefTriangle{{{1*/1981 Int4 Triangles::ConsRefTriangle(Int4* reft) const {1982 long int verbosity=0;1983 register Triangle *t0,*t;1984 register Int4 k=0, num;1985 1986 //initialize all triangles as -1 (outside)1987 for (Int4 it=0;it<nbt;it++) reft[it]=-1;1988 1989 //loop over all subdomains1990 for (Int4 i=0;i<NbSubDomains;i++){1991 1992 //first triangle of the subdomain i1993 t=t0=subdomains[i].head;1994 1995 //check that the subdomain is not empty1996 if (!t0){ // no empty sub domai{1997 throw ErrorException(__FUNCT__,exprintf("At least one subdomain is empty"));1998 }1999 2000 //loop2001 do{2002 k++;2003 2004 //get current triangle number2005 num = Number(t);2006 2007 //check that num is in [0 nbt[2008 if (num<0 || num>=nbt){2009 throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt"));2010 }2011 2012 //reft of this triangle is the subdomain number2013 reft[num]=i;2014 } while (t0 != (t=t->link));2015 //stop when all triangles of subdomains have been tagged2016 2017 }2018 return k;2019 }2020 /*}}}1*/2021 2052 /*FUNCTION Triangles::Crack{{{1*/ 2022 2053 int Triangles::Crack() { … … 2215 2246 if (verbosity > 3) printf(" number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 2216 2247 } 2217 /*}}}1*/2218 /*FUNCTION Triangles::FillHoleInMesh{{{1*/2219 void Triangles::FillHoleInMesh() {2220 2221 int verbosity=0;2222 2223 // generation of the integer coordinate2224 {2225 2226 // find extrema coordinates of vertices pmin,pmax2227 Int4 i;2228 if(verbosity>2) printf(" Filling holes in mesh of %i vertices\n",nbv);2229 2230 //initialize ordre2231 if (!ordre){2232 throw ErrorException(__FUNCT__,exprintf("!ordre"));2233 }2234 for (i=0;i<nbv;i++) ordre[i]=0;2235 2236 NbSubDomains =0;2237 2238 /* generation of the adjacence of the triangles*/2239 2240 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);2241 2242 //initialize st2243 Int4* st = new Int4[nbt*3];2244 for (i=0;i<nbt*3;i++) st[i]=-1;2245 2246 //check number of edges2247 Int4 kk =0;2248 for (i=0;i<nbe;i++){2249 kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));2250 }2251 if (kk != nbe) {2252 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe));2253 }2254 2255 //2256 for (i=0;i<nbt;i++){2257 for (int j=0;j<3;j++) {2258 Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),2259 Number(triangles[i][VerticesOfTriangularEdge[j][1]]));2260 Int4 invisible = triangles[i].Hidden(j);2261 if(st[k]==-1){2262 st[k]=3*i+j;2263 }2264 else if(st[k]>=0) {2265 if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){2266 throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))"));2267 }2268 2269 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));2270 if (invisible) triangles[i].SetHidden(j);2271 if (k<nbe) {2272 triangles[i].SetLocked(j);2273 }2274 st[k]=-2-st[k];2275 }2276 else {2277 throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",2278 Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));2279 }2280 }2281 }2282 if(verbosity>5) {2283 printf(" info on Mesh:\n");2284 printf(" - number of vertices = %i \n",nbv);2285 printf(" - number of triangles = %i \n",nbt);2286 printf(" - number of given edges = %i \n",nbe);2287 printf(" - number of all edges = %i \n" ,edge4->nb());2288 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv);2289 }2290 2291 // check the consistant of edge[].adj and the geometrical required vertex2292 Int4 k=0;2293 for (i=0;i<edge4->nb();i++){2294 if (st[i] >=0){ // edge alone2295 if (i < nbe) {2296 Int4 i0=edge4->i(i);2297 ordre[i0] = vertices+i0;2298 Int4 i1=edge4->j(i);2299 ordre[i1] = vertices+i1;2300 }2301 else {2302 k=k+1;2303 if (k <20) {2304 throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));2305 }2306 }2307 }2308 }2309 if(k != 0) {2310 throw ErrorException(__FUNCT__,exprintf("%i boundary edges are not defined as edges",k));2311 }2312 2313 /* mesh generation with boundary points*/2314 Int4 nbvb = 0;2315 for (i=0;i<nbv;i++){2316 vertices[i].t=0;2317 vertices[i].vint=0;2318 if (ordre[i]){2319 ordre[nbvb++] = ordre[i];2320 }2321 }2322 2323 Triangle* savetriangles= triangles;2324 Int4 savenbt=nbt;2325 Int4 savenbtx=nbtx;2326 SubDomain * savesubdomains = subdomains;2327 subdomains = 0;2328 2329 Int4 Nbtriafillhole = 2*nbvb;2330 Triangle* triafillhole =new Triangle[Nbtriafillhole];2331 triangles = triafillhole;2332 2333 nbt=2;2334 nbtx= Nbtriafillhole;2335 2336 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)2337 if ( ++i >= nbvb) {2338 throw ErrorException(__FUNCT__,exprintf("FillHoleInMesh: All the vertices are aligned"));2339 }2340 Exchange( ordre[2], ordre[i]);2341 2342 Vertex * v0=ordre[0], *v1=ordre[1];2343 2344 2345 triangles[0](0) = 0; // sommet pour infini2346 triangles[0](1) = v0;2347 triangles[0](2) = v1;2348 2349 triangles[1](0) = 0;// sommet pour infini2350 triangles[1](2) = v0;2351 triangles[1](1) = v1;2352 const int e0 = OppositeEdge[0];2353 const int e1 = NextEdge[e0];2354 const int e2 = PreviousEdge[e0];2355 triangles[0].SetAdj2(e0, &triangles[1] ,e0);2356 triangles[0].SetAdj2(e1, &triangles[1] ,e2);2357 triangles[0].SetAdj2(e2, &triangles[1] ,e1);2358 2359 triangles[0].det = -1; // faux triangles2360 triangles[1].det = -1; // faux triangles2361 2362 triangles[0].SetTriangleContainingTheVertex();2363 triangles[1].SetTriangleContainingTheVertex();2364 2365 triangles[0].link=&triangles[1];2366 triangles[1].link=&triangles[0];2367 2368 if (!quadtree )2369 delete quadtree; // ->ReInitialise();2370 2371 quadtree = new QuadTree(this,0);2372 quadtree->Add(*v0);2373 quadtree->Add(*v1);2374 2375 // We add the vertices one by one2376 Int4 NbSwap=0;2377 for (Int4 icount=2; icount<nbvb; icount++) {2378 Vertex *vi = ordre[icount];2379 Icoor2 dete[3];2380 Triangle *tcvi = FindTriangleContaining(vi->i,dete);2381 quadtree->Add(*vi);2382 Add(*vi,tcvi,dete);2383 NbSwap += vi->Optim(1,1);2384 }2385 2386 // inforce the boundary2387 TriangleAdjacent ta(0,0);2388 Int4 nbloss = 0,knbe=0;2389 for ( i = 0; i < nbe; i++){2390 if (st[i] >=0){ // edge alone => on border ... FH oct 20092391 Vertex & a=edges[i][0], & b = edges[i][1];2392 if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.2393 {2394 knbe++;2395 if (ForceEdge(a,b,ta)<0)2396 nbloss++;2397 }2398 }2399 }2400 if(nbloss) {2401 throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe));2402 }2403 2404 FindSubDomain(1);2405 // remove all the hole2406 // remove all the good sub domain2407 Int4 krm =0;2408 for (i=0;i<nbt;i++){2409 if (triangles[i].link){ // remove triangles2410 krm++;2411 for (int j=0;j<3;j++)2412 {2413 TriangleAdjacent ta = triangles[i].Adj(j);2414 Triangle & tta = * (Triangle *) ta;2415 if(! tta.link) // edge between remove and not remove2416 { // change the link of ta;2417 int ja = ta;2418 Vertex *v0= ta.EdgeVertex(0);2419 Vertex *v1= ta.EdgeVertex(1);2420 Int4 k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);2421 if (st[k]<0){2422 throw ErrorException(__FUNCT__,exprintf("st[k]<0"));2423 }2424 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));2425 ta.SetLock();2426 st[k]=-2-st[k];2427 }2428 }2429 }2430 }2431 Int4 NbTfillHoll =0;2432 for (i=0;i<nbt;i++){2433 if (triangles[i].link) {2434 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);2435 triangles[i].color=-1;2436 }2437 else2438 {2439 triangles[i].color= savenbt+ NbTfillHoll++;2440 }2441 }2442 if (savenbt+NbTfillHoll>savenbtx ){2443 throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx"));2444 }2445 // copy of the outside triangles in saveTriangles2446 for (i=0;i<nbt;i++){2447 if(triangles[i].color>=0) {2448 savetriangles[savenbt]=triangles[i];2449 savetriangles[savenbt].link=0;2450 savenbt++;2451 }2452 }2453 // gestion of the adj2454 k =0;2455 Triangle * tmax = triangles + nbt;2456 for (i=0;i<savenbt;i++)2457 {2458 Triangle & ti = savetriangles[i];2459 for (int j=0;j<3;j++)2460 {2461 Triangle * ta = ti.TriangleAdj(j);2462 int aa = ti.NuEdgeTriangleAdj(j);2463 int lck = ti.Locked(j);2464 if (!ta) k++; // bug2465 else if ( ta >= triangles && ta < tmax)2466 {2467 ta= savetriangles + ta->color;2468 ti.SetAdj2(j,ta,aa);2469 if(lck) ti.SetLocked(j);2470 }2471 }2472 }2473 // OutSidesTriangles = triangles;2474 // Int4 NbOutSidesTriangles = nbt;2475 2476 // restore triangles;2477 nbt=savenbt;2478 nbtx=savenbtx;2479 delete [] triangles;2480 delete [] subdomains;2481 triangles = savetriangles;2482 subdomains = savesubdomains;2483 if (k) {2484 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));2485 }2486 FindSubDomain();2487 2488 delete edge4;2489 delete [] st;2490 for (i=0;i<nbv;i++)2491 quadtree->Add(vertices[i]);2492 2493 SetVertexFieldOn();2494 2495 for (i=0;i<nbe;i++)2496 if(edges[i].onGeometry)2497 for(int j=0;j<2;j++)2498 if (!edges[i].adj[j])2499 if(!edges[i][j].onGeometry->IsRequiredVertex()) {2500 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)"));2501 }2502 }2503 }2504 2248 /*}}}1*/ 2505 2249 /*FUNCTION Triangles::FindSubDomain{{{1*/ … … 2846 2590 } 2847 2591 /*}}}1*/ 2592 /*FUNCTION Triangles::GenerateMeshProperties{{{1*/ 2593 void Triangles::GenerateMeshProperties() { 2594 2595 int verbosity=0; 2596 2597 // generation of the integer coordinate 2598 2599 // find extrema coordinates of vertices pmin,pmax 2600 Int4 i; 2601 if(verbosity>2) printf(" Filling holes in mesh of %i vertices\n",nbv); 2602 2603 //initialize ordre 2604 if (!ordre){ 2605 throw ErrorException(__FUNCT__,exprintf("!ordre")); 2606 } 2607 for (i=0;i<nbv;i++) ordre[i]=0; 2608 2609 NbSubDomains =0; 2610 2611 /* generation of the adjacence of the triangles*/ 2612 2613 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 2614 2615 //initialize st 2616 Int4* st = new Int4[nbt*3]; 2617 for (i=0;i<nbt*3;i++) st[i]=-1; 2618 2619 //check number of edges 2620 Int4 kk =0; 2621 for (i=0;i<nbe;i++){ 2622 kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); 2623 } 2624 if (kk != nbe) { 2625 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe)); 2626 } 2627 2628 // 2629 for (i=0;i<nbt;i++){ 2630 for (int j=0;j<3;j++) { 2631 Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), 2632 Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 2633 Int4 invisible = triangles[i].Hidden(j); 2634 if(st[k]==-1){ 2635 st[k]=3*i+j; 2636 } 2637 else if(st[k]>=0) { 2638 if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){ 2639 throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))")); 2640 } 2641 2642 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); 2643 if (invisible) triangles[i].SetHidden(j); 2644 if (k<nbe) { 2645 triangles[i].SetLocked(j); 2646 } 2647 st[k]=-2-st[k]; 2648 } 2649 else { 2650 throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles", 2651 Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]))); 2652 } 2653 } 2654 } 2655 if(verbosity>5) { 2656 printf(" info on Mesh:\n"); 2657 printf(" - number of vertices = %i \n",nbv); 2658 printf(" - number of triangles = %i \n",nbt); 2659 printf(" - number of given edges = %i \n",nbe); 2660 printf(" - number of all edges = %i \n" ,edge4->nb()); 2661 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); 2662 } 2663 2664 // check the consistant of edge[].adj and the geometrical required vertex 2665 Int4 k=0; 2666 for (i=0;i<edge4->nb();i++){ 2667 if (st[i] >=0){ // edge alone 2668 if (i < nbe) { 2669 Int4 i0=edge4->i(i); 2670 ordre[i0] = vertices+i0; 2671 Int4 i1=edge4->j(i); 2672 ordre[i1] = vertices+i1; 2673 } 2674 else { 2675 k=k+1; 2676 if (k <20) { 2677 throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i))); 2678 } 2679 } 2680 } 2681 } 2682 if(k != 0) { 2683 throw ErrorException(__FUNCT__,exprintf("%i boundary edges are not defined as edges",k)); 2684 } 2685 2686 /* mesh generation with boundary points*/ 2687 Int4 nbvb = 0; 2688 for (i=0;i<nbv;i++){ 2689 vertices[i].t=0; 2690 vertices[i].vint=0; 2691 if (ordre[i]){ 2692 ordre[nbvb++] = ordre[i]; 2693 } 2694 } 2695 2696 Triangle* savetriangles= triangles; 2697 Int4 savenbt=nbt; 2698 Int4 savenbtx=nbtx; 2699 SubDomain * savesubdomains = subdomains; 2700 subdomains = 0; 2701 2702 Int4 Nbtriafillhole = 2*nbvb; 2703 Triangle* triafillhole =new Triangle[Nbtriafillhole]; 2704 triangles = triafillhole; 2705 2706 nbt=2; 2707 nbtx= Nbtriafillhole; 2708 2709 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 2710 if ( ++i >= nbvb) { 2711 throw ErrorException(__FUNCT__,exprintf("GenerateMeshProperties: All the vertices are aligned")); 2712 } 2713 Exchange( ordre[2], ordre[i]); 2714 2715 Vertex * v0=ordre[0], *v1=ordre[1]; 2716 2717 2718 triangles[0](0) = 0; // sommet pour infini 2719 triangles[0](1) = v0; 2720 triangles[0](2) = v1; 2721 2722 triangles[1](0) = 0;// sommet pour infini 2723 triangles[1](2) = v0; 2724 triangles[1](1) = v1; 2725 const int e0 = OppositeEdge[0]; 2726 const int e1 = NextEdge[e0]; 2727 const int e2 = PreviousEdge[e0]; 2728 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 2729 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 2730 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 2731 2732 triangles[0].det = -1; // faux triangles 2733 triangles[1].det = -1; // faux triangles 2734 2735 triangles[0].SetTriangleContainingTheVertex(); 2736 triangles[1].SetTriangleContainingTheVertex(); 2737 2738 triangles[0].link=&triangles[1]; 2739 triangles[1].link=&triangles[0]; 2740 2741 if (!quadtree ) 2742 delete quadtree; // ->ReInitialise(); 2743 2744 quadtree = new QuadTree(this,0); 2745 quadtree->Add(*v0); 2746 quadtree->Add(*v1); 2747 2748 // We add the vertices one by one 2749 Int4 NbSwap=0; 2750 for (Int4 icount=2; icount<nbvb; icount++) { 2751 Vertex *vi = ordre[icount]; 2752 Icoor2 dete[3]; 2753 Triangle *tcvi = FindTriangleContaining(vi->i,dete); 2754 quadtree->Add(*vi); 2755 AddVertex(*vi,tcvi,dete); 2756 NbSwap += vi->Optim(1,1); 2757 } 2758 2759 // inforce the boundary 2760 TriangleAdjacent ta(0,0); 2761 Int4 nbloss = 0,knbe=0; 2762 for ( i = 0; i < nbe; i++){ 2763 if (st[i] >=0){ // edge alone => on border ... FH oct 2009 2764 Vertex & a=edges[i][0], & b = edges[i][1]; 2765 if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1. 2766 { 2767 knbe++; 2768 if (ForceEdge(a,b,ta)<0) 2769 nbloss++; 2770 } 2771 } 2772 } 2773 if(nbloss) { 2774 throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe)); 2775 } 2776 2777 FindSubDomain(1); 2778 // remove all the hole 2779 // remove all the good sub domain 2780 Int4 krm =0; 2781 for (i=0;i<nbt;i++){ 2782 if (triangles[i].link){ // remove triangles 2783 krm++; 2784 for (int j=0;j<3;j++) 2785 { 2786 TriangleAdjacent ta = triangles[i].Adj(j); 2787 Triangle & tta = * (Triangle *) ta; 2788 if(! tta.link) // edge between remove and not remove 2789 { // change the link of ta; 2790 int ja = ta; 2791 Vertex *v0= ta.EdgeVertex(0); 2792 Vertex *v1= ta.EdgeVertex(1); 2793 Int4 k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv); 2794 if (st[k]<0){ 2795 throw ErrorException(__FUNCT__,exprintf("st[k]<0")); 2796 } 2797 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); 2798 ta.SetLock(); 2799 st[k]=-2-st[k]; 2800 } 2801 } 2802 } 2803 } 2804 Int4 NbTfillHoll =0; 2805 for (i=0;i<nbt;i++){ 2806 if (triangles[i].link) { 2807 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL); 2808 triangles[i].color=-1; 2809 } 2810 else 2811 { 2812 triangles[i].color= savenbt+ NbTfillHoll++; 2813 } 2814 } 2815 if (savenbt+NbTfillHoll>savenbtx ){ 2816 throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx")); 2817 } 2818 // copy of the outside triangles in saveTriangles 2819 for (i=0;i<nbt;i++){ 2820 if(triangles[i].color>=0) { 2821 savetriangles[savenbt]=triangles[i]; 2822 savetriangles[savenbt].link=0; 2823 savenbt++; 2824 } 2825 } 2826 // gestion of the adj 2827 k =0; 2828 Triangle * tmax = triangles + nbt; 2829 for (i=0;i<savenbt;i++) 2830 { 2831 Triangle & ti = savetriangles[i]; 2832 for (int j=0;j<3;j++) 2833 { 2834 Triangle * ta = ti.TriangleAdj(j); 2835 int aa = ti.NuEdgeTriangleAdj(j); 2836 int lck = ti.Locked(j); 2837 if (!ta) k++; // bug 2838 else if ( ta >= triangles && ta < tmax) 2839 { 2840 ta= savetriangles + ta->color; 2841 ti.SetAdj2(j,ta,aa); 2842 if(lck) ti.SetLocked(j); 2843 } 2844 } 2845 } 2846 // OutSidesTriangles = triangles; 2847 // Int4 NbOutSidesTriangles = nbt; 2848 2849 // restore triangles; 2850 nbt=savenbt; 2851 nbtx=savenbtx; 2852 delete [] triangles; 2853 delete [] subdomains; 2854 triangles = savetriangles; 2855 subdomains = savesubdomains; 2856 if (k) { 2857 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k)); 2858 } 2859 FindSubDomain(); 2860 2861 delete edge4; 2862 delete [] st; 2863 for (i=0;i<nbv;i++) 2864 quadtree->Add(vertices[i]); 2865 2866 SetVertexFieldOn(); 2867 2868 for (i=0;i<nbe;i++) 2869 if(edges[i].onGeometry) 2870 for(int j=0;j<2;j++) 2871 if (!edges[i].adj[j]) 2872 if(!edges[i][j].onGeometry->IsRequiredVertex()) { 2873 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)")); 2874 } 2875 } 2876 /*}}}1*/ 2848 2877 /*FUNCTION Triangles::GeomToTriangles0{{{1*/ 2849 2878 void Triangles::GeomToTriangles0(Int4 inbvx,BamgOpts* bamgopts) { … … 3503 3532 } 3504 3533 /*}}}1*/ 3505 /*FUNCTION Triangles::IntersectConsMetric{{{1*/3506 void Triangles::IntersectConsMetric(BamgOpts* bamgopts){3507 // Hessiantype = 0 => H is computed using double P2 projection3508 // Hessiantype = 1 => H is computed with green formula3509 3510 /*Options*/3511 int Hessiantype=bamgopts->Hessiantype;3512 3513 if (Hessiantype==0){3514 BuildMetric0(bamgopts);3515 }3516 else if (Hessiantype==1){3517 BuildMetric1(bamgopts);3518 }3519 else{3520 throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (1->use Green formula, 0-> double P2 projection)",Hessiantype));3521 }3522 }3523 /*}}}1*/3524 /*FUNCTION Triangles::IntersectGeomMetric{{{1*/3525 void Triangles::IntersectGeomMetric(BamgOpts* bamgopts){3526 3527 /*Get options*/3528 int verbosity=bamgopts->verbose;3529 double anisomax =bamgopts->anisomax;3530 double errg =bamgopts->errg;3531 3532 Real8 ss[2]={0.00001,0.99999};3533 Real8 errC = 2*sqrt(2*errg);3534 Real8 hmax = Gh.MaximalHmax();3535 Real8 hmin = Gh.MinimalHmin();3536 3537 //check that hmax is positive3538 if (hmax<=0){3539 throw ErrorException(__FUNCT__,exprintf("hmax<=0"));3540 }3541 3542 //errC cannot be higher than 13543 if (errC > 1) errC = 1;3544 3545 //Set all vertices to "on"3546 SetVertexFieldOn();3547 3548 //loop over all the vertices on edges3549 for (Int4 i=0;i<nbe;i++){3550 for (int j=0;j<2;j++){3551 3552 Vertex V;3553 VertexOnGeom GV;3554 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);3555 3556 GeometricalEdge * eg = GV;3557 Real8 s = GV;3558 R2 tg;3559 Real8 R1= eg->R1tg(s,tg);3560 Real8 ht=hmax;3561 // err relative to the length of the edge3562 if (R1>1.0e-20) {3563 ht = Min(Max(errC/R1,hmin),hmax);3564 }3565 Real8 hn=Min(hmax,ht*anisomax);3566 if (ht<=0 || hn<=0){3567 throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));3568 }3569 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);3570 Metric MVp(Vp);3571 edges[i][j].m.IntersectWith(MVp);3572 }3573 }3574 // the problem is for the vertex on vertex3575 }3576 /*}}}1*/3577 3534 /*FUNCTION Triangles::Insert{{{1*/ 3578 3535 void Triangles::Insert() { … … 3692 3649 3693 3650 //Add newvertex to the existing mesh 3694 Add (*newvertex,tcvi,dete);3651 AddVertex(*newvertex,tcvi,dete); 3695 3652 3696 3653 //Make the mesh Delaunay around newvertex by swaping the edges … … 3782 3739 } 3783 3740 quadtree->Add(vj); 3784 Add (vj,tcvj,dete);3741 AddVertex(vj,tcvj,dete); 3785 3742 NbSwap += vj.Optim(1); 3786 3743 iv++; … … 5262 5219 // ReMakeTriangleContainingTheVertex(); 5263 5220 5264 FillHoleInMesh();5221 GenerateMeshProperties(); 5265 5222 5266 5223 if (verbosity>2){ … … 5339 5296 throw ErrorException(__FUNCT__,exprintf("!tcvi || tcvi->det < 0")); 5340 5297 } 5341 Add (vi,tcvi,dete);5298 AddVertex(vi,tcvi,dete); 5342 5299 NbSwap += vi.Optim(1); 5343 5300 iv++; … … 5354 5311 5355 5312 return NbSplitEdge; 5313 } 5314 /*}}}1*/ 5315 /*FUNCTION Triangles::TriangleReferenceList{{{1*/ 5316 Int4 Triangles::TriangleReferenceList(Int4* reft) const { 5317 long int verbosity=0; 5318 register Triangle *t0,*t; 5319 register Int4 k=0, num; 5320 5321 //initialize all triangles as -1 (outside) 5322 for (Int4 it=0;it<nbt;it++) reft[it]=-1; 5323 5324 //loop over all subdomains 5325 for (Int4 i=0;i<NbSubDomains;i++){ 5326 5327 //first triangle of the subdomain i 5328 t=t0=subdomains[i].head; 5329 5330 //check that the subdomain is not empty 5331 if (!t0){ throw ErrorException(__FUNCT__,exprintf("At least one subdomain is empty"));} 5332 5333 //loop 5334 do{ 5335 k++; 5336 5337 //get current triangle number 5338 num = Number(t); 5339 5340 //check that num is in [0 nbt[ 5341 if (num<0 || num>=nbt){ throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt"));} 5342 5343 //reft of this triangle is the subdomain number 5344 reft[num]=i; 5345 5346 } while (t0 != (t=t->link)); 5347 //stop when all triangles of subdomains have been tagged 5348 5349 } 5350 return k; 5356 5351 } 5357 5352 /*}}}1*/
Note:
See TracChangeset
for help on using the changeset viewer.
