![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2937
- Timestamp:
- 01/29/10 15:59:55 (15 years ago)
- Location:
- issm/trunk/src
- Files:
-
- 7 edited
-
c/Bamgx/Bamgx.cpp (modified) (2 diffs)
-
c/Bamgx/Mesh2.h (modified) (1 diff)
-
c/Bamgx/objects/Triangles.cpp (modified) (17 diffs)
-
c/Bamgx/objects/Vertex.cpp (modified) (1 diff)
-
c/objects/BamgOpts.h (modified) (1 diff)
-
m/classes/public/bamg.m (modified) (1 diff)
-
mex/Bamg/Bamg.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Bamgx.cpp
r2934 r2937 71 71 throw ErrorException(__FUNCT__,exprintf("maxsubdiv %g is should be in ]1,%g]",maxsubdiv,boundmaxsubdiv)); 72 72 } 73 // if iso, then anisomax = 1; 73 74 if (iso) anisomax=1; 74 75 // no metric -> no smoothing … … 76 77 NbSmooth=0; 77 78 } 79 80 /*If no mesh in input, generate one*/ 78 81 79 82 if(bamgmesh->NumTriangles==0){ -
issm/trunk/src/c/Bamgx/Mesh2.h
r2924 r2937 137 137 // ~Vertex(){} 138 138 Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1); 139 void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,BamgOpts* bamgopts); 139 140 int ref() const { return ReferenceNumber;} 140 141 -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2934 r2937 900 900 /*Options*/ 901 901 const int dim = 2; 902 int AbsError;903 902 double* s=NULL; 904 903 Int4 nbsol; 905 904 int* typsols=NULL; 906 Real8 hmin1;907 Real8 hmax1;908 Real8 coef;909 Real8 anisomax;910 Real8 CutOff;911 905 int NbJacobi; 912 int Rescaling;913 double power;914 906 int verbosity; 915 907 916 908 /*Recover options*/ 917 909 verbosity=bamgopts->verbose; 918 AbsError=bamgopts->AbsError;919 CutOff=bamgopts->cutoff;920 hmin1=bamgopts->hmin;921 hmax1=bamgopts->hmax;922 coef=bamgopts->coef;923 910 NbJacobi=bamgopts->nbjacobi; 924 Rescaling=bamgopts->Rescaling; //do normalization925 power=bamgopts->power;926 anisomax=bamgopts->anisomax;927 928 /*process options*/929 if (AbsError) CutOff=0.0;930 coef=sqrt(bamgopts->err)*coef;931 911 932 912 /*Get and process fields*/ … … 955 935 Int4 i,k,iA,iB,iC,iv; 956 936 R2 O(0,0); 957 int RelativeMetric = CutOff>1e-30;958 Real8 hmin = Max(hmin1,MinimalHmin());959 Real8 hmax = Min(hmax1,MaximalHmax());960 Real8 coef2 = 1/(coef*coef);961 937 double* ss=(double*)s; 962 938 double sA,sB,sC; … … 974 950 if(verbosity>1) { 975 951 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 976 printf(" coef = %g\n",coef);977 printf(" hmin = %g hmax = %g\n",hmin,hmax);978 printf(" anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power);979 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff);980 else printf(" Absolute error\n");981 952 } 982 953 … … 1049 1020 Real8 smin=ss[0],smax=ss[0]; 1050 1021 Real8 h1=1.e30,h2=1e-30,rx=0; 1051 Real8 coef = 1./(anisomax*anisomax);1052 1022 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1053 1023 int nbfield=typsols?sizeoftype[typsols[nusol]]:1; … … 1078 1048 Real8 sdelta=smax-smin; 1079 1049 Real8 absmax=Max(Abs(smin),Abs(smax)); 1080 Real8 cnorm =Rescaling ? coef2/sdelta : coef2;1081 1050 1082 1051 //display info 1083 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield);1052 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield); 1084 1053 1085 1054 //skip constant field … … 1236 1205 } 1237 1206 1238 //constuction of the metric from the Hessian dxdx. dxdy,dydy 1239 Real8 rCutOff=CutOff*absmax;// relative cut off 1240 1241 //loop over the nodes 1242 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 1243 1244 MetricIso Miso; 1245 Real8 ci ; 1246 1247 // compute norm of the solution 1248 if (RelativeMetric){ 1249 double xx =0,*sfk=sf+k; 1250 for (int ifield=0;ifield<nbfield;ifield++,sfk++){ 1251 xx += *sfk* *sfk; 1252 } 1253 xx=sqrt(xx); 1254 ci=coef2/Max(xx,rCutOff); 1255 } 1256 else ci=cnorm; 1257 1258 //initialize metric Miv with ci*H 1259 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 1260 1261 //Get eigen values and vectors of Miv 1262 MatVVP2x2 Vp(Miv); 1263 1264 //move eigen valuse to their absolute values 1265 Vp.Abs(); 1266 1267 //Allpy a power if requested by user 1268 if(power!=1.0) Vp.pow(power); 1269 1270 //get minimum and maximum eigen values 1271 h1=Min(h1,Vp.lmin()); 1272 h2=Max(h2,Vp.lmax()); 1273 1274 //modify eigen values according to hmin and hmax 1275 Vp.Maxh(hmax); 1276 Vp.Minh(hmin); 1277 1278 //multiply eigen values by coef 1279 Vp.BoundAniso2(coef); 1280 1281 //rebuild Metric from Vp 1282 Metric MVp(Vp); 1283 1284 //Apply Metric to vertex 1285 vertices[iv].m.IntersectWith(MVp); 1286 1287 //info to be displayed 1288 //rx = max(lmax/lmin) (anisotropy ratio) 1289 rx = Max(rx,Vp.Aniso2()); 1290 hn1=Min(hn1,Vp.lmin()); 1291 hn2=Max(hn2,Vp.lmax()); 1292 rnx = Max(rnx,Vp.Aniso2()); 1293 } 1294 1295 //display info 1296 if (verbosity>2) { 1297 1298 1299 printf(" Field %i of solution %i\n",nufield,nusol); 1300 printf(" before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5)); 1301 printf(" after bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); 1302 } 1207 /*Compute Metric from Hessian*/ 1208 for ( iv=0;iv<nbv;iv++){ 1209 vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts); 1210 } 1211 1303 1212 } // end of for all field 1304 1213 }// end for all solution … … 1325 1234 /*Options*/ 1326 1235 const int dim = 2; 1327 int AbsError;1328 1236 double* s=NULL; 1329 1237 Int4 nbsol; 1330 1238 int* typsols=NULL; 1331 Real8 hmin1;1332 Real8 hmax1;1333 Real8 coef;1334 Real8 anisomax;1335 Real8 CutOff;1336 1239 int NbJacobi; 1337 int Rescaling;1338 double power;1339 1240 int verbosity; 1340 1241 1341 1242 /*Recover options*/ 1342 1243 verbosity=bamgopts->verbose; 1343 AbsError=bamgopts->AbsError;1344 CutOff=bamgopts->cutoff;1345 hmin1=bamgopts->hmin;1346 hmax1=bamgopts->hmax;1347 coef=bamgopts->coef;1348 1244 NbJacobi=bamgopts->nbjacobi; 1349 Rescaling=bamgopts->Rescaling; //do normalization1350 power=bamgopts->power;1351 anisomax=bamgopts->anisomax;1352 1353 /*process options*/1354 if (AbsError) CutOff=0.0;1355 coef=sqrt(bamgopts->err)*coef;1356 1245 1357 1246 /*Get and process fields*/ … … 1380 1269 Int4 i,k,iA,iB,iC,iv; 1381 1270 R2 O(0,0); 1382 int RelativeMetric = CutOff>1e-30;1383 Real8 hmin = Max(hmin1,MinimalHmin());1384 Real8 hmax = Min(hmax1,MaximalHmax());1385 Real8 coef2 = 1/(coef*coef);1386 1271 double* ss=(double*)s; 1387 1272 double sA,sB,sC; … … 1399 1284 if(verbosity>1) { 1400 1285 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 1401 printf(" coef = %g\n",coef);1402 printf(" hmin = %g hmax = %g\n",hmin,hmax);1403 printf(" anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power);1404 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff);1405 else printf(" Absolute error\n");1406 1286 } 1407 1287 … … 1476 1356 Real8 smin=ss[0],smax=ss[0]; 1477 1357 Real8 h1=1.e30,h2=1e-30,rx=0; 1478 Real8 coef = 1./(anisomax*anisomax);1479 1358 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1480 1359 int nbfield=typsols?sizeoftype[typsols[nusol]]:1; … … 1505 1384 Real8 sdelta=smax-smin; 1506 1385 Real8 absmax=Max(Abs(smin),Abs(smax)); 1507 Real8 cnorm =Rescaling ? coef2/sdelta : coef2;1508 1386 1509 1387 //display info 1510 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield);1388 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield); 1511 1389 1512 1390 //skip constant field … … 1657 1535 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 1658 1536 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); 1659 Real8 coef = 1.0/(3*dd*det33);1660 Real8 coef2 = 2*coef;1661 Hxx *= coef2;1662 Hyy *= coef2;1663 Hxy *= coef2;1664 1537 if(nbb==0){ 1665 1538 dxdx[iA] += Hxx; … … 1739 1612 } 1740 1613 1741 //constuction of the metric from the Hessian dxdx. dxdy,dydy 1742 Real8 rCutOff=CutOff*absmax;// relative cut off 1743 1744 //loop over the nodes 1745 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 1746 1747 MetricIso Miso; 1748 Real8 ci ; 1749 1750 // compute norm of the solution 1751 if (RelativeMetric){ 1752 double xx =0,*sfk=sf+k; 1753 for (int ifield=0;ifield<nbfield;ifield++,sfk++){ 1754 xx += *sfk* *sfk; 1755 } 1756 xx=sqrt(xx); 1757 ci=coef2/Max(xx,rCutOff); 1758 } 1759 else ci=cnorm; 1760 1761 //initialize metric Miv with ci*H 1762 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 1763 1764 //Get eigen values and vectors of Miv 1765 MatVVP2x2 Vp(Miv); 1766 1767 //move eigen valuse to their absolute values 1768 Vp.Abs(); 1769 1770 //Allpy a power if requested by user 1771 if(power!=1.0) Vp.pow(power); 1772 1773 //get minimum and maximum eigen values 1774 h1=Min(h1,Vp.lmin()); 1775 h2=Max(h2,Vp.lmax()); 1776 1777 //modify eigen values according to hmin and hmax 1778 Vp.Maxh(hmax); 1779 Vp.Minh(hmin); 1780 1781 //multiply eigen values by coef 1782 Vp.BoundAniso2(coef); 1783 1784 //rebuild Metric from Vp 1785 Metric MVp(Vp); 1786 1787 //Apply Metric to vertex 1788 vertices[iv].m.IntersectWith(MVp); 1789 1790 //info to be displayed 1791 //rx = max(lmax/lmin) (anisotropy ratio) 1792 rx = Max(rx,Vp.Aniso2()); 1793 hn1=Min(hn1,Vp.lmin()); 1794 hn2=Max(hn2,Vp.lmax()); 1795 rnx = Max(rnx,Vp.Aniso2()); 1796 } 1797 1798 //display info 1799 if (verbosity>2) { 1800 1801 1802 printf(" Field %i of solution %i\n",nufield,nusol); 1803 printf(" before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5)); 1804 printf(" after bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); 1805 } 1614 /*Compute Metric from Hessian*/ 1615 for ( iv=0;iv<nbv;iv++){ 1616 vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts); 1617 } 1618 1806 1619 } // end of for all field 1807 1620 }// end for all solution … … 1828 1641 int* typsols=NULL; 1829 1642 int verbosity; 1830 Real8 hmin1;1831 Real8 hmax1;1832 double power;1833 double anisomax;1834 1643 1835 1644 int i,j,k,iA,iB,iC; 1836 1645 int iv,nbfield; 1837 double coef;1838 1646 1839 1647 /*Recover options*/ 1840 1648 verbosity=bamgopts->verbose; 1841 power=bamgopts->power;1842 hmin1=bamgopts->hmin;1843 hmax1=bamgopts->hmax;1844 1845 /*Recover options*/1846 verbosity=bamgopts->verbose;1847 hmin1=bamgopts->hmin;1848 hmax1=bamgopts->hmax;1849 anisomax=bamgopts->anisomax;1850 1649 1851 1650 /*Get and process fields*/ … … 1872 1671 1873 1672 //initialization of some variables 1874 Real8 hmin = Max(hmin1,MinimalHmin());1875 Real8 hmax = Min(hmax1,MaximalHmax());1876 1673 double* ss=(double*)s; 1877 1674 double sA,sB,sC; … … 1953 1750 for (Int4 nusol=0;nusol<nbsol;nusol++) { 1954 1751 int nbfield=typsols?sizeoftype[typsols[nusol]]:1; 1752 Real8 smin=ss[0],smax=ss[0]; 1753 1754 //only one field 1755 if (nbfield == 1) { 1756 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 1757 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){ 1758 smin=Min(smin,ss[k]); 1759 smax=Max(smax,ss[k]); 1760 } 1761 } 1762 1763 //vectorial case 1764 else{ 1765 for (iv=0,k=0;iv<nbv;iv++,k+=n ){ 1766 //compute v = √sum(s[i]^2) 1767 double v=0; 1768 for (int i=0;i<nbfield;i++){ 1769 v += ss[k+i]*ss[k+i]; 1770 } 1771 v = sqrt(v); 1772 smin=Min(smin,v); 1773 smax=Max(smax,v); 1774 } 1775 } 1776 Real8 sdelta=smax-smin; 1777 Real8 absmax=Max(Abs(smin),Abs(smax)); 1778 1779 //display info 1780 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield); 1781 1782 //skip constant field 1783 if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){ 1784 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 1785 continue; 1786 } 1787 1955 1788 1956 1789 //loop over all the fields of the solution … … 2017 1850 2018 1851 /*Compute Metric from Hessian*/ 2019 2020 //compute multiplicative coefficient (2/9 because it is 2d)2021 Real8 ci=(2.0/9.0)*(1.0/(bamgopts->err));2022 2023 1852 for ( iv=0;iv<nbv;iv++){ 2024 2025 //initialize metric Miv with ci*H 2026 Metric Miv(dxdx_vertex[iv]*ci,dxdy_vertex[iv]*ci,dydy_vertex[iv]*ci); 2027 2028 //Get eigen values and vectors of Miv 2029 MatVVP2x2 Vp(Miv); 2030 2031 //move eigen valuse to their absolute values 2032 Vp.Abs(); 2033 2034 //Apply a power if requested by user 2035 if(power!=1.0) Vp.pow(power); 2036 2037 //modify eigen values according to hmin and hmax 2038 Vp.Maxh(hmax); 2039 Vp.Minh(hmin); 2040 2041 //Bound anisotropy by 1/(anisomax)^2 2042 Vp.BoundAniso2(1/(anisomax*anisomax)); 2043 2044 //rebuild Metric from Vp 2045 Metric MVp(Vp); 2046 2047 //Apply Metric to vertex 2048 vertices[iv].m.IntersectWith(MVp); 2049 } 1853 vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[n*iv],bamgopts); 1854 } 1855 2050 1856 }//for all fields 2051 1857 }//for all solutions -
issm/trunk/src/c/Bamgx/objects/Vertex.cpp
r2865 r2937 127 127 } 128 128 /*}}}1*/ 129 /*FUNCTION Vertex::MetricFromHessian{{{1*/ 130 void Vertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,BamgOpts* bamgopts){ 131 /*Compute Metric from Hessian*/ 132 133 /*get options*/ 134 double power=(bamgopts->power); 135 double anisomax=(bamgopts->anisomax); 136 double CutOff=bamgopts->cutoff; 137 double hmin=(bamgopts->hmin); 138 double hmax=(bamgopts->hmax); 139 double coef=bamgopts->coef; 140 int Metrictype=(bamgopts->Metrictype); 141 int Rescaling=(bamgopts->Rescaling); 142 143 /*Intermediary*/ 144 double ci; 145 146 /*compute multiplicative coefficient depending on Metric Type (2/9 because it is 2d)*/ 147 148 //Absolute Error 149 /* 150 * 2 1 151 *Metric M = --- ------------ Abs(Hessian) 152 * 9 err * coeff^2 153 */ 154 if (Metrictype==0){ 155 if (Rescaling){ 156 ci= 1/(bamgopts->err*coef*coef) * 1/(smax-smin); 157 } 158 else{ 159 ci=1/(bamgopts->err*coef*coef); 160 } 161 } 162 163 //Relative Error 164 /* 165 * 2 1 Abs(Hessian) 166 *Metric M = --- ------------ ---------------------- 167 * 9 err * coeff^2 max( |s| , cutoff*max(|s|) ) 168 * 169 */ 170 else if (Metrictype==1){ 171 ci=1/(bamgopts->err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax)))); 172 } 173 174 //Interpolation Error 175 /* 176 * 2 1 177 *Metric M = --- ------------ Abs(Hessian) 178 * 9 err * coeff^2 179 */ 180 else if (Metrictype==2){ 181 ci=(2.0/9.0)*(1.0/(bamgopts->err)); 182 } 183 else{ 184 throw ErrorException(__FUNCT__,exprintf("Metrictype %i not supported yet (use 0,1 or 2(default))",Metrictype)); 185 } 186 187 //initialize metric Miv with ci*H 188 Metric Miv(Hxx*ci,Hyx*ci,Hyy*ci); 189 190 //Get eigen values and vectors of Miv 191 MatVVP2x2 Vp(Miv); 192 193 //move eigen valuse to their absolute values 194 Vp.Abs(); 195 196 //Apply a power if requested by user 197 if(power!=1.0) Vp.pow(power); 198 199 //modify eigen values according to hmin and hmax 200 Vp.Maxh(hmax); 201 Vp.Minh(hmin); 202 203 //Bound anisotropy by 1/(anisomax)^2 204 Vp.BoundAniso2(1/(anisomax*anisomax)); 205 206 //rebuild Metric from Vp 207 Metric MVp(Vp); 208 209 //Apply Metric to vertex 210 m.IntersectWith(MVp); 211 212 } 213 /*}}}1*/ 129 214 130 215 /*Intermediary*/ -
issm/trunk/src/c/objects/BamgOpts.h
r2934 r2937 16 16 int NbSmooth; 17 17 int Rescaling; 18 int Metrictype; 18 19 int nbjacobi; 19 20 int AbsError; -
issm/trunk/src/m/classes/public/bamg.m
r2934 r2937 84 84 bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1); 85 85 bamg_options.AbsError=getfieldvalue(options,'AbsError',0); 86 bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0); 86 bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',2); 87 bamg_options.Metrictype=getfieldvalue(options,'Metrictype',2); 87 88 bamg_options.NbSmooth=getfieldvalue(options,'NbSmooth',3); 88 89 bamg_options.omega=getfieldvalue(options,'omega',1.8); -
issm/trunk/src/mex/Bamg/Bamg.cpp
r2934 r2937 36 36 double err,errg,coef; 37 37 double power; 38 int Hessiantype, NbSmooth;38 int Hessiantype,Metrictype,NbSmooth; 39 39 int Rescaling,nbjacobi,AbsError; 40 40 double omega; … … 117 117 FetchData(&Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype")); 118 118 bamgopts.Hessiantype=Hessiantype; 119 FetchData(&Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype")); 120 bamgopts.Metrictype=Metrictype; 119 121 FetchData(&power,mxGetField(BAMGOPTIONS,0,"power")); 120 122 bamgopts.power=power;
Note:
See TracChangeset
for help on using the changeset viewer.
