![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2920
- Timestamp:
- 01/27/10 10:57:29 (15 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/Bamgx/objects/Triangles.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2917 r2920 889 889 printf(" output: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); 890 890 } 891 } 892 /*}}}1*/ 893 /*FUNCTION Triangles::BuildMetric0 (Green formula){{{1*/ 894 void Triangles::BuildMetric0(BamgOpts* bamgopts){ 895 // the array of solution s is store 896 // sol0,sol1,...,soln on vertex 0 897 // sol0,sol1,...,soln on vertex 1 898 // etc. 899 900 /*Options*/ 901 const int dim = 2; 902 int AbsError; 903 double* s=NULL; 904 Int4 nbsol; 905 int* typsols=NULL; 906 Real8 hmin1; 907 Real8 hmax1; 908 Real8 coef; 909 Real8 anisomax=1e30; 910 Real8 CutOff; 911 int NbJacobi; 912 int Rescaling; 913 double power; 914 int verbosity; 915 916 /*Recover options*/ 917 verbosity=bamgopts->verbose; 918 AbsError=bamgopts->AbsError; 919 CutOff=bamgopts->cutoff; 920 hmin1=bamgopts->hmin; 921 hmax1=bamgopts->hmax; 922 coef=bamgopts->coef; 923 NbJacobi=bamgopts->nbjacobi; 924 Rescaling=bamgopts->Rescaling; //do normalization 925 power=bamgopts->power; 926 927 /*process options*/ 928 if (AbsError) CutOff=0.0; 929 coef=sqrt(bamgopts->err)*coef; 930 931 /*Get and process fields*/ 932 s=bamgopts->field; 933 nbsol=1; //for now, only one field 934 typsols=(int*)xmalloc(1*sizeof(int)); 935 typsols[0]=0; // only one dof per node 936 937 //sizeoftype = {1,2,3,4} 938 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 939 940 // computation of the number of fields 941 Int4 ntmp = 0; 942 if (typsols){ 943 //if there is more than one dof per vertex for one field 944 //increase ntmp to take into account all the fields 945 //if nbsol=1 946 for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]]; 947 } 948 else ntmp = nbsol; 949 950 //n is the total number of fields 951 const Int4 n = ntmp; 952 953 //initialization of some variables 954 Int4 i,k,iA,iB,iC,iv; 955 R2 O(0,0); 956 int RelativeMetric = CutOff>1e-30; 957 Real8 hmin = Max(hmin1,MinimalHmin()); 958 Real8 hmax = Min(hmax1,MaximalHmax()); 959 Real8 coef2 = 1/(coef*coef); 960 double* ss=(double*)s; 961 double sA,sB,sC; 962 Real8* detT = new Real8[nbt]; 963 Real8* Mmass= new Real8[nbv]; 964 Real8* Mmassxx= new Real8[nbv]; 965 Real8* dxdx= new Real8[nbv]; 966 Real8* dxdy= new Real8[nbv]; 967 Real8* dydy= new Real8[nbv]; 968 Real8* workT= new Real8[nbt]; 969 Real8* workV= new Real8[nbv]; 970 int* OnBoundary = new int[nbv]; 971 972 //display infos 973 if(verbosity>1) { 974 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 975 printf(" coef = %g\n",coef); 976 printf(" hmin = %g hmax = %g\n",hmin,hmax); 977 printf(" anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power); 978 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); 979 else printf(" Absolute error\n"); 980 } 981 982 //initialize Mmass, OnBoundary and Massxx by zero 983 for (iv=0;iv<nbv;iv++){ 984 Mmass[iv]=0; 985 OnBoundary[iv]=0; 986 Mmassxx[iv]=0; 987 } 988 989 //Build detT Mmas Mmassxx workT and OnBoundary 990 for (i=0;i<nbt;i++){ 991 992 //lopp over the real triangles (no boundary elements) 993 if(triangles[i].link){ 994 995 //get current triangle t 996 const Triangle &t=triangles[i]; 997 998 // coor of 3 vertices 999 R2 A=t[0]; 1000 R2 B=t[1]; 1001 R2 C=t[2]; 1002 1003 // number of the 3 vertices 1004 iA = Number(t[0]); 1005 iB = Number(t[1]); 1006 iC = Number(t[2]); 1007 1008 //compute triangle determinant (2*Area) 1009 Real8 dett = bamg::Area2(A,B,C); 1010 detT[i]=dett; 1011 dett /= 6; 1012 1013 // construction of OnBoundary (flag=1 if on boundary, else 0) 1014 int nbb=0; 1015 for(int j=0;j<3;j++){ 1016 //get adjacent triangle 1017 Triangle *ta=t.Adj(j); 1018 //if there is no adjacent triangle, the edge of the triangle t is on boundary 1019 if ( !ta || !ta->link){ 1020 //mark the two vertices of the edge as OnBoundary 1021 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1; 1022 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1; 1023 nbb++; 1024 } 1025 } 1026 1027 //number of vertices on boundary for current triangle t 1028 workT[i] = nbb; 1029 1030 //Build Mmass Mmass[i] = Mmass[i] + Area/3 1031 Mmass[iA] += dett; 1032 Mmass[iB] += dett; 1033 Mmass[iC] += dett; 1034 1035 //Build Massxx = Mmass 1036 Mmassxx[iA] += dett; 1037 Mmassxx[iB] += dett; 1038 Mmassxx[iC] += dett; 1039 } 1040 1041 //else: the triangle is a boundary triangle -> workT=-1 1042 else workT[i]=-1; 1043 } 1044 1045 //for all Solution 1046 for (Int4 nusol=0;nusol<nbsol;nusol++) { 1047 1048 Real8 smin=ss[0],smax=ss[0]; 1049 Real8 h1=1.e30,h2=1e-30,rx=0; 1050 Real8 coef = 1./(anisomax*anisomax); 1051 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1052 int nbfield=typsols?sizeoftype[typsols[nusol]]:1; 1053 1054 //only one field 1055 if (nbfield == 1) { 1056 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 1057 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){ 1058 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1059 smin=Min(smin,ss[k]); 1060 smax=Max(smax,ss[k]); 1061 } 1062 } 1063 1064 //vectorial case 1065 else{ 1066 for (iv=0,k=0;iv<nbv;iv++,k+=n ){ 1067 //compute v = √sum(s[i]^2) 1068 double v=0; 1069 for (int i=0;i<nbfield;i++){ 1070 v += ss[k+i]*ss[k+i]; 1071 } 1072 v = sqrt(v); 1073 smin=Min(smin,v); 1074 smax=Max(smax,v); 1075 } 1076 } 1077 Real8 sdelta=smax-smin; 1078 Real8 absmax=Max(Abs(smin),Abs(smax)); 1079 Real8 cnorm =Rescaling ? coef2/sdelta : coef2; 1080 1081 //display info 1082 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); 1083 1084 //skip constant field 1085 if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){ 1086 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 1087 continue; 1088 } 1089 1090 //pointer toward ss that is also a pointer toward s (solutions) 1091 double* sf=ss; 1092 1093 //loop over all the fields of the solution 1094 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){ 1095 //ss++ so that for each iteration ss points toward the right field 1096 1097 //initialize the hessian matrix 1098 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1099 1100 //loop over the triangles 1101 for (i=0;i<nbt;i++){ 1102 1103 //for real all triangles 1104 if(triangles[i].link){ 1105 1106 // coor of 3 vertices 1107 R2 A=triangles[i][0]; 1108 R2 B=triangles[i][1]; 1109 R2 C=triangles[i][2]; 1110 1111 //warning: the normal is internal and the size is the length of the edge 1112 R2 nAB = Orthogonal(B-A); 1113 R2 nBC = Orthogonal(C-B); 1114 R2 nCA = Orthogonal(A-C); 1115 //note that : nAB + nBC + nCA == 0 1116 1117 // number of the 3 vertices 1118 iA = Number(triangles[i][0]); 1119 iB = Number(triangles[i][1]); 1120 iC = Number(triangles[i][2]); 1121 1122 // for the test of boundary edge 1123 // the 3 adj triangles 1124 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 1125 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 1126 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 1127 1128 // value of the P1 fonction on 3 vertices 1129 sA = ss[iA*n]; 1130 sB = ss[iB*n]; 1131 sC = ss[iC*n]; 1132 1133 /*The nodal functions are such that for a vertex A: 1134 N_A(x,y)=alphaA x + beta_A y +gamma_A 1135 N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 1136 solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 1137 leads to: 1138 N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 1139 and this gives: 1140 alpha_A = (yB-yC)/(2*Area(T)) 1141 beta_A = (xC-xB)/(2*Area(T)) 1142 and therefore: 1143 grad N_A = nA / detT 1144 for an interpolation of a solution s: 1145 grad(s) = s * sum_{i=A,B,C} grad(N_i) */ 1146 1147 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i]; 1148 1149 //Use Green to compute Hessian Matrix 1150 1151 // if edge on boundary no contribution => normal = 0 1152 if ( !tBC || !tBC->link ) nBC=O; 1153 if ( !tCA || !tCA->link ) nCA=O; 1154 if ( !tAB || !tAB->link ) nAB=O; 1155 1156 // remark we forgot a 1/2 because 1157 // int_{edge} w_i = 1/2 if i is in edge 1158 // 0 if not 1159 // if we don't take the boundary 1160 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; 1161 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; 1162 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; 1163 1164 //warning optimization (1) the division by 2 is done on the metric construction 1165 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; 1166 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; 1167 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 1168 1169 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; 1170 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; 1171 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; 1172 1173 } // for real all triangles 1174 } 1175 1176 Int4 kk=0; 1177 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 1178 if(Mmassxx[iv]>0){ 1179 dxdx[iv] /= 2*Mmassxx[iv]; 1180 // warning optimization (1) on term dxdy[iv]*ci/2 1181 dxdy[iv] /= 4*Mmassxx[iv]; 1182 dydy[iv] /= 2*Mmassxx[iv]; 1183 // Compute the matrix with abs(eigen value) 1184 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 1185 MatVVP2x2 Vp(M); 1186 Vp.Abs(); 1187 M = Vp; 1188 dxdx[iv] = M.a11; 1189 dxdy[iv] = M.a21; 1190 dydy[iv] = M.a22; 1191 } 1192 else kk++; 1193 } 1194 1195 // correction of second derivative 1196 // by a laplacien 1197 Real8* d2[3] = {dxdx, dxdy, dydy}; 1198 Real8* dd; 1199 for (int xy = 0;xy<3;xy++) { 1200 dd = d2[xy]; 1201 // do leat 2 iteration for boundary problem 1202 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ 1203 for (i=0;i<nbt;i++) 1204 if(triangles[i].link){// the real triangles 1205 // number of the 3 vertices 1206 iA = Number(triangles[i][0]); 1207 iB = Number(triangles[i][1]); 1208 iC = Number(triangles[i][2]); 1209 Real8 cc=3; 1210 if(ijacobi==0) 1211 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 1212 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 1213 } 1214 for (iv=0;iv<nbv;iv++) workV[iv]=0; 1215 1216 for (i=0;i<nbt;i++){ 1217 if(triangles[i].link){ // the real triangles 1218 // number of the 3 vertices 1219 iA = Number(triangles[i][0]); 1220 iB = Number(triangles[i][1]); 1221 iC = Number(triangles[i][2]); 1222 Real8 cc = workT[i]*detT[i]; 1223 workV[iA] += cc; 1224 workV[iB] += cc; 1225 workV[iC] += cc; 1226 } 1227 } 1228 1229 for (iv=0;iv<nbv;iv++){ 1230 if( ijacobi<NbJacobi || OnBoundary[iv]){ 1231 dd[iv] = workV[iv]/(Mmass[iv]*6); 1232 } 1233 } 1234 } 1235 } 1236 1237 //constuction of the metric from the Hessian dxdx. dxdy,dydy 1238 Real8 rCutOff=CutOff*absmax;// relative cut off 1239 1240 //loop over the nodes 1241 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 1242 1243 MetricIso Miso; 1244 Real8 ci ; 1245 1246 // compute norm of the solution 1247 if (RelativeMetric){ 1248 double xx =0,*sfk=sf+k; 1249 for (int ifield=0;ifield<nbfield;ifield++,sfk++){ 1250 xx += *sfk* *sfk; 1251 } 1252 xx=sqrt(xx); 1253 ci=coef2/Max(xx,rCutOff); 1254 } 1255 else ci=cnorm; 1256 1257 //initialize metric Miv with ci*H 1258 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 1259 1260 //Get eigen values and vectors of Miv 1261 MatVVP2x2 Vp(Miv); 1262 1263 //move eigen valuse to their absolute values 1264 Vp.Abs(); 1265 1266 //Allpy a power if requested by user 1267 if(power!=1.0) Vp.pow(power); 1268 1269 //get minimum and maximum eigen values 1270 h1=Min(h1,Vp.lmin()); 1271 h2=Max(h2,Vp.lmax()); 1272 1273 //modify eigen values according to hmin and hmax 1274 Vp.Maxh(hmin); 1275 Vp.Minh(hmax); 1276 1277 //multiply eigen values by coef 1278 Vp.BoundAniso2(coef); 1279 1280 //rebuild Metric from Vp 1281 Metric MVp(Vp); 1282 1283 //Apply Metric to vertex 1284 vertices[iv].m.IntersectWith(MVp); 1285 1286 //info to be displayed 1287 //rx = max(lmax/lmin) (anisotropy ratio) 1288 rx = Max(rx,Vp.Aniso2()); 1289 hn1=Min(hn1,Vp.lmin()); 1290 hn2=Max(hn2,Vp.lmax()); 1291 rnx = Max(rnx,Vp.Aniso2()); 1292 } 1293 1294 //display info 1295 if (verbosity>2) { 1296 1297 1298 printf(" Field %i of solution %i\n",nufield,nusol); 1299 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)); 1300 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)); 1301 } 1302 } // end of for all field 1303 }// end for all solution 1304 1305 delete [] detT; 1306 delete [] Mmass; 1307 delete [] dxdx; 1308 delete [] dxdy; 1309 delete [] dydy; 1310 delete [] workT; 1311 delete [] workV; 1312 delete [] Mmassxx; 1313 delete [] OnBoundary; 1314 1315 } 1316 /*}}}1*/ 1317 /*FUNCTION Triangles::BuildMetric1 (from P2 on 4T){{{1*/ 1318 void Triangles::BuildMetric1(BamgOpts* bamgopts){ 1319 // the array of solution s is store 1320 // sol0,sol1,...,soln on vertex 0 1321 // sol0,sol1,...,soln on vertex 1 1322 // etc. 1323 1324 /*Options*/ 1325 const int dim = 2; 1326 int AbsError; 1327 double* s=NULL; 1328 Int4 nbsol; 1329 int* typsols=NULL; 1330 Real8 hmin1; 1331 Real8 hmax1; 1332 Real8 coef; 1333 Real8 anisomax=1e30; 1334 Real8 CutOff; 1335 int NbJacobi; 1336 int Rescaling; 1337 double power; 1338 int verbosity; 1339 1340 /*Recover options*/ 1341 verbosity=bamgopts->verbose; 1342 AbsError=bamgopts->AbsError; 1343 CutOff=bamgopts->cutoff; 1344 hmin1=bamgopts->hmin; 1345 hmax1=bamgopts->hmax; 1346 coef=bamgopts->coef; 1347 NbJacobi=bamgopts->nbjacobi; 1348 Rescaling=bamgopts->Rescaling; //do normalization 1349 power=bamgopts->power; 1350 1351 /*process options*/ 1352 if (AbsError) CutOff=0.0; 1353 coef=sqrt(bamgopts->err)*coef; 1354 1355 /*Get and process fields*/ 1356 s=bamgopts->field; 1357 nbsol=1; //for now, only one field 1358 typsols=(int*)xmalloc(1*sizeof(int)); 1359 typsols[0]=0; // only one dof per node 1360 1361 //sizeoftype = {1,2,3,4} 1362 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 1363 1364 // computation of the number of fields 1365 Int4 ntmp = 0; 1366 if (typsols){ 1367 //if there is more than one dof per vertex for one field 1368 //increase ntmp to take into account all the fields 1369 //if nbsol=1 1370 for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]]; 1371 } 1372 else ntmp = nbsol; 1373 1374 //n is the total number of fields 1375 const Int4 n = ntmp; 1376 1377 //initialization of some variables 1378 Int4 i,k,iA,iB,iC,iv; 1379 R2 O(0,0); 1380 int RelativeMetric = CutOff>1e-30; 1381 Real8 hmin = Max(hmin1,MinimalHmin()); 1382 Real8 hmax = Min(hmax1,MaximalHmax()); 1383 Real8 coef2 = 1/(coef*coef); 1384 double* ss=(double*)s; 1385 double sA,sB,sC; 1386 Real8* detT = new Real8[nbt]; 1387 Real8* Mmass= new Real8[nbv]; 1388 Real8* Mmassxx= new Real8[nbv]; 1389 Real8* dxdx= new Real8[nbv]; 1390 Real8* dxdy= new Real8[nbv]; 1391 Real8* dydy= new Real8[nbv]; 1392 Real8* workT= new Real8[nbt]; 1393 Real8* workV= new Real8[nbv]; 1394 int* OnBoundary = new int[nbv]; 1395 1396 //display infos 1397 if(verbosity>1) { 1398 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 1399 printf(" coef = %g\n",coef); 1400 printf(" hmin = %g hmax = %g\n",hmin,hmax); 1401 printf(" anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power); 1402 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); 1403 else printf(" Absolute error\n"); 1404 } 1405 1406 //initialize Mmass, OnBoundary and Massxx by zero 1407 for (iv=0;iv<nbv;iv++){ 1408 Mmass[iv]=0; 1409 OnBoundary[iv]=0; 1410 Mmassxx[iv]=0; 1411 } 1412 1413 //Build detT Mmas Mmassxx workT and OnBoundary 1414 for (i=0;i<nbt;i++){ 1415 1416 //lopp over the real triangles (no boundary elements) 1417 if(triangles[i].link){ 1418 1419 //get current triangle t 1420 const Triangle &t=triangles[i]; 1421 1422 // coor of 3 vertices 1423 R2 A=t[0]; 1424 R2 B=t[1]; 1425 R2 C=t[2]; 1426 1427 // number of the 3 vertices 1428 iA = Number(t[0]); 1429 iB = Number(t[1]); 1430 iC = Number(t[2]); 1431 1432 //compute triangle determinant (2*Area) 1433 Real8 dett = bamg::Area2(A,B,C); 1434 detT[i]=dett; 1435 dett /= 6; 1436 1437 // construction of OnBoundary (flag=1 if on boundary, else 0) 1438 int nbb=0; 1439 for(int j=0;j<3;j++){ 1440 //get adjacent triangle 1441 Triangle *ta=t.Adj(j); 1442 //if there is no adjacent triangle, the edge of the triangle t is on boundary 1443 if ( !ta || !ta->link){ 1444 //mark the two vertices of the edge as OnBoundary 1445 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1; 1446 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1; 1447 nbb++; 1448 } 1449 } 1450 1451 //number of vertices on boundary for current triangle t 1452 workT[i] = nbb; 1453 1454 //Build Mmass Mmass[i] = Mmass[i] + Area/3 1455 Mmass[iA] += dett; 1456 Mmass[iB] += dett; 1457 Mmass[iC] += dett; 1458 1459 //Build Massxx = Mmass 1460 if(nbb==0){ 1461 Mmassxx[iA] += dett; 1462 Mmassxx[iB] += dett; 1463 Mmassxx[iC] += dett; 1464 } 1465 } 1466 1467 //else: the triangle is a boundary triangle -> workT=-1 1468 else workT[i]=-1; 1469 } 1470 1471 //for all Solution 1472 for (Int4 nusol=0;nusol<nbsol;nusol++) { 1473 1474 Real8 smin=ss[0],smax=ss[0]; 1475 Real8 h1=1.e30,h2=1e-30,rx=0; 1476 Real8 coef = 1./(anisomax*anisomax); 1477 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1478 int nbfield=typsols?sizeoftype[typsols[nusol]]:1; 1479 1480 //only one field 1481 if (nbfield == 1) { 1482 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 1483 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){ 1484 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1485 smin=Min(smin,ss[k]); 1486 smax=Max(smax,ss[k]); 1487 } 1488 } 1489 1490 //vectorial case 1491 else{ 1492 for (iv=0,k=0;iv<nbv;iv++,k+=n ){ 1493 //compute v = √sum(s[i]^2) 1494 double v=0; 1495 for (int i=0;i<nbfield;i++){ 1496 v += ss[k+i]*ss[k+i]; 1497 } 1498 v = sqrt(v); 1499 smin=Min(smin,v); 1500 smax=Max(smax,v); 1501 } 1502 } 1503 Real8 sdelta=smax-smin; 1504 Real8 absmax=Max(Abs(smin),Abs(smax)); 1505 Real8 cnorm =Rescaling ? coef2/sdelta : coef2; 1506 1507 //display info 1508 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); 1509 1510 //skip constant field 1511 if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){ 1512 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 1513 continue; 1514 } 1515 1516 //pointer toward ss that is also a pointer toward s (solutions) 1517 double* sf=ss; 1518 1519 //loop over all the fields of the solution 1520 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){ 1521 //ss++ so that for each iteration ss points toward the right field 1522 1523 //initialize the hessian matrix 1524 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1525 1526 //loop over the triangles 1527 for (i=0;i<nbt;i++){ 1528 1529 //for real all triangles 1530 if(triangles[i].link){ 1531 1532 // coor of 3 vertices 1533 R2 A=triangles[i][0]; 1534 R2 B=triangles[i][1]; 1535 R2 C=triangles[i][2]; 1536 1537 //warning: the normal is internal and the size is the length of the edge 1538 R2 nAB = Orthogonal(B-A); 1539 R2 nBC = Orthogonal(C-B); 1540 R2 nCA = Orthogonal(A-C); 1541 //note that : nAB + nBC + nCA == 0 1542 1543 // number of the 3 vertices 1544 iA = Number(triangles[i][0]); 1545 iB = Number(triangles[i][1]); 1546 iC = Number(triangles[i][2]); 1547 1548 // for the test of boundary edge 1549 // the 3 adj triangles 1550 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 1551 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 1552 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 1553 1554 // value of the P1 fonction on 3 vertices 1555 sA = ss[iA*n]; 1556 sB = ss[iB*n]; 1557 sC = ss[iC*n]; 1558 1559 /*The nodal functions are such that for a vertex A: 1560 N_A(x,y)=alphaA x + beta_A y +gamma_A 1561 N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 1562 solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 1563 leads to: 1564 N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 1565 and this gives: 1566 alpha_A = (yB-yC)/(2*Area(T)) 1567 beta_A = (xC-xB)/(2*Area(T)) 1568 and therefore: 1569 grad N_A = nA / detT 1570 for an interpolation of a solution s: 1571 grad(s) = s * sum_{i=A,B,C} grad(N_i) */ 1572 1573 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i]; 1574 1575 //from P2 on 4T to compute Hessian 1576 1577 int nbb=0; 1578 Real8 dd = detT[i]; 1579 Real8 lla,llb,llc,llf; 1580 Real8 taa[3][3],bb[3]; 1581 1582 // construction of the transpose of lin system 1583 for (int j=0;j<3;j++){ 1584 int ie = OppositeEdge[j]; 1585 TriangleAdjacent ta = triangles[i].Adj(ie); 1586 Triangle* tt = ta; 1587 1588 //if the adjacent triangle is not a boundary triangle: 1589 if (tt && tt->link){ 1590 Vertex &v = *ta.OppositeVertex(); 1591 R2 V = v; 1592 Int4 iV = Number(v); 1593 Real8 lA = bamg::Area2(V,B,C)/dd; 1594 Real8 lB = bamg::Area2(A,V,C)/dd; 1595 Real8 lC = bamg::Area2(A,B,V)/dd; 1596 taa[0][j] = lB*lC; 1597 taa[1][j] = lC*lA; 1598 taa[2][j] = lA*lB; 1599 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; 1600 bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ; 1601 } 1602 else{ 1603 nbb++; 1604 taa[0][j]=0; 1605 taa[1][j]=0; 1606 taa[2][j]=0; 1607 taa[j][j]=1; 1608 bb[j]=0; 1609 } 1610 } 1611 1612 // resolution of 3x3 linear system transpose 1613 Real8 det33 = det3x3(taa[0],taa[1],taa[2]); 1614 Real8 cBC = det3x3(bb,taa[1],taa[2]); 1615 Real8 cCA = det3x3(taa[0],bb,taa[2]); 1616 Real8 cAB = det3x3(taa[0],taa[1],bb); 1617 1618 if (!det33){ 1619 throw ErrorException(__FUNCT__,exprintf("!det33")); 1620 } 1621 // computation of the Hessian in the element 1622 1623 // H( li*lj) = grad li grad lj + grad lj grad lj 1624 // grad li = njk / detT ; with i j k =(A,B,C) 1625 Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x); 1626 Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y); 1627 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 1628 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); 1629 Real8 coef = 1.0/(3*dd*det33); 1630 Real8 coef2 = 2*coef; 1631 Hxx *= coef2; 1632 Hyy *= coef2; 1633 Hxy *= coef2; 1634 if(nbb==0){ 1635 dxdx[iA] += Hxx; 1636 dydy[iA] += Hyy; 1637 dxdy[iA] += Hxy; 1638 1639 dxdx[iB] += Hxx; 1640 dydy[iB] += Hyy; 1641 dxdy[iB] += Hxy; 1642 1643 dxdx[iC] += Hxx; 1644 dydy[iC] += Hyy; 1645 dxdy[iC] += Hxy; 1646 } 1647 } // for real all triangles 1648 } 1649 1650 Int4 kk=0; 1651 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 1652 if(Mmassxx[iv]>0){ 1653 dxdx[iv] /= 2*Mmassxx[iv]; 1654 // warning optimization (1) on term dxdy[iv]*ci/2 1655 dxdy[iv] /= 4*Mmassxx[iv]; 1656 dydy[iv] /= 2*Mmassxx[iv]; 1657 // Compute the matrix with abs(eigen value) 1658 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 1659 MatVVP2x2 Vp(M); 1660 Vp.Abs(); 1661 M = Vp; 1662 dxdx[iv] = M.a11; 1663 dxdy[iv] = M.a21; 1664 dydy[iv] = M.a22; 1665 } 1666 else kk++; 1667 } 1668 1669 // correction of second derivative 1670 // by a laplacien 1671 Real8* d2[3] = {dxdx, dxdy, dydy}; 1672 Real8* dd; 1673 for (int xy = 0;xy<3;xy++) { 1674 dd = d2[xy]; 1675 // do leat 2 iteration for boundary problem 1676 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ 1677 for (i=0;i<nbt;i++) 1678 if(triangles[i].link){// the real triangles 1679 // number of the 3 vertices 1680 iA = Number(triangles[i][0]); 1681 iB = Number(triangles[i][1]); 1682 iC = Number(triangles[i][2]); 1683 Real8 cc=3; 1684 if(ijacobi==0) 1685 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 1686 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 1687 } 1688 for (iv=0;iv<nbv;iv++) workV[iv]=0; 1689 1690 for (i=0;i<nbt;i++){ 1691 if(triangles[i].link){ // the real triangles 1692 // number of the 3 vertices 1693 iA = Number(triangles[i][0]); 1694 iB = Number(triangles[i][1]); 1695 iC = Number(triangles[i][2]); 1696 Real8 cc = workT[i]*detT[i]; 1697 workV[iA] += cc; 1698 workV[iB] += cc; 1699 workV[iC] += cc; 1700 } 1701 } 1702 1703 for (iv=0;iv<nbv;iv++){ 1704 if( ijacobi<NbJacobi || OnBoundary[iv]){ 1705 dd[iv] = workV[iv]/(Mmass[iv]*6); 1706 } 1707 } 1708 } 1709 } 1710 1711 //constuction of the metric from the Hessian dxdx. dxdy,dydy 1712 Real8 rCutOff=CutOff*absmax;// relative cut off 1713 1714 //loop over the nodes 1715 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 1716 1717 MetricIso Miso; 1718 Real8 ci ; 1719 1720 // compute norm of the solution 1721 if (RelativeMetric){ 1722 double xx =0,*sfk=sf+k; 1723 for (int ifield=0;ifield<nbfield;ifield++,sfk++){ 1724 xx += *sfk* *sfk; 1725 } 1726 xx=sqrt(xx); 1727 ci=coef2/Max(xx,rCutOff); 1728 } 1729 else ci=cnorm; 1730 1731 //initialize metric Miv with ci*H 1732 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 1733 1734 //Get eigen values and vectors of Miv 1735 MatVVP2x2 Vp(Miv); 1736 1737 //move eigen valuse to their absolute values 1738 Vp.Abs(); 1739 1740 //Allpy a power if requested by user 1741 if(power!=1.0) Vp.pow(power); 1742 1743 //get minimum and maximum eigen values 1744 h1=Min(h1,Vp.lmin()); 1745 h2=Max(h2,Vp.lmax()); 1746 1747 //modify eigen values according to hmin and hmax 1748 Vp.Maxh(hmin); 1749 Vp.Minh(hmax); 1750 1751 //multiply eigen values by coef 1752 Vp.BoundAniso2(coef); 1753 1754 //rebuild Metric from Vp 1755 Metric MVp(Vp); 1756 1757 //Apply Metric to vertex 1758 vertices[iv].m.IntersectWith(MVp); 1759 1760 //info to be displayed 1761 //rx = max(lmax/lmin) (anisotropy ratio) 1762 rx = Max(rx,Vp.Aniso2()); 1763 hn1=Min(hn1,Vp.lmin()); 1764 hn2=Max(hn2,Vp.lmax()); 1765 rnx = Max(rnx,Vp.Aniso2()); 1766 } 1767 1768 //display info 1769 if (verbosity>2) { 1770 1771 1772 printf(" Field %i of solution %i\n",nufield,nusol); 1773 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)); 1774 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)); 1775 } 1776 } // end of for all field 1777 }// end for all solution 1778 1779 delete [] detT; 1780 delete [] Mmass; 1781 delete [] dxdx; 1782 delete [] dxdy; 1783 delete [] dydy; 1784 delete [] workT; 1785 delete [] workV; 1786 delete [] Mmassxx; 1787 delete [] OnBoundary; 1788 1789 } 1790 /*}}}1*/ 1791 /*FUNCTION Triangles::BuildMetric2{{{1*/ 1792 void Triangles::BuildMetric2(BamgOpts* bamgopts){ 1793 1794 throw ErrorException(__FUNCT__,exprintf("not supported yet")); 891 1795 } 892 1796 /*}}}1*/ … … 2738 3642 2739 3643 /*Options*/ 2740 const int dim = 2; 2741 int AbsError; 2742 double* s=NULL; 2743 Int4 nbsol; 2744 int* typsols=NULL; 2745 Real8 hmin1; 2746 Real8 hmax1; 2747 Real8 coef; 2748 Real8 anisomax=1e30; 2749 Real8 CutOff; 2750 int NbJacobi; 2751 int Rescaling; 2752 double power; 2753 int Hessiantype; 2754 int verbosity; 2755 2756 /*Recover options*/ 2757 verbosity=bamgopts->verbose; 2758 AbsError=bamgopts->AbsError; 2759 CutOff=bamgopts->cutoff; 2760 hmin1=bamgopts->hmin; 2761 hmax1=bamgopts->hmax; 2762 coef=bamgopts->coef; 2763 NbJacobi=bamgopts->nbjacobi; 2764 Rescaling=bamgopts->Rescaling; //do normalization 2765 power=bamgopts->power; 2766 Hessiantype=bamgopts->Hessiantype; 2767 2768 /*process options*/ 2769 if (AbsError) CutOff=0.0; 2770 coef=sqrt(bamgopts->err)*coef; 2771 2772 /*Get and process fields*/ 2773 s=bamgopts->field; 2774 nbsol=1; //for now, only one field 2775 typsols=(int*)xmalloc(1*sizeof(int)); 2776 typsols[0]=0; // only one dof per node 2777 2778 //sizeoftype = {1,2,3,4} 2779 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 2780 2781 // computation of the number of fields 2782 Int4 ntmp = 0; 2783 if (typsols){ 2784 //if there is more than one dof per vertex for one field 2785 //increase ntmp to take into account all the fields 2786 //if nbsol=1 2787 for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]]; 3644 int Hessiantype=bamgopts->Hessiantype; 3645 3646 if (Hessiantype==0){ 3647 BuildMetric0(bamgopts); 2788 3648 } 2789 else ntmp = nbsol; 2790 2791 //n is the total number of fields 2792 const Int4 n = ntmp; 2793 2794 //initialization of some variables 2795 Int4 i,k,iA,iB,iC,iv; 2796 R2 O(0,0); 2797 int RelativeMetric = CutOff>1e-30; 2798 Real8 hmin = Max(hmin1,MinimalHmin()); 2799 Real8 hmax = Min(hmax1,MaximalHmax()); 2800 Real8 coef2 = 1/(coef*coef); 2801 double* ss=(double*)s; 2802 double sA,sB,sC; 2803 Real8* detT = new Real8[nbt]; 2804 Real8* Mmass= new Real8[nbv]; 2805 Real8* Mmassxx= new Real8[nbv]; 2806 Real8* dxdx= new Real8[nbv]; 2807 Real8* dxdy= new Real8[nbv]; 2808 Real8* dydy= new Real8[nbv]; 2809 Real8* workT= new Real8[nbt]; 2810 Real8* workV= new Real8[nbv]; 2811 int* OnBoundary = new int[nbv]; 2812 2813 //display infos 2814 if(verbosity>1) { 2815 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 2816 printf(" coef = %g\n",coef); 2817 printf(" hmin = %g hmax = %g\n",hmin,hmax); 2818 printf(" anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power); 2819 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); 2820 else printf(" Absolute error\n"); 3649 else if (Hessiantype==1){ 3650 BuildMetric1(bamgopts); 2821 3651 } 2822 2823 //initialize Mmass, OnBoundary and Massxx by zero 2824 for (iv=0;iv<nbv;iv++){ 2825 Mmass[iv]=0; 2826 OnBoundary[iv]=0; 2827 Mmassxx[iv]=0; 3652 else{ 3653 throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (0->use Green formula, 1-> from P2 on 4T, 2-> double P2 projection)",Hessiantype)); 2828 3654 } 2829 2830 //Build detT Mmas Mmassxx workT and OnBoundary2831 for (i=0;i<nbt;i++){2832 2833 //lopp over the real triangles (no boundary elements)2834 if(triangles[i].link){2835 2836 //get current triangle t2837 const Triangle &t=triangles[i];2838 2839 // coor of 3 vertices2840 R2 A=t[0];2841 R2 B=t[1];2842 R2 C=t[2];2843 2844 // number of the 3 vertices2845 iA = Number(t[0]);2846 iB = Number(t[1]);2847 iC = Number(t[2]);2848 2849 //compute triangle determinant (2*Area)2850 Real8 dett = bamg::Area2(A,B,C);2851 detT[i]=dett;2852 dett /= 6;2853 2854 // construction of OnBoundary (flag=1 if on boundary, else 0)2855 int nbb=0;2856 for(int j=0;j<3;j++){2857 //get adjacent triangle2858 Triangle *ta=t.Adj(j);2859 //if there is no adjacent triangle, the edge of the triangle t is on boundary2860 if ( !ta || !ta->link){2861 //mark the two vertices of the edge as OnBoundary2862 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1;2863 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1;2864 nbb++;2865 }2866 }2867 2868 //number of vertices on boundary for current triangle t2869 workT[i] = nbb;2870 2871 //Build Mmass Mmass[i] = Mmass[i] + Area/32872 Mmass[iA] += dett;2873 Mmass[iB] += dett;2874 Mmass[iC] += dett;2875 2876 //Build Massxx = Mmass2877 if((nbb==0) || (Hessiantype==0)){2878 Mmassxx[iA] += dett;2879 Mmassxx[iB] += dett;2880 Mmassxx[iC] += dett;2881 }2882 }2883 2884 //else: the triangle is a boundary triangle -> workT=-12885 else workT[i]=-1;2886 }2887 2888 //for all Solution2889 for (Int4 nusol=0;nusol<nbsol;nusol++) {2890 2891 Real8 smin=ss[0],smax=ss[0];2892 Real8 h1=1.e30,h2=1e-30,rx=0;2893 Real8 coef = 1./(anisomax*anisomax);2894 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;2895 int nbfield=typsols?sizeoftype[typsols[nusol]]:1;2896 2897 //only one field2898 if (nbfield == 1) {2899 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)2900 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){2901 dxdx[iv]=dxdy[iv]=dydy[iv]=0;2902 smin=Min(smin,ss[k]);2903 smax=Max(smax,ss[k]);2904 }2905 }2906 2907 //vectorial case2908 else{2909 for (iv=0,k=0;iv<nbv;iv++,k+=n ){2910 //compute v = √sum(s[i]^2)2911 double v=0;2912 for (int i=0;i<nbfield;i++){2913 v += ss[k+i]*ss[k+i];2914 }2915 v = sqrt(v);2916 smin=Min(smin,v);2917 smax=Max(smax,v);2918 }2919 }2920 Real8 sdelta=smax-smin;2921 Real8 absmax=Max(Abs(smin),Abs(smax));2922 Real8 cnorm =Rescaling ? coef2/sdelta : coef2;2923 2924 //display info2925 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);2926 2927 //skip constant field2928 if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){2929 if (verbosity>2) printf(" Solution %i is constant, skipping...\n");2930 continue;2931 }2932 2933 //pointer toward ss that is also a pointer toward s (solutions)2934 double* sf=ss;2935 2936 //loop over all the fields of the solution2937 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){2938 //ss++ so that for each iteration ss points toward the right field2939 2940 //initialize the hessian matrix2941 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0;2942 2943 //loop over the triangles2944 for (i=0;i<nbt;i++){2945 2946 //for real all triangles2947 if(triangles[i].link){2948 2949 // coor of 3 vertices2950 R2 A=triangles[i][0];2951 R2 B=triangles[i][1];2952 R2 C=triangles[i][2];2953 2954 //warning: the normal is internal and the size is the length of the edge2955 R2 nAB = Orthogonal(B-A);2956 R2 nBC = Orthogonal(C-B);2957 R2 nCA = Orthogonal(A-C);2958 //note that : nAB + nBC + nCA == 02959 2960 // number of the 3 vertices2961 iA = Number(triangles[i][0]);2962 iB = Number(triangles[i][1]);2963 iC = Number(triangles[i][2]);2964 2965 // for the test of boundary edge2966 // the 3 adj triangles2967 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);2968 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);2969 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);2970 2971 // value of the P1 fonction on 3 vertices2972 sA = ss[iA*n];2973 sB = ss[iB*n];2974 sC = ss[iC*n];2975 2976 /*The nodal functions are such that for a vertex A:2977 N_A(x,y)=alphaA x + beta_A y +gamma_A2978 N_A(A) = 1, N_A(B) = 0, N_A(C) = 02979 solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)2980 leads to:2981 N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))2982 and this gives:2983 alpha_A = (yB-yC)/(2*Area(T))2984 beta_A = (xC-xB)/(2*Area(T))2985 and therefore:2986 grad N_A = nA / detT2987 for an interpolation of a solution s:2988 grad(s) = s * sum_{i=A,B,C} grad(N_i) */2989 2990 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i];2991 2992 //from P2 on 4T to compute Hessian2993 if(Hessiantype==1){2994 int nbb=0;2995 Real8 dd = detT[i];2996 Real8 lla,llb,llc,llf;2997 Real8 taa[3][3],bb[3];2998 2999 // construction of the transpose of lin system3000 for (int j=0;j<3;j++){3001 int ie = OppositeEdge[j];3002 TriangleAdjacent ta = triangles[i].Adj(ie);3003 Triangle* tt = ta;3004 3005 //if the adjacent triangle is not a boundary triangle:3006 if (tt && tt->link){3007 Vertex &v = *ta.OppositeVertex();3008 R2 V = v;3009 Int4 iV = Number(v);3010 Real8 lA = bamg::Area2(V,B,C)/dd;3011 Real8 lB = bamg::Area2(A,V,C)/dd;3012 Real8 lC = bamg::Area2(A,B,V)/dd;3013 taa[0][j] = lB*lC;3014 taa[1][j] = lC*lA;3015 taa[2][j] = lA*lB;3016 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;3017 bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ;3018 }3019 else{3020 nbb++;3021 taa[0][j]=0;3022 taa[1][j]=0;3023 taa[2][j]=0;3024 taa[j][j]=1;3025 bb[j]=0;3026 }3027 }3028 3029 // resolution of 3x3 linear system transpose3030 Real8 det33 = det3x3(taa[0],taa[1],taa[2]);3031 Real8 cBC = det3x3(bb,taa[1],taa[2]);3032 Real8 cCA = det3x3(taa[0],bb,taa[2]);3033 Real8 cAB = det3x3(taa[0],taa[1],bb);3034 3035 if (!det33){3036 throw ErrorException(__FUNCT__,exprintf("!det33"));3037 }3038 // computation of the Hessian in the element3039 3040 // H( li*lj) = grad li grad lj + grad lj grad lj3041 // grad li = njk / detT ; with i j k =(A,B,C)3042 Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);3043 Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);3044 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x)3045 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);3046 Real8 coef = 1.0/(3*dd*det33);3047 Real8 coef2 = 2*coef;3048 Hxx *= coef2;3049 Hyy *= coef2;3050 Hxy *= coef2;3051 if(nbb==0){3052 dxdx[iA] += Hxx;3053 dydy[iA] += Hyy;3054 dxdy[iA] += Hxy;3055 3056 dxdx[iB] += Hxx;3057 dydy[iB] += Hyy;3058 dxdy[iB] += Hxy;3059 3060 dxdx[iC] += Hxx;3061 dydy[iC] += Hyy;3062 dxdy[iC] += Hxy;3063 }3064 }3065 3066 //Use Green to compute Hessian Matrix3067 else if (Hessiantype==0){3068 3069 // if edge on boundary no contribution => normal = 03070 if ( !tBC || !tBC->link ) nBC=O;3071 if ( !tCA || !tCA->link ) nCA=O;3072 if ( !tAB || !tAB->link ) nAB=O;3073 3074 // remark we forgot a 1/2 because3075 // int_{edge} w_i = 1/2 if i is in edge3076 // 0 if not3077 // if we don't take the boundary3078 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;3079 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;3080 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;3081 3082 //warning optimization (1) the division by 2 is done on the metric construction3083 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;3084 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;3085 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ;3086 3087 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;3088 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;3089 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;3090 }3091 3092 } // for real all triangles3093 }3094 3095 Int4 kk=0;3096 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){3097 if(Mmassxx[iv]>0){3098 dxdx[iv] /= 2*Mmassxx[iv];3099 // warning optimization (1) on term dxdy[iv]*ci/23100 dxdy[iv] /= 4*Mmassxx[iv];3101 dydy[iv] /= 2*Mmassxx[iv];3102 // Compute the matrix with abs(eigen value)3103 Metric M(dxdx[iv], dxdy[iv], dydy[iv]);3104 MatVVP2x2 Vp(M);3105 Vp.Abs();3106 M = Vp;3107 dxdx[iv] = M.a11;3108 dxdy[iv] = M.a21;3109 dydy[iv] = M.a22;3110 }3111 else kk++;3112 }3113 3114 // correction of second derivative3115 // by a laplacien3116 Real8* d2[3] = {dxdx, dxdy, dydy};3117 Real8* dd;3118 for (int xy = 0;xy<3;xy++) {3119 dd = d2[xy];3120 // do leat 2 iteration for boundary problem3121 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){3122 for (i=0;i<nbt;i++)3123 if(triangles[i].link){// the real triangles3124 // number of the 3 vertices3125 iA = Number(triangles[i][0]);3126 iB = Number(triangles[i][1]);3127 iC = Number(triangles[i][2]);3128 Real8 cc=3;3129 if(ijacobi==0)3130 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);3131 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;3132 }3133 for (iv=0;iv<nbv;iv++) workV[iv]=0;3134 3135 for (i=0;i<nbt;i++){3136 if(triangles[i].link){ // the real triangles3137 // number of the 3 vertices3138 iA = Number(triangles[i][0]);3139 iB = Number(triangles[i][1]);3140 iC = Number(triangles[i][2]);3141 Real8 cc = workT[i]*detT[i];3142 workV[iA] += cc;3143 workV[iB] += cc;3144 workV[iC] += cc;3145 }3146 }3147 3148 for (iv=0;iv<nbv;iv++){3149 if( ijacobi<NbJacobi || OnBoundary[iv]){3150 dd[iv] = workV[iv]/(Mmass[iv]*6);3151 }3152 }3153 }3154 }3155 3156 //constuction of the metric from the Hessian dxdx. dxdy,dydy3157 Real8 rCutOff=CutOff*absmax;// relative cut off3158 3159 //loop over the nodes3160 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){3161 3162 MetricIso Miso;3163 Real8 ci ;3164 3165 // compute norm of the solution3166 if (RelativeMetric){3167 double xx =0,*sfk=sf+k;3168 for (int ifield=0;ifield<nbfield;ifield++,sfk++){3169 xx += *sfk* *sfk;3170 }3171 xx=sqrt(xx);3172 ci=coef2/Max(xx,rCutOff);3173 }3174 else ci=cnorm;3175 3176 //initialize metric Miv with ci*H3177 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci);3178 3179 //Get eigen values and vectors of Miv3180 MatVVP2x2 Vp(Miv);3181 3182 //move eigen valuse to their absolute values3183 Vp.Abs();3184 3185 //Allpy a power if requested by user3186 if(power!=1.0) Vp.pow(power);3187 3188 //get minimum and maximum eigen values3189 h1=Min(h1,Vp.lmin());3190 h2=Max(h2,Vp.lmax());3191 3192 //modify eigen values according to hmin and hmax3193 Vp.Maxh(hmin);3194 Vp.Minh(hmax);3195 3196 //multiply eigen values by coef3197 Vp.BoundAniso2(coef);3198 3199 //rebuild Metric from Vp3200 Metric MVp(Vp);3201 3202 //Apply Metric to vertex3203 vertices[iv].m.IntersectWith(MVp);3204 3205 //info to be displayed3206 //rx = max(lmax/lmin) (anisotropy ratio)3207 rx = Max(rx,Vp.Aniso2());3208 hn1=Min(hn1,Vp.lmin());3209 hn2=Max(hn2,Vp.lmax());3210 rnx = Max(rnx,Vp.Aniso2());3211 }3212 3213 //display info3214 if (verbosity>2) {3215 3216 3217 printf(" Field %i of solution %i\n",nufield,nusol);3218 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));3219 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));3220 }3221 } // end of for all field3222 }// end for all solution3223 3224 delete [] detT;3225 delete [] Mmass;3226 delete [] dxdx;3227 delete [] dxdy;3228 delete [] dydy;3229 delete [] workT;3230 delete [] workV;3231 delete [] Mmassxx;3232 delete [] OnBoundary;3233 3234 3655 } 3235 3656 /*}}}1*/
Note:
See TracChangeset
for help on using the changeset viewer.
