![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5401 for issm/trunk/src/c/objects/Bamg/Mesh.cpp
- Timestamp:
- 08/19/10 09:38:48 (15 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/objects/Bamg/Mesh.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Bamg/Mesh.cpp
r5400 r5401 5110 5110 if (b->Required() ) break; 5111 5111 int kprev=k; 5112 k = e->AdjVertex Number[kprev];// next vertices5112 k = e->AdjVertexIndex[kprev];// next vertices 5113 5113 e = e->Adj[kprev]; 5114 5114 ISSMASSERT(e); … … 5329 5329 for(int phase=0;phase<=step;phase++){ 5330 5330 5331 /*Loop over all curves from icurve till the last curve*/ 5332 for(Curve *curve= Gh.curves+icurve;curve;curve= curve->next){ 5333 5334 /*Get index of current curve*/ 5335 int icurveequi= Gh.GetId(curve); 5336 5337 /*For phase 0, check that we are at the begining of the curve only*/ 5338 if(phase==0 && icurveequi!=icurve) continue; 5339 5340 int k0=jedge,k1; 5341 Edge* pe= BTh.edges+iedge; 5342 int iedgeequi=bcurve[icurveequi]/2; 5343 int jedgeequi=bcurve[icurveequi]%2; 5344 5345 int k0equi=jedgeequi,k1equi; 5346 Edge * peequi= BTh.edges+iedgeequi; 5347 GeometricalEdge *ongequi = peequi->GeometricalEdgeHook; 5348 5349 double sNew=Lstep;// abscisse of the new points (phase==1) 5350 L=0;// length of the curve 5351 long i=0;// index of new points on the curve 5352 register GeometricalVertex * GA0 = *(*peequi)[k0equi].GeometricalEdgeHook; 5353 BamgVertex *A0; 5354 A0 = GA0->to; // the vertex in new mesh 5355 BamgVertex *A1; 5356 VertexOnGeom *GA1; 5357 Edge* PreviousNewEdge = 0; 5358 5359 // New Curve phase 5360 ISSMASSERT(A0-vertices>=0 && A0-vertices<nbv); 5361 if(ongequi->Required()){ 5362 GeometricalVertex *GA1 = *(*peequi)[1-k0equi].GeometricalEdgeHook; 5363 A1 = GA1->to; // 5364 } 5365 else { 5366 for(;;){ 5367 Edge &ee=*pe; 5368 Edge &eeequi=*peequi; 5369 k1 = 1-k0; // next vertex of the edge 5370 k1equi= 1 - k0equi; 5371 ISSMASSERT(pe && ee.GeometricalEdgeHook); 5372 ee.GeometricalEdgeHook->SetMark(); 5373 BamgVertex & v0=ee[0], & v1=ee[1]; 5374 R2 AB=(R2)v1-(R2)v0; 5375 double L0=L,LAB; 5376 LAB=LengthInterpole(v0.m,v1.m,AB); 5377 L+= LAB; 5378 5379 if (phase){ 5380 // computation of the new points for the given curve 5381 while ((i!=NbCreatePointOnCurve) && sNew<=L) { 5382 5383 //some checks 5384 ISSMASSERT(sNew>=L0); 5385 ISSMASSERT(LAB); 5386 ISSMASSERT(vertices && nbv<maxnbv); 5387 ISSMASSERT(edges && nbe<nbex); 5388 ISSMASSERT(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex); 5389 5390 // new vertex on edge 5391 A1=vertices+nbv++; 5392 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge; 5393 Edge* e = edges + nbe++; 5394 double se= (sNew-L0)/LAB; 5395 if (se<0 || se>=1.000000001){ 5396 ISSMERROR("Problem creating point on a boundary: se=%g should be in [0 1]",se); 5397 } 5398 se = abscisseInterpole(v0.m,v1.m,AB,se,1); 5399 if (se<0 || se>1){ 5400 ISSMERROR("Problem creating point on a boundary: se=%g should be in [0 1]",se); 5401 } 5402 se = k1 ? se : 1. - se; 5403 se = k1==k1equi ? se : 1. - se; 5404 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 5405 ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 5406 A1->ReferenceNumber = eeequi.ReferenceNumber; 5407 A1->DirOfSearch =NoDirOfSearch; 5408 e->GeometricalEdgeHook = ongequi; 5409 e->v[0]=A0; 5410 e->v[1]=A1; 5411 e->ReferenceNumber = eeequi.ReferenceNumber; 5412 e->adj[0]=PreviousNewEdge; 5413 5414 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; 5415 PreviousNewEdge=e; 5416 A0=A1; 5417 sNew += Lstep; 5418 if (++i== NbCreatePointOnCurve) break; 5331 /*Current curve pointer*/ 5332 Curve *curve= Gh.curves+icurve; 5333 5334 /*Get index of current curve*/ 5335 int icurveequi= Gh.GetId(curve); 5336 5337 /*For phase 0, check that we are at the begining of the curve only*/ 5338 if(phase==0 && icurveequi!=icurve) continue; 5339 5340 int k0=jedge,k1; 5341 Edge* pe= BTh.edges+iedge; 5342 int iedgeequi=bcurve[icurveequi]/2; 5343 int jedgeequi=bcurve[icurveequi]%2; 5344 5345 int k0equi=jedgeequi,k1equi; 5346 Edge * peequi= BTh.edges+iedgeequi; 5347 GeometricalEdge *ongequi = peequi->GeometricalEdgeHook; 5348 5349 double sNew=Lstep;// abscisse of the new points (phase==1) 5350 L=0;// length of the curve 5351 long i=0;// index of new points on the curve 5352 register GeometricalVertex * GA0 = *(*peequi)[k0equi].GeometricalEdgeHook; 5353 BamgVertex *A0; 5354 A0 = GA0->to; // the vertex in new mesh 5355 BamgVertex *A1; 5356 VertexOnGeom *GA1; 5357 Edge* PreviousNewEdge = 0; 5358 5359 // New Curve phase 5360 ISSMASSERT(A0-vertices>=0 && A0-vertices<nbv); 5361 if(ongequi->Required()){ 5362 GeometricalVertex *GA1 = *(*peequi)[1-k0equi].GeometricalEdgeHook; 5363 A1 = GA1->to; // 5364 } 5365 else { 5366 for(;;){ 5367 Edge &ee=*pe; 5368 Edge &eeequi=*peequi; 5369 k1 = 1-k0; // next vertex of the edge 5370 k1equi= 1 - k0equi; 5371 ISSMASSERT(pe && ee.GeometricalEdgeHook); 5372 ee.GeometricalEdgeHook->SetMark(); 5373 BamgVertex & v0=ee[0], & v1=ee[1]; 5374 R2 AB=(R2)v1-(R2)v0; 5375 double L0=L,LAB; 5376 LAB=LengthInterpole(v0.m,v1.m,AB); 5377 L+= LAB; 5378 5379 if (phase){ 5380 // computation of the new points for the given curve 5381 while ((i!=NbCreatePointOnCurve) && sNew<=L) { 5382 5383 //some checks 5384 ISSMASSERT(sNew>=L0); 5385 ISSMASSERT(LAB); 5386 ISSMASSERT(vertices && nbv<maxnbv); 5387 ISSMASSERT(edges && nbe<nbex); 5388 ISSMASSERT(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex); 5389 5390 // new vertex on edge 5391 A1=vertices+nbv++; 5392 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge; 5393 Edge* e = edges + nbe++; 5394 double se= (sNew-L0)/LAB; 5395 if (se<0 || se>=1.000000001){ 5396 ISSMERROR("Problem creating point on a boundary: se=%g should be in [0 1]",se); 5419 5397 } 5398 se = abscisseInterpole(v0.m,v1.m,AB,se,1); 5399 if (se<0 || se>1){ 5400 ISSMERROR("Problem creating point on a boundary: se=%g should be in [0 1]",se); 5401 } 5402 se = k1 ? se : 1. - se; 5403 se = k1==k1equi ? se : 1. - se; 5404 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 5405 ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 5406 A1->ReferenceNumber = eeequi.ReferenceNumber; 5407 A1->DirOfSearch =NoDirOfSearch; 5408 e->GeometricalEdgeHook = ongequi; 5409 e->v[0]=A0; 5410 e->v[1]=A1; 5411 e->ReferenceNumber = eeequi.ReferenceNumber; 5412 e->adj[0]=PreviousNewEdge; 5413 5414 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; 5415 PreviousNewEdge=e; 5416 A0=A1; 5417 sNew += Lstep; 5418 if (++i== NbCreatePointOnCurve) break; 5420 5419 } 5421 5422 //some checks 5423 ISSMASSERT(ee.GeometricalEdgeHook->CurveNumber==ei.GeometricalEdgeHook->CurveNumber);5424 if (ee[k1].GeometricalEdgeHook->IsRequiredVertex()) {5425 ISSMASSERT(eeequi[k1equi].GeometricalEdgeHook->IsRequiredVertex());5426 register GeometricalVertex * GA1 = *eeequi[k1equi].GeometricalEdgeHook;5427 A1=GA1->to;// the vertex in new mesh5428 ISSMASSERT(A1-vertices>=0 && A1-vertices<nbv);5429 break;5430 }5431 if (!ee.adj[k1]) {5432 ISSMERROR(" adj edge %i, nbe=%i, Gh.vertices=%i",BTh.GetId(ee),nbe,Gh.vertices);5433 }5434 pe = ee.adj[k1]; // next edge5435 k0 = pe->Intersection(ee);5436 peequi= eeequi.adj[k1equi]; // next edge5437 k0equi=peequi->Intersection(eeequi);5438 }// for(;;) end of the curve5439 } 5440 5441 5442 if (phase){ // construction of the last edge 5443 Edge* e=edges + nbe++;5444 e->GeometricalEdgeHook = ongequi;5445 e->v[0]=A0;5446 e->v[1]=A1;5447 e->ReferenceNumber = peequi->ReferenceNumber;5448 e->adj[0]=PreviousNewEdge;5449 e->adj[1]=0;5450 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;5451 PreviousNewEdge =e;5452 5453 ISSMASSERT(i==NbCreatePointOnCurve); 5454 }5455 } // end loop on equi curve5420 } 5421 5422 //some checks 5423 ISSMASSERT(ee.GeometricalEdgeHook->CurveNumber==ei.GeometricalEdgeHook->CurveNumber); 5424 if (ee[k1].GeometricalEdgeHook->IsRequiredVertex()) { 5425 ISSMASSERT(eeequi[k1equi].GeometricalEdgeHook->IsRequiredVertex()); 5426 register GeometricalVertex * GA1 = *eeequi[k1equi].GeometricalEdgeHook; 5427 A1=GA1->to;// the vertex in new mesh 5428 ISSMASSERT(A1-vertices>=0 && A1-vertices<nbv); 5429 break; 5430 } 5431 if (!ee.adj[k1]) { 5432 ISSMERROR(" adj edge %i, nbe=%i, Gh.vertices=%i",BTh.GetId(ee),nbe,Gh.vertices); 5433 } 5434 pe = ee.adj[k1]; // next edge 5435 k0 = pe->Intersection(ee); 5436 peequi= eeequi.adj[k1equi]; // next edge 5437 k0equi=peequi->Intersection(eeequi); 5438 }// for(;;) end of the curve 5439 } 5440 5441 5442 if (phase){ // construction of the last edge 5443 Edge* e=edges + nbe++; 5444 e->GeometricalEdgeHook = ongequi; 5445 e->v[0]=A0; 5446 e->v[1]=A1; 5447 e->ReferenceNumber = peequi->ReferenceNumber; 5448 e->adj[0]=PreviousNewEdge; 5449 e->adj[1]=0; 5450 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; 5451 PreviousNewEdge = e; 5452 5453 ISSMASSERT(i==NbCreatePointOnCurve); 5454 } 5456 5455 5457 5456 if (!phase) { // … … 5460 5459 Lcurve = L; 5461 5460 NbCreatePointOnCurve = NbSegOnCurve-1; 5462 5463 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next){ 5464 NbOfNewEdge += NbSegOnCurve; 5465 NbOfNewPoints += NbCreatePointOnCurve; 5466 } 5461 NbOfNewEdge += NbSegOnCurve; 5462 NbOfNewPoints += NbCreatePointOnCurve; 5467 5463 } 5468 5464 }
Note:
See TracChangeset
for help on using the changeset viewer.
