Ignore:
Timestamp:
06/10/11 07:54:45 (14 years ago)
Author:
Mathieu Morlighem
Message:

Removed linearization of Cost function in CM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r8592 r8596  
    16551655        double     dux,duy,meanvel,epsvel;
    16561656        double     scalex=0,scaley=0,scale=0,S=0;
     1657        double     vx,vy,vxobs,vyobs,weight;
    16571658        double     xyz_list[NUMVERTICES][3];
    1658         double     vx_list[NUMVERTICES];
    1659         double     vy_list[NUMVERTICES];
    1660         double     obs_vx_list[NUMVERTICES];
    1661         double     obs_vy_list[NUMVERTICES];
    1662         double     dux_list[NUMVERTICES];
    1663         double     duy_list[NUMVERTICES];
    1664         double     weights_list[NUMVERTICES];
    16651659        double     l1l2l3[3];
    16661660        GaussTria* gauss=NULL;
    … …  
    16731667        this->parameters->FindParam(&meanvel,MeanVelEnum);
    16741668        this->parameters->FindParam(&epsvel,EpsVelEnum);
    1675         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    1676         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    1677         GetParameterListOnVertices(&vx_list[0],VxEnum);
    1678         GetParameterListOnVertices(&vy_list[0],VyEnum);
    1679         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
     1669        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     1670        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     1671        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     1672        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     1673        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
     1674
    16801675        inputs->GetParameterValue(&response,CmResponseEnum);
    1681         if(response==SurfaceAverageVelMisfitEnum){
    1682                 inputs->GetParameterValue(&S,SurfaceAreaEnum);
    1683         }
    1684 
    1685         /*Get Du at the 3 nodes (integration of the linearized function)
    1686          * Here we integrate linearized functions:
    1687          *               
    1688          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    1689          *
    1690          *       d J                  dJ_i
    1691          * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
    1692          *       d u                  du_i
    1693          *
    1694          * where J_i are the misfits at the 3 nodes of the triangle
    1695          *       Phi_i is the nodal function (P1) with respect to
    1696          *       the vertex i
    1697          */
    1698         if(response==SurfaceAbsVelMisfitEnum){
    1699                 /*We are using an absolute misfit:
    1700                  *
    1701                  *      1  [           2              2 ]
    1702                  * J = --- | (u - u   )  +  (v - v   )  |
    1703                  *      2  [       obs            obs   ]
    1704                  *
    1705                  *        dJ
    1706                  * DU = - -- = (u   - u )
    1707                  *        du     obs
    1708                  */
    1709                 for (i=0;i<NUMVERTICES;i++){
    1710                         dux_list[i]=obs_vx_list[i]-vx_list[i];
    1711                         duy_list[i]=obs_vy_list[i]-vy_list[i];
    1712                 }
    1713         }
    1714         else if(response==SurfaceRelVelMisfitEnum){
    1715                 /*We are using a relative misfit:
    1716                  *                       
    1717                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    1718                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    1719                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    1720                  *              obs                        obs                     
    1721                  *
    1722                  *        dJ     \bar{v}^2
    1723                  * DU = - -- = ------------- (u   - u )
    1724                  *        du   (u   + eps)^2    obs
    1725                  *               obs
    1726                  */
    1727                 for (i=0;i<NUMVERTICES;i++){
    1728                         scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
    1729                         scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
    1730                         if(obs_vx_list[i]==0)scalex=0;
    1731                         if(obs_vy_list[i]==0)scaley=0;
    1732                         dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
    1733                         duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
    1734                 }
    1735         }
    1736         else if(response==SurfaceLogVelMisfitEnum){
    1737                 /*We are using a logarithmic misfit:
    1738                  *                       
    1739                  *                 [        vel + eps     ] 2
    1740                  * J = 4 \bar{v}^2 | log ( -----------  ) | 
    1741                  *                 [       vel   + eps    ]
    1742                  *                            obs
    1743                  *
    1744                  *        dJ                 2 * log(...)
    1745                  * DU = - -- = - 4 \bar{v}^2 -------------  u
    1746                  *        du                 vel^2 + eps
    1747                  *           
    1748                  */
    1749                 for (i=0;i<NUMVERTICES;i++){
    1750                         velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
    1751                         obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
    1752                         scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    1753                         dux_list[i]=scale*vx_list[i];
    1754                         duy_list[i]=scale*vy_list[i];
    1755                 }
    1756         }
    1757         else if(response==SurfaceAverageVelMisfitEnum){
    1758                 /*We are using an spacially average absolute misfit:
    1759                  *
    1760                  *      1                    2              2
    1761                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    1762                  *      S                obs            obs
    1763                  *
    1764                  *        dJ      1       1
    1765                  * DU = - -- = - --- ----------- * 2 (u - u   )
    1766                  *        du      S  2 sqrt(...)           obs
    1767                  */
    1768                 for (i=0;i<NUMVERTICES;i++){
    1769                         scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
    1770                         dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
    1771                         duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
    1772                 }
    1773         }
    1774         else if(response==SurfaceLogVxVyMisfitEnum){
    1775                 /*We are using an logarithmic 2 misfit:
    1776                  *
    1777                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    1778                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    1779                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    1780                  *                              obs                       obs
    1781                  *        dJ                              1      u                             1
    1782                  * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    1783                  *        du                         |u| + eps  |u|                           u + eps
    1784                  */
    1785                 for (i=0;i<NUMVERTICES;i++){
    1786                         dux_list[i] = - pow(meanvel,(double)2)*(
    1787                                                 log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
    1788                         duy_list[i] = - pow(meanvel,(double)2)*(
    1789                                                 log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
    1790                 }
    1791         }
    1792         else{
    1793                 /*Not supported yet! : */
    1794                 _error_("response %s not supported yet",EnumToStringx(response));
    1795         }
    1796 
    1797         /*Apply weights to DU*/
    1798         for (i=0;i<NUMVERTICES;i++){
    1799                 dux_list[i]=weights_list[i]*dux_list[i];
    1800                 duy_list[i]=weights_list[i]*duy_list[i];
    1801         }
     1676        if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum);
    18021677
    18031678        /* Start  looping on the number of gaussian points: */
    1804         gauss=new GaussTria(2);
    1805         for(ig=gauss->begin();ig<gauss->end();ig++){
     1679        gauss=new GaussTria(4);
     1680        for (ig=gauss->begin();ig<gauss->end();ig++){
    18061681
    18071682                gauss->GaussPoint(ig);
    18081683
     1684                /* Get Jacobian determinant: */
    18091685                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1686
     1687                /*Get all parameters at gaussian point*/
     1688                weights_input->GetParameterValue(&weight,gauss,0);
     1689                vx_input->GetParameterValue(&vx,gauss);
     1690                vy_input->GetParameterValue(&vy,gauss);
     1691                vxobs_input->GetParameterValue(&vxobs,gauss);
     1692                vyobs_input->GetParameterValue(&vyobs,gauss);
    18101693                GetNodalFunctions(l1l2l3, gauss);
    18111694
    1812                 TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
    1813                 TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
    1814 
    1815                 for (i=0;i<NUMVERTICES;i++){
    1816                         pe->values[i*NDOF2+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
    1817                         pe->values[i*NDOF2+1]+=duy*Jdet*gauss->weight*l1l2l3[i];
     1695                if(response==SurfaceAbsVelMisfitEnum){
     1696                        /*
     1697                         *      1  [           2              2 ]
     1698                         * J = --- | (u - u   )  +  (v - v   )  |
     1699                         *      2  [       obs            obs   ]
     1700                         *
     1701                         *        dJ
     1702                         * DU = - -- = (u   - u )
     1703                         *        du     obs
     1704                         */
     1705                        for (i=0;i<NUMVERTICES;i++){
     1706                                dux=vxobs-vx;
     1707                                duy=vyobs-vy;
     1708                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1709                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1710                        }
     1711                }
     1712                else if(response==SurfaceRelVelMisfitEnum){
     1713                        /*
     1714                         *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     1715                         * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     1716                         *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     1717                         *              obs                        obs                     
     1718                         *
     1719                         *        dJ     \bar{v}^2
     1720                         * DU = - -- = ------------- (u   - u )
     1721                         *        du   (u   + eps)^2    obs
     1722                         *               obs
     1723                         */
     1724                        for (i=0;i<NUMVERTICES;i++){
     1725                                scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
     1726                                scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     1727                                dux=scalex*(vxobs-vx);
     1728                                duy=scaley*(vyobs-vy);
     1729                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1730                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1731                        }
     1732                }
     1733                else if(response==SurfaceLogVelMisfitEnum){
     1734                        /*
     1735                         *                 [        vel + eps     ] 2
     1736                         * J = 4 \bar{v}^2 | log ( -----------  ) | 
     1737                         *                 [       vel   + eps    ]
     1738                         *                            obs
     1739                         *
     1740                         *        dJ                 2 * log(...)
     1741                         * DU = - -- = - 4 \bar{v}^2 -------------  u
     1742                         *        du                 vel^2 + eps
     1743                         *           
     1744                         */
     1745                        for (i=0;i<NUMVERTICES;i++){
     1746                                velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
     1747                                obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
     1748                                scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
     1749                                dux=scale*vx;
     1750                                duy=scale*vy;
     1751                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1752                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1753                        }
     1754                }
     1755                else if(response==SurfaceAverageVelMisfitEnum){
     1756                        /*
     1757                         *      1                    2              2
     1758                         * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     1759                         *      S                obs            obs
     1760                         *
     1761                         *        dJ      1       1
     1762                         * DU = - -- = - --- ----------- * 2 (u - u   )
     1763                         *        du      S  2 sqrt(...)           obs
     1764                         */
     1765                        for (i=0;i<NUMVERTICES;i++){
     1766                                scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
     1767                                dux=scale*(vxobs-vx);
     1768                                duy=scale*(vyobs-vy);
     1769                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1770                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1771                        }
     1772                }
     1773                else if(response==SurfaceLogVxVyMisfitEnum){
     1774                        /*
     1775                         *      1            [        |u| + eps     2          |v| + eps     2  ]
     1776                         * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     1777                         *      2            [       |u    |+ eps              |v    |+ eps     ]
     1778                         *                              obs                       obs
     1779                         *        dJ                              1      u                             1
     1780                         * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     1781                         *        du                         |u| + eps  |u|                           u + eps
     1782                         */
     1783                        for (i=0;i<NUMVERTICES;i++){
     1784                                dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     1785                                duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     1786                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1787                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1788                        }
     1789                }
     1790                else{
     1791                        /*Not supported yet! : */
     1792                        _error_("response %s not supported yet",EnumToStringx(response));
    18181793                }
    18191794        }
    … …  
    18281803
    18291804        /*Intermediaries */
    1830         int        i,ig;
    1831         int        fit=-1;
    1832         int        response;
     1805        int        i,ig,response;
    18331806        double     Jdet;
    18341807        double     obs_velocity_mag,velocity_mag;
    18351808        double     dux,duy,meanvel,epsvel;
    18361809        double     scalex=0,scaley=0,scale=0,S=0;
     1810        double     vx,vy,vxobs,vyobs,weight;
    18371811        double     xyz_list[NUMVERTICES][3];
    1838         double     vx_list[NUMVERTICES];
    1839         double     vy_list[NUMVERTICES];
    1840         double     obs_vx_list[NUMVERTICES];
    1841         double     obs_vy_list[NUMVERTICES];
    1842         double     dux_list[NUMVERTICES];
    1843         double     duy_list[NUMVERTICES];
    1844         double     weights_list[NUMVERTICES];
    18451812        double     l1l2l3[3];
    18461813        GaussTria* gauss=NULL;
    … …  
    18531820        this->parameters->FindParam(&meanvel,MeanVelEnum);
    18541821        this->parameters->FindParam(&epsvel,EpsVelEnum);
    1855         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    1856         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    1857         GetParameterListOnVertices(&vx_list[0],VxEnum);
    1858         GetParameterListOnVertices(&vy_list[0],VyEnum);
    1859         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
     1822        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     1823        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     1824        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     1825        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     1826        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
     1827
    18601828        inputs->GetParameterValue(&response,CmResponseEnum);
    1861         if(response==SurfaceAverageVelMisfitEnum){
    1862                 inputs->GetParameterValue(&S,SurfaceAreaEnum);
    1863         }
    1864 
    1865         /*Get Du at the 3 nodes (integration of the linearized function)
    1866          * Here we integrate linearized functions:
    1867          *               
    1868          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    1869          *
    1870          *       d J                  dJ_i
    1871          * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
    1872          *       d u                  du_i
    1873          *
    1874          * where J_i are the misfits at the 3 nodes of the triangle
    1875          *       Phi_i is the nodal function (P1) with respect to
    1876          *       the vertex i
    1877          */
    1878         if(response==SurfaceAbsVelMisfitEnum){
    1879                 /*We are using an absolute misfit:
    1880                  *
    1881                  *      1  [           2              2 ]
    1882                  * J = --- | (u - u   )  +  (v - v   )  |
    1883                  *      2  [       obs            obs   ]
    1884                  *
    1885                  *        dJ             2
    1886                  * DU = - -- = (u   - u )
    1887                  *        du     obs
    1888                  */
    1889                 for (i=0;i<NUMVERTICES;i++){
    1890                         dux_list[i]=obs_vx_list[i]-vx_list[i];
    1891                         duy_list[i]=obs_vy_list[i]-vy_list[i];
    1892                 }
    1893         }
    1894         else if(response==SurfaceRelVelMisfitEnum){
    1895                 /*We are using a relative misfit:
    1896                  *                       
    1897                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    1898                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    1899                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    1900                  *              obs                        obs                     
    1901                  *
    1902                  *        dJ     \bar{v}^2
    1903                  * DU = - -- = ------------- (u   - u )
    1904                  *        du   (u   + eps)^2    obs
    1905                  *               obs
    1906                  */
    1907                 for (i=0;i<NUMVERTICES;i++){
    1908                         scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
    1909                         scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
    1910                         if(obs_vx_list[i]==0)scalex=0;
    1911                         if(obs_vy_list[i]==0)scaley=0;
    1912                         dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
    1913                         duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
    1914                 }
    1915         }
    1916         else if(response==SurfaceLogVelMisfitEnum){
    1917                 /*We are using a logarithmic misfit:
    1918                  *                       
    1919                  *                 [        vel + eps     ] 2
    1920                  * J = 4 \bar{v}^2 | log ( -----------  ) | 
    1921                  *                 [       vel   + eps    ]
    1922                  *                            obs
    1923                  *
    1924                  *        dJ                 2 * log(...)
    1925                  * DU = - -- = - 4 \bar{v}^2 -------------  u
    1926                  *        du                 vel^2 + eps
    1927                  *           
    1928                  */
    1929                 for (i=0;i<NUMVERTICES;i++){
    1930                         velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
    1931                         obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
    1932                         scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    1933                         dux_list[i]=scale*vx_list[i];
    1934                         duy_list[i]=scale*vy_list[i];
    1935                 }
    1936         }
    1937         else if(response==SurfaceAverageVelMisfitEnum){
    1938                 /*We are using an spacially average absolute misfit:
    1939                  *
    1940                  *      1                    2              2
    1941                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    1942                  *      S                obs            obs
    1943                  *
    1944                  *        dJ      1       1
    1945                  * DU = - -- = - --- ----------- * 2 (u - u   )
    1946                  *        du      S  2 sqrt(...)           obs
    1947                  */
    1948                 for (i=0;i<NUMVERTICES;i++){
    1949                         scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
    1950                         dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
    1951                         duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
    1952                 }
    1953         }
    1954         else if(response==SurfaceLogVxVyMisfitEnum){
    1955                 /*We are using an logarithmic 2 misfit:
    1956                  *
    1957                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    1958                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    1959                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    1960                  *                              obs                       obs
    1961                  *        dJ                              1      u                             1
    1962                  * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    1963                  *        du                         |u| + eps  |u|                           u + eps
    1964                  */
    1965                 for (i=0;i<NUMVERTICES;i++){
    1966                         dux_list[i] = - pow(meanvel,(double)2)*(
    1967                                                 log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
    1968                         duy_list[i] = - pow(meanvel,(double)2)*(
    1969                                                 log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
    1970                 }
    1971         }
    1972         else{
    1973                 /*Not supported yet! : */
    1974                 _error_("response %s not supported yet",EnumToStringx(response));
    1975         }
    1976 
    1977         /*Apply weights to DU*/
    1978         for (i=0;i<NUMVERTICES;i++){
    1979                 dux_list[i]=weights_list[i]*dux_list[i];
    1980                 duy_list[i]=weights_list[i]*duy_list[i];
    1981         }
     1829        if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum);
    19821830
    19831831        /* Start  looping on the number of gaussian points: */
    1984         gauss=new GaussTria(2);
    1985         for(ig=gauss->begin();ig<gauss->end();ig++){
     1832        gauss=new GaussTria(4);
     1833        for (ig=gauss->begin();ig<gauss->end();ig++){
    19861834
    19871835                gauss->GaussPoint(ig);
    19881836
     1837                /* Get Jacobian determinant: */
    19891838                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1839
     1840                /*Get all parameters at gaussian point*/
     1841                weights_input->GetParameterValue(&weight,gauss,0);
     1842                vx_input->GetParameterValue(&vx,gauss);
     1843                vy_input->GetParameterValue(&vy,gauss);
     1844                vxobs_input->GetParameterValue(&vxobs,gauss);
     1845                vyobs_input->GetParameterValue(&vyobs,gauss);
    19901846                GetNodalFunctions(l1l2l3, gauss);
    19911847
    1992                 TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
    1993                 TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
    1994 
    1995                 for (i=0;i<NUMVERTICES;i++){
    1996                         pe->values[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
    1997                         pe->values[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i];
     1848                if(response==SurfaceAbsVelMisfitEnum){
     1849                        /*
     1850                         *      1  [           2              2 ]
     1851                         * J = --- | (u - u   )  +  (v - v   )  |
     1852                         *      2  [       obs            obs   ]
     1853                         *
     1854                         *        dJ
     1855                         * DU = - -- = (u   - u )
     1856                         *        du     obs
     1857                         */
     1858                        for (i=0;i<NUMVERTICES;i++){
     1859                                dux=vxobs-vx;
     1860                                duy=vyobs-vy;
     1861                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1862                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1863                        }
     1864                }
     1865                else if(response==SurfaceRelVelMisfitEnum){
     1866                        /*
     1867                         *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     1868                         * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     1869                         *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     1870                         *              obs                        obs                     
     1871                         *
     1872                         *        dJ     \bar{v}^2
     1873                         * DU = - -- = ------------- (u   - u )
     1874                         *        du   (u   + eps)^2    obs
     1875                         *               obs
     1876                         */
     1877                        for (i=0;i<NUMVERTICES;i++){
     1878                                scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
     1879                                scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     1880                                dux=scalex*(vxobs-vx);
     1881                                duy=scaley*(vyobs-vy);
     1882                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1883                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1884                        }
     1885                }
     1886                else if(response==SurfaceLogVelMisfitEnum){
     1887                        /*
     1888                         *                 [        vel + eps     ] 2
     1889                         * J = 4 \bar{v}^2 | log ( -----------  ) | 
     1890                         *                 [       vel   + eps    ]
     1891                         *                            obs
     1892                         *
     1893                         *        dJ                 2 * log(...)
     1894                         * DU = - -- = - 4 \bar{v}^2 -------------  u
     1895                         *        du                 vel^2 + eps
     1896                         *           
     1897                         */
     1898                        for (i=0;i<NUMVERTICES;i++){
     1899                                velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
     1900                                obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
     1901                                scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
     1902                                dux=scale*vx;
     1903                                duy=scale*vy;
     1904                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1905                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1906                        }
     1907                }
     1908                else if(response==SurfaceAverageVelMisfitEnum){
     1909                        /*
     1910                         *      1                    2              2
     1911                         * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     1912                         *      S                obs            obs
     1913                         *
     1914                         *        dJ      1       1
     1915                         * DU = - -- = - --- ----------- * 2 (u - u   )
     1916                         *        du      S  2 sqrt(...)           obs
     1917                         */
     1918                        for (i=0;i<NUMVERTICES;i++){
     1919                                scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
     1920                                dux=scale*(vxobs-vx);
     1921                                duy=scale*(vyobs-vy);
     1922                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1923                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1924                        }
     1925                }
     1926                else if(response==SurfaceLogVxVyMisfitEnum){
     1927                        /*
     1928                         *      1            [        |u| + eps     2          |v| + eps     2  ]
     1929                         * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     1930                         *      2            [       |u    |+ eps              |v    |+ eps     ]
     1931                         *                              obs                       obs
     1932                         *        dJ                              1      u                             1
     1933                         * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     1934                         *        du                         |u| + eps  |u|                           u + eps
     1935                         */
     1936                        for (i=0;i<NUMVERTICES;i++){
     1937                                dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     1938                                duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     1939                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1940                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1941                        }
     1942                }
     1943                else{
     1944                        /*Not supported yet! : */
     1945                        _error_("response %s not supported yet",EnumToStringx(response));
    19981946                }
    19991947        }
    … …  
    48454793        int        i,ig;
    48464794        double     Jelem=0;
    4847         double     velocity_mag,obs_velocity_mag;
    4848         double     meanvel, epsvel,misfit,Jdet;
     4795        double     misfit,Jdet;
     4796        double     vx,vy,vxobs,vyobs,weight;
    48494797        double     xyz_list[NUMVERTICES][3];
    4850         double     vx_list[NUMVERTICES];
    4851         double     vy_list[NUMVERTICES];
    4852         double     obs_vx_list[NUMVERTICES];
    4853         double     obs_vy_list[NUMVERTICES];
    4854         double     misfit_square_list[NUMVERTICES];
    4855         double     misfit_list[NUMVERTICES];
    4856         double     weights_list[NUMVERTICES];
    48574798        GaussTria *gauss=NULL;
    48584799
    … …  
    48634804        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    48644805
    4865         /* Recover input data: */
    4866         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    4867         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    4868         GetParameterListOnVertices(&vx_list[0],VxEnum);
    4869         GetParameterListOnVertices(&vy_list[0],VyEnum);
    4870         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
    4871 
    4872         /*retrieve some parameters: */
    4873         this->parameters->FindParam(&meanvel,MeanVelEnum);
    4874         this->parameters->FindParam(&epsvel,EpsVelEnum);
    4875        
    4876         /* Compute SurfaceAbsVelMisfit at the 3 nodes
    4877          * Here we integrate linearized functions:
    4878          *               
    4879          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    4880          *
    4881          * where J_i are the misfits at the 3 nodes of the triangle
    4882          *       Phi_i is the nodal function (P1) with respect to
    4883          *       the vertex i
    4884          */
    4885 
    4886         /*We are using an absolute misfit:
    4887          *
    4888          *      1  [           2              2 ]
    4889          * J = --- | (u - u   )  +  (v - v   )  |
    4890          *      2  [       obs            obs   ]
    4891          *
    4892          */
    4893         for (i=0;i<NUMVERTICES;i++){
    4894                 misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
    4895         }
    4896         /*Process units: */
    4897         if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters);
    4898 
    4899         /*Apply weights to misfits*/
    4900         for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
     4806        /*Retrieve all inputs we will be needing: */
     4807        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     4808        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     4809        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     4810        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     4811        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
    49014812
    49024813        /* Start  looping on the number of gaussian points: */
    49034814        gauss=new GaussTria(2);
    4904         for (ig=gauss->begin();ig<gauss->end(); ig++){
     4815        for (ig=gauss->begin();ig<gauss->end();ig++){
    49054816
    49064817                gauss->GaussPoint(ig);
    49074818
     4819                /* Get Jacobian determinant: */
    49084820                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    4909                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
    4910                 Jelem+=misfit*Jdet*gauss->weight;
     4821
     4822                /*Get all parameters at gaussian point*/
     4823                weights_input->GetParameterValue(&weight,gauss,0);
     4824                vx_input->GetParameterValue(&vx,gauss);
     4825                vy_input->GetParameterValue(&vy,gauss);
     4826                vxobs_input->GetParameterValue(&vxobs,gauss);
     4827                vyobs_input->GetParameterValue(&vyobs,gauss);
     4828
     4829                /*Compute SurfaceAbsVelMisfitEnum:
     4830                 *
     4831                 *      1  [           2              2 ]
     4832                 * J = --- | (u - u   )  +  (v - v   )  |
     4833                 *      2  [       obs            obs   ]
     4834                 *
     4835                 */
     4836                misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );
     4837
     4838                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
     4839
     4840                /*Add to cost function*/
     4841                Jelem+=misfit*weight*Jdet*gauss->weight;
    49114842        }
    49124843
    … …  
    49514882
    49524883        int        i,ig;
    4953         int        fit=-1;
    4954         double     Jelem=0,S=0;
    4955         double     scalex=1, scaley=1;
    4956         double     meanvel, epsvel,Jdet;
    4957         double     velocity_mag,obs_velocity_mag,misfit;
     4884        double     Jelem=0,S,Jdet;
     4885        double     misfit;
     4886        double     vx,vy,vxobs,vyobs,weight;
    49584887        double     xyz_list[NUMVERTICES][3];
    4959         double     vx_list[NUMVERTICES];
    4960         double     vy_list[NUMVERTICES];
    4961         double     obs_vx_list[NUMVERTICES];
    4962         double     obs_vy_list[NUMVERTICES];
    4963         double     misfit_square_list[NUMVERTICES];
    4964         double     misfit_list[NUMVERTICES];
    4965         double     weights_list[NUMVERTICES];
    49664888        GaussTria *gauss=NULL;
    49674889
    … …  
    49724894        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    49734895
    4974         /* Recover input data: */
     4896        /*Retrieve all inputs we will be needing: */
    49754897        inputs->GetParameterValue(&S,SurfaceAreaEnum);
    4976         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    4977         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    4978         GetParameterListOnVertices(&vx_list[0],VxEnum);
    4979         GetParameterListOnVertices(&vy_list[0],VyEnum);
    4980         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
    4981 
    4982         /*retrieve some parameters: */
    4983         this->parameters->FindParam(&meanvel,MeanVelEnum);
    4984         this->parameters->FindParam(&epsvel,EpsVelEnum);
    4985 
    4986         /* Compute SurfaceAverageVelMisfit at the 3 nodes
    4987          * Here we integrate linearized functions:
    4988          *               
    4989          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    4990          *
    4991          * where J_i are the misfits at the 3 nodes of the triangle
    4992          *       Phi_i is the nodal function (P1) with respect to
    4993          *       the vertex i
    4994          */
    4995 
    4996         /*We are using a spacially average absolute misfit:
    4997          *
    4998          *      1                    2              2
    4999          * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    5000          *      S                obs            obs
    5001          */
    5002         for (i=0;i<NUMVERTICES;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
    5003 
    5004         /*Process units: */
    5005         if(process_units)UnitConversion(&misfit_square_list[0],NUMVERTICES,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
    5006 
    5007         /*Take the square root, and scale by surface: */
    5008         for (i=0;i<NUMVERTICES;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
    5009 
    5010         /*Apply weights to misfits*/
    5011         for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
     4898        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     4899        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     4900        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     4901        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     4902        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
    50124903
    50134904        /* Start  looping on the number of gaussian points: */
    5014         gauss=new GaussTria(2);
     4905        gauss=new GaussTria(3);
    50154906        for (ig=gauss->begin();ig<gauss->end();ig++){
    50164907
    50174908                gauss->GaussPoint(ig);
    50184909
     4910                /* Get Jacobian determinant: */
    50194911                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    5020                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
    5021                 Jelem+=misfit*Jdet*gauss->weight;
     4912
     4913                /*Get all parameters at gaussian point*/
     4914                weights_input->GetParameterValue(&weight,gauss,0);
     4915                vx_input->GetParameterValue(&vx,gauss);
     4916                vy_input->GetParameterValue(&vy,gauss);
     4917                vxobs_input->GetParameterValue(&vxobs,gauss);
     4918                vyobs_input->GetParameterValue(&vyobs,gauss);
     4919
     4920                /*Compute SurfaceAverageVelMisfitEnum:
     4921                 *
     4922                 *      1                    2              2
     4923                 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     4924                 *      S                obs            obs
     4925                 */
     4926                misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);
     4927
     4928                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
     4929
     4930                /*Add to cost function*/
     4931                Jelem+=misfit*weight*Jdet*gauss->weight;
    50224932        }
    50234933
    … …  
    50314941
    50324942        const int    numdof=NDOF2*NUMVERTICES;
     4943
     4944        int        i,ig;
     4945        double     Jelem=0;
     4946        double     meanvel, epsvel,misfit,Jdet;
     4947        double     velocity_mag,obs_velocity_mag;
     4948        double     xyz_list[NUMVERTICES][3];
     4949        double     vx,vy,vxobs,vyobs,weight;
     4950        GaussTria *gauss=NULL;
     4951
     4952        /*If on water, return 0: */
     4953        if(IsOnWater())return 0;
     4954
     4955        /* Get node coordinates and dof list: */
     4956        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4957
     4958        /*Retrieve all inputs we will be needing: */
     4959        this->parameters->FindParam(&meanvel,MeanVelEnum);
     4960        this->parameters->FindParam(&epsvel,EpsVelEnum);
     4961        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     4962        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     4963        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     4964        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     4965        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
     4966
     4967        /* Start  looping on the number of gaussian points: */
     4968        gauss=new GaussTria(4);
     4969        for (ig=gauss->begin();ig<gauss->end();ig++){
     4970
     4971                gauss->GaussPoint(ig);
     4972
     4973                /* Get Jacobian determinant: */
     4974                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     4975
     4976                /*Get all parameters at gaussian point*/
     4977                weights_input->GetParameterValue(&weight,gauss,0);
     4978                vx_input->GetParameterValue(&vx,gauss);
     4979                vy_input->GetParameterValue(&vy,gauss);
     4980                vxobs_input->GetParameterValue(&vxobs,gauss);
     4981                vyobs_input->GetParameterValue(&vyobs,gauss);
     4982
     4983                /*Compute SurfaceLogVelMisfit:
     4984                 *                 [        vel + eps     ] 2
     4985                 * J = 4 \bar{v}^2 | log ( -----------  ) | 
     4986                 *                 [       vel   + eps    ]
     4987                 *                            obs
     4988                 */
     4989                velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
     4990                obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
     4991                misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);
     4992
     4993                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
     4994
     4995                /*Add to cost function*/
     4996                Jelem+=misfit*weight*Jdet*gauss->weight;
     4997        }
     4998
     4999        /*clean-up and Return: */
     5000        delete gauss;
     5001        return Jelem;
     5002}
     5003/*}}}*/
     5004/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
     5005double Tria::SurfaceLogVxVyMisfit(bool process_units){
     5006
     5007        const int    numdof=NDOF2*NUMVERTICES;
     5008
     5009        int        i,ig;
     5010        int        fit=-1;
     5011        double     Jelem=0, S=0;
     5012        double     meanvel, epsvel, misfit, Jdet;
     5013        double     vx,vy,vxobs,vyobs,weight;
     5014        double     xyz_list[NUMVERTICES][3];
     5015        GaussTria *gauss=NULL;
     5016
     5017        /*If on water, return 0: */
     5018        if(IsOnWater())return 0;
     5019
     5020        /* Get node coordinates and dof list: */
     5021        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     5022
     5023        /*Retrieve all inputs we will be needing: */
     5024        this->parameters->FindParam(&meanvel,MeanVelEnum);
     5025        this->parameters->FindParam(&epsvel,EpsVelEnum);
     5026        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     5027        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     5028        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     5029        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     5030        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
     5031       
     5032        /* Start  looping on the number of gaussian points: */
     5033        gauss=new GaussTria(4);
     5034        for (ig=gauss->begin();ig<gauss->end();ig++){
     5035
     5036                gauss->GaussPoint(ig);
     5037
     5038                /* Get Jacobian determinant: */
     5039                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     5040
     5041                /*Get all parameters at gaussian point*/
     5042                weights_input->GetParameterValue(&weight,gauss,0);
     5043                vx_input->GetParameterValue(&vx,gauss);
     5044                vy_input->GetParameterValue(&vy,gauss);
     5045                vxobs_input->GetParameterValue(&vxobs,gauss);
     5046                vyobs_input->GetParameterValue(&vyobs,gauss);
     5047
     5048                /*Compute SurfaceRelVelMisfit:
     5049                 *
     5050                 *      1            [        |u| + eps     2          |v| + eps     2  ]
     5051                 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     5052                 *      2            [       |u    |+ eps              |v    |+ eps     ]
     5053                 *                              obs                       obs
     5054                 */
     5055                misfit=0.5*pow(meanvel,2.)*(
     5056                                        pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
     5057                                        pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );
     5058
     5059                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
     5060
     5061                /*Add to cost function*/
     5062                Jelem+=misfit*weight*Jdet*gauss->weight;
     5063        }
     5064
     5065        /*clean-up and Return: */
     5066        delete gauss;
     5067        return Jelem;
     5068}
     5069/*}}}*/
     5070/*FUNCTION Tria::SurfaceNormal{{{1*/
     5071void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
     5072
     5073        int i;
     5074        double v13[3],v23[3];
     5075        double normal[3];
     5076        double normal_norm;
     5077
     5078        for (i=0;i<3;i++){
     5079                v13[i]=xyz_list[0][i]-xyz_list[2][i];
     5080                v23[i]=xyz_list[1][i]-xyz_list[2][i];
     5081        }
     5082
     5083        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
     5084        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
     5085        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
     5086
     5087        normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
     5088
     5089        *(surface_normal)=normal[0]/normal_norm;
     5090        *(surface_normal+1)=normal[1]/normal_norm;
     5091        *(surface_normal+2)=normal[2]/normal_norm;
     5092}
     5093/*}}}*/
     5094/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
     5095double Tria::SurfaceRelVelMisfit(bool process_units){
     5096        const int  numdof=2*NUMVERTICES;
    50335097
    50345098        int        i,ig;
    … …  
    50365100        double     scalex=1,scaley=1;
    50375101        double     meanvel, epsvel,misfit,Jdet;
    5038         double     velocity_mag,obs_velocity_mag;
     5102        double     vx,vy,vxobs,vyobs,weight;
    50395103        double     xyz_list[NUMVERTICES][3];
    5040         double     vx_list[NUMVERTICES];
    5041         double     vy_list[NUMVERTICES];
    5042         double     obs_vx_list[NUMVERTICES];
    5043         double     obs_vy_list[NUMVERTICES];
    5044         double     misfit_square_list[NUMVERTICES];
    5045         double     misfit_list[NUMVERTICES];
    5046         double     weights_list[NUMVERTICES];
    50475104        GaussTria *gauss=NULL;
    50485105
    … …  
    50535110        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    50545111
    5055         /* Recover input data: */
    5056         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    5057         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    5058         GetParameterListOnVertices(&vx_list[0],VxEnum);
    5059         GetParameterListOnVertices(&vy_list[0],VyEnum);
    5060         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
    5061 
    5062         /*retrieve some parameters: */
     5112        /*Retrieve all inputs we will be needing: */
    50635113        this->parameters->FindParam(&meanvel,MeanVelEnum);
    50645114        this->parameters->FindParam(&epsvel,EpsVelEnum);
    5065        
    5066         /* Compute SurfaceLogVelMisfit at the 3 nodes
    5067          * Here we integrate linearized functions:
    5068          *               
    5069          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    5070          *
    5071          * where J_i are the misfits at the 3 nodes of the triangle
    5072          *       Phi_i is the nodal function (P1) with respect to
    5073          *       the vertex i
    5074          */
    5075 
    5076         /*We are using a logarithmic misfit:
    5077          *                       
    5078          *                 [        vel + eps     ] 2
    5079          * J = 4 \bar{v}^2 | log ( -----------  ) | 
    5080          *                 [       vel   + eps    ]
    5081          *                            obs
    5082          */
    5083         for (i=0;i<NUMVERTICES;i++){
    5084                 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
    5085                 obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
    5086                 misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
    5087         }
    5088 
    5089         /*Process units: */
    5090         if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
    5091 
    5092         /*Apply weights to misfits*/
    5093         for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
     5115        Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
     5116        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
     5117        Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
     5118        Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
     5119        Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
    50945120
    50955121        /* Start  looping on the number of gaussian points: */
    5096         gauss=new GaussTria(2);
     5122        gauss=new GaussTria(4);
    50975123        for (ig=gauss->begin();ig<gauss->end();ig++){
    50985124
    50995125                gauss->GaussPoint(ig);
    51005126
     5127                /* Get Jacobian determinant: */
    51015128                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    5102                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
    5103                 Jelem+=misfit*Jdet*gauss->weight;
    5104         }
    5105 
    5106         /*clean-up and Return: */
    5107         delete gauss;
    5108         return Jelem;
    5109 }
    5110 /*}}}*/
    5111 /*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
    5112 double Tria::SurfaceLogVxVyMisfit(bool process_units){
    5113 
    5114         const int    numdof=NDOF2*NUMVERTICES;
    5115 
    5116         int        i,ig;
    5117         int        fit=-1;
    5118         double     Jelem=0, S=0;
    5119         double     scalex=1,scaley=1;
    5120         double     meanvel, epsvel, misfit, Jdet;
    5121         double     velocity_mag,obs_velocity_mag;
    5122         double     xyz_list[NUMVERTICES][3];
    5123         double     vx_list[NUMVERTICES];
    5124         double     vy_list[NUMVERTICES];
    5125         double     obs_vx_list[NUMVERTICES];
    5126         double     obs_vy_list[NUMVERTICES];
    5127         double     misfit_square_list[NUMVERTICES];
    5128         double     misfit_list[NUMVERTICES];
    5129         double     weights_list[NUMVERTICES];
    5130         GaussTria *gauss=NULL;
    5131 
    5132         /*If on water, return 0: */
    5133         if(IsOnWater())return 0;
    5134 
    5135         /* Get node coordinates and dof list: */
    5136         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    5137 
    5138         /* Recover input data: */
    5139         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    5140         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    5141         GetParameterListOnVertices(&vx_list[0],VxEnum);
    5142         GetParameterListOnVertices(&vy_list[0],VyEnum);
    5143         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
    5144 
    5145         /*retrieve some parameters: */
    5146         this->parameters->FindParam(&meanvel,MeanVelEnum);
    5147         this->parameters->FindParam(&epsvel,EpsVelEnum);
    5148        
    5149         /* Compute SurfaceLogVxVyMisfit at the 3 nodes
    5150          * Here we integrate linearized functions:
    5151          *               
    5152          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    5153          *
    5154          * where J_i are the misfits at the 3 nodes of the triangle
    5155          *       Phi_i is the nodal function (P1) with respect to
    5156          *       the vertex i
    5157          */
    5158 
    5159         /*We are using an logarithmic 2 misfit:
    5160          *
    5161          *      1            [        |u| + eps     2          |v| + eps     2  ]
    5162          * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    5163          *      2            [       |u    |+ eps              |v    |+ eps     ]
    5164          *                              obs                       obs
    5165          */
    5166         for (i=0;i<NUMVERTICES;i++){
    5167                 misfit_list[i]=0.5*pow(meanvel,(double)2)*(
    5168                                         pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
    5169                                         pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) );
    5170         }
    5171 
    5172         /*Process units: */
    5173         if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
    5174 
    5175         /*Apply weights to misfits*/
    5176         for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    5177 
    5178         /* Start  looping on the number of gaussian points: */
    5179         gauss=new GaussTria(2);
    5180         for (ig=gauss->begin();ig<gauss->end();ig++){
    5181 
    5182                 gauss->GaussPoint(ig);
    5183 
    5184                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    5185                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
    5186                 Jelem+=misfit*Jdet*gauss->weight;
    5187         }
    5188 
    5189         /*clean-up and Return: */
    5190         delete gauss;
    5191         return Jelem;
    5192 }
    5193 /*}}}*/
    5194 /*FUNCTION Tria::SurfaceNormal{{{1*/
    5195 void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
    5196 
    5197         int i;
    5198         double v13[3],v23[3];
    5199         double normal[3];
    5200         double normal_norm;
    5201 
    5202         for (i=0;i<3;i++){
    5203                 v13[i]=xyz_list[0][i]-xyz_list[2][i];
    5204                 v23[i]=xyz_list[1][i]-xyz_list[2][i];
    5205         }
    5206 
    5207         normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
    5208         normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    5209         normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    5210 
    5211         normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
    5212 
    5213         *(surface_normal)=normal[0]/normal_norm;
    5214         *(surface_normal+1)=normal[1]/normal_norm;
    5215         *(surface_normal+2)=normal[2]/normal_norm;
    5216 }
    5217 /*}}}*/
    5218 /*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
    5219 double Tria::SurfaceRelVelMisfit(bool process_units){
    5220 
    5221         const int    numdof=2*NUMVERTICES;
    5222 
    5223         int        i,ig;
    5224         double     Jelem=0;
    5225         double     scalex=1,scaley=1;
    5226         double     meanvel, epsvel,misfit,Jdet;
    5227         double     velocity_mag,obs_velocity_mag;
    5228         double     xyz_list[NUMVERTICES][3];
    5229         double     vx_list[NUMVERTICES];
    5230         double     vy_list[NUMVERTICES];
    5231         double     obs_vx_list[NUMVERTICES];
    5232         double     obs_vy_list[NUMVERTICES];
    5233         double     misfit_square_list[NUMVERTICES];
    5234         double     misfit_list[NUMVERTICES];
    5235         double     weights_list[NUMVERTICES];
    5236         GaussTria *gauss=NULL;
    5237 
    5238         /*If on water, return 0: */
    5239         if(IsOnWater())return 0;
    5240 
    5241         /* Get node coordinates and dof list: */
    5242         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    5243 
    5244         /* Recover input data: */
    5245         GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
    5246         GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    5247         GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    5248         GetParameterListOnVertices(&vx_list[0],VxEnum);
    5249         GetParameterListOnVertices(&vy_list[0],VyEnum);
    5250 
    5251         /*retrieve some parameters: */
    5252         this->parameters->FindParam(&meanvel,MeanVelEnum);
    5253         this->parameters->FindParam(&epsvel,EpsVelEnum);
    5254 
    5255         /* Compute SurfaceRelVelMisfit at the 3 nodes
    5256          * Here we integrate linearized functions:
    5257          *               
    5258          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    5259          *
    5260          * where J_i are the misfits at the 3 nodes of the triangle
    5261          *       Phi_i is the nodal function (P1) with respect to
    5262          *       the vertex i
    5263          */
    5264 
    5265         /*We are using a relative misfit:
    5266          *                       
    5267          *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    5268          * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    5269          *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    5270          *              obs                        obs                     
    5271          */
    5272         for (i=0;i<NUMVERTICES;i++){
    5273                 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
    5274                 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
    5275                 if(obs_vx_list[i]==0)scalex=0;
    5276                 if(obs_vy_list[i]==0)scaley=0;
    5277                 misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2));
    5278         }
    5279 
    5280         /*Process units: */
    5281         if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
    5282 
    5283         /*Apply weights to misfits*/
    5284         for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    5285 
    5286         /* Start  looping on the number of gaussian points: */
    5287         gauss=new GaussTria(2);
    5288         for (ig=gauss->begin();ig<gauss->end();ig++){
    5289 
    5290                 gauss->GaussPoint(ig);
    5291 
    5292                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    5293                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
    5294                 Jelem+=misfit*Jdet*gauss->weight;
     5129
     5130                /*Get all parameters at gaussian point*/
     5131                weights_input->GetParameterValue(&weight,gauss,0);
     5132                vx_input->GetParameterValue(&vx,gauss);
     5133                vy_input->GetParameterValue(&vy,gauss);
     5134                vxobs_input->GetParameterValue(&vxobs,gauss);
     5135                vyobs_input->GetParameterValue(&vyobs,gauss);
     5136
     5137                /*Compute SurfaceRelVelMisfit:
     5138                 *                       
     5139                 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     5140                 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     5141                 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     5142                 *              obs                        obs                     
     5143                 */
     5144                scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
     5145                scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     5146                misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
     5147                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
     5148
     5149                /*Add to cost function*/
     5150                Jelem+=misfit*weight*Jdet*gauss->weight;
    52955151        }
    52965152
Note: See TracChangeset for help on using the changeset viewer.