![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2813
- Timestamp:
- 01/12/10 17:20:44 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 4 edited
-
Mesh2.cpp (modified) (1 diff)
-
Mesh2.h (modified) (1 diff)
-
Metric.cpp (modified) (2 diffs)
-
Triangles.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.cpp
r2811 r2813 1074 1074 } 1075 1075 1076 Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j)1077 { return quadtree->NearestVertex(i,j); }1078 1079 1080 void Triangles::PreInit(Int4 inbvx,char *fname)1081 {1082 long int verbosity=0;1083 1084 srand(19999999);1085 OnDisk =0;1086 NbRef=0;1087 // allocGeometry=0;1088 identity=0;1089 NbOfTriangleSearchFind =0;1090 NbOfSwapTriangle =0;1091 nbiv=0;1092 nbv=0;1093 nbvx=inbvx;1094 nbt=0;1095 NbOfQuad = 0;1096 nbtx=2*inbvx-2;1097 NbSubDomains=0;1098 NbVertexOnBThVertex=0;1099 NbVertexOnBThEdge=0;1100 VertexOnBThVertex=0;1101 VertexOnBThEdge=0;1102 1103 NbCrackedVertices=0;1104 NbCrackedEdges =0;1105 CrackedEdges =0;1106 nbe = 0;1107 name = fname ;1108 1109 if (inbvx) {1110 vertices=new Vertex[nbvx];1111 assert(vertices);1112 ordre=new (Vertex* [nbvx]);1113 assert(ordre);1114 triangles=new Triangle[nbtx];1115 assert(triangles);}1116 else {1117 vertices=0;1118 ordre=0;1119 triangles=0;1120 nbtx=0;1121 }1122 if ( name || inbvx)1123 {1124 time_t timer =time(0);1125 char buf[70];1126 strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer));1127 counter++;1128 char countbuf[30];1129 sprintf(countbuf,"%d",counter);1130 int lg =0 ;1131 if (&BTh != this && BTh.name)1132 lg = strlen(BTh.name)+4;1133 identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2 + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)];1134 identity[0]=0;1135 if (lg)1136 strcat(strcat(strcat(identity,"B="),BTh.name),", ");1137 1138 if (Gh.name)1139 strcat(strcat(identity,"G="),Gh.name);1140 strcat(strcat(identity,";"),countbuf);1141 strcat(identity,buf);1142 // cout << "New MAILLAGE "<< identity << endl;1143 }1144 1145 quadtree=0;1146 // edgescomponante=0;1147 edges=0;1148 VerticesOnGeomVertex=0;1149 VerticesOnGeomEdge=0;1150 NbVerticesOnGeomVertex=0;1151 NbVerticesOnGeomEdge=0;1152 // nbMaxIntersectionTriangles=0;1153 // lIntTria;1154 subdomains=0;1155 NbSubDomains=0;1156 // Meshbegin = vertices;1157 // Meshend = vertices + nbvx;1158 if (verbosity>98)1159 cout << "Triangles::PreInit() " << nbvx << " " << nbtx1160 << " " << vertices1161 << " " << ordre << " " << triangles <<endl;1162 }1163 1164 1076 1165 1077 double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) { -
issm/trunk/src/c/Bamgx/Mesh2.h
r2810 r2813 961 961 }; // End Class Geometry 962 962 963 //From Metric.cpp 964 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]) 965 { return A[0] * ( B[1]*C[2]-B[2]*C[1]) 966 - A[1] * ( B[0]*C[2]-B[2]*C[0]) 967 + A[2] * ( B[0]*C[1]-B[1]*C[0]); 968 } 969 963 970 ///////////////////////////////////////////////////////////////////////////////////// 964 971 ///////////////////////////////////////////////////////////////////////////////////// -
issm/trunk/src/c/Bamgx/Metric.cpp
r2805 r2813 31 31 32 32 namespace bamg { 33 34 35 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3])36 { return A[0] * ( B[1]*C[2]-B[2]*C[1])37 - A[1] * ( B[0]*C[2]-B[2]*C[0])38 + A[2] * ( B[0]*C[1]-B[1]*C[0]);39 }40 33 41 34 SaveMetricInterpole LastMetricInterpole; … … 206 199 return r; 207 200 } 208 void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0)209 210 {211 long int verbosity=0;212 213 if(verbosity>1)214 cout << " -- IntersectGeomMetric geometric err=" << err << (iso ? " iso " : " aniso " ) << endl;215 Real8 ss[2]={0.00001,0.99999};216 Real8 errC = 2*sqrt(2*err);217 Real8 hmax = Gh.MaximalHmax();218 Real8 hmin = Gh.MinimalHmin();219 Real8 maxaniso = 1e6;220 assert(hmax>0);221 SetVertexFieldOn();222 if (errC > 1) errC = 1;223 for (Int4 i=0;i<nbe;i++)224 for (int j=0;j<2;j++)225 {226 227 Vertex V;228 VertexOnGeom GV;229 // cerr << Number(edges[i]) << " " << ss[j] << endl;230 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);231 {232 GeometricalEdge * eg = GV;233 Real8 s = GV;234 R2 tg;235 // cerr << i << " " << j << " " << Number(V) << " on = "236 // << Gh.Number(eg) << " at s = " << s << " " << endl;237 Real8 R1= eg->R1tg(s,tg);238 // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x "239 // << V.r << errC/ Max(R1,1e-20) << " hold=" <<V.m(tg) << " " << endl;240 Real8 ht = hmax;241 if (R1>1.0e-20)242 { // err relative to the length of the edge243 ht = Min(Max(errC/R1,hmin),hmax);244 }245 Real8 hn = iso? ht : Min(hmax,ht*maxaniso);246 //cerr << ht << " " << hn << "m=" << edges[i][j].m << endl;247 assert(ht>0 && hn>0);248 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);249 //cerr << " : " ;250 Metric MVp(Vp);251 // cerr << " : " << MVp << endl;252 edges[i][j].m.IntersectWith(MVp);253 //cerr << " . " << endl;254 }255 256 }257 // the problem is for the vertex on vertex258 259 }260 /*261 void Triangles::BoundAnisotropy(Real8 anisomax)262 {263 if (verbosity > 1)264 cout << " -- BoundAnisotropy by " << anisomax << endl;265 Real8 h1=1.e30,h2=1e-30,rx=0;266 Real8 coef = 1./(anisomax*anisomax);267 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;268 for (Int4 i=0;i<nbv;i++)269 {270 271 MatVVP2x2 Vp(vertices[i]);272 273 h1=Min(h1,Vp.lmin());274 h2=Max(h2,Vp.lmax());275 rx = Max(rx,Vp.Aniso2());276 277 Vp.BoundAniso2(coef);278 279 hn1=Min(hn1,Vp.lmin());280 hn2=Max(hn2,Vp.lmax());281 rnx = Max(rnx,Vp.Aniso2());282 283 284 vertices[i].m = Vp;285 286 }287 288 if (verbosity>2)289 {290 cout << " input : Hmin = " << sqrt(1/h2) << " Hmax = " << sqrt(1/h1)291 << " factor of anisotropy max = " << sqrt(rx) << endl;292 cout << " output: Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1)293 << " factor of anisotropy max = " << sqrt(rnx) << endl;294 }295 }296 */297 void Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso)298 {299 long int verbosity=0;300 301 double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));302 if (verbosity > 1)303 cout << " -- BoundAnisotropy by " << anisomax << endl;304 Real8 h1=1.e30,h2=1e-30,rx=0;305 Real8 coef = 1./(anisomax*anisomax);306 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;307 for (Int4 i=0;i<nbv;i++)308 {309 310 MatVVP2x2 Vp(vertices[i]);311 double lmax=Vp.lmax();312 h1=Min(h1,Vp.lmin());313 h2=Max(h2,Vp.lmax());314 rx = Max(rx,Vp.Aniso2());315 316 Vp *= Min(lminaniso,lmax)/lmax;317 318 Vp.BoundAniso2(coef);319 320 hn1=Min(hn1,Vp.lmin());321 hn2=Max(hn2,Vp.lmax());322 rnx = Max(rnx,Vp.Aniso2());323 324 325 vertices[i].m = Vp;326 327 }328 329 if (verbosity>2)330 {331 cout << " input : Hmin = " << sqrt(1/h2) << " Hmax = " << sqrt(1/h1)332 << " factor of anisotropy max = " << sqrt(rx) << endl;333 cout << " output: Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1)334 << " factor of anisotropy max = " << sqrt(rnx) << endl;335 }336 }337 void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols,338 const Real8 hmin1,const Real8 hmax1,const Real8 coef,339 const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,340 const int DoNormalisation,const double power,const int choice)341 { // the array of solution s is store342 // sol0,sol1,...,soln on vertex 0343 // sol0,sol1,...,soln on vertex 1344 // etc.345 // choise = 0 => H is computed with green formule346 // otherwise => H is computed from P2 on 4T347 const int dim = 2;348 349 long int verbosity=0;350 351 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ;352 353 // computation of the nb of field354 Int4 ntmp = 0;355 if (typsols)356 {357 for (Int4 i=0;i<nbsol;i++)358 ntmp += sizeoftype[typsols[i]];359 }360 else361 ntmp = nbsol;362 363 // n is the total number of fields364 365 const Int4 n = ntmp;366 367 Int4 i,k,iA,iB,iC,iv;368 R2 O(0,0);369 int RelativeMetric = CutOff>1e-30;370 Real8 hmin = Max(hmin1,MinimalHmin());371 Real8 hmax = Min(hmax1,MaximalHmax());372 Real8 coef2 = 1/(coef*coef);373 374 if(verbosity>1)375 {376 cout << " -- Construction of Metric: Nb of field. " << n << " nbt = "377 << nbt << " nbv= " << nbv378 << " coef = " << coef << endl379 << " hmin = " << hmin << " hmax=" << hmax380 << " anisomax = " << anisomax << " Nb Jacobi " << NbJacobi << " Power = " << power ;381 if (RelativeMetric)382 cout << " RelativeErr with CutOff= " << CutOff << endl;383 else384 cout << " Absolute Err" <<endl;385 }386 double *ss=(double*)s;//, *ssiii = ss;387 388 double sA,sB,sC;389 390 Real8 *detT = new Real8[nbt];391 Real8 *Mmass= new Real8[nbv];392 Real8 *Mmassxx= new Real8[nbv];393 Real8 *dxdx= new Real8[nbv];394 Real8 *dxdy= new Real8[nbv];395 Real8 *dydy= new Real8[nbv];396 Real8 *workT= new Real8[nbt];397 Real8 *workV= new Real8[nbv];398 int *OnBoundary = new int[nbv];399 for (iv=0;iv<nbv;iv++)400 {401 Mmass[iv]=0;402 OnBoundary[iv]=0;403 Mmassxx[iv]=0;404 }405 406 for (i=0;i<nbt;i++)407 if(triangles[i].link) // the real triangles408 {409 const Triangle &t=triangles[i];410 // coor of 3 vertices411 R2 A=t[0];412 R2 B=t[1];413 R2 C=t[2];414 415 416 // number of the 3 vertices417 iA = Number(t[0]);418 iB = Number(t[1]);419 iC = Number(t[2]);420 421 Real8 dett = bamg::Area2(A,B,C);422 detT[i]=dett;423 dett /= 6;424 425 // construction of on boundary426 int nbb =0;427 for(int j=0;j<3;j++)428 {429 Triangle *ta=t.Adj(j);430 if ( ! ta || !ta->link) // no adj triangle => edge on boundary431 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1,432 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1,433 nbb++;434 }435 436 workT[i] = nbb;437 Mmass[iA] += dett;438 Mmass[iB] += dett;439 Mmass[iC] += dett;440 441 if((nbb==0)|| !choice)442 {443 Mmassxx[iA] += dett;444 Mmassxx[iB] += dett;445 Mmassxx[iC] += dett;446 }447 }448 else449 workT[i]=-1;450 451 // for (Int4 kcount=0;kcount<n;kcount++,ss++)452 for (Int4 nusol=0;nusol<nbsol;nusol++)453 { //for all Solution454 455 Real8 smin=ss[0],smax=ss[0];456 457 Real8 h1=1.e30,h2=1e-30,rx=0;458 Real8 coef = 1./(anisomax*anisomax);459 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;460 int nbfield = typsols? sizeoftype[typsols[nusol]] : 1;461 if (nbfield == 1)462 for ( iv=0,k=0; iv<nbv; iv++,k+=n )463 {464 dxdx[iv]=dxdy[iv]=dydy[iv]=0;465 smin=Min(smin,ss[k]);466 smax=Max(smax,ss[k]);467 }468 else469 {470 // cas vectoriel471 for ( iv=0,k=0; iv<nbv; iv++,k+=n )472 {473 double v=0;474 for (int i=0;i<nbfield;i++)475 v += ss[k+i]*ss[k+i];476 v = sqrt(v);477 smin=Min(smin,v);478 smax=Max(smax,v);479 }480 }481 Real8 sdelta = smax-smin;482 Real8 absmax=Max(Abs(smin),Abs(smax));483 Real8 cnorm = DoNormalisation ? coef2/sdelta : coef2;484 485 if(verbosity>2)486 cout << " Solution " << nusol << " Min = " << smin << " Max = "487 << smax << " Delta =" << sdelta << " cnorm = " << cnorm << " Nb of fields =" << nbfield << endl;488 489 490 if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1))491 {492 if (verbosity>2)493 cout << " Solution " << nusol << " is constant. We skip. "494 << " Min = " << smin << " Max = " << smax << endl;495 continue;496 }497 498 double *sf = ss;499 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++)500 {501 for ( iv=0,k=0; iv<nbv; iv++,k+=n )502 dxdx[iv]=dxdy[iv]=dydy[iv]=0;503 for (i=0;i<nbt;i++)504 if(triangles[i].link)505 {// for real all triangles506 // coor of 3 vertices507 R2 A=triangles[i][0];508 R2 B=triangles[i][1];509 R2 C=triangles[i][2];510 511 512 // warning the normal is internal and the513 // size is the length of the edge514 R2 nAB = Orthogonal(B-A);515 R2 nBC = Orthogonal(C-B);516 R2 nCA = Orthogonal(A-C);517 // remark : nAB + nBC + nCA == 0518 519 // number of the 3 vertices520 iA = Number(triangles[i][0]);521 iB = Number(triangles[i][1]);522 iC = Number(triangles[i][2]);523 524 // for the test of boundary edge525 // the 3 adj triangles526 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);527 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);528 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);529 530 // value of the P1 fonction on 3 vertices531 sA = ss[iA*n];532 sB = ss[iB*n];533 sC = ss[iC*n];534 535 R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ;536 if(choice)537 {538 int nbb = 0;539 Real8 dd = detT[i];540 Real8 lla,llb,llc,llf;541 Real8 taa[3][3],bb[3];542 // construction of the trans of lin system543 for (int j=0;j<3;j++)544 {545 int ie = OppositeEdge[j];546 TriangleAdjacent ta = triangles[i].Adj(ie);547 Triangle *tt = ta;548 if (tt && tt->link)549 {550 Vertex &v = *ta.OppositeVertex();551 R2 V = v;552 Int4 iV = Number(v);553 Real8 lA = bamg::Area2(V,B,C)/dd;554 Real8 lB = bamg::Area2(A,V,C)/dd;555 Real8 lC = bamg::Area2(A,B,V)/dd;556 taa[0][j] = lB*lC;557 taa[1][j] = lC*lA;558 taa[2][j] = lA*lB;559 //Real8 xx = V.x-V.y;560 //Real8 yy = V.x + V.y;561 //cout << " iv " << ss[iV*n] << " == " << (8*xx*xx+yy*yy)562 // << " l = " << lA << " " << lB << " " << lC563 // << " = " << lA+lB+lC << " " << V << " == " << A*lA+B*lB+C*lC << endl;564 565 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;566 567 bb[j] = ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ;568 }569 else570 {571 nbb++;572 taa[0][j]=0;573 taa[1][j]=0;574 taa[2][j]=0;575 taa[j][j]=1;576 bb[j]=0;577 }578 }579 580 // resolution of 3x3 lineaire system transpose581 Real8 det33 = det3x3(taa[0],taa[1],taa[2]);582 Real8 cBC = det3x3(bb,taa[1],taa[2]);583 Real8 cCA = det3x3(taa[0],bb,taa[2]);584 Real8 cAB = det3x3(taa[0],taa[1],bb);585 586 assert(det33);587 // det33=1;588 // verif589 // cout << " " << (taa[0][0]*cBC + taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ;590 // cout << " " << (taa[0][1]*cBC + taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1];591 // cout << " " << (taa[0][2]*cBC + taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2]592 // << " -- " ;593 //cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB + llb*llc* cBC + llc*lla*cCA)/det33594 // << " == " << llf << endl;595 // computation of the gradient in the element596 597 // H( li*lj) = grad li grad lj + grad lj grad lj598 // grad li = njk / detT ; with i j k ={A,B,C)599 Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);600 Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);601 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x)602 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);603 Real8 coef = 1.0/(3*dd*det33);604 Real8 coef2 = 2*coef;605 // cout << " H = " << Hxx << " " << Hyy << " " << Hxy/2 << " coef2 = " << coef2 << endl;606 Hxx *= coef2;607 Hyy *= coef2;608 Hxy *= coef2;609 //cout << i << " H = " << 3*Hxx/dd << " " << 3*Hyy/dd << " " << 3*Hxy/(dd*2) << " nbb = " << nbb << endl;610 if(nbb==0)611 {612 dxdx[iA] += Hxx;613 dydy[iA] += Hyy;614 dxdy[iA] += Hxy;615 616 dxdx[iB] += Hxx;617 dydy[iB] += Hyy;618 dxdy[iB] += Hxy;619 620 dxdx[iC] += Hxx;621 dydy[iC] += Hyy;622 dxdy[iC] += Hxy;623 }624 625 }626 else627 {628 629 // if edge on boundary no contribution => normal = 0630 if ( ! tBC || ! tBC->link ) nBC = O;631 if ( ! tCA || ! tCA->link ) nCA = O;632 if ( ! tAB || ! tAB->link ) nAB = O;633 634 // remark we forgot a 1/2 because635 // $\\int_{edge} w_i = 1/2 $ if $i$ is in edge636 // 0 if not637 // if we don't take the boundary638 // dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;639 640 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;641 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;642 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;643 644 // warning optimization (1) the divide by 2 is done on the metrix construction645 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;646 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;647 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ;648 649 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;650 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;651 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;652 }653 654 } // for real all triangles655 Int4 kk=0;656 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )657 if(Mmassxx[iv]>0)658 {659 dxdx[iv] /= 2*Mmassxx[iv];660 // warning optimization (1) on term dxdy[iv]*ci/2661 dxdy[iv] /= 4*Mmassxx[iv];662 dydy[iv] /= 2*Mmassxx[iv];663 // Compute the matrix with abs(eigen value)664 Metric M(dxdx[iv], dxdy[iv], dydy[iv]);665 MatVVP2x2 Vp(M);666 //cout <<iv << " M = " << M << " aniso= " << Vp.Aniso() ;667 Vp.Abs();668 M = Vp;669 dxdx[iv] = M.a11;670 dxdy[iv] = M.a21;671 dydy[iv] = M.a22;672 // cout << " (abs) iv M = " << M << " aniso= " << Vp.Aniso() <<endl;673 }674 else kk++;675 676 677 // correction of second derivate678 // by a laplacien679 680 Real8 *d2[3] = { dxdx, dxdy, dydy};681 Real8 *dd;682 for (int xy = 0;xy<3;xy++)683 {684 dd = d2[xy];685 // do leat 2 iteration for boundary problem686 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++)687 {688 for (i=0;i<nbt;i++)689 if(triangles[i].link) // the real triangles690 {691 // number of the 3 vertices692 iA = Number(triangles[i][0]);693 iB = Number(triangles[i][1]);694 iC = Number(triangles[i][2]);695 Real8 cc=3;696 if(ijacobi==0)697 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);698 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;699 }700 for (iv=0;iv<nbv;iv++)701 workV[iv]=0;702 703 for (i=0;i<nbt;i++)704 if(triangles[i].link) // the real triangles705 {706 // number of the 3 vertices707 iA = Number(triangles[i][0]);708 iB = Number(triangles[i][1]);709 iC = Number(triangles[i][2]);710 Real8 cc = workT[i]*detT[i];711 workV[iA] += cc;712 workV[iB] += cc;713 workV[iC] += cc;714 }715 716 for (iv=0;iv<nbv;iv++)717 if( ijacobi<NbJacobi || OnBoundary[iv])718 dd[iv] = workV[iv]/(Mmass[iv]*6);719 720 721 }722 723 724 }725 726 // constuction of the metrix from the Hessian dxdx. dxdy,dydy727 728 Real8 rCutOff=CutOff*absmax;// relative cut off729 730 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )731 { // for all vertices732 //{733 //Metric M(dxdx[iv], dxdy[iv], dydy[iv]);734 // MatVVP2x2 Vp(M);735 //cout << " iv M="<< M << " Vp = " << Vp << " aniso " << Vp.Aniso() << endl;736 //}737 MetricIso Miso;738 // new code to compute ci ---739 Real8 ci ;740 if (RelativeMetric)741 { // compute the norm of the solution742 double xx =0,*sfk=sf+k;743 for (int ifield=0;ifield<nbfield;ifield++,sfk++)744 xx += *sfk* *sfk;745 xx=sqrt(xx);746 ci = coef2/Max(xx,rCutOff);747 }748 else ci = cnorm;749 750 // old751 // Real8 ci = RelativeMetric ? coef2/(Max(Abs(ss[k]),rCutOff)) : cnorm ;752 // modif F Hecht 101099753 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci);754 MatVVP2x2 Vp(Miv);755 756 Vp.Abs();757 if(power!=1.0)758 Vp.pow(power);759 760 761 762 h1=Min(h1,Vp.lmin());763 h2=Max(h2,Vp.lmax());764 765 Vp.Maxh(hmin);766 Vp.Minh(hmax);767 768 rx = Max(rx,Vp.Aniso2());769 770 Vp.BoundAniso2(coef);771 772 hn1=Min(hn1,Vp.lmin());773 hn2=Max(hn2,Vp.lmax());774 rnx = Max(rnx,Vp.Aniso2());775 776 Metric MVp(Vp);777 vertices[iv].m.IntersectWith(MVp);778 }// for all vertices779 if (verbosity>2)780 {781 cout << " Field " << nufield << " of solution " << nusol << endl;782 cout << " before bounding : Hmin = " << sqrt(1/h2) << " Hmax = "783 << sqrt(1/h1) << " factor of anisotropy max = " << sqrt(rx) << endl;784 cout << " after bounding : Hmin = " << sqrt(1/hn2) << " Hmax = "785 << sqrt(1/hn1) << " factor of anisotropy max = " << sqrt(rnx) << endl;786 }787 } // end of for all field788 }// end for all solution789 790 delete [] detT;791 delete [] Mmass;792 delete [] dxdx;793 delete [] dxdy;794 delete [] dydy;795 delete [] workT;796 delete [] workV;797 delete [] Mmassxx;798 delete [] OnBoundary;799 800 }801 802 803 void Triangles::ReadMetric(BamgOpts* bamgopts,const Real8 hmin1=1.0e-30,const Real8 hmax1=1.0e30,const Real8 coef=1)804 {805 int i,j;806 807 if(bamgopts->verbose>3) printf(" processing metric\n");808 809 Real8 hmin = Max(hmin1,MinimalHmin());810 Real8 hmax = Min(hmax1,MaximalHmax());811 812 //for now we only use j==3813 j=3;814 815 for (i=0;i<nbv;i++){816 Real8 h;817 if (j == 1){818 h=bamgopts->metric[i];819 vertices[i].m=Metric(Max(hmin,Min(hmax, h*coef)));820 }821 else if (j==3){822 Real8 a,b,c;823 a=bamgopts->metric[i*3+0];824 b=bamgopts->metric[i*3+1];825 c=bamgopts->metric[i*3+2];826 MetricAnIso M(a,b,c);827 MatVVP2x2 Vp(M/coef);828 829 Vp.Maxh(hmin);830 Vp.Minh(hmax);831 vertices[i].m = Vp;832 }833 }834 }835 836 void Triangles::WriteMetric(ostream & f,int iso)837 {838 if (iso)839 {840 f << nbv <<" " << 1 << endl ;841 for (Int4 iv=0;iv<nbv;iv++)842 {843 MatVVP2x2 V=vertices[iv].m;844 f << V.hmin() << endl;845 }846 }847 else848 {849 f << nbv <<" " << 3 << endl ;850 for (Int4 iv=0;iv<nbv;iv++)851 f << vertices[iv].m.a11 << " "852 << vertices[iv].m.a21 << " "853 << vertices[iv].m.a22 << endl;854 }855 }856 void Triangles::MaxSubDivision(Real8 maxsubdiv)857 {858 long int verbosity=0;859 860 const Real8 maxsubdiv2 = maxsubdiv*maxsubdiv;861 if(verbosity>1)862 cout << " -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv << endl ;863 // for all the edges864 // if the len of the edge is to long865 Int4 it,nbchange=0;866 Real8 lmax=0;867 for (it=0;it<nbt;it++)868 {869 Triangle &t=triangles[it];870 for (int j=0;j<3;j++)871 {872 Triangle &tt = *t.TriangleAdj(j);873 if ( ! &tt || it < Number(tt) && ( tt.link || t.link))874 {875 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];876 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];877 R2 AB= (R2) v1-(R2) v0;878 Metric M = v0;879 Real8 l = M(AB,AB);880 lmax = Max(lmax,l);881 if(l> maxsubdiv2)882 { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M883 Real8 lc = M(AC,AC);884 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;885 D2xD2 Rt1(Rt.inv());886 D2xD2 D(maxsubdiv2,0,0,lc);887 D2xD2 MM = Rt1*D*Rt1.t();888 v0.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);889 nbchange++;890 }891 M = v1;892 l = M(AB,AB);893 lmax = Max(lmax,l);894 if(l> maxsubdiv2)895 { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M896 Real8 lc = M(AC,AC);897 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;898 D2xD2 Rt1(Rt.inv());899 D2xD2 D(maxsubdiv2,0,0,lc);900 D2xD2 MM = Rt1*D*Rt1.t();901 v1.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);902 nbchange++;903 }904 905 906 }907 }908 }909 if(verbosity>3)910 cout << " Nb of metric change = " << nbchange911 << " Max of the subdivision of a edges before change = " << sqrt(lmax) << endl;912 913 }914 915 void Triangles::SmoothMetric(Real8 raisonmax)916 {917 long int verbosity=0;918 919 if(raisonmax<1.1) return;920 if(verbosity > 1)921 cout << " -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl;922 ReMakeTriangleContainingTheVertex();923 Int4 i,j,kch,kk,ip;924 Int4 *first_np_or_next_t0 = new Int4[nbv];925 Int4 *first_np_or_next_t1 = new Int4[nbv];926 Int4 Head0 =0,Head1=-1;927 Real8 logseuil= log(raisonmax);928 929 for(i=0;i<nbv-1;i++)930 first_np_or_next_t0[i]=i+1;931 first_np_or_next_t0[nbv-1]=-1;// end;932 for(i=0;i<nbv;i++)933 first_np_or_next_t1[i]=-1;934 kk=0;935 while (Head0>=0&& kk++<100) {936 kch=0;937 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1)938 { // pour tous les triangles autour du sommet s939 // cout << kk << " i = " << i << " " << ip << endl;940 register Triangle * t= vertices[i].t;941 assert(t);942 Vertex & vi = vertices[i];943 TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);944 Vertex *pvj0 = ta.EdgeVertex(0);945 while (1) {946 // cout << i << " " << Number(ta.EdgeVertex(0)) << " "947 // << Number(ta.EdgeVertex(1)) << " ---> " ;948 ta=Previous(Adj(ta));949 // cout << Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl;950 assert(vertices+i == ta.EdgeVertex(1));951 Vertex & vj = *(ta.EdgeVertex(0));952 if ( &vj ) {953 j= &vj-vertices;954 assert(j>=0 && j < nbv);955 R2 Aij = (R2) vj - (R2) vi;956 Real8 ll = Norme2(Aij);957 if (0) {958 Real8 hi = ll/vi.m(Aij);959 Real8 hj = ll/vj.m(Aij);960 if(hi < hj)961 {962 Real8 dh=(hj-hi)/ll;963 //cout << " dh = " << dh << endl;964 if (dh>logseuil) {965 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));966 if(first_np_or_next_t1[j]<0)967 kch++,first_np_or_next_t1[j]=Head1,Head1=j;968 }969 }970 }971 else972 {973 Real8 li = vi.m(Aij);974 //Real8 lj = vj.m(Aij);975 // if ( i == 2 || j == 2)976 // cout << " inter " << i << " " << j << " " << ((1 +logseuil*li)) << endl;977 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )978 //if( vj.m.IntersectWith(vi.m*(lj/li/(1 +logseuil*lj))) )979 if(first_np_or_next_t1[j]<0) // if the metrix change980 kch++,first_np_or_next_t1[j]=Head1,Head1=j;981 }982 }983 if ( &vj == pvj0 ) break;984 }985 }986 Head0 = Head1;987 Head1 = -1;988 Exchange(first_np_or_next_t0,first_np_or_next_t1);989 if(verbosity>5)990 cout << " Iteration " << kk << " Nb de vertices with change " << kch << endl;991 }992 if(verbosity>2 && verbosity < 5)993 cout << " Nb of Loop " << kch << endl;994 delete [] first_np_or_next_t0;995 delete [] first_np_or_next_t1;996 }997 201 998 202 Real8 LengthInterpole(const MetricAnIso Ma,const MetricAnIso Mb, R2 AB) -
issm/trunk/src/c/Bamgx/Triangles.cpp
r2811 r2813 233 233 234 234 /*IO*/ 235 /*FUNCTION Triangles::WriteMesh {{{1*/236 void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){237 238 int i,j;239 int verbose;240 241 verbose=bamgopts->verbose;242 Int4 *reft = new Int4[nbt];243 Int4 nbInT = ConsRefTriangle(reft);244 245 //Vertices246 if(verbose>3) printf(" writing Vertices\n");247 bamgmesh->NumVertices=nbv;248 xfree((void**)&bamgmesh->Vertices);249 if (nbv){250 bamgmesh->Vertices=(double*)xmalloc(3*nbv*sizeof(double));251 for (i=0;i<nbv;i++){252 bamgmesh->Vertices[i*3+0]=vertices[i].r.x;253 bamgmesh->Vertices[i*3+1]=vertices[i].r.y;254 bamgmesh->Vertices[i*3+2]=vertices[i].ref();255 }256 }257 else{258 bamgmesh->Vertices=NULL;259 }260 261 //Edges262 if(verbose>3) printf(" writing Edges\n");263 bamgmesh->NumEdges=nbe;264 xfree((void**)&bamgmesh->Edges);265 if (nbe){266 bamgmesh->Edges=(double*)xmalloc(3*nbe*sizeof(double));267 for (i=0;i<nbe;i++){268 bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing269 bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing270 bamgmesh->Edges[i*3+2]=edges[i].ref;271 }272 }273 else{274 bamgmesh->Edges=NULL;275 }276 277 //CrackedEdges278 if(verbose>3) printf(" writing CrackedEdges\n");279 bamgmesh->NumCrackedEdges=NbCrackedEdges;280 xfree((void**)&bamgmesh->CrackedEdges);281 if (NbCrackedEdges){282 bamgmesh->CrackedEdges=(double*)xmalloc(2*NbCrackedEdges*sizeof(double));283 for (i=0;i<NbCrackedEdges;i++){284 bamgmesh->CrackedEdges[i*2+0]=Number(CrackedEdges[i].a.edge)+1; //back to M indexing285 bamgmesh->CrackedEdges[i*2+1]=Number(CrackedEdges[i].b.edge)+1; //back to M indexing286 }287 }288 else{289 bamgmesh->CrackedEdges=NULL;290 }291 292 //Triangles293 if(verbose>3) printf(" writing Triangles\n");294 Int4 k=nbInT-NbOfQuad*2;295 Int4 num=0;296 bamgmesh->NumTriangles=k;297 xfree((void**)&bamgmesh->Triangles);298 if (k){299 bamgmesh->Triangles=(double*)xmalloc(4*k*sizeof(double));300 for (i=0;i<nbt;i++){301 Triangle &t=triangles[i];302 if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){303 bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing304 bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing305 bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing306 bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;307 num=num+1;308 }309 }310 }311 else{312 bamgmesh->Triangles=NULL;313 }314 315 //Quadrilaterals316 if(verbose>3) printf(" writing Quadrilaterals\n");317 bamgmesh->NumQuadrilaterals=NbOfQuad;318 xfree((void**)&bamgmesh->Quadrilaterals);319 if (NbOfQuad){320 bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double));321 for (i=0;i<nbt;i++){322 Triangle &t =triangles[i];323 Triangle* ta;324 Vertex *v0,*v1,*v2,*v3;325 if (reft[i]<0) continue;326 if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) {327 bamgmesh->Quadrilaterals[i*5+0]=Number(v0)+1; //back to M indexing328 bamgmesh->Quadrilaterals[i*5+1]=Number(v1)+1; //back to M indexing329 bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing330 bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing331 bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref;332 }333 }334 }335 else{336 bamgmesh->Quadrilaterals=NULL;337 }338 339 //SubDomains340 if(verbose>3) printf(" writing SubDomains\n");341 bamgmesh->NumSubDomains=NbSubDomains;342 xfree((void**)&bamgmesh->SubDomains);343 if (NbSubDomains){344 bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));345 for (i=0;i<NbSubDomains;i++){346 bamgmesh->SubDomains[i*4+0]=3;347 bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)];348 bamgmesh->SubDomains[i*4+2]=1;349 bamgmesh->SubDomains[i*4+3]=subdomains[i].ref;350 }351 }352 else{353 bamgmesh->SubDomains=NULL;354 }355 356 //SubDomainsFromGeom357 if(verbose>3) printf(" writing SubDomainsFromGeom\n");358 bamgmesh->NumSubDomainsFromGeom=Gh.NbSubDomains;359 xfree((void**)&bamgmesh->SubDomainsFromGeom);360 if (Gh.NbSubDomains){361 bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double));362 for (i=0;i<Gh.NbSubDomains;i++){363 bamgmesh->SubDomainsFromGeom[i*4+0]=2;364 bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing365 bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens;366 bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref;367 }368 }369 else{370 bamgmesh->SubDomainsFromGeom=NULL;371 }372 373 //VerticesOnGeomVertex374 if(verbose>3) printf(" writing VerticesOnGeometricVertex\n");375 bamgmesh->NumVerticesOnGeometricVertex=NbVerticesOnGeomVertex;376 xfree((void**)&bamgmesh->VerticesOnGeometricVertex);377 if (NbVerticesOnGeomVertex){378 bamgmesh->VerticesOnGeometricVertex=(double*)xmalloc(2*NbVerticesOnGeomVertex*sizeof(double));379 for (i=0;i<NbVerticesOnGeomVertex;i++){380 VertexOnGeom &v=VerticesOnGeomVertex[i];381 if (!v.OnGeomVertex()){382 throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not"));383 }384 bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing385 bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing386 }387 }388 else{389 bamgmesh->VerticesOnGeometricVertex=NULL;390 }391 392 //VertexOnGeometricEdge393 if(verbose>3) printf(" writing VerticesOnGeometricEdge\n");394 bamgmesh->NumVerticesOnGeometricEdge=NbVerticesOnGeomEdge;395 xfree((void**)&bamgmesh->VerticesOnGeometricEdge);396 if (NbVerticesOnGeomEdge){397 bamgmesh->VerticesOnGeometricEdge=(double*)xmalloc(3*NbVerticesOnGeomEdge*sizeof(double));398 for (i=0;i<NbVerticesOnGeomEdge;i++){399 const VertexOnGeom &v=VerticesOnGeomEdge[i];400 if (!v.OnGeomEdge()){401 throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricEdge is actually not"));402 }403 bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing404 bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing405 bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v;406 }407 }408 else{409 bamgmesh->VerticesOnGeometricEdge=NULL;410 }411 412 //EdgesOnGeometricEdge413 if(verbose>3) printf(" writing EdgesOnGeometricEdge\n");414 k=0;415 for (i=0;i<nbe;i++){416 if (edges[i].on) k=k+1;417 }418 bamgmesh->NumEdgesOnGeometricEdge=k;419 xfree((void**)&bamgmesh->EdgesOnGeometricEdge);420 if (k){421 bamgmesh->EdgesOnGeometricEdge=(double*)xmalloc(2*(int)k*sizeof(double));422 int count=0;423 for (i=0;i<nbe;i++){424 if (edges[i].on){425 bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing426 bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].on)+1; //back to Matlab indexing427 count=count+1;428 }429 }430 }431 else{432 bamgmesh->EdgesOnGeometricEdge=NULL;433 }434 }435 /*}}}1*/436 235 /*FUNCTION Triangles::ReadMesh{{{1*/ 437 236 void Triangles::ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){ … … 662 461 Gh.AfterRead(); 663 462 } 463 } 464 /*}}}1*/ 465 /*FUNCTION Triangles::WriteMesh {{{1*/ 466 void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){ 467 468 int i,j; 469 int verbose; 470 471 verbose=bamgopts->verbose; 472 Int4 *reft = new Int4[nbt]; 473 Int4 nbInT = ConsRefTriangle(reft); 474 475 //Vertices 476 if(verbose>3) printf(" writing Vertices\n"); 477 bamgmesh->NumVertices=nbv; 478 xfree((void**)&bamgmesh->Vertices); 479 if (nbv){ 480 bamgmesh->Vertices=(double*)xmalloc(3*nbv*sizeof(double)); 481 for (i=0;i<nbv;i++){ 482 bamgmesh->Vertices[i*3+0]=vertices[i].r.x; 483 bamgmesh->Vertices[i*3+1]=vertices[i].r.y; 484 bamgmesh->Vertices[i*3+2]=vertices[i].ref(); 485 } 486 } 487 else{ 488 bamgmesh->Vertices=NULL; 489 } 490 491 //Edges 492 if(verbose>3) printf(" writing Edges\n"); 493 bamgmesh->NumEdges=nbe; 494 xfree((void**)&bamgmesh->Edges); 495 if (nbe){ 496 bamgmesh->Edges=(double*)xmalloc(3*nbe*sizeof(double)); 497 for (i=0;i<nbe;i++){ 498 bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing 499 bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing 500 bamgmesh->Edges[i*3+2]=edges[i].ref; 501 } 502 } 503 else{ 504 bamgmesh->Edges=NULL; 505 } 506 507 //CrackedEdges 508 if(verbose>3) printf(" writing CrackedEdges\n"); 509 bamgmesh->NumCrackedEdges=NbCrackedEdges; 510 xfree((void**)&bamgmesh->CrackedEdges); 511 if (NbCrackedEdges){ 512 bamgmesh->CrackedEdges=(double*)xmalloc(2*NbCrackedEdges*sizeof(double)); 513 for (i=0;i<NbCrackedEdges;i++){ 514 bamgmesh->CrackedEdges[i*2+0]=Number(CrackedEdges[i].a.edge)+1; //back to M indexing 515 bamgmesh->CrackedEdges[i*2+1]=Number(CrackedEdges[i].b.edge)+1; //back to M indexing 516 } 517 } 518 else{ 519 bamgmesh->CrackedEdges=NULL; 520 } 521 522 //Triangles 523 if(verbose>3) printf(" writing Triangles\n"); 524 Int4 k=nbInT-NbOfQuad*2; 525 Int4 num=0; 526 bamgmesh->NumTriangles=k; 527 xfree((void**)&bamgmesh->Triangles); 528 if (k){ 529 bamgmesh->Triangles=(double*)xmalloc(4*k*sizeof(double)); 530 for (i=0;i<nbt;i++){ 531 Triangle &t=triangles[i]; 532 if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){ 533 bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing 534 bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing 535 bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing 536 bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref; 537 num=num+1; 538 } 539 } 540 } 541 else{ 542 bamgmesh->Triangles=NULL; 543 } 544 545 //Quadrilaterals 546 if(verbose>3) printf(" writing Quadrilaterals\n"); 547 bamgmesh->NumQuadrilaterals=NbOfQuad; 548 xfree((void**)&bamgmesh->Quadrilaterals); 549 if (NbOfQuad){ 550 bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double)); 551 for (i=0;i<nbt;i++){ 552 Triangle &t =triangles[i]; 553 Triangle* ta; 554 Vertex *v0,*v1,*v2,*v3; 555 if (reft[i]<0) continue; 556 if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { 557 bamgmesh->Quadrilaterals[i*5+0]=Number(v0)+1; //back to M indexing 558 bamgmesh->Quadrilaterals[i*5+1]=Number(v1)+1; //back to M indexing 559 bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing 560 bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing 561 bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref; 562 } 563 } 564 } 565 else{ 566 bamgmesh->Quadrilaterals=NULL; 567 } 568 569 //SubDomains 570 if(verbose>3) printf(" writing SubDomains\n"); 571 bamgmesh->NumSubDomains=NbSubDomains; 572 xfree((void**)&bamgmesh->SubDomains); 573 if (NbSubDomains){ 574 bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double)); 575 for (i=0;i<NbSubDomains;i++){ 576 bamgmesh->SubDomains[i*4+0]=3; 577 bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)]; 578 bamgmesh->SubDomains[i*4+2]=1; 579 bamgmesh->SubDomains[i*4+3]=subdomains[i].ref; 580 } 581 } 582 else{ 583 bamgmesh->SubDomains=NULL; 584 } 585 586 //SubDomainsFromGeom 587 if(verbose>3) printf(" writing SubDomainsFromGeom\n"); 588 bamgmesh->NumSubDomainsFromGeom=Gh.NbSubDomains; 589 xfree((void**)&bamgmesh->SubDomainsFromGeom); 590 if (Gh.NbSubDomains){ 591 bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double)); 592 for (i=0;i<Gh.NbSubDomains;i++){ 593 bamgmesh->SubDomainsFromGeom[i*4+0]=2; 594 bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing 595 bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens; 596 bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref; 597 } 598 } 599 else{ 600 bamgmesh->SubDomainsFromGeom=NULL; 601 } 602 603 //VerticesOnGeomVertex 604 if(verbose>3) printf(" writing VerticesOnGeometricVertex\n"); 605 bamgmesh->NumVerticesOnGeometricVertex=NbVerticesOnGeomVertex; 606 xfree((void**)&bamgmesh->VerticesOnGeometricVertex); 607 if (NbVerticesOnGeomVertex){ 608 bamgmesh->VerticesOnGeometricVertex=(double*)xmalloc(2*NbVerticesOnGeomVertex*sizeof(double)); 609 for (i=0;i<NbVerticesOnGeomVertex;i++){ 610 VertexOnGeom &v=VerticesOnGeomVertex[i]; 611 if (!v.OnGeomVertex()){ 612 throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not")); 613 } 614 bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing 615 bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing 616 } 617 } 618 else{ 619 bamgmesh->VerticesOnGeometricVertex=NULL; 620 } 621 622 //VertexOnGeometricEdge 623 if(verbose>3) printf(" writing VerticesOnGeometricEdge\n"); 624 bamgmesh->NumVerticesOnGeometricEdge=NbVerticesOnGeomEdge; 625 xfree((void**)&bamgmesh->VerticesOnGeometricEdge); 626 if (NbVerticesOnGeomEdge){ 627 bamgmesh->VerticesOnGeometricEdge=(double*)xmalloc(3*NbVerticesOnGeomEdge*sizeof(double)); 628 for (i=0;i<NbVerticesOnGeomEdge;i++){ 629 const VertexOnGeom &v=VerticesOnGeomEdge[i]; 630 if (!v.OnGeomEdge()){ 631 throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricEdge is actually not")); 632 } 633 bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing 634 bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing 635 bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; 636 } 637 } 638 else{ 639 bamgmesh->VerticesOnGeometricEdge=NULL; 640 } 641 642 //EdgesOnGeometricEdge 643 if(verbose>3) printf(" writing EdgesOnGeometricEdge\n"); 644 k=0; 645 for (i=0;i<nbe;i++){ 646 if (edges[i].on) k=k+1; 647 } 648 bamgmesh->NumEdgesOnGeometricEdge=k; 649 xfree((void**)&bamgmesh->EdgesOnGeometricEdge); 650 if (k){ 651 bamgmesh->EdgesOnGeometricEdge=(double*)xmalloc(2*(int)k*sizeof(double)); 652 int count=0; 653 for (i=0;i<nbe;i++){ 654 if (edges[i].on){ 655 bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing 656 bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].on)+1; //back to Matlab indexing 657 count=count+1; 658 } 659 } 660 } 661 else{ 662 bamgmesh->EdgesOnGeometricEdge=NULL; 663 } 664 } 665 /*}}}1*/ 666 /*FUNCTION Triangles::ReadMetric{{{1*/ 667 void Triangles::ReadMetric(BamgOpts* bamgopts,const Real8 hmin1=1.0e-30,const Real8 hmax1=1.0e30,const Real8 coef=1) { 668 int i,j; 669 670 if(bamgopts->verbose>3) printf(" processing metric\n"); 671 672 Real8 hmin = Max(hmin1,MinimalHmin()); 673 Real8 hmax = Min(hmax1,MaximalHmax()); 674 675 //for now we only use j==3 676 j=3; 677 678 for (i=0;i<nbv;i++){ 679 Real8 h; 680 if (j == 1){ 681 h=bamgopts->metric[i]; 682 vertices[i].m=Metric(Max(hmin,Min(hmax, h*coef))); 683 } 684 else if (j==3){ 685 Real8 a,b,c; 686 a=bamgopts->metric[i*3+0]; 687 b=bamgopts->metric[i*3+1]; 688 c=bamgopts->metric[i*3+2]; 689 MetricAnIso M(a,b,c); 690 MatVVP2x2 Vp(M/coef); 691 692 Vp.Maxh(hmin); 693 Vp.Minh(hmax); 694 vertices[i].m = Vp; 695 } 696 } 697 } 698 /*}}}1*/ 699 /*FUNCTION Triangles::WriteMetric TO UPDATE{{{1*/ 700 void Triangles::WriteMetric(ostream & f,int iso) { 701 if (iso) 702 { 703 f << nbv <<" " << 1 << endl ; 704 for (Int4 iv=0;iv<nbv;iv++) 705 { 706 MatVVP2x2 V=vertices[iv].m; 707 f << V.hmin() << endl; 708 } 709 } 710 else 711 { 712 f << nbv <<" " << 3 << endl ; 713 for (Int4 iv=0;iv<nbv;iv++) 714 f << vertices[iv].m.a11 << " " 715 << vertices[iv].m.a21 << " " 716 << vertices[iv].m.a22 << endl; 717 } 664 718 } 665 719 /*}}}1*/ … … 4918 4972 } 4919 4973 /*}}}1*/ 4974 /*FUNCTION Triangles::IntersectGeomMetric{{{1*/ 4975 void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0){ 4976 long int verbosity=0; 4977 4978 if(verbosity>1) 4979 cout << " -- IntersectGeomMetric geometric err=" << err << (iso ? " iso " : " aniso " ) << endl; 4980 Real8 ss[2]={0.00001,0.99999}; 4981 Real8 errC = 2*sqrt(2*err); 4982 Real8 hmax = Gh.MaximalHmax(); 4983 Real8 hmin = Gh.MinimalHmin(); 4984 Real8 maxaniso = 1e6; 4985 assert(hmax>0); 4986 SetVertexFieldOn(); 4987 if (errC > 1) errC = 1; 4988 for (Int4 i=0;i<nbe;i++) 4989 for (int j=0;j<2;j++) 4990 { 4991 4992 Vertex V; 4993 VertexOnGeom GV; 4994 // cerr << Number(edges[i]) << " " << ss[j] << endl; 4995 Gh.ProjectOnCurve(edges[i],ss[j],V,GV); 4996 { 4997 GeometricalEdge * eg = GV; 4998 Real8 s = GV; 4999 R2 tg; 5000 // cerr << i << " " << j << " " << Number(V) << " on = " 5001 // << Gh.Number(eg) << " at s = " << s << " " << endl; 5002 Real8 R1= eg->R1tg(s,tg); 5003 // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x " 5004 // << V.r << errC/ Max(R1,1e-20) << " hold=" <<V.m(tg) << " " << endl; 5005 Real8 ht = hmax; 5006 if (R1>1.0e-20) 5007 { // err relative to the length of the edge 5008 ht = Min(Max(errC/R1,hmin),hmax); 5009 } 5010 Real8 hn = iso? ht : Min(hmax,ht*maxaniso); 5011 //cerr << ht << " " << hn << "m=" << edges[i][j].m << endl; 5012 assert(ht>0 && hn>0); 5013 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg); 5014 //cerr << " : " ; 5015 Metric MVp(Vp); 5016 // cerr << " : " << MVp << endl; 5017 edges[i][j].m.IntersectWith(MVp); 5018 //cerr << " . " << endl; 5019 } 5020 5021 } 5022 // the problem is for the vertex on vertex 5023 } 5024 /*}}}1*/ 5025 /*FUNCTION Triangles::BoundAnisotropy{{{1*/ 5026 void Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso) { 5027 long int verbosity=0; 5028 5029 double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100)); 5030 if (verbosity > 1) 5031 cout << " -- BoundAnisotropy by " << anisomax << endl; 5032 Real8 h1=1.e30,h2=1e-30,rx=0; 5033 Real8 coef = 1./(anisomax*anisomax); 5034 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 5035 for (Int4 i=0;i<nbv;i++) 5036 { 5037 5038 MatVVP2x2 Vp(vertices[i]); 5039 double lmax=Vp.lmax(); 5040 h1=Min(h1,Vp.lmin()); 5041 h2=Max(h2,Vp.lmax()); 5042 rx = Max(rx,Vp.Aniso2()); 5043 5044 Vp *= Min(lminaniso,lmax)/lmax; 5045 5046 Vp.BoundAniso2(coef); 5047 5048 hn1=Min(hn1,Vp.lmin()); 5049 hn2=Max(hn2,Vp.lmax()); 5050 rnx = Max(rnx,Vp.Aniso2()); 5051 5052 5053 vertices[i].m = Vp; 5054 5055 } 5056 5057 if (verbosity>2) 5058 { 5059 cout << " input : Hmin = " << sqrt(1/h2) << " Hmax = " << sqrt(1/h1) 5060 << " factor of anisotropy max = " << sqrt(rx) << endl; 5061 cout << " output: Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 5062 << " factor of anisotropy max = " << sqrt(rnx) << endl; 5063 } 5064 } 5065 /*}}}1*/ 5066 /*FUNCTION Triangles::IntersectConsMetric{{{1*/ 5067 void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols, 5068 const Real8 hmin1,const Real8 hmax1,const Real8 coef, 5069 const Real8 anisomax ,const Real8 CutOff,const int NbJacobi, 5070 const int DoNormalisation,const double power,const int choice) 5071 { // the array of solution s is store 5072 // sol0,sol1,...,soln on vertex 0 5073 // sol0,sol1,...,soln on vertex 1 5074 // etc. 5075 // choise = 0 => H is computed with green formule 5076 // otherwise => H is computed from P2 on 4T 5077 const int dim = 2; 5078 5079 long int verbosity=0; 5080 5081 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 5082 5083 // computation of the nb of field 5084 Int4 ntmp = 0; 5085 if (typsols) 5086 { 5087 for (Int4 i=0;i<nbsol;i++) 5088 ntmp += sizeoftype[typsols[i]]; 5089 } 5090 else 5091 ntmp = nbsol; 5092 5093 // n is the total number of fields 5094 5095 const Int4 n = ntmp; 5096 5097 Int4 i,k,iA,iB,iC,iv; 5098 R2 O(0,0); 5099 int RelativeMetric = CutOff>1e-30; 5100 Real8 hmin = Max(hmin1,MinimalHmin()); 5101 Real8 hmax = Min(hmax1,MaximalHmax()); 5102 Real8 coef2 = 1/(coef*coef); 5103 5104 if(verbosity>1) 5105 { 5106 cout << " -- Construction of Metric: Nb of field. " << n << " nbt = " 5107 << nbt << " nbv= " << nbv 5108 << " coef = " << coef << endl 5109 << " hmin = " << hmin << " hmax=" << hmax 5110 << " anisomax = " << anisomax << " Nb Jacobi " << NbJacobi << " Power = " << power ; 5111 if (RelativeMetric) 5112 cout << " RelativeErr with CutOff= " << CutOff << endl; 5113 else 5114 cout << " Absolute Err" <<endl; 5115 } 5116 double *ss=(double*)s;//, *ssiii = ss; 5117 5118 double sA,sB,sC; 5119 5120 Real8 *detT = new Real8[nbt]; 5121 Real8 *Mmass= new Real8[nbv]; 5122 Real8 *Mmassxx= new Real8[nbv]; 5123 Real8 *dxdx= new Real8[nbv]; 5124 Real8 *dxdy= new Real8[nbv]; 5125 Real8 *dydy= new Real8[nbv]; 5126 Real8 *workT= new Real8[nbt]; 5127 Real8 *workV= new Real8[nbv]; 5128 int *OnBoundary = new int[nbv]; 5129 for (iv=0;iv<nbv;iv++) 5130 { 5131 Mmass[iv]=0; 5132 OnBoundary[iv]=0; 5133 Mmassxx[iv]=0; 5134 } 5135 5136 for (i=0;i<nbt;i++) 5137 if(triangles[i].link) // the real triangles 5138 { 5139 const Triangle &t=triangles[i]; 5140 // coor of 3 vertices 5141 R2 A=t[0]; 5142 R2 B=t[1]; 5143 R2 C=t[2]; 5144 5145 5146 // number of the 3 vertices 5147 iA = Number(t[0]); 5148 iB = Number(t[1]); 5149 iC = Number(t[2]); 5150 5151 Real8 dett = bamg::Area2(A,B,C); 5152 detT[i]=dett; 5153 dett /= 6; 5154 5155 // construction of on boundary 5156 int nbb =0; 5157 for(int j=0;j<3;j++) 5158 { 5159 Triangle *ta=t.Adj(j); 5160 if ( ! ta || !ta->link) // no adj triangle => edge on boundary 5161 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1, 5162 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1, 5163 nbb++; 5164 } 5165 5166 workT[i] = nbb; 5167 Mmass[iA] += dett; 5168 Mmass[iB] += dett; 5169 Mmass[iC] += dett; 5170 5171 if((nbb==0)|| !choice) 5172 { 5173 Mmassxx[iA] += dett; 5174 Mmassxx[iB] += dett; 5175 Mmassxx[iC] += dett; 5176 } 5177 } 5178 else 5179 workT[i]=-1; 5180 5181 for (Int4 nusol=0;nusol<nbsol;nusol++) 5182 { //for all Solution 5183 5184 Real8 smin=ss[0],smax=ss[0]; 5185 5186 Real8 h1=1.e30,h2=1e-30,rx=0; 5187 Real8 coef = 1./(anisomax*anisomax); 5188 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 5189 int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 5190 if (nbfield == 1) 5191 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) 5192 { 5193 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 5194 smin=Min(smin,ss[k]); 5195 smax=Max(smax,ss[k]); 5196 } 5197 else 5198 { 5199 // cas vectoriel 5200 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) 5201 { 5202 double v=0; 5203 for (int i=0;i<nbfield;i++) 5204 v += ss[k+i]*ss[k+i]; 5205 v = sqrt(v); 5206 smin=Min(smin,v); 5207 smax=Max(smax,v); 5208 } 5209 } 5210 Real8 sdelta = smax-smin; 5211 Real8 absmax=Max(Abs(smin),Abs(smax)); 5212 Real8 cnorm = DoNormalisation ? coef2/sdelta : coef2; 5213 5214 if(verbosity>2) 5215 cout << " Solution " << nusol << " Min = " << smin << " Max = " 5216 << smax << " Delta =" << sdelta << " cnorm = " << cnorm << " Nb of fields =" << nbfield << endl; 5217 5218 5219 if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) 5220 { 5221 if (verbosity>2) 5222 cout << " Solution " << nusol << " is constant. We skip. " 5223 << " Min = " << smin << " Max = " << smax << endl; 5224 continue; 5225 } 5226 5227 double *sf = ss; 5228 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++) 5229 { 5230 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) 5231 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 5232 for (i=0;i<nbt;i++) 5233 if(triangles[i].link) 5234 {// for real all triangles 5235 // coor of 3 vertices 5236 R2 A=triangles[i][0]; 5237 R2 B=triangles[i][1]; 5238 R2 C=triangles[i][2]; 5239 5240 5241 // warning the normal is internal and the 5242 // size is the length of the edge 5243 R2 nAB = Orthogonal(B-A); 5244 R2 nBC = Orthogonal(C-B); 5245 R2 nCA = Orthogonal(A-C); 5246 // remark : nAB + nBC + nCA == 0 5247 5248 // number of the 3 vertices 5249 iA = Number(triangles[i][0]); 5250 iB = Number(triangles[i][1]); 5251 iC = Number(triangles[i][2]); 5252 5253 // for the test of boundary edge 5254 // the 3 adj triangles 5255 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 5256 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 5257 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 5258 5259 // value of the P1 fonction on 3 vertices 5260 sA = ss[iA*n]; 5261 sB = ss[iB*n]; 5262 sC = ss[iC*n]; 5263 5264 R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ; 5265 if(choice) 5266 { 5267 int nbb = 0; 5268 Real8 dd = detT[i]; 5269 Real8 lla,llb,llc,llf; 5270 Real8 taa[3][3],bb[3]; 5271 // construction of the trans of lin system 5272 for (int j=0;j<3;j++) 5273 { 5274 int ie = OppositeEdge[j]; 5275 TriangleAdjacent ta = triangles[i].Adj(ie); 5276 Triangle *tt = ta; 5277 if (tt && tt->link) 5278 { 5279 Vertex &v = *ta.OppositeVertex(); 5280 R2 V = v; 5281 Int4 iV = Number(v); 5282 Real8 lA = bamg::Area2(V,B,C)/dd; 5283 Real8 lB = bamg::Area2(A,V,C)/dd; 5284 Real8 lC = bamg::Area2(A,B,V)/dd; 5285 taa[0][j] = lB*lC; 5286 taa[1][j] = lC*lA; 5287 taa[2][j] = lA*lB; 5288 //Real8 xx = V.x-V.y; 5289 //Real8 yy = V.x + V.y; 5290 //cout << " iv " << ss[iV*n] << " == " << (8*xx*xx+yy*yy) 5291 // << " l = " << lA << " " << lB << " " << lC 5292 // << " = " << lA+lB+lC << " " << V << " == " << A*lA+B*lB+C*lC << endl; 5293 5294 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; 5295 5296 bb[j] = ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ; 5297 } 5298 else 5299 { 5300 nbb++; 5301 taa[0][j]=0; 5302 taa[1][j]=0; 5303 taa[2][j]=0; 5304 taa[j][j]=1; 5305 bb[j]=0; 5306 } 5307 } 5308 5309 // resolution of 3x3 lineaire system transpose 5310 Real8 det33 = det3x3(taa[0],taa[1],taa[2]); 5311 Real8 cBC = det3x3(bb,taa[1],taa[2]); 5312 Real8 cCA = det3x3(taa[0],bb,taa[2]); 5313 Real8 cAB = det3x3(taa[0],taa[1],bb); 5314 5315 assert(det33); 5316 // det33=1; 5317 // verif 5318 // cout << " " << (taa[0][0]*cBC + taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ; 5319 // cout << " " << (taa[0][1]*cBC + taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1]; 5320 // cout << " " << (taa[0][2]*cBC + taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2] 5321 // << " -- " ; 5322 //cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB + llb*llc* cBC + llc*lla*cCA)/det33 5323 // << " == " << llf << endl; 5324 // computation of the gradient in the element 5325 5326 // H( li*lj) = grad li grad lj + grad lj grad lj 5327 // grad li = njk / detT ; with i j k =(A,B,C) 5328 Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x); 5329 Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y); 5330 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 5331 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); 5332 Real8 coef = 1.0/(3*dd*det33); 5333 Real8 coef2 = 2*coef; 5334 // cout << " H = " << Hxx << " " << Hyy << " " << Hxy/2 << " coef2 = " << coef2 << endl; 5335 Hxx *= coef2; 5336 Hyy *= coef2; 5337 Hxy *= coef2; 5338 //cout << i << " H = " << 3*Hxx/dd << " " << 3*Hyy/dd << " " << 3*Hxy/(dd*2) << " nbb = " << nbb << endl; 5339 if(nbb==0) 5340 { 5341 dxdx[iA] += Hxx; 5342 dydy[iA] += Hyy; 5343 dxdy[iA] += Hxy; 5344 5345 dxdx[iB] += Hxx; 5346 dydy[iB] += Hyy; 5347 dxdy[iB] += Hxy; 5348 5349 dxdx[iC] += Hxx; 5350 dydy[iC] += Hyy; 5351 dxdy[iC] += Hxy; 5352 } 5353 5354 } 5355 else 5356 { 5357 5358 // if edge on boundary no contribution => normal = 0 5359 if ( ! tBC || ! tBC->link ) nBC = O; 5360 if ( ! tCA || ! tCA->link ) nCA = O; 5361 if ( ! tAB || ! tAB->link ) nAB = O; 5362 5363 // remark we forgot a 1/2 because 5364 // $\\int_{edge} w_i = 1/2 $ if $i$ is in edge 5365 // 0 if not 5366 // if we don't take the boundary 5367 // dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; 5368 5369 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; 5370 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; 5371 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; 5372 5373 // warning optimization (1) the divide by 2 is done on the metrix construction 5374 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; 5375 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; 5376 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 5377 5378 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; 5379 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; 5380 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; 5381 } 5382 5383 } // for real all triangles 5384 Int4 kk=0; 5385 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ) 5386 if(Mmassxx[iv]>0) 5387 { 5388 dxdx[iv] /= 2*Mmassxx[iv]; 5389 // warning optimization (1) on term dxdy[iv]*ci/2 5390 dxdy[iv] /= 4*Mmassxx[iv]; 5391 dydy[iv] /= 2*Mmassxx[iv]; 5392 // Compute the matrix with abs(eigen value) 5393 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 5394 MatVVP2x2 Vp(M); 5395 //cout <<iv << " M = " << M << " aniso= " << Vp.Aniso() ; 5396 Vp.Abs(); 5397 M = Vp; 5398 dxdx[iv] = M.a11; 5399 dxdy[iv] = M.a21; 5400 dydy[iv] = M.a22; 5401 // cout << " (abs) iv M = " << M << " aniso= " << Vp.Aniso() <<endl; 5402 } 5403 else kk++; 5404 5405 5406 // correction of second derivate 5407 // by a laplacien 5408 5409 Real8 *d2[3] = { dxdx, dxdy, dydy}; 5410 Real8 *dd; 5411 for (int xy = 0;xy<3;xy++) 5412 { 5413 dd = d2[xy]; 5414 // do leat 2 iteration for boundary problem 5415 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++) 5416 { 5417 for (i=0;i<nbt;i++) 5418 if(triangles[i].link) // the real triangles 5419 { 5420 // number of the 3 vertices 5421 iA = Number(triangles[i][0]); 5422 iB = Number(triangles[i][1]); 5423 iC = Number(triangles[i][2]); 5424 Real8 cc=3; 5425 if(ijacobi==0) 5426 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 5427 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 5428 } 5429 for (iv=0;iv<nbv;iv++) 5430 workV[iv]=0; 5431 5432 for (i=0;i<nbt;i++) 5433 if(triangles[i].link) // the real triangles 5434 { 5435 // number of the 3 vertices 5436 iA = Number(triangles[i][0]); 5437 iB = Number(triangles[i][1]); 5438 iC = Number(triangles[i][2]); 5439 Real8 cc = workT[i]*detT[i]; 5440 workV[iA] += cc; 5441 workV[iB] += cc; 5442 workV[iC] += cc; 5443 } 5444 5445 for (iv=0;iv<nbv;iv++) 5446 if( ijacobi<NbJacobi || OnBoundary[iv]) 5447 dd[iv] = workV[iv]/(Mmass[iv]*6); 5448 5449 5450 } 5451 5452 5453 } 5454 5455 // constuction of the metrix from the Hessian dxdx. dxdy,dydy 5456 5457 Real8 rCutOff=CutOff*absmax;// relative cut off 5458 5459 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ) 5460 { // for all vertices 5461 //{ 5462 //Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 5463 // MatVVP2x2 Vp(M); 5464 //cout << " iv M="<< M << " Vp = " << Vp << " aniso " << Vp.Aniso() << endl; 5465 //} 5466 MetricIso Miso; 5467 // new code to compute ci --- 5468 Real8 ci ; 5469 if (RelativeMetric) 5470 { // compute the norm of the solution 5471 double xx =0,*sfk=sf+k; 5472 for (int ifield=0;ifield<nbfield;ifield++,sfk++) 5473 xx += *sfk* *sfk; 5474 xx=sqrt(xx); 5475 ci = coef2/Max(xx,rCutOff); 5476 } 5477 else ci = cnorm; 5478 5479 // old 5480 // Real8 ci = RelativeMetric ? coef2/(Max(Abs(ss[k]),rCutOff)) : cnorm ; 5481 // modif F Hecht 101099 5482 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 5483 MatVVP2x2 Vp(Miv); 5484 5485 Vp.Abs(); 5486 if(power!=1.0) 5487 Vp.pow(power); 5488 5489 5490 5491 h1=Min(h1,Vp.lmin()); 5492 h2=Max(h2,Vp.lmax()); 5493 5494 Vp.Maxh(hmin); 5495 Vp.Minh(hmax); 5496 5497 rx = Max(rx,Vp.Aniso2()); 5498 5499 Vp.BoundAniso2(coef); 5500 5501 hn1=Min(hn1,Vp.lmin()); 5502 hn2=Max(hn2,Vp.lmax()); 5503 rnx = Max(rnx,Vp.Aniso2()); 5504 5505 Metric MVp(Vp); 5506 vertices[iv].m.IntersectWith(MVp); 5507 }// for all vertices 5508 if (verbosity>2) 5509 { 5510 cout << " Field " << nufield << " of solution " << nusol << endl; 5511 cout << " before bounding : Hmin = " << sqrt(1/h2) << " Hmax = " 5512 << sqrt(1/h1) << " factor of anisotropy max = " << sqrt(rx) << endl; 5513 cout << " after bounding : Hmin = " << sqrt(1/hn2) << " Hmax = " 5514 << sqrt(1/hn1) << " factor of anisotropy max = " << sqrt(rnx) << endl; 5515 } 5516 } // end of for all field 5517 }// end for all solution 5518 5519 delete [] detT; 5520 delete [] Mmass; 5521 delete [] dxdx; 5522 delete [] dxdy; 5523 delete [] dydy; 5524 delete [] workT; 5525 delete [] workV; 5526 delete [] Mmassxx; 5527 delete [] OnBoundary; 5528 5529 } 5530 /*}}}1*/ 5531 /*FUNCTION Triangles::MaxSubDivision{{{1*/ 5532 void Triangles::MaxSubDivision(Real8 maxsubdiv) { 5533 long int verbosity=0; 5534 5535 const Real8 maxsubdiv2 = maxsubdiv*maxsubdiv; 5536 if(verbosity>1) 5537 cout << " -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv << endl ; 5538 // for all the edges 5539 // if the len of the edge is to long 5540 Int4 it,nbchange=0; 5541 Real8 lmax=0; 5542 for (it=0;it<nbt;it++) 5543 { 5544 Triangle &t=triangles[it]; 5545 for (int j=0;j<3;j++) 5546 { 5547 Triangle &tt = *t.TriangleAdj(j); 5548 if ( ! &tt || it < Number(tt) && ( tt.link || t.link)) 5549 { 5550 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 5551 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 5552 R2 AB= (R2) v1-(R2) v0; 5553 Metric M = v0; 5554 Real8 l = M(AB,AB); 5555 lmax = Max(lmax,l); 5556 if(l> maxsubdiv2) 5557 { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 5558 Real8 lc = M(AC,AC); 5559 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; 5560 D2xD2 Rt1(Rt.inv()); 5561 D2xD2 D(maxsubdiv2,0,0,lc); 5562 D2xD2 MM = Rt1*D*Rt1.t(); 5563 v0.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y); 5564 nbchange++; 5565 } 5566 M = v1; 5567 l = M(AB,AB); 5568 lmax = Max(lmax,l); 5569 if(l> maxsubdiv2) 5570 { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 5571 Real8 lc = M(AC,AC); 5572 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; 5573 D2xD2 Rt1(Rt.inv()); 5574 D2xD2 D(maxsubdiv2,0,0,lc); 5575 D2xD2 MM = Rt1*D*Rt1.t(); 5576 v1.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y); 5577 nbchange++; 5578 } 5579 5580 5581 } 5582 } 5583 } 5584 if(verbosity>3) 5585 cout << " Nb of metric change = " << nbchange 5586 << " Max of the subdivision of a edges before change = " << sqrt(lmax) << endl; 5587 5588 } 5589 /*}}}1*/ 5590 /*FUNCTION Triangles::SmoothMetric{{{1*/ 5591 void Triangles::SmoothMetric(Real8 raisonmax) { 5592 long int verbosity=0; 5593 5594 if(raisonmax<1.1) return; 5595 if(verbosity > 1) 5596 cout << " -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl; 5597 ReMakeTriangleContainingTheVertex(); 5598 Int4 i,j,kch,kk,ip; 5599 Int4 *first_np_or_next_t0 = new Int4[nbv]; 5600 Int4 *first_np_or_next_t1 = new Int4[nbv]; 5601 Int4 Head0 =0,Head1=-1; 5602 Real8 logseuil= log(raisonmax); 5603 5604 for(i=0;i<nbv-1;i++) 5605 first_np_or_next_t0[i]=i+1; 5606 first_np_or_next_t0[nbv-1]=-1;// end; 5607 for(i=0;i<nbv;i++) 5608 first_np_or_next_t1[i]=-1; 5609 kk=0; 5610 while (Head0>=0&& kk++<100) { 5611 kch=0; 5612 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) 5613 { // pour tous les triangles autour du sommet s 5614 // cout << kk << " i = " << i << " " << ip << endl; 5615 register Triangle * t= vertices[i].t; 5616 assert(t); 5617 Vertex & vi = vertices[i]; 5618 TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]); 5619 Vertex *pvj0 = ta.EdgeVertex(0); 5620 while (1) { 5621 // cout << i << " " << Number(ta.EdgeVertex(0)) << " " 5622 // << Number(ta.EdgeVertex(1)) << " ---> " ; 5623 ta=Previous(Adj(ta)); 5624 // cout << Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl; 5625 assert(vertices+i == ta.EdgeVertex(1)); 5626 Vertex & vj = *(ta.EdgeVertex(0)); 5627 if ( &vj ) { 5628 j= &vj-vertices; 5629 assert(j>=0 && j < nbv); 5630 R2 Aij = (R2) vj - (R2) vi; 5631 Real8 ll = Norme2(Aij); 5632 if (0) { 5633 Real8 hi = ll/vi.m(Aij); 5634 Real8 hj = ll/vj.m(Aij); 5635 if(hi < hj) 5636 { 5637 Real8 dh=(hj-hi)/ll; 5638 //cout << " dh = " << dh << endl; 5639 if (dh>logseuil) { 5640 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi)); 5641 if(first_np_or_next_t1[j]<0) 5642 kch++,first_np_or_next_t1[j]=Head1,Head1=j; 5643 } 5644 } 5645 } 5646 else 5647 { 5648 Real8 li = vi.m(Aij); 5649 //Real8 lj = vj.m(Aij); 5650 // if ( i == 2 || j == 2) 5651 // cout << " inter " << i << " " << j << " " << ((1 +logseuil*li)) << endl; 5652 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) 5653 //if( vj.m.IntersectWith(vi.m*(lj/li/(1 +logseuil*lj))) ) 5654 if(first_np_or_next_t1[j]<0) // if the metrix change 5655 kch++,first_np_or_next_t1[j]=Head1,Head1=j; 5656 } 5657 } 5658 if ( &vj == pvj0 ) break; 5659 } 5660 } 5661 Head0 = Head1; 5662 Head1 = -1; 5663 Exchange(first_np_or_next_t0,first_np_or_next_t1); 5664 if(verbosity>5) 5665 cout << " Iteration " << kk << " Nb de vertices with change " << kch << endl; 5666 } 5667 if(verbosity>2 && verbosity < 5) 5668 cout << " Nb of Loop " << kch << endl; 5669 delete [] first_np_or_next_t0; 5670 delete [] first_np_or_next_t1; 5671 } 5672 /*}}}1*/ 5673 /*FUNCTION Triangles::NearestVertex{{{1*/ 5674 Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j) { 5675 return quadtree->NearestVertex(i,j); 5676 } 5677 /*}}}1*/ 5678 /*FUNCTION Triangles::PreInit{{{1*/ 5679 void Triangles::PreInit(Int4 inbvx,char *fname) { 5680 long int verbosity=0; 5681 5682 srand(19999999); 5683 OnDisk =0; 5684 NbRef=0; 5685 // allocGeometry=0; 5686 identity=0; 5687 NbOfTriangleSearchFind =0; 5688 NbOfSwapTriangle =0; 5689 nbiv=0; 5690 nbv=0; 5691 nbvx=inbvx; 5692 nbt=0; 5693 NbOfQuad = 0; 5694 nbtx=2*inbvx-2; 5695 NbSubDomains=0; 5696 NbVertexOnBThVertex=0; 5697 NbVertexOnBThEdge=0; 5698 VertexOnBThVertex=0; 5699 VertexOnBThEdge=0; 5700 5701 NbCrackedVertices=0; 5702 NbCrackedEdges =0; 5703 CrackedEdges =0; 5704 nbe = 0; 5705 name = fname ; 5706 5707 if (inbvx) { 5708 vertices=new Vertex[nbvx]; 5709 assert(vertices); 5710 ordre=new (Vertex* [nbvx]); 5711 assert(ordre); 5712 triangles=new Triangle[nbtx]; 5713 assert(triangles);} 5714 else { 5715 vertices=0; 5716 ordre=0; 5717 triangles=0; 5718 nbtx=0; 5719 } 5720 if ( name || inbvx) 5721 { 5722 time_t timer =time(0); 5723 char buf[70]; 5724 strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer)); 5725 counter++; 5726 char countbuf[30]; 5727 sprintf(countbuf,"%d",counter); 5728 int lg =0 ; 5729 if (&BTh != this && BTh.name) 5730 lg = strlen(BTh.name)+4; 5731 identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2 + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)]; 5732 identity[0]=0; 5733 if (lg) 5734 strcat(strcat(strcat(identity,"B="),BTh.name),", "); 5735 5736 if (Gh.name) 5737 strcat(strcat(identity,"G="),Gh.name); 5738 strcat(strcat(identity,";"),countbuf); 5739 strcat(identity,buf); 5740 // cout << "New MAILLAGE "<< identity << endl; 5741 } 5742 5743 quadtree=0; 5744 // edgescomponante=0; 5745 edges=0; 5746 VerticesOnGeomVertex=0; 5747 VerticesOnGeomEdge=0; 5748 NbVerticesOnGeomVertex=0; 5749 NbVerticesOnGeomEdge=0; 5750 // nbMaxIntersectionTriangles=0; 5751 // lIntTria; 5752 subdomains=0; 5753 NbSubDomains=0; 5754 // Meshbegin = vertices; 5755 // Meshend = vertices + nbvx; 5756 if (verbosity>98) 5757 cout << "Triangles::PreInit() " << nbvx << " " << nbtx 5758 << " " << vertices 5759 << " " << ordre << " " << triangles <<endl; 5760 } 5761 /*}}}1*/ 4920 5762 4921 5763 } // end of namespace bamg
Note:
See TracChangeset
for help on using the changeset viewer.
