![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 8596
- Timestamp:
- 06/10/11 07:54:45 (14 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r8592 r8596 1655 1655 double dux,duy,meanvel,epsvel; 1656 1656 double scalex=0,scaley=0,scale=0,S=0; 1657 double vx,vy,vxobs,vyobs,weight; 1657 1658 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];1665 1659 double l1l2l3[3]; 1666 1660 GaussTria* gauss=NULL; … … 1673 1667 this->parameters->FindParam(&meanvel,MeanVelEnum); 1674 1668 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 1680 1675 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); 1802 1677 1803 1678 /* 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++){ 1806 1681 1807 1682 gauss->GaussPoint(ig); 1808 1683 1684 /* Get Jacobian determinant: */ 1809 1685 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); 1810 1693 GetNodalFunctions(l1l2l3, gauss); 1811 1694 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)); 1818 1793 } 1819 1794 } … … 1828 1803 1829 1804 /*Intermediaries */ 1830 int i,ig; 1831 int fit=-1; 1832 int response; 1805 int i,ig,response; 1833 1806 double Jdet; 1834 1807 double obs_velocity_mag,velocity_mag; 1835 1808 double dux,duy,meanvel,epsvel; 1836 1809 double scalex=0,scaley=0,scale=0,S=0; 1810 double vx,vy,vxobs,vyobs,weight; 1837 1811 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];1845 1812 double l1l2l3[3]; 1846 1813 GaussTria* gauss=NULL; … … 1853 1820 this->parameters->FindParam(&meanvel,MeanVelEnum); 1854 1821 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 1860 1828 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); 1982 1830 1983 1831 /* 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++){ 1986 1834 1987 1835 gauss->GaussPoint(ig); 1988 1836 1837 /* Get Jacobian determinant: */ 1989 1838 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); 1990 1846 GetNodalFunctions(l1l2l3, gauss); 1991 1847 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)); 1998 1946 } 1999 1947 } … … 4845 4793 int i,ig; 4846 4794 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; 4849 4797 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];4857 4798 GaussTria *gauss=NULL; 4858 4799 … … 4863 4804 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 4864 4805 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); 4901 4812 4902 4813 /* Start looping on the number of gaussian points: */ 4903 4814 gauss=new GaussTria(2); 4904 for (ig=gauss->begin();ig<gauss->end(); 4815 for (ig=gauss->begin();ig<gauss->end();ig++){ 4905 4816 4906 4817 gauss->GaussPoint(ig); 4907 4818 4819 /* Get Jacobian determinant: */ 4908 4820 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; 4911 4842 } 4912 4843 … … 4951 4882 4952 4883 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; 4958 4887 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];4966 4888 GaussTria *gauss=NULL; 4967 4889 … … 4972 4894 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 4973 4895 4974 /* Recover input data: */4896 /*Retrieve all inputs we will be needing: */ 4975 4897 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); 5012 4903 5013 4904 /* Start looping on the number of gaussian points: */ 5014 gauss=new GaussTria( 2);4905 gauss=new GaussTria(3); 5015 4906 for (ig=gauss->begin();ig<gauss->end();ig++){ 5016 4907 5017 4908 gauss->GaussPoint(ig); 5018 4909 4910 /* Get Jacobian determinant: */ 5019 4911 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; 5022 4932 } 5023 4933 … … 5031 4941 5032 4942 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*/ 5005 double 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*/ 5071 void 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*/ 5095 double Tria::SurfaceRelVelMisfit(bool process_units){ 5096 const int numdof=2*NUMVERTICES; 5033 5097 5034 5098 int i,ig; … … 5036 5100 double scalex=1,scaley=1; 5037 5101 double meanvel, epsvel,misfit,Jdet; 5038 double v elocity_mag,obs_velocity_mag;5102 double vx,vy,vxobs,vyobs,weight; 5039 5103 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];5047 5104 GaussTria *gauss=NULL; 5048 5105 … … 5053 5110 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 5054 5111 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: */ 5063 5113 this->parameters->FindParam(&meanvel,MeanVelEnum); 5064 5114 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); 5094 5120 5095 5121 /* Start looping on the number of gaussian points: */ 5096 gauss=new GaussTria( 2);5122 gauss=new GaussTria(4); 5097 5123 for (ig=gauss->begin();ig<gauss->end();ig++){ 5098 5124 5099 5125 gauss->GaussPoint(ig); 5100 5126 5127 /* Get Jacobian determinant: */ 5101 5128 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; 5295 5151 } 5296 5152
Note:
See TracChangeset
for help on using the changeset viewer.
