![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 8603
- Timestamp:
- 06/10/11 14:58:43 (14 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r8601 r8603 1599 1599 1600 1600 /*Intermediaries */ 1601 int i,ig ;1601 int i,ig,resp; 1602 1602 double Jdet; 1603 1603 double thickness,thicknessobs,weight; 1604 int *responses = NULL; 1605 int num_responses; 1604 1606 double xyz_list[NUMVERTICES][3]; 1605 1607 double l1l2l3[3]; … … 1613 1615 /*Retrieve all inputs and parameters*/ 1614 1616 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1615 Input* thickness_input =inputs->GetInput(ThicknessEnum); _assert_(thickness_input); 1616 Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);_assert_(thicknessobs_input); 1617 Input* weights_input =inputs->GetInput(WeightsEnum); _assert_(weights_input); 1617 this->parameters->FindParam(&num_responses,NumResponsesEnum); 1618 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 1619 Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); 1620 Input* thicknessobs_input = inputs->GetInput(ThicknessObsEnum);_assert_(thicknessobs_input); 1621 Input* weights_input = inputs->GetInput(WeightsEnum); _assert_(weights_input); 1618 1622 1619 1623 /* Start looping on the number of gaussian points: */ … … 1632 1636 weights_input->GetParameterValue(&weight, gauss,0); 1633 1637 1634 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i]; 1635 /*Regularization of the constraint: 2000000 79 N*/ 1636 //for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[0]*dbasis[0][i]*Jdet*gauss->weight; 1637 //for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[1]*dbasis[1][i]*Jdet*gauss->weight; 1638 /*Loop over all requested responses*/ 1639 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 1640 1641 case ThicknessAbsMisfitEnum: 1642 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i]; 1643 /*Regularization of the constraint: 2000000 79 N*/ 1644 //for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[0]*dbasis[0][i]*Jdet*gauss->weight; 1645 //for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[1]*dbasis[1][i]*Jdet*gauss->weight; 1646 break; 1647 default: 1648 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1649 } 1638 1650 } 1639 1651 … … 1650 1662 1651 1663 /*Intermediaries */ 1652 int i, ig;1664 int i,resp,ig; 1653 1665 int *responses=NULL; 1654 int response;1666 int num_responses; 1655 1667 double Jdet; 1656 1668 double obs_velocity_mag,velocity_mag; … … 1667 1679 /*Retrieve all inputs and parameters*/ 1668 1680 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1681 this->parameters->FindParam(&num_responses,NumResponsesEnum); 1682 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 1669 1683 this->parameters->FindParam(&meanvel,MeanVelEnum); 1670 1684 this->parameters->FindParam(&epsvel,EpsVelEnum); 1671 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); response=responses[0];1672 1685 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 1673 1686 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); … … 1676 1689 Input* vyobs_input =inputs->GetInput(VyObsEnum); _assert_(vyobs_input); 1677 1690 1678 if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum); 1691 /*Get Surface if required by one response*/ 1692 for(resp=0;resp<num_responses;resp++){ 1693 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 1694 inputs->GetParameterValue(&S,SurfaceAreaEnum); break; 1695 } 1696 } 1679 1697 1680 1698 /* Start looping on the number of gaussian points: */ … … 1695 1713 GetNodalFunctions(l1l2l3, gauss); 1696 1714 1697 if(response==SurfaceAbsVelMisfitEnum){ 1698 /* 1699 * 1 [ 2 2 ] 1700 * J = --- | (u - u ) + (v - v ) | 1701 * 2 [ obs obs ] 1702 * 1703 * dJ 1704 * DU = - -- = (u - u ) 1705 * du obs 1706 */ 1707 for (i=0;i<NUMVERTICES;i++){ 1708 dux=vxobs-vx; 1709 duy=vyobs-vy; 1710 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1711 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1712 } 1713 } 1714 else if(response==SurfaceRelVelMisfitEnum){ 1715 /* 1716 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1717 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1718 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1719 * obs obs 1720 * 1721 * dJ \bar{v}^2 1722 * DU = - -- = ------------- (u - u ) 1723 * du (u + eps)^2 obs 1724 * obs 1725 */ 1726 for (i=0;i<NUMVERTICES;i++){ 1727 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1728 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1729 dux=scalex*(vxobs-vx); 1730 duy=scaley*(vyobs-vy); 1731 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1732 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1733 } 1734 } 1735 else if(response==SurfaceLogVelMisfitEnum){ 1736 /* 1737 * [ vel + eps ] 2 1738 * J = 4 \bar{v}^2 | log ( ----------- ) | 1739 * [ vel + eps ] 1740 * obs 1741 * 1742 * dJ 2 * log(...) 1743 * DU = - -- = - 4 \bar{v}^2 ------------- u 1744 * du vel^2 + eps 1745 * 1746 */ 1747 for (i=0;i<NUMVERTICES;i++){ 1748 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1749 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1750 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1751 dux=scale*vx; 1752 duy=scale*vy; 1753 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1754 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1755 } 1756 } 1757 else if(response==SurfaceAverageVelMisfitEnum){ 1758 /* 1759 * 1 2 2 1760 * J = --- sqrt( (u - u ) + (v - v ) ) 1761 * S obs obs 1762 * 1763 * dJ 1 1 1764 * DU = - -- = - --- ----------- * 2 (u - u ) 1765 * du S 2 sqrt(...) obs 1766 */ 1767 for (i=0;i<NUMVERTICES;i++){ 1768 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1769 dux=scale*(vxobs-vx); 1770 duy=scale*(vyobs-vy); 1771 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1772 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1773 } 1774 } 1775 else if(response==SurfaceLogVxVyMisfitEnum){ 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 = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1787 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1788 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1789 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1790 } 1791 } 1792 else{ 1793 /*Not supported yet! : */ 1794 _error_("response %s not supported yet",EnumToStringx(response)); 1715 /*Loop over all requested responses*/ 1716 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 1717 1718 case SurfaceAbsVelMisfitEnum: 1719 /* 1720 * 1 [ 2 2 ] 1721 * J = --- | (u - u ) + (v - v ) | 1722 * 2 [ obs obs ] 1723 * 1724 * dJ 1725 * DU = - -- = (u - u ) 1726 * du obs 1727 */ 1728 for (i=0;i<NUMVERTICES;i++){ 1729 dux=vxobs-vx; 1730 duy=vyobs-vy; 1731 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1732 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1733 } 1734 break; 1735 case SurfaceRelVelMisfitEnum: 1736 /* 1737 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1738 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1739 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1740 * obs obs 1741 * 1742 * dJ \bar{v}^2 1743 * DU = - -- = ------------- (u - u ) 1744 * du (u + eps)^2 obs 1745 * obs 1746 */ 1747 for (i=0;i<NUMVERTICES;i++){ 1748 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1749 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1750 dux=scalex*(vxobs-vx); 1751 duy=scaley*(vyobs-vy); 1752 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1753 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1754 } 1755 break; 1756 case SurfaceLogVelMisfitEnum: 1757 /* 1758 * [ vel + eps ] 2 1759 * J = 4 \bar{v}^2 | log ( ----------- ) | 1760 * [ vel + eps ] 1761 * obs 1762 * 1763 * dJ 2 * log(...) 1764 * DU = - -- = - 4 \bar{v}^2 ------------- u 1765 * du vel^2 + eps 1766 * 1767 */ 1768 for (i=0;i<NUMVERTICES;i++){ 1769 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1770 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1771 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1772 dux=scale*vx; 1773 duy=scale*vy; 1774 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1775 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1776 } 1777 break; 1778 case SurfaceAverageVelMisfitEnum: 1779 /* 1780 * 1 2 2 1781 * J = --- sqrt( (u - u ) + (v - v ) ) 1782 * S obs obs 1783 * 1784 * dJ 1 1 1785 * DU = - -- = - --- ----------- * 2 (u - u ) 1786 * du S 2 sqrt(...) obs 1787 */ 1788 for (i=0;i<NUMVERTICES;i++){ 1789 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1790 dux=scale*(vxobs-vx); 1791 duy=scale*(vyobs-vy); 1792 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1793 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1794 } 1795 break; 1796 case SurfaceLogVxVyMisfitEnum: 1797 /* 1798 * 1 [ |u| + eps 2 |v| + eps 2 ] 1799 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1800 * 2 [ |u |+ eps |v |+ eps ] 1801 * obs obs 1802 * dJ 1 u 1 1803 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1804 * du |u| + eps |u| u + eps 1805 */ 1806 for (i=0;i<NUMVERTICES;i++){ 1807 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1808 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1809 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1810 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1811 } 1812 break; 1813 default: 1814 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1795 1815 } 1796 1816 } … … 1806 1826 1807 1827 /*Intermediaries */ 1808 int i, ig;1828 int i,resp,ig; 1809 1829 int *responses=NULL; 1810 int response;1830 int num_responses; 1811 1831 double Jdet; 1812 1832 double obs_velocity_mag,velocity_mag; … … 1823 1843 /*Retrieve all inputs and parameters*/ 1824 1844 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1845 this->parameters->FindParam(&num_responses,NumResponsesEnum); 1846 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 1825 1847 this->parameters->FindParam(&meanvel,MeanVelEnum); 1826 1848 this->parameters->FindParam(&epsvel,EpsVelEnum); 1827 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); response=responses[0]; 1828 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 1829 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); 1830 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); 1831 Input* vxobs_input =inputs->GetInput(VxObsEnum); _assert_(vxobs_input); 1832 Input* vyobs_input =inputs->GetInput(VyObsEnum); _assert_(vyobs_input); 1833 1834 if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum); 1849 Input* weights_input = inputs->GetInput(WeightsEnum); _assert_(weights_input); 1850 Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); 1851 Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); 1852 Input* vxobs_input = inputs->GetInput(VxObsEnum); _assert_(vxobs_input); 1853 Input* vyobs_input = inputs->GetInput(VyObsEnum); _assert_(vyobs_input); 1854 1855 /*Get Surface if required by one response*/ 1856 for(resp=0;resp<num_responses;resp++){ 1857 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 1858 inputs->GetParameterValue(&S,SurfaceAreaEnum); break; 1859 } 1860 } 1835 1861 1836 1862 /* Start looping on the number of gaussian points: */ … … 1851 1877 GetNodalFunctions(l1l2l3, gauss); 1852 1878 1853 if(response==SurfaceAbsVelMisfitEnum){ 1854 /* 1855 * 1 [ 2 2 ] 1856 * J = --- | (u - u ) + (v - v ) | 1857 * 2 [ obs obs ] 1858 * 1859 * dJ 1860 * DU = - -- = (u - u ) 1861 * du obs 1862 */ 1863 for (i=0;i<NUMVERTICES;i++){ 1864 dux=vxobs-vx; 1865 duy=vyobs-vy; 1866 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1867 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1868 } 1869 } 1870 else if(response==SurfaceRelVelMisfitEnum){ 1871 /* 1872 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1873 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1874 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1875 * obs obs 1876 * 1877 * dJ \bar{v}^2 1878 * DU = - -- = ------------- (u - u ) 1879 * du (u + eps)^2 obs 1880 * obs 1881 */ 1882 for (i=0;i<NUMVERTICES;i++){ 1883 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1884 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1885 dux=scalex*(vxobs-vx); 1886 duy=scaley*(vyobs-vy); 1887 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1888 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1889 } 1890 } 1891 else if(response==SurfaceLogVelMisfitEnum){ 1892 /* 1893 * [ vel + eps ] 2 1894 * J = 4 \bar{v}^2 | log ( ----------- ) | 1895 * [ vel + eps ] 1896 * obs 1897 * 1898 * dJ 2 * log(...) 1899 * DU = - -- = - 4 \bar{v}^2 ------------- u 1900 * du vel^2 + eps 1901 * 1902 */ 1903 for (i=0;i<NUMVERTICES;i++){ 1904 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1905 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1906 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1907 dux=scale*vx; 1908 duy=scale*vy; 1909 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1910 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1911 } 1912 } 1913 else if(response==SurfaceAverageVelMisfitEnum){ 1914 /* 1915 * 1 2 2 1916 * J = --- sqrt( (u - u ) + (v - v ) ) 1917 * S obs obs 1918 * 1919 * dJ 1 1 1920 * DU = - -- = - --- ----------- * 2 (u - u ) 1921 * du S 2 sqrt(...) obs 1922 */ 1923 for (i=0;i<NUMVERTICES;i++){ 1924 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1925 dux=scale*(vxobs-vx); 1926 duy=scale*(vyobs-vy); 1927 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1928 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1929 } 1930 } 1931 else if(response==SurfaceLogVxVyMisfitEnum){ 1932 /* 1933 * 1 [ |u| + eps 2 |v| + eps 2 ] 1934 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1935 * 2 [ |u |+ eps |v |+ eps ] 1936 * obs obs 1937 * dJ 1 u 1 1938 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1939 * du |u| + eps |u| u + eps 1940 */ 1941 for (i=0;i<NUMVERTICES;i++){ 1942 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1943 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1944 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1945 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1946 } 1947 } 1948 else{ 1949 /*Not supported yet! : */ 1950 _error_("response %s not supported yet",EnumToStringx(response)); 1879 /*Loop over all requested responses*/ 1880 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 1881 1882 case SurfaceAbsVelMisfitEnum: 1883 /* 1884 * 1 [ 2 2 ] 1885 * J = --- | (u - u ) + (v - v ) | 1886 * 2 [ obs obs ] 1887 * 1888 * dJ 1889 * DU = - -- = (u - u ) 1890 * du obs 1891 */ 1892 for (i=0;i<NUMVERTICES;i++){ 1893 dux=vxobs-vx; 1894 duy=vyobs-vy; 1895 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1896 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1897 } 1898 break; 1899 case SurfaceRelVelMisfitEnum: 1900 /* 1901 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1902 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1903 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1904 * obs obs 1905 * 1906 * dJ \bar{v}^2 1907 * DU = - -- = ------------- (u - u ) 1908 * du (u + eps)^2 obs 1909 * obs 1910 */ 1911 for (i=0;i<NUMVERTICES;i++){ 1912 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1913 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1914 dux=scalex*(vxobs-vx); 1915 duy=scaley*(vyobs-vy); 1916 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1917 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1918 } 1919 break; 1920 case SurfaceLogVelMisfitEnum: 1921 /* 1922 * [ vel + eps ] 2 1923 * J = 4 \bar{v}^2 | log ( ----------- ) | 1924 * [ vel + eps ] 1925 * obs 1926 * 1927 * dJ 2 * log(...) 1928 * DU = - -- = - 4 \bar{v}^2 ------------- u 1929 * du vel^2 + eps 1930 * 1931 */ 1932 for (i=0;i<NUMVERTICES;i++){ 1933 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1934 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1935 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1936 dux=scale*vx; 1937 duy=scale*vy; 1938 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1939 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1940 } 1941 break; 1942 case SurfaceAverageVelMisfitEnum: 1943 /* 1944 * 1 2 2 1945 * J = --- sqrt( (u - u ) + (v - v ) ) 1946 * S obs obs 1947 * 1948 * dJ 1 1 1949 * DU = - -- = - --- ----------- * 2 (u - u ) 1950 * du S 2 sqrt(...) obs 1951 */ 1952 for (i=0;i<NUMVERTICES;i++){ 1953 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1954 dux=scale*(vxobs-vx); 1955 duy=scale*(vyobs-vy); 1956 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1957 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1958 } 1959 break; 1960 case SurfaceLogVxVyMisfitEnum: 1961 /* 1962 * 1 [ |u| + eps 2 |v| + eps 2 ] 1963 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1964 * 2 [ |u |+ eps |v |+ eps ] 1965 * obs obs 1966 * dJ 1 u 1 1967 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1968 * du |u| + eps |u| u + eps 1969 */ 1970 for (i=0;i<NUMVERTICES;i++){ 1971 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1972 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1973 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1974 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1975 } 1976 break; 1977 default: 1978 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1951 1979 } 1952 1980 }
Note:
See TracChangeset
for help on using the changeset viewer.
