![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2932
- Timestamp:
- 01/29/10 09:49:13 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 5 edited
-
Bamgx.cpp (modified) (1 diff)
-
Metric.h (modified) (2 diffs)
-
R2.h (modified) (2 diffs)
-
objects/MetricAnIso.cpp (modified) (5 diffs)
-
objects/Triangles.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Bamgx.cpp
r2928 r2932 149 149 150 150 BTh.IntersectGeomMetric(errg,iso); 151 if (verbosity>1) printf(" Smoothing Metric..."); 151 152 if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation); 153 if (verbosity>1) printf(" done\n"); 152 154 BTh.MaxSubDivision(maxsubdiv); 153 155 BTh.BoundAnisotropy(anisomax,hminaniso); -
issm/trunk/src/c/Bamgx/Metric.h
r2928 r2932 64 64 Real8 b,const MetricAnIso mb); 65 65 int IntersectWith(const MetricAnIso M2); 66 int IntersectWith_new(const MetricAnIso M2);67 66 MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return MetricAnIso(a11*c2,a21*c2,a22*c2);} 68 67 MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return MetricAnIso(a11*c2,a21*c2,a22*c2);} … … 123 122 } 124 123 125 void SimultaneousMatrixReduction(const MetricAnIso M1,const MetricAnIso M2, MetricAnIso* M); 126 void SimultaneousMatrixReduction_bamg( MetricAnIso M1, MetricAnIso M2,double & l1,double & l2, D2xD2 & V); 127 void eigen2(double* M,double lambda[2],double vp[2][2]);; 124 void SimultaneousMatrixReduction( MetricAnIso M1, MetricAnIso M2,double & l1,double & l2, D2xD2 & V); 128 125 129 126 inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) { -
issm/trunk/src/c/Bamgx/R2.h
r2924 r2932 8 8 R x,y; 9 9 void Echo(){ 10 printf("Member of P2:\n"); 10 11 printf(" x: %g\n",x); 11 12 printf(" y: %g\n",y); … … 35 36 public: 36 37 P2<R,RR> x,y; 38 void Echo(){ 39 printf("Member of P2xP2:\n"); 40 printf(" x.x: %g x.y: %g\n",x.x,x.y); 41 printf(" y.x: %g y.x: %g\n",y.x,y.y); 42 } 37 43 P2xP2 (): x(),y() {} 38 44 P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {} -
issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
r2928 r2932 69 69 } 70 70 /*}}}*/ 71 /*FUNCTION MetricAnIso::IntersectWith_new{{{1*/72 int MetricAnIso::IntersectWith_new(MetricAnIso M2) {73 /*Get a new metric from an existing metric (M1=this)74 * and a new metric given in input M2 using a75 * Simultaneous Matrix Reduction76 * */77 78 MetricAnIso M;79 80 //Compute M intersection of M1 and M281 SimultaneousMatrixReduction(*this,M2,&M);82 83 //check wether something has been done84 if(M.a11!=a11 || M.a21!=a21 || M.a22!=a22){85 86 //update this with M87 a11=M.a11;88 a21=M.a21;89 a22=M.a22;90 91 //return 1 -> something has been done92 return 1;93 }94 else{95 96 //nothing has changed, return 097 return 0;98 }99 }100 /*}}}1*/101 71 /*FUNCTION MetricAnIso::IntersectWith{{{1*/ 102 72 int MetricAnIso::IntersectWith(const MetricAnIso M2) { … … 111 81 double l1,l2; 112 82 113 SimultaneousMatrixReduction _bamg(*this,M2,l1,l2,M);83 SimultaneousMatrixReduction(*this,M2,l1,l2,M); 114 84 115 85 R2 v1(M.x.x,M.y.x); … … 127 97 a11=Mi.x.x; 128 98 a21=0.5*(Mi.x.y+Mi.y.x); 129 a22=Mi.y.y; } 130 return r; 99 a22=Mi.y.y; 100 } 101 return r; 131 102 } 132 103 /*}}}1*/ 133 104 134 105 /*Intermediary*/ 135 /*FUNCTION eigen2{{{1*/136 void eigen2(double* M,double* lambda,double vp[2][2]) {137 /*m=[a b]138 * [c d] */139 140 double normM,detM,trM;141 double norm;142 const double a=M[2*0+0];143 const double b=M[2*0+1];144 const double c=M[2*1+0];145 const double d=M[2*1+1];146 147 //compute norm(M)148 normM=pow(pow(a,2.0)+pow(b,2.0)+pow(c,2.0)+pow(d,2.0),0.5);149 150 //if normM<10-30, M=zeros(2,2)151 if (normM<1.0e-30){152 lambda[0]=0.0;153 lambda[1]=0.0;154 vp[0][0]=vp[0][1]=vp[1][0]=vp[1][1]=0.0;155 return;156 }157 158 //check that matrix is already diagonal if det(M)=0159 if(Abs(b)<(1.0e-12*normM) & Abs(c)<(1.0e-12*normM)){160 lambda[0]=a;161 lambda[1]=d;162 vp[0][0]=1;163 vp[0][1]=0;164 vp[1][0]=0;165 vp[1][1]=1;166 return;167 }168 169 //compute det(M)170 detM=a*d-c*b;171 172 //not symetrical and still det<10^-30173 if (Abs(detM)<1.0e-30){174 throw ErrorException(__FUNCT__,exprintf("Cannot find eigen values of a non invertible matrix"));175 }176 177 //See http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html178 trM=a+d;179 180 //general case181 lambda[0] = 0.5*(trM + pow( trM*trM - 4.0*detM ,0.5));182 lambda[1] = 0.5*(trM - pow( trM*trM - 4.0*detM ,0.5));183 184 if(Abs(b)>1.0e-30){ //if b!=0185 norm=pow(b*b+pow(lambda[0]-a,2.0),0.5);186 vp[0][0]=b/norm;187 vp[1][0]=(lambda[0]-a)/norm;188 norm=pow(b*b+pow(lambda[1]-a,2.0),0.5);189 vp[0][1]=b/norm;190 vp[1][1]=(lambda[1]-a)/norm;191 }192 else{193 norm=pow(c*c+pow(lambda[0]-d,2.0),0.5);194 vp[0][0]=(lambda[0]-d)/norm;195 vp[1][0]=c/norm;196 norm=pow(c*c+pow(lambda[1]-d,2.0),0.5);197 vp[0][1]=(lambda[1]-d)/norm;198 vp[1][1]=c/norm;199 }200 }201 /*}}}1*/202 106 /*FUNCTION LengthInterpole{{{1*/ 203 107 Real8 LengthInterpole(const MetricAnIso Ma,const MetricAnIso Mb, R2 AB) { … … 267 171 /*}}}1*/ 268 172 /*FUNCTION SimultaneousMatrixReduction{{{1*/ 269 void SimultaneousMatrixReduction(MetricAnIso M1,MetricAnIso M2,MetricAnIso* pM) { 270 271 /*intermediary*/ 272 double a11,a21,a22; 273 double b11,b21,b22; 274 double mu1,mu2; 275 double lambda1,lambda2; 276 double detM1,detP; 277 double N[2][2]; 278 double lambda[2]; 279 double vp[2][2]; 280 double invP[2][2]; 281 MetricAnIso M; 282 283 //get components of both metrics 284 a11=M1.a11; a21=M1.a21; a22=M1.a22; 285 b11=M2.a11; b21=M2.a21; b22=M2.a22; 286 287 //check diagnonal matrices (easy) 288 if (Abs(a21)< 1.0e-20 && Abs(b21)< 1.0e-20){ 289 M.a11=Max(a11,b11); 290 M.a22=Max(a22,b22); 291 M.a21=0; 292 *pM=M; 293 return; 294 } 295 296 //Build N=M1^-1 M2 (Alauzet2003 p16) 297 detM1=a11*a22-a21*a21; 298 if (!detM1){ 299 printf("a11 a21 a22 = %g %g %g\n",a11,a21,a22); 300 throw ErrorException(__FUNCT__,exprintf("One of the metric has a zero eigen value or cannot be inverted... (detP=%g)",detP)); 301 } 302 N[0][0]= 1.0/detM1 * ( a22*b11-a21*b21); 303 N[0][1]= 1.0/detM1 * ( a22*b21-a21*b22); 304 N[1][0]= 1.0/detM1 * (-a21*b11+a11*b21); 305 N[1][1]= 1.0/detM1 * (-a21*b21+a11*b22); 306 307 //now, we must find the eigen values and vectors of N 308 eigen2(&N[0][0],lambda,vp); 309 310 //Now that we have P=vp, we can recompute M = M1 inter M2 311 // 312 // -T [ max(lambda1, mu1) 0 ] -1 313 // M = P [ ] P 314 // [ 0 max(lambda2, mu2)] 315 316 //get det(P) 317 detP=vp[0][0]*vp[1][1]-vp[0][1]*vp[1][0]; 318 if (!detP){ 319 printf("vp= [%g %g ; %g %g]\n",vp[0][0],vp[0][1],vp[1][0],vp[1][1]); 320 throw ErrorException(__FUNCT__,exprintf("One of the metric has a zero eigen value or cannot be inverted... (detP=%g)",detM1)); 321 } 322 invP[0][0]= 1.0/detP * ( vp[1][1]); 323 invP[0][1]= 1.0/detP * (-vp[0][1]); 324 invP[1][0]= 1.0/detP * (-vp[1][0]); 325 invP[1][1]= 1.0/detP * ( vp[0][0]); 326 327 //Compute lambdai and mui using Rayleigh formula: lambdai = ei' M1 ei 328 lambda1=vp[0][0]*(a11*vp[0][0]+a21*vp[1][0]) + vp[1][0]*(a21*vp[0][0]+a22*vp[1][0]); 329 lambda2=vp[0][1]*(a11*vp[0][1]+a21*vp[1][1]) + vp[1][1]*(a21*vp[0][1]+a22*vp[1][1]); 330 mu1 =vp[0][0]*(b11*vp[0][0]+b21*vp[1][0]) + vp[1][0]*(b21*vp[0][0]+b22*vp[1][0]); 331 mu2 =vp[0][1]*(b11*vp[0][1]+b21*vp[1][1]) + vp[1][1]*(b21*vp[0][1]+b22*vp[1][1]); 332 333 //finally compute M 334 M.a11=invP[0][0]*Max(lambda1,mu1)*invP[0][0] + invP[1][0]*Max(lambda2,mu2)*invP[1][0]; 335 M.a21=invP[0][0]*Max(lambda1,mu1)*invP[0][1] + invP[1][0]*Max(lambda2,mu2)*invP[1][1]; 336 M.a22=invP[0][1]*Max(lambda1,mu1)*invP[0][1] + invP[1][1]*Max(lambda2,mu2)*invP[1][1]; 337 *pM=M; 338 } 339 /*}}}1*/ 340 /*FUNCTION SimultaneousMatrixReduction_bamg{{{1*/ 341 void SimultaneousMatrixReduction_bamg( MetricAnIso M1, MetricAnIso M2,double & l1,double & l2, D2xD2 & V) { 173 void SimultaneousMatrixReduction( MetricAnIso M1, MetricAnIso M2,double & l1,double & l2, D2xD2 & V) { 342 174 double a11=M1.a11,a21=M1.a21,a22=M1.a22; 343 175 double b11=M2.a11,b21=M2.a21,b22=M2.a22; … … 353 185 const double c=-c21+a11*a22; 354 186 const double bb = b*b,ac= a*c; 187 const double b2a = b/(2*a); 355 188 const double delta = bb - 4 * ac; 356 if ( bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb) ){189 if ( (bb + Abs(ac) < 1.0e-34 ) || (delta < 1.0e-6*bb) ){ 357 190 // racine double; 358 191 if (Abs(a) < 1.e-30 ) -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2929 r2932 2008 2008 } 2009 2009 2010 2011 2010 /*Compute Metric from Hessian*/ 2012 2011 … … 2033 2032 2034 2033 //Apply Metric to vertex 2035 vertices[iv].m.IntersectWith _new(MVp);2034 vertices[iv].m.IntersectWith(MVp); 2036 2035 } 2037 2036 }//for all fields
Note:
See TracChangeset
for help on using the changeset viewer.
