Changeset 2924


Ignore:
Timestamp:
01/28/10 08:38:17 (15 years ago)
Author:
Mathieu Morlighem
Message:

minor changes

Location:
issm/trunk/src/c/Bamgx
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r2921 r2924  
    148148
    149149                public:
    150                         Triangle * t; // le triangle
    151                         int  a; // le numero de l arete
     150                        Triangle* t; // le triangle
     151                        int  a;      // le numero de l arete
    152152
    153153                        TriangleAdjacent(Triangle  * tt,int  aa): t(tt),a(aa &3) {};

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

    r2917 r2924  
    9191                void  Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;}
    9292
    93                 void Minh(double h) {Max(1.0/(h*h));}
    94                 void Maxh(double h) {Min(1.0/(h*h));}
     93                void Minh(double h) {Min(1.0/(h*h));}
     94                void Maxh(double h) {Max(1.0/(h*h));}
    9595                void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
    9696                MatVVP2x2(const MetricAnIso );
    … …  
    122122        }
    123123
    124         void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
     124        void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
    125125        MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
    126126

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

    r2863 r2924  
    77                  public: 
    88                          R x,y;
     9                          void Echo(){
     10                                  printf("   x: %g\n",x);
     11                                  printf("   y: %g\n",y);
     12                          }
    913                          P2 () :x(0),y(0) {};
    1014                          P2 (R a,R b)  :x(a),y(b)  {}

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

    r2899 r2924  
    7171        /*FUNCTION MetricAnIso::IntersectWith{{{1*/
    7272        int MetricAnIso::IntersectWith(const MetricAnIso M2) {
     73                /*Get a new metric from an existing metric (M1=this)
     74                 * and a new metric given in input M2 using a
     75                 * Simultaneous Matrix Reduction
     76                 * */
     77
    7378                int r=0;
    7479                MetricAnIso & M1 = *this;
    … …  
    7681                double l1,l2;
    7782
    78                 ReductionSimultanee(*this,M2,l1,l2,M);
     83                SimultaneousMatrixReduction(*this,M2,l1,l2,M);
     84
    7985                R2 v1(M.x.x,M.y.x);
    8086                R2 v2(M.x.y,M.y.y);
    … …  
    97103
    98104        /*Intermediary*/
    99         /*FUNCTION ReductionSimultanee{{{1*/
    100         void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
     105        /*FUNCTION SimultaneousMatrixReduction{{{1*/
     106        void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
    101107                double a11=M1.a11,a21=M1.a21,a22=M1.a22;
    102108                double b11=M2.a11,b21=M2.a21,b22=M2.a22;
    … …  
    113119                const double bb = b*b,ac= a*c;
    114120                const double delta = bb - 4 * ac;
    115                 if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) )
    116                   {
     121                if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
    117122                        // racine double;
    118123                        if (Abs(a) < 1.e-30 )
    … …  
    153158                D2xD2 M;
    154159                double l1,l2;
    155                 ReductionSimultanee(M1,M2,l1,l2,M);
     160                SimultaneousMatrixReduction(M1,M2,l1,l2,M);
    156161                R2 v0(M.x.x,M.y.x);
    157162                R2 v1(M.x.y,M.y.y);

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

    r2920 r2924  
    700700                                        MatVVP2x2 Vp(M/coef);
    701701
    702                                         Vp.Maxh(hmin);
    703                                         Vp.Minh(hmax);
     702                                        Vp.Maxh(hmax);
     703                                        Vp.Minh(hmin);
    704704                                        vertices[i].m = Vp;
    705705                                }
    … …  
    12721272
    12731273                                        //modify eigen values according to hmin and hmax
    1274                                         Vp.Maxh(hmin);
    1275                                         Vp.Minh(hmax);
     1274                                        Vp.Maxh(hmax);
     1275                                        Vp.Minh(hmin);
    12761276
    12771277                                        //multiply eigen values by coef
    … …  
    15771577                                                int   nbb=0;
    15781578                                                Real8 dd = detT[i];
    1579                                                 Real8 lla,llb,llc,llf;
    15801579                                                Real8 taa[3][3],bb[3];
    15811580
    … …  
    15861585                                                        Triangle*        tt = ta;
    15871586
     1587
     1588                                                        /* V is the vertex opposed to the edge shared by the
     1589                                                         * current triangle i and its neighbor ta
     1590                                                         *         C       
     1591                                                         *       / | \     
     1592                                                         *      /  |  \   
     1593                                                         *     /   |   \   
     1594                                                         * V  / ta | i  \ B
     1595                                                         *    \    |    / 
     1596                                                         *     \   |   /   
     1597                                                         *      \  |  /
     1598                                                         *       \ | /     
     1599                                                         *         A
     1600                                                         * lA=area(VBC)/Area(ABC)
     1601                                                         * lB=area(AVC)/Area(ABC)
     1602                                                         * lC=area(ABV)/Area(ABC)
     1603                                                         */
     1604
     1605                                                        //first, get V
     1606                                                        Vertex &v = *ta.OppositeVertex();
     1607                                                        Int4   iV = Number(v);
     1608
    15881609                                                        //if the adjacent triangle is not a boundary triangle:
    15891610                                                        if (tt && tt->link){
    15901611                                                                Vertex &v = *ta.OppositeVertex();
    1591                                                                 R2     V = v;
     1612                                                                R2      V = v;
    15921613                                                                Int4   iV = Number(v);
     1614
     1615                                                                //get lA lB and lC
    15931616                                                                Real8  lA = bamg::Area2(V,B,C)/dd;
    15941617                                                                Real8  lB = bamg::Area2(A,V,C)/dd;
    15951618                                                                Real8  lC = bamg::Area2(A,B,V)/dd;
     1619
     1620                                                                //fill taa
    15961621                                                                taa[0][j] = lB*lC;
    15971622                                                                taa[1][j] = lC*lA;
    15981623                                                                taa[2][j] = lA*lB;
    1599                                                                 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;
    1600                                                                 bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ;
     1624
     1625                                                                //fill bb
     1626                                                                bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ; //value of the solution on V minus...
    16011627                                                        }
    16021628                                                        else{
     1629                                                                //there might be a problem here: if 2 nodes are inline lA, lB or lC
     1630                                                                //is 0. This means that one column has only one term.
     1631                                                                //if instead of 0.1 we use 0, and taa[j][j]=1, we might have det(taa)=0
    16031632                                                                nbb++;
    1604                                                                 taa[0][j]=0;
    1605                                                                 taa[1][j]=0;
    1606                                                                 taa[2][j]=0;
     1633                                                                taa[0][j]=0.1;
     1634                                                                taa[1][j]=0.1;
     1635                                                                taa[2][j]=0.1;
    16071636                                                                taa[j][j]=1;
    16081637                                                                bb[j]=0;
    … …  
    16111640
    16121641                                                // resolution of 3x3 linear system transpose
    1613                                                 Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);           
     1642                                                Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);   
    16141643                                                Real8 cBC   =  det3x3(bb,taa[1],taa[2]);
    16151644                                                Real8 cCA   =  det3x3(taa[0],bb,taa[2]);
    16161645                                                Real8 cAB   =  det3x3(taa[0],taa[1],bb);
    1617 
    16181646                                                if (!det33){
    1619                                                         throw ErrorException(__FUNCT__,exprintf("!det33"));
     1647                                                        throw ErrorException(__FUNCT__,exprintf("for triangle %i, det33==0! cannot compute theHessian matrix using Hessiantype=1",i+1));
    16201648                                                }
     1649
    16211650                                                // computation of the Hessian in the element
    1622 
    16231651                                                // H( li*lj) = grad li grad lj + grad lj grad lj
    16241652                                                // grad li = njk  / detT ; with i j k =(A,B,C)
    … …  
    16731701                                for (int xy = 0;xy<3;xy++) {
    16741702                                        dd = d2[xy];
    1675                                         // do leat 2 iteration for boundary problem
     1703                                        // do least 2 iteration for boundary problem
    16761704                                        for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
    16771705                                                for (i=0;i<nbt;i++)
    … …  
    17461774
    17471775                                        //modify eigen values according to hmin and hmax
    1748                                         Vp.Maxh(hmin);
    1749                                         Vp.Minh(hmax);
     1776                                        Vp.Maxh(hmax);
     1777                                        Vp.Minh(hmin);
    17501778
    17511779                                        //multiply eigen values by coef
    … …  
    17891817        }
    17901818        /*}}}1*/
    1791         /*FUNCTION Triangles::BuildMetric2{{{1*/
     1819        /*FUNCTION Triangles::BuildMetric2 (double P2 projection){{{1*/
    17921820        void Triangles::BuildMetric2(BamgOpts* bamgopts){
    17931821
    1794                 throw ErrorException(__FUNCT__,exprintf("not supported yet"));
     1822                /*Options*/
     1823                const int dim = 2;
     1824                double* s=NULL;
     1825                Int4 nbsol;
     1826                int* typsols=NULL;
     1827                int verbosity;
     1828                Real8 hmin1;
     1829                Real8 hmax1;
     1830                int   i,j,k,iA,iB,iC;
     1831                int   iv,nbfield;
     1832
     1833                /*Recover options*/
     1834                verbosity=bamgopts->verbose;
     1835                hmin1=bamgopts->hmin;
     1836                hmax1=bamgopts->hmax;
     1837
     1838                /*Recover options*/
     1839                verbosity=bamgopts->verbose;
     1840                hmin1=bamgopts->hmin;
     1841                hmax1=bamgopts->hmax;
     1842
     1843                /*Get and process fields*/
     1844                s=bamgopts->field;
     1845                nbsol=1;       //for now, only one field
     1846                typsols=(int*)xmalloc(1*sizeof(int));
     1847                typsols[0]=0; // only one dof per node
     1848
     1849                //sizeoftype = {1,2,3,4}
     1850                int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ;
     1851
     1852                // computation of the number of fields
     1853                Int4 ntmp = 0;
     1854                if (typsols){
     1855                        //if there is more than one dof per vertex for one field
     1856                        //increase ntmp to take into account all the fields
     1857                        //if nbsol=1
     1858                        for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
     1859                }
     1860                else ntmp = nbsol;
     1861
     1862                //n is the total number of fields
     1863                const Int4 n = ntmp;
     1864
     1865                //initialization of some variables
     1866                Real8   hmin = Max(hmin1,MinimalHmin());
     1867                Real8   hmax = Min(hmax1,MaximalHmax());
     1868                double* ss=(double*)s;
     1869                double  sA,sB,sC;
     1870                Real8*  detT = new Real8[nbt];
     1871                Real8*  sumareas = new Real8[nbv];
     1872                Real8*  alpha= new Real8[nbt*3];
     1873                Real8*  beta = new Real8[nbt*3];
     1874                Real8*  dx_elem    = new Real8[nbt];
     1875                Real8*  dy_elem    = new Real8[nbt];
     1876                Real8*  dx_vertex  = new Real8[nbv];
     1877                Real8*  dy_vertex  = new Real8[nbv];
     1878                Real8*  dxdx_elem  = new Real8[nbt];
     1879                Real8*  dxdy_elem  = new Real8[nbt];
     1880                Real8*  dydy_elem  = new Real8[nbt];
     1881                Real8*  dxdx_vertex= new Real8[nbv];
     1882                Real8*  dxdy_vertex= new Real8[nbv];
     1883                Real8*  dydy_vertex= new Real8[nbv];
     1884
     1885                //display infos
     1886                if(verbosity>1) {
     1887                        printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
     1888                }
     1889
     1890                //first, build the chains that will be used for the Hessian computation, as weel as the area of each element
     1891                int head_s[nbv];
     1892                int next_p[3*nbt];
     1893                int  p=0;
     1894                //initialization
     1895                for(i=0;i<nbv;i++){
     1896                        sumareas[i]=0;
     1897                        head_s[i]=-1;
     1898                }
     1899                for(i=0;i<nbt;i++){
     1900
     1901                        //lopp over the real triangles (no boundary elements)
     1902                        if(triangles[i].link){
     1903
     1904                                //get current triangle t
     1905                                const Triangle &t=triangles[i];
     1906
     1907                                // coor of 3 vertices
     1908                                R2 A=t[0];
     1909                                R2 B=t[1];
     1910                                R2 C=t[2];
     1911
     1912                                //compute triangle determinant (2*Area)
     1913                                Real8 dett = bamg::Area2(A,B,C);
     1914                                detT[i]=dett;
     1915
     1916                                /*The nodal functions are such that for a vertex A:
     1917                                 *    N_A(x,y)=alphaA x + beta_A y +gamma_A
     1918                                 *    N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
     1919                                 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
     1920                                 * leads to:
     1921                                 *    N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
     1922                                 * and this gives:
     1923                                 *    alpha_A = (yB-yC)/(2*Area(T))*/
     1924                                alpha[i*3+0]=(B.y-C.y)/dett;
     1925                                alpha[i*3+1]=(C.y-A.y)/dett;
     1926                                alpha[i*3+2]=(A.y-B.y)/dett;
     1927                                beta[ i*3+0]=(C.x-B.x)/dett;
     1928                                beta[ i*3+1]=(A.x-C.x)/dett;
     1929                                beta[ i*3+2]=(B.x-A.x)/dett;
     1930
     1931                                //compute chains
     1932                                for(j=0;j<3;j++){
     1933                                        k=Number(triangles[i][j]);
     1934                                        next_p[p]=head_s[k];
     1935                                        head_s[k]=p++;
     1936
     1937                                        //add area to sumareas
     1938                                        sumareas[k]+=dett;
     1939                                }
     1940
     1941                        }
     1942                }
     1943
     1944                //for all Solutions
     1945                for (Int4 nusol=0;nusol<nbsol;nusol++) {
     1946                        int   nbfield=typsols?sizeoftype[typsols[nusol]]:1;
     1947
     1948                        //loop over all the fields of the solution
     1949                        for (Int4 nufield=0;nufield<nbfield;nufield++,s++){
     1950                                //ss++ so that for each iteration ss points toward the right field
     1951
     1952                                //initialize the hessian and gradient matrices
     1953                                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;
     1954
     1955                                //1: Compute gradient for each element (exact)
     1956                                for (i=0;i<nbt;i++){
     1957                                        if(triangles[i].link){
     1958                                                // number of the 3 vertices
     1959                                                iA = Number(triangles[i][0]);
     1960                                                iB = Number(triangles[i][1]);
     1961                                                iC = Number(triangles[i][2]);
     1962
     1963                                                // value of the P1 fonction on 3 vertices
     1964                                                sA = ss[iA*n];
     1965                                                sB = ss[iB*n];
     1966                                                sC = ss[iC*n];
     1967
     1968                                                //gradient = (sum alpha_i s_i, sum_i beta_i s_i)
     1969                                                dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
     1970                                                dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
     1971                                        }
     1972                                }
     1973
     1974                                //2: then compute a gradient for each vertex using a P2 projection
     1975                                for(i=0;i<nbv;i++){
     1976                                        for(p=head_s[i];p!=-1;p=next_p[p]){
     1977                                                //Get triangle number
     1978                                                k=(Int4)(p/3);
     1979                                                dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
     1980                                                dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
     1981                                        }
     1982                                }
     1983
     1984                                //3: compute Hessian matrix on each element
     1985                                for (i=0;i<nbt;i++){
     1986                                        if(triangles[i].link){
     1987                                                // number of the 3 vertices
     1988                                                iA = Number(triangles[i][0]);
     1989                                                iB = Number(triangles[i][1]);
     1990                                                iC = Number(triangles[i][2]);
     1991
     1992                                                //Hessian
     1993                                                dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
     1994                                                dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
     1995                                                dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
     1996                                        }
     1997                                }
     1998
     1999                                //4: finaly compute Hessian on each vertex using the second P2 projection
     2000                                for(i=0;i<nbv;i++){
     2001                                        for(p=head_s[i];p!=-1;p=next_p[p]){
     2002                                                //Get triangle number
     2003                                                k=(Int4)(p/3);
     2004                                                dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
     2005                                                dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
     2006                                                dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
     2007                                        }
     2008                                }
     2009
     2010
     2011                                /*Compute Metric from Hessian*/
     2012
     2013                                //compute multiplicative coefficient (2/9 because it is 2d)
     2014                                Real8 ci=(2.0/9.0)*(1.0/(bamgopts->err));
     2015
     2016                                for ( iv=0;iv<nbv;iv++){
     2017
     2018                                        //initialize metric Miv with ci*H
     2019                                        Metric Miv(dxdx_vertex[iv]*ci,dxdy_vertex[iv]*ci,dydy_vertex[iv]*ci);
     2020
     2021                                        //Get eigen values and vectors of Miv
     2022                                        MatVVP2x2 Vp(Miv);
     2023                                        if(iv==2){
     2024                                                printf("Hessien' = [%g %g %g]\n",dxdx_vertex[iv]*ci,dxdy_vertex[iv]*ci,dydy_vertex[iv]*ci);
     2025                                                Vp.Echo();
     2026                                        }
     2027
     2028                                        //move eigen valuse to their absolute values
     2029                                        Vp.Abs();
     2030
     2031                                        //modify eigen values according to hmin and hmax
     2032                                        Vp.Maxh(hmax);
     2033                                        Vp.Minh(hmin);
     2034
     2035                                        //rebuild Metric from Vp
     2036                                        Metric MVp(Vp);
     2037                                        if(iv==2){
     2038                                                MVp.Echo();
     2039                                                vertices[iv].m.Echo();
     2040                                        }
     2041
     2042                                        //Apply Metric to vertex
     2043                                        vertices[iv].m.IntersectWith(MVp);
     2044                                        if(iv==2){
     2045                                                vertices[iv].m.Echo();
     2046                                        }
     2047                                }
     2048
     2049                                for(i=0;i<3;i++){
     2050                                        vertices[i].m.Echo();
     2051                                }
     2052                        }//for all fields
     2053                }//for all solutions
     2054                for(i=0;i<5;i++){
     2055                        vertices[i].m.Echo();
     2056                }
     2057
     2058                //clean up
     2059                delete [] detT;
     2060                delete [] alpha;
     2061                delete [] beta;
     2062                delete [] sumareas;
     2063                delete [] dx_elem;
     2064                delete [] dy_elem;
     2065                delete [] dx_vertex;
     2066                delete [] dy_vertex;
     2067                delete [] dxdx_elem;
     2068                delete [] dxdy_elem;
     2069                delete [] dydy_elem;
     2070                delete [] dxdx_vertex;
     2071                delete [] dxdy_vertex;
     2072                delete [] dydy_vertex;
    17952073        }
    17962074        /*}}}1*/
    … …  
    36503928                BuildMetric1(bamgopts);
    36513929        }
     3930        else if (Hessiantype==2){
     3931                BuildMetric2(bamgopts);
     3932        }
    36523933        else{
    36533934                throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (0->use Green formula, 1-> from P2 on 4T, 2-> double P2 projection)",Hessiantype));
Note: See TracChangeset for help on using the changeset viewer.