Changeset 2940


Ignore:
Timestamp:
02/01/10 11:03:26 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added multifield adaptation

Location:
issm/trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Bamgx/Mesh2.h

    r2939 r2940  
    138138                Real8  Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,Real8 =1);
    139139                void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,BamgOpts* bamgopts);
     140                void Echo();
    140141                int ref() const { return ReferenceNumber;}
    141142
    … …  
    304305
    305306                private: // les arete sont opposes a un sommet
    306                 Vertex * ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
    307                 Triangle * at [3]; // nu triangle adjacent 
    308                 Int1  aa[3];  // les nu des arete dans le triangles (mod 4)
     307                Vertex* ns [3];    // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
     308                Triangle* at [3];  // 3 adjacent triangles (nu)
     309                Int1  aa[3];       // les nu des arete dans le triangles (mod 4)
    309310                public:
    310311                Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
    … …  
    313314                        Int4 color;
    314315                };
     316                void Echo();
    315317                void SetDet() {
    316318                        if(ns[0] && ns[1] && ns[2])    det = bamg::det(*ns[0],*ns[1],*ns[2]);
    … …  
    323325                TriangleAdjacent FindBoundaryEdge(int ) const;
    324326
    325                 void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu)
    326                   {
     327                void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
    327328                        if (link  >=tb && link  <te) link  = tb + renu[link -tb];
    328329                        if (at[0] >=tb && at[0] <te) at[0] = tb + renu[at[0]-tb];
    329330                        if (at[1] >=tb && at[1] <te) at[1] = tb + renu[at[1]-tb];
    330331                        if (at[2] >=tb && at[2] <te) at[2] = tb + renu[at[2]-tb];   
    331                   }
    332                 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu)
    333                   {
     332                }
     333                void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
    334334                        if (ns[0] >=vb && ns[0] <ve) ns[0] = vb + renu[ns[0]-vb];
    335335                        if (ns[1] >=vb && ns[1] <ve) ns[1] = vb + renu[ns[1]-vb];
    336336                        if (ns[2] >=vb && ns[2] <ve) ns[2] = vb + renu[ns[2]-vb];   
    337                   }
    338 
     337                }
    339338
    340339                const Vertex & operator[](int i) const {return *ns[i];};
    … …  
    353352
    354353                inline Real4 qualite() ;
    355 
    356354
    357355                void SetAdjAdj(Int1 a)
    … …  
    390388                        register Triangle * t = at[a];
    391389                        if(t) t->aa[aa[a] & 3] |=16;
    392                         aa[a] |= 16;}
    393                         void SetCracked(int a){
    394                                 register Triangle * t = at[a];
    395                                 if(t) t->aa[aa[a] & 3] |=32;
    396                                 aa[a] |= 32;}
    397 
    398                                 double   QualityQuad(int a,int option=1) const;
    399                                 Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
    400 
    401                                 void SetLocked(int a){
    402                                         register Triangle * t = at[a];
    403                                         t->aa[aa[a] & 3] |=4;
    404                                         aa[a] |= 4;}
    405 
    406                                         void SetMarkUnSwap(int a){
    407                                                 register Triangle * t = at[a];
    408                                                 t->aa[aa[a] & 3] |=8;
    409                                                 aa[a] |=8 ;}
    410 
    411 
    412                                                 void SetUnMarkUnSwap(int a){
    413                                                         register Triangle * t = at[a];
    414                                                         t->aa[aa[a] & 3] &=55; // 23 + 32
    415                                                         aa[a] &=55 ;}
     390                        aa[a] |= 16;
     391                }
     392                void SetCracked(int a){
     393                        register Triangle * t = at[a];
     394                        if(t) t->aa[aa[a] & 3] |=32;
     395                        aa[a] |= 32;
     396                }
     397
     398                double   QualityQuad(int a,int option=1) const;
     399                Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
     400
     401                void SetLocked(int a){
     402                        register Triangle * t = at[a];
     403                        t->aa[aa[a] & 3] |=4;
     404                        aa[a] |= 4;
     405                }
     406
     407                void SetMarkUnSwap(int a){
     408                        register Triangle * t = at[a];
     409                        t->aa[aa[a] & 3] |=8;
     410                        aa[a] |=8 ;
     411                }
     412
     413
     414                void SetUnMarkUnSwap(int a){
     415                        register Triangle * t = at[a];
     416                        t->aa[aa[a] & 3] &=55; // 23 + 32
     417                        aa[a] &=55 ;
     418                }
    416419
    417420        };  // end of Triangle class
    … …  
    799802                                                 void BuildMetric0(BamgOpts* bamgopts);
    800803                                                 void BuildMetric1(BamgOpts* bamgopts);
    801                                                  void BuildMetric2(BamgOpts* bamgopts);
    802804                                                 void IntersectGeomMetric(const Real8 err,const int iso);
    803805

  • issm/trunk/src/c/Bamgx/objects/Triangle.cpp

    r2865 r2940  
    4444          }
    4545        /*}}}1*/
     46        /*FUNCTION Triangle::Echo {{{1*/
     47
     48        void Triangle::Echo(void){
     49
     50                int i;
     51
     52                printf("Triangle:\n");
     53                printf("   ns pointer towards three vertices\n");
     54                printf("      ns[0] ns[1] ns[2] = %p %p %p\n",ns[0],ns[1],ns[2]);
     55                printf("   at pointer towards three adjacent triangles\n");
     56                printf("      at[0] at[1] at[2] = %p %p %p\n",at[0],at[1],at[2]);
     57                printf("   det (integer triangle determinant) = %i\n",det);
     58                if (link){
     59                        printf("   link (pointer toward duplicate triangle)= %p\n",link);
     60                }
     61                else{
     62                        printf("   color = %i\n",color);
     63                }
     64
     65                printf("\nThree vertices:\n");
     66                for(i=0;i<3;i++){
     67                        if (ns[i]){
     68                                ns[i]->Echo();
     69                        }
     70                        else{
     71                                printf("   vertex %i does not exist\n",i+1);
     72                        }
     73                }
     74
     75                return;
     76        }
     77        /*}}}*/
    4678
    4779}

  • issm/trunk/src/c/Bamgx/objects/Triangles.cpp

    r2939 r2940  
    891891        }
    892892        /*}}}1*/
    893         /*FUNCTION Triangles::BuildMetric0 (Green formula){{{1*/
     893        /*FUNCTION Triangles::BuildMetric0 (double P2 projection){{{1*/
    894894        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.
     895
     896                /*Options*/
     897                const int dim = 2;
     898                double* s=NULL;
     899                Int4    nbsol;
     900                int     verbosity;
     901
     902                int   i,j,k,iA,iB,iC;
     903                int   iv;
     904
     905                /*Recover options*/
     906                verbosity=bamgopts->verbose;
     907
     908                /*Get and process fields*/
     909                s=bamgopts->field;
     910                nbsol=bamgopts->numfields;
     911
     912                //initialization of some variables
     913                double* ss=(double*)s;
     914                double  sA,sB,sC;
     915                Real8*  detT = new Real8[nbt];
     916                Real8*  sumareas = new Real8[nbv];
     917                Real8*  alpha= new Real8[nbt*3];
     918                Real8*  beta = new Real8[nbt*3];
     919                Real8*  dx_elem    = new Real8[nbt];
     920                Real8*  dy_elem    = new Real8[nbt];
     921                Real8*  dx_vertex  = new Real8[nbv];
     922                Real8*  dy_vertex  = new Real8[nbv];
     923                Real8*  dxdx_elem  = new Real8[nbt];
     924                Real8*  dxdy_elem  = new Real8[nbt];
     925                Real8*  dydy_elem  = new Real8[nbt];
     926                Real8*  dxdx_vertex= new Real8[nbv];
     927                Real8*  dxdy_vertex= new Real8[nbv];
     928                Real8*  dydy_vertex= new Real8[nbv];
     929
     930                //display infos
     931                if(verbosity>1) {
     932                        printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv);
     933                }
     934
     935                //first, build the chains that will be used for the Hessian computation, as weel as the area of each element
     936                int head_s[nbv];
     937                int next_p[3*nbt];
     938                int  p=0;
     939                //initialization
     940                for(i=0;i<nbv;i++){
     941                        sumareas[i]=0;
     942                        head_s[i]=-1;
     943                }
     944                for(i=0;i<nbt;i++){
     945
     946                        //lopp over the real triangles (no boundary elements)
     947                        if(triangles[i].link){
     948
     949                                //get current triangle t
     950                                const Triangle &t=triangles[i];
     951
     952                                // coor of 3 vertices
     953                                R2 A=t[0];
     954                                R2 B=t[1];
     955                                R2 C=t[2];
     956
     957                                //compute triangle determinant (2*Area)
     958                                Real8 dett = bamg::Area2(A,B,C);
     959                                detT[i]=dett;
     960
     961                                /*The nodal functions are such that for a vertex A:
     962                                 *    N_A(x,y)=alphaA x + beta_A y +gamma_A
     963                                 *    N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
     964                                 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
     965                                 * leads to:
     966                                 *    N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
     967                                 * and this gives:
     968                                 *    alpha_A = (yB-yC)/(2*Area(T))*/
     969                                alpha[i*3+0]=(B.y-C.y)/dett;
     970                                alpha[i*3+1]=(C.y-A.y)/dett;
     971                                alpha[i*3+2]=(A.y-B.y)/dett;
     972                                beta[ i*3+0]=(C.x-B.x)/dett;
     973                                beta[ i*3+1]=(A.x-C.x)/dett;
     974                                beta[ i*3+2]=(B.x-A.x)/dett;
     975
     976                                //compute chains
     977                                for(j=0;j<3;j++){
     978                                        k=Number(triangles[i][j]);
     979                                        next_p[p]=head_s[k];
     980                                        head_s[k]=p++;
     981
     982                                        //add area to sumareas
     983                                        sumareas[k]+=dett;
     984                                }
     985
     986                        }
     987                }
     988
     989                //for all Solutions
     990                for (Int4 nusol=0;nusol<nbsol;nusol++) {
     991                        Real8 smin=ss[nusol],smax=ss[nusol];
     992
     993                        //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
     994                        for ( iv=0,k=0; iv<nbv; iv++){
     995                                smin=Min(smin,ss[iv*nbsol+nusol]);
     996                                smax=Max(smax,ss[iv*nbsol+nusol]);
     997                        }
     998                        Real8 sdelta=smax-smin;
     999                        Real8 absmax=Max(Abs(smin),Abs(smax));
     1000
     1001                        //display info
     1002                        if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g\n",nusol,smin,smax,sdelta);
     1003
     1004                        //skip constant field
     1005                        if (sdelta < 1.0e-10*Max(absmax,1e-20)){
     1006                                printf("      Solution %i is constant, skipping...\n",nusol);
     1007                                continue;
     1008                        }
     1009
     1010                        //initialize the hessian and gradient matrices
     1011                        for ( iv=0,k=0; iv<nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
     1012
     1013                        //1: Compute gradient for each element (exact)
     1014                        for (i=0;i<nbt;i++){
     1015                                if(triangles[i].link){
     1016                                        // number of the 3 vertices
     1017                                        iA = Number(triangles[i][0]);
     1018                                        iB = Number(triangles[i][1]);
     1019                                        iC = Number(triangles[i][2]);
     1020
     1021                                        // value of the P1 fonction on 3 vertices
     1022                                        sA = ss[iA*nbsol+nusol];
     1023                                        sB = ss[iB*nbsol+nusol];
     1024                                        sC = ss[iC*nbsol+nusol];
     1025
     1026                                        //gradient = (sum alpha_i s_i, sum_i beta_i s_i)
     1027                                        dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
     1028                                        dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
     1029                                }
     1030                        }
     1031
     1032                        //2: then compute a gradient for each vertex using a P2 projection
     1033                        for(i=0;i<nbv;i++){
     1034                                for(p=head_s[i];p!=-1;p=next_p[p]){
     1035                                        //Get triangle number
     1036                                        k=(Int4)(p/3);
     1037                                        dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
     1038                                        dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
     1039                                }
     1040                        }
     1041
     1042                        //3: compute Hessian matrix on each element
     1043                        for (i=0;i<nbt;i++){
     1044                                if(triangles[i].link){
     1045                                        // number of the 3 vertices
     1046                                        iA = Number(triangles[i][0]);
     1047                                        iB = Number(triangles[i][1]);
     1048                                        iC = Number(triangles[i][2]);
     1049
     1050                                        //Hessian
     1051                                        dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
     1052                                        dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
     1053                                        dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
     1054                                }
     1055                        }
     1056
     1057                        //4: finaly compute Hessian on each vertex using the second P2 projection
     1058                        for(i=0;i<nbv;i++){
     1059                                for(p=head_s[i];p!=-1;p=next_p[p]){
     1060                                        //Get triangle number
     1061                                        k=(Int4)(p/3);
     1062                                        dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
     1063                                        dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
     1064                                        dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
     1065                                }
     1066                        }
     1067
     1068                        /*Compute Metric from Hessian*/
     1069                        for ( iv=0;iv<nbv;iv++){
     1070                                vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts);
     1071                        }
     1072
     1073                }//for all solutions
     1074
     1075                //clean up
     1076                delete [] detT;
     1077                delete [] alpha;
     1078                delete [] beta;
     1079                delete [] sumareas;
     1080                delete [] dx_elem;
     1081                delete [] dy_elem;
     1082                delete [] dx_vertex;
     1083                delete [] dy_vertex;
     1084                delete [] dxdx_elem;
     1085                delete [] dxdy_elem;
     1086                delete [] dydy_elem;
     1087                delete [] dxdx_vertex;
     1088                delete [] dxdy_vertex;
     1089                delete [] dydy_vertex;
     1090        }
     1091        /*}}}1*/
     1092        /*FUNCTION Triangles::BuildMetric1 (Green formula){{{1*/
     1093        void Triangles::BuildMetric1(BamgOpts* bamgopts){
    8991094
    9001095                /*Options*/
    … …  
    9021097                double* s=NULL;
    9031098                Int4 nbsol;
    904                 int* typsols=NULL;
    9051099                int NbJacobi;
    9061100                int verbosity;
    … …  
    9131107                s=bamgopts->field;
    9141108                nbsol=bamgopts->numfields;
    915                 typsols=(int*)xmalloc(1*sizeof(int));
    916                 typsols[0]=0; // only one dof per node
    917 
    918                 //sizeoftype = {1,2,3,4}
    919                 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ;
    920 
    921                 // computation of the number of fields
    922                 Int4 ntmp = 0;
    923                 if (typsols){
    924                         //if there is more than one dof per vertex for one field
    925                         //increase ntmp to take into account all the fields
    926                         //if nbsol=1
    927                         for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
    928                 }
    929                 else ntmp = nbsol;
    930 
    931                 //n is the total number of fields
    932                 const Int4 n = ntmp;
    9331109
    9341110                //initialization of some variables
    … …  
    9491125                //display infos
    9501126                if(verbosity>1) {
    951                         printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
     1127                        printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv);
    9521128                }
    9531129
    … …  
    10181194                for (Int4 nusol=0;nusol<nbsol;nusol++) {
    10191195
    1020                         Real8 smin=ss[0],smax=ss[0];
     1196                        Real8 smin=ss[nusol],smax=ss[nusol];
    10211197                        Real8 h1=1.e30,h2=1e-30,rx=0;
    10221198                        Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 
    1023                         int   nbfield=typsols?sizeoftype[typsols[nusol]]:1;
    1024 
    1025                         //only one field
    1026                         if (nbfield == 1) {
    1027                                 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
    1028                                 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
    1029                                         dxdx[iv]=dxdy[iv]=dydy[iv]=0;
    1030                                         smin=Min(smin,ss[k]);
    1031                                         smax=Max(smax,ss[k]);
    1032                                 }
    1033                         }
    1034 
    1035                         //vectorial case
    1036                         else{
    1037                                 for (iv=0,k=0;iv<nbv;iv++,k+=n ){
    1038                                         //compute v = √sum(s[i]^2)
    1039                                         double v=0;                 
    1040                                         for (int i=0;i<nbfield;i++){
    1041                                                 v += ss[k+i]*ss[k+i];
    1042                                         }
    1043                                         v = sqrt(v);
    1044                                         smin=Min(smin,v);
    1045                                         smax=Max(smax,v);
    1046                                 }
     1199
     1200                        //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
     1201                        for ( iv=0,k=0; iv<nbv; iv++ ){
     1202                                dxdx[iv]=dxdy[iv]=dydy[iv]=0;
     1203                                smin=Min(smin,ss[iv*nbsol+nusol]);
     1204                                smax=Max(smax,ss[iv*nbsol+nusol]);
    10471205                        }
    10481206                        Real8 sdelta=smax-smin;
    … …  
    10501208
    10511209                        //display info
    1052                         if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
     1210                        if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbsol);
    10531211
    10541212                        //skip constant field
    1055                         if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){
     1213                        if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
    10561214                                if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
    10571215                                continue;
    … …  
    10611219                        double* sf=ss;
    10621220
    1063                         //loop over all the fields of the solution
    1064                         for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){
    1065                                 //ss++ so that for each iteration ss points toward the right field
    1066 
    10671221                                //initialize the hessian matrix
    1068                                 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
     1222                                for ( iv=0,k=0; iv<nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
    10691223
    10701224                                //loop over the triangles
    … …  
    10971251
    10981252                                                // value of the P1 fonction on 3 vertices
    1099                                                 sA = ss[iA*n];
    1100                                                 sB = ss[iB*n];
    1101                                                 sC = ss[iC*n];
     1253                                                sA = ss[iA*nbsol+nusol];
     1254                                                sB = ss[iB*nbsol+nusol];
     1255                                                sC = ss[iC*nbsol+nusol];
    11021256
    11031257                                                /*The nodal functions are such that for a vertex A:
    … …  
    11451299
    11461300                                Int4 kk=0;
    1147                                 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
     1301                                for ( iv=0,k=0 ; iv<nbv; iv++){
    11481302                                        if(Mmassxx[iv]>0){
    11491303                                                dxdx[iv] /= 2*Mmassxx[iv];
    … …  
    12071361                                /*Compute Metric from Hessian*/
    12081362                                for ( iv=0;iv<nbv;iv++){
    1209                                         vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts);
    1210                                 }
    1211 
    1212                         } //  end of for all field
     1363                                        vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts);
     1364                                }
     1365
    12131366                }// end for all solution
    12141367
    … …  
    12231376                delete []  OnBoundary;
    12241377
    1225         }
    1226         /*}}}1*/
    1227         /*FUNCTION Triangles::BuildMetric1 (from P2 on 4T){{{1*/
    1228         void Triangles::BuildMetric1(BamgOpts* bamgopts){
    1229                 //  the array of solution s is store   
    1230                 // sol0,sol1,...,soln    on vertex 0
    1231                 //  sol0,sol1,...,soln   on vertex 1
    1232                 //  etc.
    1233 
    1234                 /*Options*/
    1235                 const int dim = 2;
    1236                 double* s=NULL;
    1237                 Int4 nbsol;
    1238                 int* typsols=NULL;
    1239                 int NbJacobi;
    1240                 int verbosity;
    1241 
    1242                 /*Recover options*/
    1243                 verbosity=bamgopts->verbose;
    1244                 NbJacobi=bamgopts->nbjacobi;
    1245 
    1246                 /*Get and process fields*/
    1247                 s=bamgopts->field;
    1248                 nbsol=bamgopts->numfields;
    1249                 typsols=(int*)xmalloc(1*sizeof(int));
    1250                 typsols[0]=0; // only one dof per node
    1251 
    1252                 //sizeoftype = {1,2,3,4}
    1253                 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ;
    1254 
    1255                 // computation of the number of fields
    1256                 Int4 ntmp = 0;
    1257                 if (typsols){
    1258                         //if there is more than one dof per vertex for one field
    1259                         //increase ntmp to take into account all the fields
    1260                         //if nbsol=1
    1261                         for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
    1262                 }
    1263                 else ntmp = nbsol;
    1264 
    1265                 //n is the total number of fields
    1266                 const Int4 n = ntmp;
    1267 
    1268                 //initialization of some variables
    1269                 Int4    i,k,iA,iB,iC,iv;
    1270                 R2      O(0,0);
    1271                 double* ss=(double*)s;
    1272                 double  sA,sB,sC;
    1273                 Real8*  detT = new Real8[nbt];
    1274                 Real8*  Mmass= new Real8[nbv];
    1275                 Real8*  Mmassxx= new Real8[nbv];
    1276                 Real8*  dxdx= new Real8[nbv];
    1277                 Real8*  dxdy= new Real8[nbv];
    1278                 Real8*  dydy= new Real8[nbv];
    1279                 Real8*  workT= new Real8[nbt];
    1280                 Real8*  workV= new Real8[nbv];
    1281                 int*    OnBoundary = new int[nbv];
    1282 
    1283                 //display infos
    1284                 if(verbosity>1) {
    1285                         printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
    1286                 }
    1287 
    1288                 //initialize Mmass, OnBoundary and Massxx by zero
    1289                 for (iv=0;iv<nbv;iv++){
    1290                         Mmass[iv]=0;
    1291                         OnBoundary[iv]=0;
    1292                         Mmassxx[iv]=0;
    1293                 }
    1294 
    1295                 //Build detT Mmas Mmassxx workT and OnBoundary
    1296                 for (i=0;i<nbt;i++){
    1297 
    1298                         //lopp over the real triangles (no boundary elements)
    1299                         if(triangles[i].link){
    1300 
    1301                                 //get current triangle t
    1302                                 const Triangle &t=triangles[i];
    1303 
    1304                                 // coor of 3 vertices
    1305                                 R2 A=t[0];
    1306                                 R2 B=t[1];
    1307                                 R2 C=t[2];
    1308 
    1309                                 // number of the 3 vertices
    1310                                 iA = Number(t[0]);
    1311                                 iB = Number(t[1]);
    1312                                 iC = Number(t[2]);
    1313 
    1314                                 //compute triangle determinant (2*Area)
    1315                                 Real8 dett = bamg::Area2(A,B,C);
    1316                                 detT[i]=dett;
    1317                                 dett /= 6;
    1318 
    1319                                 // construction of OnBoundary (flag=1 if on boundary, else 0)
    1320                                 int nbb=0;
    1321                                 for(int j=0;j<3;j++){
    1322                                         //get adjacent triangle
    1323                                         Triangle *ta=t.Adj(j);
    1324                                         //if there is no adjacent triangle, the edge of the triangle t is on boundary
    1325                                         if ( !ta || !ta->link){
    1326                                                 //mark the two vertices of the edge as OnBoundary
    1327                                                 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1;
    1328                                                 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1;
    1329                                                 nbb++;
    1330                                         }
    1331                                 }
    1332 
    1333                                 //number of vertices on boundary for current triangle t
    1334                                 workT[i] = nbb;
    1335 
    1336                                 //Build Mmass Mmass[i] = Mmass[i] + Area/3
    1337                                 Mmass[iA] += dett;
    1338                                 Mmass[iB] += dett;
    1339                                 Mmass[iC] += dett;
    1340 
    1341                                 //Build Massxx = Mmass
    1342                                 if(nbb==0){
    1343                                         Mmassxx[iA] += dett;
    1344                                         Mmassxx[iB] += dett;
    1345                                         Mmassxx[iC] += dett;
    1346                                 }
    1347                         }
    1348 
    1349                         //else: the triangle is a boundary triangle -> workT=-1
    1350                         else workT[i]=-1;
    1351                 }
    1352 
    1353                 //for all Solution 
    1354                 for (Int4 nusol=0;nusol<nbsol;nusol++) {
    1355 
    1356                         Real8 smin=ss[0],smax=ss[0];
    1357                         Real8 h1=1.e30,h2=1e-30,rx=0;
    1358                         Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 
    1359                         int   nbfield=typsols?sizeoftype[typsols[nusol]]:1;
    1360 
    1361                         //only one field
    1362                         if (nbfield == 1) {
    1363                                 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
    1364                                 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
    1365                                         dxdx[iv]=dxdy[iv]=dydy[iv]=0;
    1366                                         smin=Min(smin,ss[k]);
    1367                                         smax=Max(smax,ss[k]);
    1368                                 }
    1369                         }
    1370 
    1371                         //vectorial case
    1372                         else{
    1373                                 for (iv=0,k=0;iv<nbv;iv++,k+=n ){
    1374                                         //compute v = √sum(s[i]^2)
    1375                                         double v=0;                 
    1376                                         for (int i=0;i<nbfield;i++){
    1377                                                 v += ss[k+i]*ss[k+i];
    1378                                         }
    1379                                         v = sqrt(v);
    1380                                         smin=Min(smin,v);
    1381                                         smax=Max(smax,v);
    1382                                 }
    1383                         }
    1384                         Real8 sdelta=smax-smin;
    1385                         Real8 absmax=Max(Abs(smin),Abs(smax));
    1386 
    1387                         //display info
    1388                         if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
    1389 
    1390                         //skip constant field
    1391                         if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){
    1392                                 if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
    1393                                 continue;
    1394                         }
    1395 
    1396                         //pointer toward ss that is also a pointer toward s (solutions)
    1397                         double* sf=ss;
    1398 
    1399                         //loop over all the fields of the solution
    1400                         for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){
    1401                                 //ss++ so that for each iteration ss points toward the right field
    1402 
    1403                                 //initialize the hessian matrix
    1404                                 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
    1405 
    1406                                 //loop over the triangles
    1407                                 for (i=0;i<nbt;i++){
    1408 
    1409                                         //for real all triangles
    1410                                         if(triangles[i].link){
    1411 
    1412                                                 // coor of 3 vertices
    1413                                                 R2 A=triangles[i][0];
    1414                                                 R2 B=triangles[i][1];
    1415                                                 R2 C=triangles[i][2];
    1416 
    1417                                                 //warning: the normal is internal and the size is the length of the edge
    1418                                                 R2 nAB = Orthogonal(B-A);
    1419                                                 R2 nBC = Orthogonal(C-B);
    1420                                                 R2 nCA = Orthogonal(A-C);
    1421                                                 //note that :  nAB + nBC + nCA == 0
    1422 
    1423                                                 // number of the 3 vertices
    1424                                                 iA = Number(triangles[i][0]);
    1425                                                 iB = Number(triangles[i][1]);
    1426                                                 iC = Number(triangles[i][2]);
    1427 
    1428                                                 // for the test of  boundary edge
    1429                                                 // the 3 adj triangles
    1430                                                 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
    1431                                                 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
    1432                                                 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);
    1433 
    1434                                                 // value of the P1 fonction on 3 vertices
    1435                                                 sA = ss[iA*n];
    1436                                                 sB = ss[iB*n];
    1437                                                 sC = ss[iC*n];
    1438 
    1439                                                 /*The nodal functions are such that for a vertex A:
    1440                                                   N_A(x,y)=alphaA x + beta_A y +gamma_A
    1441                                                   N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
    1442                                                   solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
    1443                                                   leads to:
    1444                                                   N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
    1445                                                   and this gives:
    1446                                                   alpha_A = (yB-yC)/(2*Area(T))
    1447                                                   beta_A = (xC-xB)/(2*Area(T))
    1448                                                   and therefore:
    1449                                                   grad N_A = nA / detT
    1450                                                   for an interpolation of a solution s:
    1451                                                   grad(s) = s * sum_{i=A,B,C} grad(N_i) */
    1452 
    1453                                                 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i];
    1454 
    1455                                                 //from P2 on 4T to compute Hessian
    1456 
    1457                                                 int   nbb=0;
    1458                                                 Real8 dd = detT[i];
    1459                                                 Real8 taa[3][3],bb[3];
    1460 
    1461                                                 // construction of the transpose of lin system
    1462                                                 for (int j=0;j<3;j++){
    1463                                                         int              ie = OppositeEdge[j];
    1464                                                         TriangleAdjacent ta = triangles[i].Adj(ie);
    1465                                                         Triangle*        tt = ta;
    1466 
    1467 
    1468                                                         /* V is the vertex opposed to the edge shared by the
    1469                                                          * current triangle i and its neighbor ta
    1470                                                          *         C       
    1471                                                          *       / | \     
    1472                                                          *      /  |  \   
    1473                                                          *     /   |   \   
    1474                                                          * V  / ta | i  \ B
    1475                                                          *    \    |    / 
    1476                                                          *     \   |   /   
    1477                                                          *      \  |  /
    1478                                                          *       \ | /     
    1479                                                          *         A
    1480                                                          * lA=area(VBC)/Area(ABC)
    1481                                                          * lB=area(AVC)/Area(ABC)
    1482                                                          * lC=area(ABV)/Area(ABC)
    1483                                                          */
    1484 
    1485                                                         //first, get V
    1486                                                         Vertex &v = *ta.OppositeVertex();
    1487                                                         Int4   iV = Number(v);
    1488 
    1489                                                         //if the adjacent triangle is not a boundary triangle:
    1490                                                         if (tt && tt->link){
    1491                                                                 Vertex &v = *ta.OppositeVertex();
    1492                                                                 R2      V = v;
    1493                                                                 Int4   iV = Number(v);
    1494 
    1495                                                                 //get lA lB and lC
    1496                                                                 Real8  lA = bamg::Area2(V,B,C)/dd;
    1497                                                                 Real8  lB = bamg::Area2(A,V,C)/dd;
    1498                                                                 Real8  lC = bamg::Area2(A,B,V)/dd;
    1499 
    1500                                                                 //fill taa
    1501                                                                 taa[0][j] = lB*lC;
    1502                                                                 taa[1][j] = lC*lA;
    1503                                                                 taa[2][j] = lA*lB;
    1504 
    1505                                                                 //fill bb
    1506                                                                 bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ; //value of the solution on V minus...
    1507                                                         }
    1508                                                         else{
    1509                                                                 //there might be a problem here: if 2 nodes are inline lA, lB or lC
    1510                                                                 //is 0. This means that one column has only one term.
    1511                                                                 //if instead of 0.1 we use 0, and taa[j][j]=1, we might have det(taa)=0
    1512                                                                 nbb++;
    1513                                                                 taa[0][j]=0.1;
    1514                                                                 taa[1][j]=0.1;
    1515                                                                 taa[2][j]=0.1;
    1516                                                                 taa[j][j]=1;
    1517                                                                 bb[j]=0;
    1518                                                         }
    1519                                                 }
    1520 
    1521                                                 // resolution of 3x3 linear system transpose
    1522                                                 Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);   
    1523                                                 Real8 cBC   =  det3x3(bb,taa[1],taa[2]);
    1524                                                 Real8 cCA   =  det3x3(taa[0],bb,taa[2]);
    1525                                                 Real8 cAB   =  det3x3(taa[0],taa[1],bb);
    1526                                                 if (!det33){
    1527                                                         throw ErrorException(__FUNCT__,exprintf("for triangle %i, det33==0! cannot compute theHessian matrix using Hessiantype=1",i+1));
    1528                                                 }
    1529 
    1530                                                 // computation of the Hessian in the element
    1531                                                 // H( li*lj) = grad li grad lj + grad lj grad lj
    1532                                                 // grad li = njk  / detT ; with i j k =(A,B,C)
    1533                                                 Real8 Hxx = cAB * ( nBC.x*nCA.x) +  cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);
    1534                                                 Real8 Hyy = cAB * ( nBC.y*nCA.y) +  cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);
    1535                                                 Real8 Hxy = cAB * ( nBC.y*nCA.x) +  cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x)
    1536                                                   + cAB * ( nBC.x*nCA.y) +  cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);
    1537                                                 if(nbb==0){
    1538                                                         dxdx[iA] += Hxx;
    1539                                                         dydy[iA] += Hyy;
    1540                                                         dxdy[iA] += Hxy;
    1541 
    1542                                                         dxdx[iB] += Hxx;
    1543                                                         dydy[iB] += Hyy;
    1544                                                         dxdy[iB] += Hxy;
    1545 
    1546                                                         dxdx[iC] += Hxx;
    1547                                                         dydy[iC] += Hyy;
    1548                                                         dxdy[iC] += Hxy;
    1549                                                 }
    1550                                         } // for real all triangles
    1551                                 }
    1552 
    1553                                 Int4 kk=0;
    1554                                 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
    1555                                         if(Mmassxx[iv]>0){
    1556                                                 dxdx[iv] /= 2*Mmassxx[iv];
    1557                                                 // warning optimization (1) on term dxdy[iv]*ci/2
    1558                                                 dxdy[iv] /= 4*Mmassxx[iv];
    1559                                                 dydy[iv] /= 2*Mmassxx[iv];
    1560                                                 // Compute the matrix with abs(eigen value)
    1561                                                 Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
    1562                                                 MatVVP2x2 Vp(M);
    1563                                                 Vp.Abs();
    1564                                                 M = Vp;
    1565                                                 dxdx[iv] = M.a11;
    1566                                                 dxdy[iv] = M.a21;
    1567                                                 dydy[iv] = M.a22;
    1568                                         }
    1569                                         else kk++;
    1570                                 }
    1571 
    1572                                 // correction of second derivative
    1573                                 // by a laplacien
    1574                                 Real8* d2[3] = {dxdx, dxdy, dydy};
    1575                                 Real8* dd;
    1576                                 for (int xy = 0;xy<3;xy++) {
    1577                                         dd = d2[xy];
    1578                                         // do least 2 iteration for boundary problem
    1579                                         for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
    1580                                                 for (i=0;i<nbt;i++)
    1581                                                  if(triangles[i].link){// the real triangles
    1582                                                          // number of the 3 vertices
    1583                                                          iA = Number(triangles[i][0]);
    1584                                                          iB = Number(triangles[i][1]);
    1585                                                          iC = Number(triangles[i][2]);
    1586                                                          Real8 cc=3;
    1587                                                          if(ijacobi==0)
    1588                                                           cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
    1589                                                          workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
    1590                                                  }
    1591                                                 for (iv=0;iv<nbv;iv++) workV[iv]=0;
    1592 
    1593                                                 for (i=0;i<nbt;i++){
    1594                                                         if(triangles[i].link){ // the real triangles
    1595                                                                 // number of the 3 vertices
    1596                                                                 iA = Number(triangles[i][0]);
    1597                                                                 iB = Number(triangles[i][1]);
    1598                                                                 iC = Number(triangles[i][2]);
    1599                                                                 Real8 cc =  workT[i]*detT[i];
    1600                                                                 workV[iA] += cc;
    1601                                                                 workV[iB] += cc;
    1602                                                                 workV[iC] += cc;
    1603                                                         }
    1604                                                 }
    1605 
    1606                                                 for (iv=0;iv<nbv;iv++){
    1607                                                         if( ijacobi<NbJacobi || OnBoundary[iv]){
    1608                                                                 dd[iv] = workV[iv]/(Mmass[iv]*6);
    1609                                                         }
    1610                                                 }
    1611                                         }
    1612                                 }
    1613 
    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 
    1619                         } //  end of for all field
    1620                 }// end for all solution
    1621 
    1622                 delete [] detT;
    1623                 delete [] Mmass;
    1624                 delete [] dxdx;
    1625                 delete [] dxdy;
    1626                 delete [] dydy;
    1627                 delete []  workT;
    1628                 delete [] workV;
    1629                 delete [] Mmassxx;
    1630                 delete []  OnBoundary;
    1631 
    1632         }
    1633         /*}}}1*/
    1634         /*FUNCTION Triangles::BuildMetric2 (double P2 projection){{{1*/
    1635         void Triangles::BuildMetric2(BamgOpts* bamgopts){
    1636 
    1637                 /*Options*/
    1638                 const int dim = 2;
    1639                 double* s=NULL;
    1640                 Int4   nbsol;
    1641                 int*   typsols=NULL;
    1642                 int    verbosity;
    1643 
    1644                 int   i,j,k,iA,iB,iC;
    1645                 int   iv,nbfield;
    1646 
    1647                 /*Recover options*/
    1648                 verbosity=bamgopts->verbose;
    1649 
    1650                 /*Get and process fields*/
    1651                 s=bamgopts->field;
    1652                 nbsol=bamgopts->numfields;
    1653                 typsols=(int*)xmalloc(1*sizeof(int));
    1654                 typsols[0]=0; // only one dof per node
    1655 
    1656                 //sizeoftype = {1,2,3,4}
    1657                 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ;
    1658 
    1659                 // computation of the number of fields
    1660                 Int4 ntmp = 0;
    1661                 if (typsols){
    1662                         //if there is more than one dof per vertex for one field
    1663                         //increase ntmp to take into account all the fields
    1664                         //if nbsol=1
    1665                         for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
    1666                 }
    1667                 else ntmp = nbsol;
    1668 
    1669                 //n is the total number of fields
    1670                 const Int4 n = ntmp;
    1671 
    1672                 //initialization of some variables
    1673                 double* ss=(double*)s;
    1674                 double  sA,sB,sC;
    1675                 Real8*  detT = new Real8[nbt];
    1676                 Real8*  sumareas = new Real8[nbv];
    1677                 Real8*  alpha= new Real8[nbt*3];
    1678                 Real8*  beta = new Real8[nbt*3];
    1679                 Real8*  dx_elem    = new Real8[nbt];
    1680                 Real8*  dy_elem    = new Real8[nbt];
    1681                 Real8*  dx_vertex  = new Real8[nbv];
    1682                 Real8*  dy_vertex  = new Real8[nbv];
    1683                 Real8*  dxdx_elem  = new Real8[nbt];
    1684                 Real8*  dxdy_elem  = new Real8[nbt];
    1685                 Real8*  dydy_elem  = new Real8[nbt];
    1686                 Real8*  dxdx_vertex= new Real8[nbv];
    1687                 Real8*  dxdy_vertex= new Real8[nbv];
    1688                 Real8*  dydy_vertex= new Real8[nbv];
    1689 
    1690                 //display infos
    1691                 if(verbosity>1) {
    1692                         printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
    1693                 }
    1694 
    1695                 //first, build the chains that will be used for the Hessian computation, as weel as the area of each element
    1696                 int head_s[nbv];
    1697                 int next_p[3*nbt];
    1698                 int  p=0;
    1699                 //initialization
    1700                 for(i=0;i<nbv;i++){
    1701                         sumareas[i]=0;
    1702                         head_s[i]=-1;
    1703                 }
    1704                 for(i=0;i<nbt;i++){
    1705 
    1706                         //lopp over the real triangles (no boundary elements)
    1707                         if(triangles[i].link){
    1708 
    1709                                 //get current triangle t
    1710                                 const Triangle &t=triangles[i];
    1711 
    1712                                 // coor of 3 vertices
    1713                                 R2 A=t[0];
    1714                                 R2 B=t[1];
    1715                                 R2 C=t[2];
    1716 
    1717                                 //compute triangle determinant (2*Area)
    1718                                 Real8 dett = bamg::Area2(A,B,C);
    1719                                 detT[i]=dett;
    1720 
    1721                                 /*The nodal functions are such that for a vertex A:
    1722                                  *    N_A(x,y)=alphaA x + beta_A y +gamma_A
    1723                                  *    N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
    1724                                  * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
    1725                                  * leads to:
    1726                                  *    N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
    1727                                  * and this gives:
    1728                                  *    alpha_A = (yB-yC)/(2*Area(T))*/
    1729                                 alpha[i*3+0]=(B.y-C.y)/dett;
    1730                                 alpha[i*3+1]=(C.y-A.y)/dett;
    1731                                 alpha[i*3+2]=(A.y-B.y)/dett;
    1732                                 beta[ i*3+0]=(C.x-B.x)/dett;
    1733                                 beta[ i*3+1]=(A.x-C.x)/dett;
    1734                                 beta[ i*3+2]=(B.x-A.x)/dett;
    1735 
    1736                                 //compute chains
    1737                                 for(j=0;j<3;j++){
    1738                                         k=Number(triangles[i][j]);
    1739                                         next_p[p]=head_s[k];
    1740                                         head_s[k]=p++;
    1741 
    1742                                         //add area to sumareas
    1743                                         sumareas[k]+=dett;
    1744                                 }
    1745 
    1746                         }
    1747                 }
    1748 
    1749                 //for all Solutions
    1750                 for (Int4 nusol=0;nusol<nbsol;nusol++) {
    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 
    1788 
    1789                         //loop over all the fields of the solution
    1790                         for (Int4 nufield=0;nufield<nbfield;nufield++,s++){
    1791                                 //ss++ so that for each iteration ss points toward the right field
    1792 
    1793                                 //initialize the hessian and gradient matrices
    1794                                 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
    1795 
    1796                                 //1: Compute gradient for each element (exact)
    1797                                 for (i=0;i<nbt;i++){
    1798                                         if(triangles[i].link){
    1799                                                 // number of the 3 vertices
    1800                                                 iA = Number(triangles[i][0]);
    1801                                                 iB = Number(triangles[i][1]);
    1802                                                 iC = Number(triangles[i][2]);
    1803 
    1804                                                 // value of the P1 fonction on 3 vertices
    1805                                                 sA = ss[iA*n];
    1806                                                 sB = ss[iB*n];
    1807                                                 sC = ss[iC*n];
    1808 
    1809                                                 //gradient = (sum alpha_i s_i, sum_i beta_i s_i)
    1810                                                 dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
    1811                                                 dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
    1812                                         }
    1813                                 }
    1814 
    1815                                 //2: then compute a gradient for each vertex using a P2 projection
    1816                                 for(i=0;i<nbv;i++){
    1817                                         for(p=head_s[i];p!=-1;p=next_p[p]){
    1818                                                 //Get triangle number
    1819                                                 k=(Int4)(p/3);
    1820                                                 dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
    1821                                                 dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
    1822                                         }
    1823                                 }
    1824 
    1825                                 //3: compute Hessian matrix on each element
    1826                                 for (i=0;i<nbt;i++){
    1827                                         if(triangles[i].link){
    1828                                                 // number of the 3 vertices
    1829                                                 iA = Number(triangles[i][0]);
    1830                                                 iB = Number(triangles[i][1]);
    1831                                                 iC = Number(triangles[i][2]);
    1832 
    1833                                                 //Hessian
    1834                                                 dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
    1835                                                 dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
    1836                                                 dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
    1837                                         }
    1838                                 }
    1839 
    1840                                 //4: finaly compute Hessian on each vertex using the second P2 projection
    1841                                 for(i=0;i<nbv;i++){
    1842                                         for(p=head_s[i];p!=-1;p=next_p[p]){
    1843                                                 //Get triangle number
    1844                                                 k=(Int4)(p/3);
    1845                                                 dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
    1846                                                 dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
    1847                                                 dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
    1848                                         }
    1849                                 }
    1850 
    1851                                 /*Compute Metric from Hessian*/
    1852                                 for ( iv=0;iv<nbv;iv++){
    1853                                         vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[n*iv],bamgopts);
    1854                                 }
    1855 
    1856                         }//for all fields
    1857                 }//for all solutions
    1858 
    1859                 //clean up
    1860                 delete [] detT;
    1861                 delete [] alpha;
    1862                 delete [] beta;
    1863                 delete [] sumareas;
    1864                 delete [] dx_elem;
    1865                 delete [] dy_elem;
    1866                 delete [] dx_vertex;
    1867                 delete [] dy_vertex;
    1868                 delete [] dxdx_elem;
    1869                 delete [] dxdy_elem;
    1870                 delete [] dydy_elem;
    1871                 delete [] dxdx_vertex;
    1872                 delete [] dxdy_vertex;
    1873                 delete [] dydy_vertex;
    18741378        }
    18751379        /*}}}1*/
    … …  
    30322536                else {
    30332537                        if (!quadtree){
    3034                                 throw ErrorException(__FUNCT__,exprintf("!quadtree"));
     2538                                throw ErrorException(__FUNCT__,exprintf("no starting triangle provided and no quadtree available"));
    30352539                        }
    30362540                        Vertex *a = quadtree->NearestVertex(B.x,B.y) ;
    30372541
    3038                         if (! a || !a->t ) {
     2542                        if (!a || !a->t ) {
    30392543                                if (a) {
    30402544                                        printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a));
    … …  
    30512555                }
    30522556                Icoor2  detop ;
    3053                 int kkkk =0; // number of test triangle
    3054 
    3055                 while ( t->det < 0){ // the initial triangles is outside 
     2557
     2558                // number of test triangle
     2559                int kkkk =0;
     2560
     2561                // if the initial triangles is outside 
     2562                while ( t->det < 0){
    30562563                        int k0=(*t)(0) ?  ((  (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1  )) : 0;
    3057                         if (k0<0){ // k0 the NULL  vertex
     2564                        if (k0<0){ // k0 the NULL vertex
    30582565                                throw ErrorException(__FUNCT__,exprintf("k0<0"));
    30592566                        }
    … …  
    31512658                nbv=0;
    31522659                for (i=0;i<Gh.nbv;i++){
    3153                         /*Add vertex only if required and not duplicate
     2660                        /* Add vertex only if required and not duplicate
    31542661                         * (IsThe is a method of GeometricalVertex that checks
    31552662                         * that we are taking into acount only one vertex is
    31562663                         * 2 vertices are superimposed) */
    31572664                        if (Gh[i].Required() && Gh[i].IsThe()) {//Gh  vertices Required
    3158                                 if (nbv<nbvx) vertices[nbv]=Gh[i];
    3159                                 Gh[i].to = vertices + nbv;// save Geom -> Th
     2665
     2666                                //Add the vertex (provided that nbv<nbvx)
     2667                                if (nbv<nbvx){
     2668                                        vertices[nbv]=Gh[i];
     2669                                }
     2670                                else{
     2671                                        throw ErrorException(__FUNCT__,exprintf("Maximum number of vertices (nbvx = %i) too small",nbvx));
     2672                                }
     2673                               
     2674                                //Add pointer from geometry (Gh) to vertex from mesh (Th)
     2675                                Gh[i].to=vertices+nbv;
     2676
     2677
     2678                                //Build VerticesOnGeomVertex for current point
    31602679                                VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);
     2680
     2681                                //nbv increment
    31612682                                nbv++;
    31622683                        }
    … …  
    31652686                //check that edges is empty
    31662687                if (edges){
    3167                         throw ErrorException(__FUNCT__,exprintf("edges"));
     2688                        throw ErrorException(__FUNCT__,exprintf("edges is empty"));
    31682689                }
    31692690
    … …  
    31842705                        //go through the edges of the geometry
    31852706                        for (i=0;i<Gh.nbe;i++) {
     2707
    31862708                                GeometricalEdge &ei=Gh.edges[i];   
    3187                                 if (!ei.Dup()){ // a good curve (not dup)
     2709
     2710                                // a good curve (not dup)
     2711                                if (!ei.Dup()){
    31882712                                        for(int j=0;j<2;j++) {
    31892713                                                if (!ei.Mark() && ei[j].Required()) {
    … …  
    37343258                BuildMetric1(bamgopts);
    37353259        }
    3736         else if (Hessiantype==2){
    3737                 BuildMetric2(bamgopts);
    3738         }
    37393260        else{
    3740                 throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (0->use Green formula, 1-> from P2 on 4T, 2-> double P2 projection)",Hessiantype));
     3261                throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (1->use Green formula, 0-> double P2 projection)",Hessiantype));
    37413262        }
    37423263}
    … …  
    38983419        /*FUNCTION Triangles::InsertNewPoints{{{1*/
    38993420        Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) {
    3900                 long int verbosity=2;
     3421
     3422                long int verbosity=0;
    39013423                Real8 seuil= 1.414/2 ;// for two close point
    39023424                Int4 i;
    3903                 // insertion part ---
    3904 
    3905                 const Int4 nbvnew = nbv-nbvold;
    3906                 if (verbosity>5) printf("      Try to Insert the %i new points\n",nbvnew);
    39073425                Int4 NbSwap=0;
    39083426                Icoor2 dete[3];
    39093427
    3910                 // construction d'un ordre aleatoire
    3911                 if (! nbvnew)
    3912                  return 0;
    3913                 if (nbvnew) {
    3914                         const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv)  ;
    3915                         Int4 k3 = rand()%nbvnew ;
    3916                         for (Int4 is3=0; is3<nbvnew; is3++) {
    3917                                 register Int4 j = nbvold +(k3 = (k3 + PrimeNumber)% nbvnew);
    3918                                 register Int4 i = nbvold+is3;
    3919                                 ordre[i]= vertices + j;
    3920                                 ordre[i]->ReferenceNumber=i;
    3921                         }
    3922                         // be carefull
    3923                         Int4  iv = nbvold;
    3924                         for (i=nbvold;i<nbv;i++)
    3925                           {// for all the new point
    3926                                 Vertex & vi = *ordre[i];
    3927                                 vi.i = toI2(vi.r);
    3928                                 vi.r = toR2(vi.i);
    3929                                 Real4 hx,hy;
    3930                                 vi.m.Box(hx,hy);
    3931                                 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
    3932                                 if (!quadtree->ToClose(vi,seuil,hi,hj))
    3933                                   {
    3934                                         // a good new point
    3935                                         Vertex & vj = vertices[iv];
    3936                                         Int4 j = vj.ReferenceNumber;
    3937                                         if ( &vj!= ordre[j]){
    3938                                                 throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]"));
    3939                                         }
    3940                                         if(i!=j)
    3941                                           { //  for valgring
    3942                                                 Exchange(vi,vj);
    3943                                                 Exchange(ordre[j],ordre[i]);
    3944                                           }
    3945                                         vj.ReferenceNumber=0;
    3946                                         Triangle *tcvj = FindTriangleContening(vj.i,dete);
    3947                                         if (tcvj && !tcvj->link){
    3948                                                 throw ErrorException(__FUNCT__,exprintf("problem inserting point"));
    3949                                         }
    3950                                         quadtree->Add(vj);
    3951                                         Add(vj,tcvj,dete);
    3952                                         NbSwap += vj.Optim(1);         
    3953                                         iv++;
    3954                                   }
    3955                           }
    3956                         if (verbosity>3) {
    3957                                 printf("   number of new points: %i\n",iv);
    3958                                 printf("   number of to close (?) points: %i\n",nbv-iv);
    3959                                 printf("   number of swap after: %i\n",NbSwap);
    3960                         }
    3961                         nbv = iv;
    3962                 }
     3428                //number of new points
     3429                const Int4 nbvnew=nbv-nbvold;
     3430
     3431                //display info if required
     3432                if (verbosity>5) printf("      Try to Insert %i new points\n",nbvnew);
     3433
     3434                //return if no new points
     3435                if (!nbvnew) return 0;
     3436
     3437                /*construction of a random order*/
     3438                const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv)  ;
     3439                //remainder of the division of rand() by nbvnew
     3440                Int4 k3 = rand()%nbvnew;
     3441                //loop over the new points
     3442                for (Int4 is3=0; is3<nbvnew; is3++){
     3443                        register Int4 j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
     3444                        register Int4 i=nbvold+is3;
     3445                        ordre[i]= vertices + j;
     3446                        ordre[i]->ReferenceNumber=i;
     3447                }
     3448
     3449                // for all the new point
     3450                Int4 iv=nbvold;
     3451                for (i=nbvold;i<nbv;i++){
     3452                        Vertex &vi=*ordre[i];
     3453                        vi.i=toI2(vi.r);
     3454                        vi.r=toR2(vi.i);
     3455                        Real4 hx,hy;
     3456                        vi.m.Box(hx,hy);
     3457                        Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
     3458                        if (!quadtree->ToClose(vi,seuil,hi,hj)){
     3459                                // a good new point
     3460                                Vertex &vj = vertices[iv];
     3461                                Int4  j=vj.ReferenceNumber;
     3462                                if (&vj!=ordre[j]){
     3463                                        throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]"));
     3464                                }
     3465                                if(i!=j){
     3466                                        Exchange(vi,vj);
     3467                                        Exchange(ordre[j],ordre[i]);
     3468                                }
     3469                                vj.ReferenceNumber=0;
     3470                                Triangle *tcvj=FindTriangleContening(vj.i,dete);
     3471                                if (tcvj && !tcvj->link){
     3472                                        tcvj->Echo();
     3473                                        throw ErrorException(__FUNCT__,exprintf("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link));
     3474                                }
     3475                                quadtree->Add(vj);
     3476                                Add(vj,tcvj,dete);
     3477                                NbSwap += vj.Optim(1);         
     3478                                iv++;
     3479                        }
     3480                }
     3481                if (verbosity>3) {
     3482                        printf("         number of new points: %i\n",iv);
     3483                        printf("         number of to close (?) points: %i\n",nbv-iv);
     3484                        printf("         number of swap after: %i\n",NbSwap);
     3485                }
     3486                nbv = iv;
    39633487
    39643488                for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1); 
    … …  
    54885012                        Triangle *tcvi = FindTriangleContening(vi.i,dete);
    54895013                        if (tcvi && !tcvi->link) {
    5490                                 printf("problem inserting point\n");
     5014                                printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
    54915015                        }
    54925016
    … …  
    59095433        /*FUNCTION AGoodNumberPrimeWith{{{1*/
    59105434        Int4 AGoodNumberPrimeWith(Int4 n){
     5435
     5436                //list of big prime numbers
    59115437                const Int4 BigPrimeNumber[] ={ 567890359L,
    59125438                        567890431L,  567890437L,  567890461L,  567890471L,
    … …  
    59145440                        567890591L,  567890599L,  567890621L,  567890629L , 0};
    59155441
    5916                 Int4 o = 0;
    5917                 Int4 pi = BigPrimeNumber[1];
    5918                 for (int i=0; BigPrimeNumber[i]; i++) {
     5442                //initialize o and pi
     5443                Int4 o =0;
     5444                Int4 pi=BigPrimeNumber[1];
     5445
     5446                //loop until BigPrimeNumber[i]==0 (end of BigPrimeNumber)
     5447                for (int i=0; BigPrimeNumber[i]; i++){
     5448
     5449                        //compute r, rest of the remainder of the division of BigPrimeNumber[i] by n
    59195450                        Int4 r = BigPrimeNumber[i] % n;
     5451
     5452                        /*compute oo = min ( r , n-r , |n - 2r|, |n-3r|)*/
    59205453                        Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));
    5921                         if ( o < oo)
    5922                          o=oo,pi=BigPrimeNumber[i];}
    5923                         return pi;
     5454                        if ( o < oo){
     5455                                o=oo;
     5456                                pi=BigPrimeNumber[i];
     5457                        }
     5458                }
     5459                return pi;
    59245460        }
    59255461        /*}}}1*/

  • issm/trunk/src/c/Bamgx/objects/Vertex.cpp

    r2938 r2940  
    206206        }
    207207        /*}}}1*/
     208        /*FUNCTION Vertex::Echo {{{1*/
     209
     210        void Vertex::Echo(void){
     211
     212                printf("Vertex:\n");
     213                printf("  integer   coordinates i.x: %i, i.y: %i\n",i.x,i.y);
     214                printf("  Euclidean coordinates r.x: %g, r.y: %g\n",r.x,r.y);
     215                printf("  ReferenceNumber = %i\n",ReferenceNumber);
     216
     217                return;
     218        }
     219        /*}}}*/
    208220
    209221        /*Intermediary*/
    210         /*FUNCTION QuadQuality{{{{1*/
     222        /*FUNCTION QuadQuality{{{1*/
    211223        double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) {
    212224                // calcul de 4 angles --

  • issm/trunk/src/m/classes/public/bamg.m

    r2938 r2940  
    8383bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1);
    8484bamg_options.AbsError=getfieldvalue(options,'AbsError',0);
    85 bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',2);
     85bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
    8686bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
    8787bamg_options.NbSmooth=getfieldvalue(options,'NbSmooth',3);

  • issm/trunk/src/mex/Bamg/Bamg.cpp

    r2938 r2940  
    152152        FetchData(&field,NULL,&numfields,mxGetField(BAMGOPTIONS,0,"field"));
    153153        bamgopts.numfields=numfields;
    154         printf("numfields = %i\n",numfields);
    155154        bamgopts.field=field;
    156155
Note: See TracChangeset for help on using the changeset viewer.