![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2957
- Timestamp:
- 02/04/10 08:32:31 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 3 edited
-
Mesh2.h (modified) (5 diffs)
-
objects/Geometry.cpp (modified) (4 diffs)
-
objects/Triangles.cpp (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2955 r2957 660 660 Int4 NbOfTriangleSearchFind; 661 661 Int4 NbOfSwapTriangle; 662 char * name, *identity;663 662 Vertex * vertices; // data of vertices des sommets 664 663 … … 710 709 711 710 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR 712 Triangles(const Triangles &,const int *flag,const int *bb ); // truncature711 Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature 713 712 714 713 void SetIntCoor(const char * from =0); … … 786 785 int UnCrack(); 787 786 788 void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo787 void BuildGeometryFromMesh(BamgOpts* bamgopts); 789 788 void FillHoleInMesh() ; 790 789 int CrackMesh(); … … 792 791 void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption 793 792 void GeomToTriangles0(Int4 nbvx); // the real constructor mesh generator 794 void PreInit(Int4 ,char* =NULL);793 void PreInit(Int4); 795 794 796 795 }; … … 801 800 Int4 NbRef; // counter of ref on the this class if 0 we can delete 802 801 803 char *name;804 802 Int4 nbvx,nbtx; // nombre max de sommets , de Geometry 805 803 Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux, -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2945 r2957 26 26 NbRef =0; 27 27 quadtree=0; 28 name = new char[strlen(Gh.name)+4];29 strcpy(name,"cp:");30 strcat(name,Gh.name);31 28 vertices = nbv ? new GeometricalVertex[nbv] : NULL; 32 29 triangles = nbt ? new Triangle[nbt]:NULL; … … 60 57 if(quadtree) delete quadtree;quadtree=0; 61 58 if(curves) delete []curves;curves=0;NbOfCurves=0; 62 if(name) delete [] name;name=0;63 59 if(subdomains) delete [] subdomains;subdomains=0; 64 60 // if(ordre) delete [] ordre; … … 877 873 878 874 printf("Geometry:\n"); 879 printf(" name: %s\n",name);880 875 printf(" nbvx (maximum number of vertices) : %i\n",nbvx); 881 876 printf(" nbtx (maximum number of triangles): %i\n",nbtx); … … 905 900 void Geometry::EmptyGeometry() { 906 901 NbRef=0; 907 name =0;908 902 quadtree=0; 909 903 curves=0; -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2955 r2957 26 26 Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 27 27 28 PreInit(0 ,"none");28 PreInit(0); 29 29 ReadMesh(bamgmesh,bamgopts); 30 30 SetIntCoor(); … … 33 33 /*}}}1*/ 34 34 /*FUNCTION Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb){{{1*/ 35 Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb) : Gh(*(new Geometry())), BTh(*this) { 36 37 char cname[] = "trunc"; 35 Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb,BamgOpts* bamgopts) : Gh(*(new Geometry())), BTh(*this) { 38 36 39 37 int i,k,itadj; … … 83 81 printf(" number of New boundary edge %i\n",nbNewBedge); 84 82 Int4 inbvx =k; 85 PreInit(inbvx ,cname);83 PreInit(inbvx); 86 84 for (i=0;i<Tho.nbv;i++) 87 85 if (kk[i]>=0) … … 121 119 delete [] reft; 122 120 delete [] refv; 123 double cutoffradian = 10.0/180.0*Pi;124 ConsGeometry(cutoffradian);121 //double cutoffradian = 10.0/180.0*Pi; 122 BuildGeometryFromMesh(bamgopts); 125 123 Gh.AfterRead(); 126 124 SetIntCoor(); … … 148 146 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; 149 147 if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex; 150 //if (name) delete [] name; //TEST crash if not commented151 if (identity) delete [] identity;152 148 if (VertexOnBThVertex) delete [] VertexOnBThVertex; 153 149 if (VertexOnBThEdge) delete [] VertexOnBThEdge; … … 172 168 // do all the allocation to be sure all the pointer existe 173 169 174 char * cname = 0; 175 if (Th.name) 176 { 177 cname = new char[strlen(Th.name)+1]; 178 strcpy(cname,Th.name); 179 } 180 PreInit(nbvxx,cname);// to make the allocation 170 PreInit(nbvxx);// to make the allocation 181 171 // copy of triangles 182 172 nt=Th.nt; … … 302 292 triangles =new Triangle[nbt]; 303 293 for (i=0;i<bamgmesh->NumQuadrilaterals;i++){ 294 //divide the quad into two triangles 304 295 Triangle & t1 = triangles[2*i]; 305 296 Triangle & t2 = triangles[2*i+1]; … … 466 457 printf("WARNING: mesh present but no geometry found. Reconstructing...\n"); 467 458 /*Recreate geometry: */ 468 ConsGeometry(bamgopts->MaximalAngleOfCorner*Pi/180);459 BuildGeometryFromMesh(bamgopts); 469 460 Gh.AfterRead(); 470 461 } … … 885 876 } 886 877 /*}}}1*/ 887 /*FUNCTION Triangles::BuildMetric0 (double P2 projection){{{1*/ 888 void Triangles::BuildMetric0(BamgOpts* bamgopts){ 889 890 /*Options*/ 891 const int dim = 2; 892 double* s=NULL; 893 Int4 nbsol; 894 int verbosity; 895 896 int i,j,k,iA,iB,iC; 897 int iv; 878 /*FUNCTION Triangles::BuildGeometryFromMesh{{{1*/ 879 void Triangles::BuildGeometryFromMesh(BamgOpts* bamgopts){ 880 /*Reconstruct Geometry from Mesh*/ 881 882 /*Intermediary*/ 883 int i,j,k; 898 884 899 885 /*Recover options*/ 900 verbosity=bamgopts->verbose; 901 902 /*Get and process fields*/ 903 s=bamgopts->field; 904 nbsol=bamgopts->numfields; 905 906 //initialization of some variables 907 double* ss=(double*)s; 908 double sA,sB,sC; 909 Real8* detT = new Real8[nbt]; 910 Real8* sumareas = new Real8[nbv]; 911 Real8* alpha= new Real8[nbt*3]; 912 Real8* beta = new Real8[nbt*3]; 913 Real8* dx_elem = new Real8[nbt]; 914 Real8* dy_elem = new Real8[nbt]; 915 Real8* dx_vertex = new Real8[nbv]; 916 Real8* dy_vertex = new Real8[nbv]; 917 Real8* dxdx_elem = new Real8[nbt]; 918 Real8* dxdy_elem = new Real8[nbt]; 919 Real8* dydy_elem = new Real8[nbt]; 920 Real8* dxdx_vertex= new Real8[nbv]; 921 Real8* dxdy_vertex= new Real8[nbv]; 922 Real8* dydy_vertex= new Real8[nbv]; 923 924 //display infos 925 if(verbosity>1) { 926 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv); 927 } 928 929 //first, build the chains that will be used for the Hessian computation, as weel as the area of each element 930 int head_s[nbv]; 931 int next_p[3*nbt]; 932 int p=0; 933 //initialization 934 for(i=0;i<nbv;i++){ 935 sumareas[i]=0; 936 head_s[i]=-1; 937 } 938 for(i=0;i<nbt;i++){ 939 940 //lopp over the real triangles (no boundary elements) 941 if(triangles[i].link){ 942 943 //get current triangle t 944 const Triangle &t=triangles[i]; 945 946 // coor of 3 vertices 947 R2 A=t[0]; 948 R2 B=t[1]; 949 R2 C=t[2]; 950 951 //compute triangle determinant (2*Area) 952 Real8 dett = bamg::Area2(A,B,C); 953 detT[i]=dett; 954 955 /*The nodal functions are such that for a vertex A: 956 * N_A(x,y)=alphaA x + beta_A y +gamma_A 957 * N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 958 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 959 * leads to: 960 * N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 961 * and this gives: 962 * alpha_A = (yB-yC)/(2*Area(T))*/ 963 alpha[i*3+0]=(B.y-C.y)/dett; 964 alpha[i*3+1]=(C.y-A.y)/dett; 965 alpha[i*3+2]=(A.y-B.y)/dett; 966 beta[ i*3+0]=(C.x-B.x)/dett; 967 beta[ i*3+1]=(A.x-C.x)/dett; 968 beta[ i*3+2]=(B.x-A.x)/dett; 969 970 //compute chains 971 for(j=0;j<3;j++){ 972 k=Number(triangles[i][j]); 973 next_p[p]=head_s[k]; 974 head_s[k]=p++; 975 976 //add area to sumareas 977 sumareas[k]+=dett; 978 } 979 980 } 981 } 982 983 //for all Solutions 984 for (Int4 nusol=0;nusol<nbsol;nusol++) { 985 Real8 smin=ss[nusol],smax=ss[nusol]; 986 987 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 988 for ( iv=0,k=0; iv<nbv; iv++){ 989 smin=Min(smin,ss[iv*nbsol+nusol]); 990 smax=Max(smax,ss[iv*nbsol+nusol]); 991 } 992 Real8 sdelta=smax-smin; 993 Real8 absmax=Max(Abs(smin),Abs(smax)); 994 995 //display info 996 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g\n",nusol,smin,smax,sdelta); 997 998 //skip constant field 999 if (sdelta < 1.0e-10*Max(absmax,1e-20)){ 1000 printf(" Solution %i is constant, skipping...\n",nusol); 1001 continue; 1002 } 1003 1004 //initialize the hessian and gradient matrices 1005 for ( iv=0,k=0; iv<nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0; 1006 1007 //1: Compute gradient for each element (exact) 1008 for (i=0;i<nbt;i++){ 1009 if(triangles[i].link){ 1010 // number of the 3 vertices 1011 iA = Number(triangles[i][0]); 1012 iB = Number(triangles[i][1]); 1013 iC = Number(triangles[i][2]); 1014 1015 // value of the P1 fonction on 3 vertices 1016 sA = ss[iA*nbsol+nusol]; 1017 sB = ss[iB*nbsol+nusol]; 1018 sC = ss[iC*nbsol+nusol]; 1019 1020 //gradient = (sum alpha_i s_i, sum_i beta_i s_i) 1021 dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2]; 1022 dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2]; 1023 } 1024 } 1025 1026 //2: then compute a gradient for each vertex using a P2 projection 1027 for(i=0;i<nbv;i++){ 1028 for(p=head_s[i];p!=-1;p=next_p[p]){ 1029 //Get triangle number 1030 k=(Int4)(p/3); 1031 dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i]; 1032 dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i]; 1033 } 1034 } 1035 1036 //3: compute Hessian matrix on each element 1037 for (i=0;i<nbt;i++){ 1038 if(triangles[i].link){ 1039 // number of the 3 vertices 1040 iA = Number(triangles[i][0]); 1041 iB = Number(triangles[i][1]); 1042 iC = Number(triangles[i][2]); 1043 1044 //Hessian 1045 dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2]; 1046 dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2]; 1047 dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2]; 1048 } 1049 } 1050 1051 //4: finaly compute Hessian on each vertex using the second P2 projection 1052 for(i=0;i<nbv;i++){ 1053 for(p=head_s[i];p!=-1;p=next_p[p]){ 1054 //Get triangle number 1055 k=(Int4)(p/3); 1056 dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i]; 1057 dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i]; 1058 dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i]; 1059 } 1060 } 1061 1062 /*Compute Metric from Hessian*/ 1063 for ( iv=0;iv<nbv;iv++){ 1064 vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts); 1065 } 1066 1067 }//for all solutions 1068 1069 //clean up 1070 delete [] detT; 1071 delete [] alpha; 1072 delete [] beta; 1073 delete [] sumareas; 1074 delete [] dx_elem; 1075 delete [] dy_elem; 1076 delete [] dx_vertex; 1077 delete [] dy_vertex; 1078 delete [] dxdx_elem; 1079 delete [] dxdy_elem; 1080 delete [] dydy_elem; 1081 delete [] dxdx_vertex; 1082 delete [] dxdy_vertex; 1083 delete [] dydy_vertex; 1084 } 1085 /*}}}1*/ 1086 /*FUNCTION Triangles::BuildMetric1 (Green formula){{{1*/ 1087 void Triangles::BuildMetric1(BamgOpts* bamgopts){ 1088 1089 /*Options*/ 1090 const int dim = 2; 1091 double* s=NULL; 1092 Int4 nbsol; 1093 int NbJacobi; 1094 int verbosity; 1095 1096 /*Recover options*/ 1097 verbosity=bamgopts->verbose; 1098 NbJacobi=bamgopts->nbjacobi; 1099 1100 /*Get and process fields*/ 1101 s=bamgopts->field; 1102 nbsol=bamgopts->numfields; 1103 1104 //initialization of some variables 1105 Int4 i,k,iA,iB,iC,iv; 1106 R2 O(0,0); 1107 double* ss=(double*)s; 1108 double sA,sB,sC; 1109 Real8* detT = new Real8[nbt]; 1110 Real8* Mmass= new Real8[nbv]; 1111 Real8* Mmassxx= new Real8[nbv]; 1112 Real8* dxdx= new Real8[nbv]; 1113 Real8* dxdy= new Real8[nbv]; 1114 Real8* dydy= new Real8[nbv]; 1115 Real8* workT= new Real8[nbt]; 1116 Real8* workV= new Real8[nbv]; 1117 int* OnBoundary = new int[nbv]; 1118 1119 //display infos 1120 if(verbosity>1) { 1121 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv); 1122 } 1123 1124 //initialize Mmass, OnBoundary and Massxx by zero 1125 for (iv=0;iv<nbv;iv++){ 1126 Mmass[iv]=0; 1127 OnBoundary[iv]=0; 1128 Mmassxx[iv]=0; 1129 } 1130 1131 //Build detT Mmas Mmassxx workT and OnBoundary 1132 for (i=0;i<nbt;i++){ 1133 1134 //lopp over the real triangles (no boundary elements) 1135 if(triangles[i].link){ 1136 1137 //get current triangle t 1138 const Triangle &t=triangles[i]; 1139 1140 // coor of 3 vertices 1141 R2 A=t[0]; 1142 R2 B=t[1]; 1143 R2 C=t[2]; 1144 1145 // number of the 3 vertices 1146 iA = Number(t[0]); 1147 iB = Number(t[1]); 1148 iC = Number(t[2]); 1149 1150 //compute triangle determinant (2*Area) 1151 Real8 dett = bamg::Area2(A,B,C); 1152 detT[i]=dett; 1153 dett /= 6; 1154 1155 // construction of OnBoundary (flag=1 if on boundary, else 0) 1156 int nbb=0; 1157 for(int j=0;j<3;j++){ 1158 //get adjacent triangle 1159 Triangle *ta=t.Adj(j); 1160 //if there is no adjacent triangle, the edge of the triangle t is on boundary 1161 if ( !ta || !ta->link){ 1162 //mark the two vertices of the edge as OnBoundary 1163 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1; 1164 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1; 1165 nbb++; 1166 } 1167 } 1168 1169 //number of vertices on boundary for current triangle t 1170 workT[i] = nbb; 1171 1172 //Build Mmass Mmass[i] = Mmass[i] + Area/3 1173 Mmass[iA] += dett; 1174 Mmass[iB] += dett; 1175 Mmass[iC] += dett; 1176 1177 //Build Massxx = Mmass 1178 Mmassxx[iA] += dett; 1179 Mmassxx[iB] += dett; 1180 Mmassxx[iC] += dett; 1181 } 1182 1183 //else: the triangle is a boundary triangle -> workT=-1 1184 else workT[i]=-1; 1185 } 1186 1187 //for all Solution 1188 for (Int4 nusol=0;nusol<nbsol;nusol++) { 1189 1190 Real8 smin=ss[nusol],smax=ss[nusol]; 1191 Real8 h1=1.e30,h2=1e-30,rx=0; 1192 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1193 1194 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 1195 for ( iv=0,k=0; iv<nbv; iv++ ){ 1196 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1197 smin=Min(smin,ss[iv*nbsol+nusol]); 1198 smax=Max(smax,ss[iv*nbsol+nusol]); 1199 } 1200 Real8 sdelta=smax-smin; 1201 Real8 absmax=Max(Abs(smin),Abs(smax)); 1202 1203 //display info 1204 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbsol); 1205 1206 //skip constant field 1207 if (sdelta < 1.0e-10*Max(absmax,1e-20) ){ 1208 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 1209 continue; 1210 } 1211 1212 //pointer toward ss that is also a pointer toward s (solutions) 1213 double* sf=ss; 1214 1215 //initialize the hessian matrix 1216 for ( iv=0,k=0; iv<nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1217 1218 //loop over the triangles 1219 for (i=0;i<nbt;i++){ 1220 1221 //for real all triangles 1222 if(triangles[i].link){ 1223 1224 // coor of 3 vertices 1225 R2 A=triangles[i][0]; 1226 R2 B=triangles[i][1]; 1227 R2 C=triangles[i][2]; 1228 1229 //warning: the normal is internal and the size is the length of the edge 1230 R2 nAB = Orthogonal(B-A); 1231 R2 nBC = Orthogonal(C-B); 1232 R2 nCA = Orthogonal(A-C); 1233 //note that : nAB + nBC + nCA == 0 1234 1235 // number of the 3 vertices 1236 iA = Number(triangles[i][0]); 1237 iB = Number(triangles[i][1]); 1238 iC = Number(triangles[i][2]); 1239 1240 // for the test of boundary edge 1241 // the 3 adj triangles 1242 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 1243 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 1244 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 1245 1246 // value of the P1 fonction on 3 vertices 1247 sA = ss[iA*nbsol+nusol]; 1248 sB = ss[iB*nbsol+nusol]; 1249 sC = ss[iC*nbsol+nusol]; 1250 1251 /*The nodal functions are such that for a vertex A: 1252 N_A(x,y)=alphaA x + beta_A y +gamma_A 1253 N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 1254 solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 1255 leads to: 1256 N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 1257 and this gives: 1258 alpha_A = (yB-yC)/(2*Area(T)) 1259 beta_A = (xC-xB)/(2*Area(T)) 1260 and therefore: 1261 grad N_A = nA / detT 1262 for an interpolation of a solution s: 1263 grad(s) = s * sum_{i=A,B,C} grad(N_i) */ 1264 1265 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i]; 1266 1267 //Use Green to compute Hessian Matrix 1268 1269 // if edge on boundary no contribution => normal = 0 1270 if ( !tBC || !tBC->link ) nBC=O; 1271 if ( !tCA || !tCA->link ) nCA=O; 1272 if ( !tAB || !tAB->link ) nAB=O; 1273 1274 // remark we forgot a 1/2 because 1275 // int_{edge} w_i = 1/2 if i is in edge 1276 // 0 if not 1277 // if we don't take the boundary 1278 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; 1279 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; 1280 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; 1281 1282 //warning optimization (1) the division by 2 is done on the metric construction 1283 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; 1284 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; 1285 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 1286 1287 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; 1288 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; 1289 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; 1290 1291 } // for real all triangles 1292 } 1293 1294 Int4 kk=0; 1295 for ( iv=0,k=0 ; iv<nbv; iv++){ 1296 if(Mmassxx[iv]>0){ 1297 dxdx[iv] /= 2*Mmassxx[iv]; 1298 // warning optimization (1) on term dxdy[iv]*ci/2 1299 dxdy[iv] /= 4*Mmassxx[iv]; 1300 dydy[iv] /= 2*Mmassxx[iv]; 1301 // Compute the matrix with abs(eigen value) 1302 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 1303 MatVVP2x2 Vp(M); 1304 Vp.Abs(); 1305 M = Vp; 1306 dxdx[iv] = M.a11; 1307 dxdy[iv] = M.a21; 1308 dydy[iv] = M.a22; 1309 } 1310 else kk++; 1311 } 1312 1313 // correction of second derivative 1314 // by a laplacien 1315 Real8* d2[3] = {dxdx, dxdy, dydy}; 1316 Real8* dd; 1317 for (int xy = 0;xy<3;xy++) { 1318 dd = d2[xy]; 1319 // do leat 2 iteration for boundary problem 1320 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ 1321 for (i=0;i<nbt;i++) 1322 if(triangles[i].link){// the real triangles 1323 // number of the 3 vertices 1324 iA = Number(triangles[i][0]); 1325 iB = Number(triangles[i][1]); 1326 iC = Number(triangles[i][2]); 1327 Real8 cc=3; 1328 if(ijacobi==0) 1329 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 1330 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 1331 } 1332 for (iv=0;iv<nbv;iv++) workV[iv]=0; 1333 1334 for (i=0;i<nbt;i++){ 1335 if(triangles[i].link){ // the real triangles 1336 // number of the 3 vertices 1337 iA = Number(triangles[i][0]); 1338 iB = Number(triangles[i][1]); 1339 iC = Number(triangles[i][2]); 1340 Real8 cc = workT[i]*detT[i]; 1341 workV[iA] += cc; 1342 workV[iB] += cc; 1343 workV[iC] += cc; 1344 } 1345 } 1346 1347 for (iv=0;iv<nbv;iv++){ 1348 if( ijacobi<NbJacobi || OnBoundary[iv]){ 1349 dd[iv] = workV[iv]/(Mmass[iv]*6); 1350 } 1351 } 1352 } 1353 } 1354 1355 /*Compute Metric from Hessian*/ 1356 for ( iv=0;iv<nbv;iv++){ 1357 vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts); 1358 } 1359 1360 }// end for all solution 1361 1362 delete [] detT; 1363 delete [] Mmass; 1364 delete [] dxdx; 1365 delete [] dxdy; 1366 delete [] dydy; 1367 delete [] workT; 1368 delete [] workV; 1369 delete [] Mmassxx; 1370 delete [] OnBoundary; 1371 1372 } 1373 /*}}}1*/ 1374 /*FUNCTION Triangles::ConsGeometry{{{1*/ 1375 void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges){ 1376 // if equiedges existe taille nbe 1377 // equiedges[i]/2 == i original 1378 // equiedges[i]/2 = j => equivalence entre i et j => meme maillage 1379 // equiedges[i]%2 : 0 meme sens , 1 pas meme sens 1380 // 1381 // -------------------------- 1382 long int verbosity=0; 886 int verbosity=bamgopts->verbose; 887 double cutoffradian=bamgopts->MaximalAngleOfCorner*Pi/180; 888 CurrentTh=this; 889 890 //display info 1383 891 if (verbosity>1) printf(" construction of the geometry from the 2d mesh\n"); 892 893 //check that the mesh is not empty 1384 894 if (nbt<=0 || nbv <=0 ) { 1385 throw ErrorException(__FUNCT__,exprintf("nbt or nbv is negative")); 1386 } 1387 1388 // construction of the edges 1389 // Triangles * OldCurrentTh =CurrentTh; 1390 CurrentTh=this; 1391 // Int4 NbTold = nbt; 1392 // generation of the integer coor 1393 // generation of the adjacence of the triangles 1394 if (cutoffradian>=0) 1395 Gh.MaximalAngleOfCorner = cutoffradian; 1396 SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv); 1397 Int4 * st = new Int4[nbt*3]; 1398 Int4 i,k; 1399 int j; 1400 if (Gh.name) delete Gh.name; 1401 Gh.name = new char [ name ? strlen(name) + 15 : 50 ]; 1402 Gh.name[0]=0; 1403 strcat(Gh.name,"cons from: "); 1404 if (name) strcat(Gh.name,name); 1405 else strcat(Gh.name," a mesh with no name"); 895 throw ErrorException(__FUNCT__,exprintf("nbt or nbv is negative (Mesh empty?)")); 896 } 897 898 //Gh is the geometry of the mesh (this), initialize MaximalAngleOfCorner 899 if (cutoffradian>=0) Gh.MaximalAngleOfCorner = cutoffradian; 900 901 //Construction of the edges 902 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 903 Int4* st = new Int4[nbt*3]; 1406 904 for (i=0;i<nbt*3;i++) 1407 905 st[i]=-1; … … 1445 943 1446 944 if(verbosity>5) { 1447 printf(" info on Mesh %s:\n",name);945 printf(" info on Meshs:\n"); 1448 946 printf(" - number of vertices = %i \n",nbv); 1449 947 printf(" - number of triangles = %i \n",nbt); … … 1697 1195 if(requis) kreq++; 1698 1196 edges[i].onGeometry = Gh.edges + i; 1699 if(equiedges && i < nbeold ) {1700 int j=equiedges[i]/2;1701 int sens=equiedges[i]%2;1702 if(i!=j && equiedges[i]>=0) {1703 if( sens==0)1704 Gh.edges[i].SetEqui();1705 else1706 Gh.edges[i].SetReverseEqui();1707 Gh.edges[i].link= & Gh.edges[j];1708 }1709 1710 }1711 1197 if(requis) { // correction fevr 2009 JYU ... 1712 1198 Gh.edges[i].v[0]->SetRequired(); … … 1766 1252 triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j)); 1767 1253 1768 } 1254 } 1255 /*}}}1*/ 1256 /*FUNCTION Triangles::BuildMetric0 (double P2 projection){{{1*/ 1257 void Triangles::BuildMetric0(BamgOpts* bamgopts){ 1258 1259 /*Options*/ 1260 const int dim = 2; 1261 double* s=NULL; 1262 Int4 nbsol; 1263 int verbosity; 1264 1265 int i,j,k,iA,iB,iC; 1266 int iv; 1267 1268 /*Recover options*/ 1269 verbosity=bamgopts->verbose; 1270 1271 /*Get and process fields*/ 1272 s=bamgopts->field; 1273 nbsol=bamgopts->numfields; 1274 1275 //initialization of some variables 1276 double* ss=(double*)s; 1277 double sA,sB,sC; 1278 Real8* detT = new Real8[nbt]; 1279 Real8* sumareas = new Real8[nbv]; 1280 Real8* alpha= new Real8[nbt*3]; 1281 Real8* beta = new Real8[nbt*3]; 1282 Real8* dx_elem = new Real8[nbt]; 1283 Real8* dy_elem = new Real8[nbt]; 1284 Real8* dx_vertex = new Real8[nbv]; 1285 Real8* dy_vertex = new Real8[nbv]; 1286 Real8* dxdx_elem = new Real8[nbt]; 1287 Real8* dxdy_elem = new Real8[nbt]; 1288 Real8* dydy_elem = new Real8[nbt]; 1289 Real8* dxdx_vertex= new Real8[nbv]; 1290 Real8* dxdy_vertex= new Real8[nbv]; 1291 Real8* dydy_vertex= new Real8[nbv]; 1292 1293 //display infos 1294 if(verbosity>1) { 1295 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv); 1296 } 1297 1298 //first, build the chains that will be used for the Hessian computation, as weel as the area of each element 1299 int head_s[nbv]; 1300 int next_p[3*nbt]; 1301 int p=0; 1302 //initialization 1303 for(i=0;i<nbv;i++){ 1304 sumareas[i]=0; 1305 head_s[i]=-1; 1306 } 1307 for(i=0;i<nbt;i++){ 1308 1309 //lopp over the real triangles (no boundary elements) 1310 if(triangles[i].link){ 1311 1312 //get current triangle t 1313 const Triangle &t=triangles[i]; 1314 1315 // coor of 3 vertices 1316 R2 A=t[0]; 1317 R2 B=t[1]; 1318 R2 C=t[2]; 1319 1320 //compute triangle determinant (2*Area) 1321 Real8 dett = bamg::Area2(A,B,C); 1322 detT[i]=dett; 1323 1324 /*The nodal functions are such that for a vertex A: 1325 * N_A(x,y)=alphaA x + beta_A y +gamma_A 1326 * N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 1327 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 1328 * leads to: 1329 * N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 1330 * and this gives: 1331 * alpha_A = (yB-yC)/(2*Area(T))*/ 1332 alpha[i*3+0]=(B.y-C.y)/dett; 1333 alpha[i*3+1]=(C.y-A.y)/dett; 1334 alpha[i*3+2]=(A.y-B.y)/dett; 1335 beta[ i*3+0]=(C.x-B.x)/dett; 1336 beta[ i*3+1]=(A.x-C.x)/dett; 1337 beta[ i*3+2]=(B.x-A.x)/dett; 1338 1339 //compute chains 1340 for(j=0;j<3;j++){ 1341 k=Number(triangles[i][j]); 1342 next_p[p]=head_s[k]; 1343 head_s[k]=p++; 1344 1345 //add area to sumareas 1346 sumareas[k]+=dett; 1347 } 1348 1349 } 1350 } 1351 1352 //for all Solutions 1353 for (Int4 nusol=0;nusol<nbsol;nusol++) { 1354 Real8 smin=ss[nusol],smax=ss[nusol]; 1355 1356 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 1357 for ( iv=0,k=0; iv<nbv; iv++){ 1358 smin=Min(smin,ss[iv*nbsol+nusol]); 1359 smax=Max(smax,ss[iv*nbsol+nusol]); 1360 } 1361 Real8 sdelta=smax-smin; 1362 Real8 absmax=Max(Abs(smin),Abs(smax)); 1363 1364 //display info 1365 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g\n",nusol,smin,smax,sdelta); 1366 1367 //skip constant field 1368 if (sdelta < 1.0e-10*Max(absmax,1e-20)){ 1369 printf(" Solution %i is constant, skipping...\n",nusol); 1370 continue; 1371 } 1372 1373 //initialize the hessian and gradient matrices 1374 for ( iv=0,k=0; iv<nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0; 1375 1376 //1: Compute gradient for each element (exact) 1377 for (i=0;i<nbt;i++){ 1378 if(triangles[i].link){ 1379 // number of the 3 vertices 1380 iA = Number(triangles[i][0]); 1381 iB = Number(triangles[i][1]); 1382 iC = Number(triangles[i][2]); 1383 1384 // value of the P1 fonction on 3 vertices 1385 sA = ss[iA*nbsol+nusol]; 1386 sB = ss[iB*nbsol+nusol]; 1387 sC = ss[iC*nbsol+nusol]; 1388 1389 //gradient = (sum alpha_i s_i, sum_i beta_i s_i) 1390 dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2]; 1391 dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2]; 1392 } 1393 } 1394 1395 //2: then compute a gradient for each vertex using a P2 projection 1396 for(i=0;i<nbv;i++){ 1397 for(p=head_s[i];p!=-1;p=next_p[p]){ 1398 //Get triangle number 1399 k=(Int4)(p/3); 1400 dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i]; 1401 dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i]; 1402 } 1403 } 1404 1405 //3: compute Hessian matrix on each element 1406 for (i=0;i<nbt;i++){ 1407 if(triangles[i].link){ 1408 // number of the 3 vertices 1409 iA = Number(triangles[i][0]); 1410 iB = Number(triangles[i][1]); 1411 iC = Number(triangles[i][2]); 1412 1413 //Hessian 1414 dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2]; 1415 dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2]; 1416 dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2]; 1417 } 1418 } 1419 1420 //4: finaly compute Hessian on each vertex using the second P2 projection 1421 for(i=0;i<nbv;i++){ 1422 for(p=head_s[i];p!=-1;p=next_p[p]){ 1423 //Get triangle number 1424 k=(Int4)(p/3); 1425 dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i]; 1426 dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i]; 1427 dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i]; 1428 } 1429 } 1430 1431 /*Compute Metric from Hessian*/ 1432 for ( iv=0;iv<nbv;iv++){ 1433 vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts); 1434 } 1435 1436 }//for all solutions 1437 1438 //clean up 1439 delete [] detT; 1440 delete [] alpha; 1441 delete [] beta; 1442 delete [] sumareas; 1443 delete [] dx_elem; 1444 delete [] dy_elem; 1445 delete [] dx_vertex; 1446 delete [] dy_vertex; 1447 delete [] dxdx_elem; 1448 delete [] dxdy_elem; 1449 delete [] dydy_elem; 1450 delete [] dxdx_vertex; 1451 delete [] dxdy_vertex; 1452 delete [] dydy_vertex; 1453 } 1454 /*}}}1*/ 1455 /*FUNCTION Triangles::BuildMetric1 (Green formula){{{1*/ 1456 void Triangles::BuildMetric1(BamgOpts* bamgopts){ 1457 1458 /*Options*/ 1459 const int dim = 2; 1460 double* s=NULL; 1461 Int4 nbsol; 1462 int NbJacobi; 1463 int verbosity; 1464 1465 /*Recover options*/ 1466 verbosity=bamgopts->verbose; 1467 NbJacobi=bamgopts->nbjacobi; 1468 1469 /*Get and process fields*/ 1470 s=bamgopts->field; 1471 nbsol=bamgopts->numfields; 1472 1473 //initialization of some variables 1474 Int4 i,k,iA,iB,iC,iv; 1475 R2 O(0,0); 1476 double* ss=(double*)s; 1477 double sA,sB,sC; 1478 Real8* detT = new Real8[nbt]; 1479 Real8* Mmass= new Real8[nbv]; 1480 Real8* Mmassxx= new Real8[nbv]; 1481 Real8* dxdx= new Real8[nbv]; 1482 Real8* dxdy= new Real8[nbv]; 1483 Real8* dydy= new Real8[nbv]; 1484 Real8* workT= new Real8[nbt]; 1485 Real8* workV= new Real8[nbv]; 1486 int* OnBoundary = new int[nbv]; 1487 1488 //display infos 1489 if(verbosity>1) { 1490 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv); 1491 } 1492 1493 //initialize Mmass, OnBoundary and Massxx by zero 1494 for (iv=0;iv<nbv;iv++){ 1495 Mmass[iv]=0; 1496 OnBoundary[iv]=0; 1497 Mmassxx[iv]=0; 1498 } 1499 1500 //Build detT Mmas Mmassxx workT and OnBoundary 1501 for (i=0;i<nbt;i++){ 1502 1503 //lopp over the real triangles (no boundary elements) 1504 if(triangles[i].link){ 1505 1506 //get current triangle t 1507 const Triangle &t=triangles[i]; 1508 1509 // coor of 3 vertices 1510 R2 A=t[0]; 1511 R2 B=t[1]; 1512 R2 C=t[2]; 1513 1514 // number of the 3 vertices 1515 iA = Number(t[0]); 1516 iB = Number(t[1]); 1517 iC = Number(t[2]); 1518 1519 //compute triangle determinant (2*Area) 1520 Real8 dett = bamg::Area2(A,B,C); 1521 detT[i]=dett; 1522 dett /= 6; 1523 1524 // construction of OnBoundary (flag=1 if on boundary, else 0) 1525 int nbb=0; 1526 for(int j=0;j<3;j++){ 1527 //get adjacent triangle 1528 Triangle *ta=t.Adj(j); 1529 //if there is no adjacent triangle, the edge of the triangle t is on boundary 1530 if ( !ta || !ta->link){ 1531 //mark the two vertices of the edge as OnBoundary 1532 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1; 1533 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1; 1534 nbb++; 1535 } 1536 } 1537 1538 //number of vertices on boundary for current triangle t 1539 workT[i] = nbb; 1540 1541 //Build Mmass Mmass[i] = Mmass[i] + Area/3 1542 Mmass[iA] += dett; 1543 Mmass[iB] += dett; 1544 Mmass[iC] += dett; 1545 1546 //Build Massxx = Mmass 1547 Mmassxx[iA] += dett; 1548 Mmassxx[iB] += dett; 1549 Mmassxx[iC] += dett; 1550 } 1551 1552 //else: the triangle is a boundary triangle -> workT=-1 1553 else workT[i]=-1; 1554 } 1555 1556 //for all Solution 1557 for (Int4 nusol=0;nusol<nbsol;nusol++) { 1558 1559 Real8 smin=ss[nusol],smax=ss[nusol]; 1560 Real8 h1=1.e30,h2=1e-30,rx=0; 1561 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1562 1563 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 1564 for ( iv=0,k=0; iv<nbv; iv++ ){ 1565 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1566 smin=Min(smin,ss[iv*nbsol+nusol]); 1567 smax=Max(smax,ss[iv*nbsol+nusol]); 1568 } 1569 Real8 sdelta=smax-smin; 1570 Real8 absmax=Max(Abs(smin),Abs(smax)); 1571 1572 //display info 1573 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbsol); 1574 1575 //skip constant field 1576 if (sdelta < 1.0e-10*Max(absmax,1e-20) ){ 1577 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 1578 continue; 1579 } 1580 1581 //pointer toward ss that is also a pointer toward s (solutions) 1582 double* sf=ss; 1583 1584 //initialize the hessian matrix 1585 for ( iv=0,k=0; iv<nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0; 1586 1587 //loop over the triangles 1588 for (i=0;i<nbt;i++){ 1589 1590 //for real all triangles 1591 if(triangles[i].link){ 1592 1593 // coor of 3 vertices 1594 R2 A=triangles[i][0]; 1595 R2 B=triangles[i][1]; 1596 R2 C=triangles[i][2]; 1597 1598 //warning: the normal is internal and the size is the length of the edge 1599 R2 nAB = Orthogonal(B-A); 1600 R2 nBC = Orthogonal(C-B); 1601 R2 nCA = Orthogonal(A-C); 1602 //note that : nAB + nBC + nCA == 0 1603 1604 // number of the 3 vertices 1605 iA = Number(triangles[i][0]); 1606 iB = Number(triangles[i][1]); 1607 iC = Number(triangles[i][2]); 1608 1609 // for the test of boundary edge 1610 // the 3 adj triangles 1611 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 1612 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 1613 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 1614 1615 // value of the P1 fonction on 3 vertices 1616 sA = ss[iA*nbsol+nusol]; 1617 sB = ss[iB*nbsol+nusol]; 1618 sC = ss[iC*nbsol+nusol]; 1619 1620 /*The nodal functions are such that for a vertex A: 1621 N_A(x,y)=alphaA x + beta_A y +gamma_A 1622 N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 1623 solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 1624 leads to: 1625 N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 1626 and this gives: 1627 alpha_A = (yB-yC)/(2*Area(T)) 1628 beta_A = (xC-xB)/(2*Area(T)) 1629 and therefore: 1630 grad N_A = nA / detT 1631 for an interpolation of a solution s: 1632 grad(s) = s * sum_{i=A,B,C} grad(N_i) */ 1633 1634 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i]; 1635 1636 //Use Green to compute Hessian Matrix 1637 1638 // if edge on boundary no contribution => normal = 0 1639 if ( !tBC || !tBC->link ) nBC=O; 1640 if ( !tCA || !tCA->link ) nCA=O; 1641 if ( !tAB || !tAB->link ) nAB=O; 1642 1643 // remark we forgot a 1/2 because 1644 // int_{edge} w_i = 1/2 if i is in edge 1645 // 0 if not 1646 // if we don't take the boundary 1647 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; 1648 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; 1649 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; 1650 1651 //warning optimization (1) the division by 2 is done on the metric construction 1652 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; 1653 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; 1654 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 1655 1656 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; 1657 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; 1658 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; 1659 1660 } // for real all triangles 1661 } 1662 1663 Int4 kk=0; 1664 for ( iv=0,k=0 ; iv<nbv; iv++){ 1665 if(Mmassxx[iv]>0){ 1666 dxdx[iv] /= 2*Mmassxx[iv]; 1667 // warning optimization (1) on term dxdy[iv]*ci/2 1668 dxdy[iv] /= 4*Mmassxx[iv]; 1669 dydy[iv] /= 2*Mmassxx[iv]; 1670 // Compute the matrix with abs(eigen value) 1671 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 1672 MatVVP2x2 Vp(M); 1673 Vp.Abs(); 1674 M = Vp; 1675 dxdx[iv] = M.a11; 1676 dxdy[iv] = M.a21; 1677 dydy[iv] = M.a22; 1678 } 1679 else kk++; 1680 } 1681 1682 // correction of second derivative 1683 // by a laplacien 1684 Real8* d2[3] = {dxdx, dxdy, dydy}; 1685 Real8* dd; 1686 for (int xy = 0;xy<3;xy++) { 1687 dd = d2[xy]; 1688 // do leat 2 iteration for boundary problem 1689 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ 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=3; 1697 if(ijacobi==0) 1698 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 1699 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 1700 } 1701 for (iv=0;iv<nbv;iv++) workV[iv]=0; 1702 1703 for (i=0;i<nbt;i++){ 1704 if(triangles[i].link){ // the real triangles 1705 // number of the 3 vertices 1706 iA = Number(triangles[i][0]); 1707 iB = Number(triangles[i][1]); 1708 iC = Number(triangles[i][2]); 1709 Real8 cc = workT[i]*detT[i]; 1710 workV[iA] += cc; 1711 workV[iB] += cc; 1712 workV[iC] += cc; 1713 } 1714 } 1715 1716 for (iv=0;iv<nbv;iv++){ 1717 if( ijacobi<NbJacobi || OnBoundary[iv]){ 1718 dd[iv] = workV[iv]/(Mmass[iv]*6); 1719 } 1720 } 1721 } 1722 } 1723 1724 /*Compute Metric from Hessian*/ 1725 for ( iv=0;iv<nbv;iv++){ 1726 vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts); 1727 } 1728 1729 }// end for all solution 1730 1731 delete [] detT; 1732 delete [] Mmass; 1733 delete [] dxdx; 1734 delete [] dxdy; 1735 delete [] dydy; 1736 delete [] workT; 1737 delete [] workV; 1738 delete [] Mmassxx; 1739 delete [] OnBoundary; 1740 1741 } 1769 1742 /*}}}1*/ 1770 1743 /*FUNCTION Triangles::ConsRefTriangle{{{1*/ … … 2067 2040 } 2068 2041 if(verbosity>5) { 2069 printf(" info on Mesh %s:\n",name);2042 printf(" info on Mesh\n"); 2070 2043 printf(" - number of vertices = %i \n",nbv); 2071 2044 printf(" - number of triangles = %i \n",nbt); … … 3851 3824 /*}}}1*/ 3852 3825 /*FUNCTION Triangles::PreInit{{{1*/ 3853 void Triangles::PreInit(Int4 inbvx ,char *fname) {3826 void Triangles::PreInit(Int4 inbvx) { 3854 3827 long int verbosity=0; 3855 3828 … … 3857 3830 NbRef=0; 3858 3831 // allocGeometry=0; 3859 identity=0;3860 3832 NbOfTriangleSearchFind =0; 3861 3833 NbOfSwapTriangle =0; … … 3876 3848 CrackedEdges =0; 3877 3849 nbe = 0; 3878 name = fname ;3879 3850 3880 3851 if (inbvx) {
Note:
See TracChangeset
for help on using the changeset viewer.
