Changeset 6216
- Timestamp:
- 10/08/10 17:03:53 (14 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
r6213 r6216 58 58 if(iomodel->my_elements[i]){ 59 59 element=(Element*)elements->GetObjectByOffset(counter); 60 element-> Update(i,iomodel); //we need i to index into elements.60 element->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements. 61 61 62 62 material=(Material*)materials->GetObjectByOffset(counter); 63 material-> Update(i,iomodel); //we need i to index into elements.63 material->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements. 64 64 counter++; 65 65 } -
issm/trunk/src/c/objects/Elements/Element.h
r6213 r6216 56 56 virtual void DeleteResults(void)=0; 57 57 virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; 58 virtual void Update(int index, IoModel* iomodel)=0;59 58 virtual void UpdateGeometry(void)=0; 60 59 virtual void InputToResult(int enum_type,int step,double time)=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r6213 r6216 470 470 void Penta::InputUpdateFromVectorDakota(bool* vector, int name, int type){ 471 471 ISSMERROR(" not supported yet!"); 472 } 473 /*}}}*/ 474 /*FUNCTION Penta::InputUpdateFromIoModel(int index,IoModel* iomodel) {{{1*/ 475 void Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ 476 477 /*Intermediaries*/ 478 IssmInt i,j; 479 int penta_vertex_ids[6]; 480 double nodeinputs[6]; 481 482 /*Checks if debuging*/ 483 /*{{{2*/ 484 ISSMASSERT(iomodel->elements); 485 /*}}}*/ 486 487 /*Recover vertices ids needed to initialize inputs*/ 488 for(i=0;i<6;i++){ 489 penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab 490 } 491 492 //add as many inputs per element as requested: 493 if (iomodel->thickness) { 494 for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1]; 495 this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); 496 } 497 if (iomodel->surface) { 498 for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1]; 499 this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); 500 } 501 if (iomodel->bed) { 502 for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1]; 503 this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs)); 504 } 505 if (iomodel->drag_coefficient) { 506 for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1]; 507 this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs)); 508 509 if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index])); 510 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); 511 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); 512 513 } 514 if (iomodel->melting_rate) { 515 for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts; 516 this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); 517 } 518 if (iomodel->accumulation_rate) { 519 for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts; 520 this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); 521 } 522 if (iomodel->geothermalflux) { 523 for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1]; 524 this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs)); 525 } 526 if (iomodel->pressure) { 527 for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1]; 528 this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); 529 } 530 if (iomodel->temperature) { 531 for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1]; 532 this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs)); 533 } 534 if (iomodel->dhdt) { 535 for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts; 536 this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs)); 537 } 538 /*vx,vy and vz: */ 539 if (iomodel->vx) { 540 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts; 541 this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); 542 this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); 543 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs)); 544 } 545 if (iomodel->vy) { 546 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts; 547 this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); 548 this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); 549 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs)); 550 } 551 if (iomodel->vz) { 552 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts; 553 this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); 554 this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs)); 555 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs)); 556 } 557 if (iomodel->vx_obs) { 558 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; 559 this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs)); 560 } 561 if (iomodel->vy_obs) { 562 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; 563 this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs)); 564 } 565 if (iomodel->vz_obs) { 566 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; 567 this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs)); 568 } 569 if (iomodel->weights) { 570 for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1]; 571 this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs)); 572 } 573 if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); 574 if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); 575 if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); 576 if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); 577 578 /*Control Inputs*/ 579 if (iomodel->control_analysis && iomodel->control_type){ 580 for(i=0;i<iomodel->num_control_type;i++){ 581 switch((int)iomodel->control_type[i]){ 582 case DhDtEnum: 583 if (iomodel->dhdt){ 584 for(j=0;j<6;j++)nodeinputs[j]=iomodel->dhdt[penta_vertex_ids[j]-1]/iomodel->yts; 585 this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs)); 586 } 587 break; 588 case VxEnum: 589 if (iomodel->vx){ 590 for(j=0;j<6;j++)nodeinputs[j]=iomodel->vx[penta_vertex_ids[j]-1]/iomodel->yts; 591 this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs)); 592 } 593 break; 594 case VyEnum: 595 if (iomodel->vy){ 596 for(j=0;j<6;j++)nodeinputs[j]=iomodel->vy[penta_vertex_ids[j]-1]/iomodel->yts; 597 this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs)); 598 } 599 break; 600 case DragCoefficientEnum: 601 if (iomodel->drag_coefficient){ 602 for(j=0;j<6;j++)nodeinputs[j]=iomodel->drag_coefficient[penta_vertex_ids[j]-1]; 603 this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs)); 604 } 605 break; 606 case RheologyBbarEnum: 607 /*Matice will take care of it*/ break; 608 default: 609 ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i])); 610 } 611 } 612 } 613 614 //Need to know the type of approximation for this element 615 if(iomodel->elements_type){ 616 if (*(iomodel->elements_type+index)==MacAyealApproximationEnum){ 617 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum)); 618 } 619 else if (*(iomodel->elements_type+index)==PattynApproximationEnum){ 620 this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); 621 } 622 else if (*(iomodel->elements_type+index)==MacAyealPattynApproximationEnum){ 623 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); 624 } 625 else if (*(iomodel->elements_type+index)==HutterApproximationEnum){ 626 this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); 627 } 628 else if (*(iomodel->elements_type+index)==StokesApproximationEnum){ 629 this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); 630 } 631 else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){ 632 this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum)); 633 } 634 else if (*(iomodel->elements_type+index)==NoneApproximationEnum){ 635 this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); 636 } 637 else{ 638 ISSMERROR("Approximation type %s not supported yet",EnumToString((int)*(iomodel->elements_type+index))); 639 } 640 } 641 472 642 } 473 643 /*}}}*/ … … 1740 1910 1741 1911 /*Fill with IoModel*/ 1742 this-> Update(index,iomodel);1912 this->InputUpdateFromIoModel(index,iomodel); 1743 1913 1744 1914 /*Defaults if not provided in iomodel*/ … … 1817 1987 } 1818 1988 /*}}}*/ 1819 /*FUNCTION Penta::Update(int index,IoModel* iomodel) {{{1*/1820 void Penta::Update(int index,IoModel* iomodel){1821 1822 /*Intermediaries*/1823 IssmInt i,j;1824 int penta_vertex_ids[6];1825 double nodeinputs[6];1826 1827 /*Checks if debuging*/1828 /*{{{2*/1829 ISSMASSERT(iomodel->elements);1830 /*}}}*/1831 1832 /*Recover vertices ids needed to initialize inputs*/1833 for(i=0;i<6;i++){1834 penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab1835 }1836 1837 //add as many inputs per element as requested:1838 if (iomodel->thickness) {1839 for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1];1840 this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));1841 }1842 if (iomodel->surface) {1843 for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1];1844 this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));1845 }1846 if (iomodel->bed) {1847 for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1];1848 this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs));1849 }1850 if (iomodel->drag_coefficient) {1851 for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1];1852 this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs));1853 1854 if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));1855 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));1856 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));1857 1858 }1859 if (iomodel->melting_rate) {1860 for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts;1861 this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));1862 }1863 if (iomodel->accumulation_rate) {1864 for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts;1865 this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));1866 }1867 if (iomodel->geothermalflux) {1868 for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1];1869 this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs));1870 }1871 if (iomodel->pressure) {1872 for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1];1873 this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));1874 }1875 if (iomodel->temperature) {1876 for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1];1877 this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs));1878 }1879 if (iomodel->dhdt) {1880 for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts;1881 this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs));1882 }1883 /*vx,vy and vz: */1884 if (iomodel->vx) {1885 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts;1886 this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));1887 this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));1888 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs));1889 }1890 if (iomodel->vy) {1891 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts;1892 this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));1893 this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));1894 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs));1895 }1896 if (iomodel->vz) {1897 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts;1898 this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));1899 this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));1900 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs));1901 }1902 if (iomodel->vx_obs) {1903 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;1904 this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs));1905 }1906 if (iomodel->vy_obs) {1907 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;1908 this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs));1909 }1910 if (iomodel->vz_obs) {1911 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;1912 this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs));1913 }1914 if (iomodel->weights) {1915 for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];1916 this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs));1917 }1918 if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));1919 if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));1920 if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));1921 if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));1922 1923 /*Control Inputs*/1924 if (iomodel->control_analysis && iomodel->control_type){1925 for(i=0;i<iomodel->num_control_type;i++){1926 switch((int)iomodel->control_type[i]){1927 case DhDtEnum:1928 if (iomodel->dhdt){1929 for(j=0;j<6;j++)nodeinputs[j]=iomodel->dhdt[penta_vertex_ids[j]-1]/iomodel->yts;1930 this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs));1931 }1932 break;1933 case VxEnum:1934 if (iomodel->vx){1935 for(j=0;j<6;j++)nodeinputs[j]=iomodel->vx[penta_vertex_ids[j]-1]/iomodel->yts;1936 this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs));1937 }1938 break;1939 case VyEnum:1940 if (iomodel->vy){1941 for(j=0;j<6;j++)nodeinputs[j]=iomodel->vy[penta_vertex_ids[j]-1]/iomodel->yts;1942 this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs));1943 }1944 break;1945 case DragCoefficientEnum:1946 if (iomodel->drag_coefficient){1947 for(j=0;j<6;j++)nodeinputs[j]=iomodel->drag_coefficient[penta_vertex_ids[j]-1];1948 this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs));1949 }1950 break;1951 case RheologyBbarEnum:1952 /*Matice will take care of it*/ break;1953 default:1954 ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i]));1955 }1956 }1957 }1958 1959 //Need to know the type of approximation for this element1960 if(iomodel->elements_type){1961 if (*(iomodel->elements_type+index)==MacAyealApproximationEnum){1962 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum));1963 }1964 else if (*(iomodel->elements_type+index)==PattynApproximationEnum){1965 this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum));1966 }1967 else if (*(iomodel->elements_type+index)==MacAyealPattynApproximationEnum){1968 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum));1969 }1970 else if (*(iomodel->elements_type+index)==HutterApproximationEnum){1971 this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum));1972 }1973 else if (*(iomodel->elements_type+index)==StokesApproximationEnum){1974 this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum));1975 }1976 else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){1977 this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum));1978 }1979 else if (*(iomodel->elements_type+index)==NoneApproximationEnum){1980 this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum));1981 }1982 else{1983 ISSMERROR("Approximation type %s not supported yet",EnumToString((int)*(iomodel->elements_type+index)));1984 }1985 }1986 1987 }1988 /*}}}*/1989 1989 /*FUNCTION Penta::UpdateGeometry{{{1*/ 1990 1990 void Penta::UpdateGeometry(void){ -
issm/trunk/src/c/objects/Elements/Penta.h
r6213 r6216 69 69 void InputUpdateFromVectorDakota(double* vector, int name, int type); 70 70 void InputUpdateFromVectorDakota(int* vector, int name, int type); 71 void InputUpdateFromIoModel(int index, IoModel* iomodel); 71 72 /*}}}*/ 72 73 /*Element virtual functions definitions: {{{1*/ … … 119 120 double SurfaceArea(void); 120 121 void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); 121 void Update(int index, IoModel* iomodel);122 122 void UpdateGeometry(void); 123 123 double TimeAdapt(); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r6213 r6216 535 535 } 536 536 /*}}}*/ 537 538 /*Element virtual functions definitions: */ 539 /*FUNCTION Tria::AverageOntoPartition {{{1*/ 540 void Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){ 541 542 bool already=false; 543 int i,j; 544 int partition[NUMVERTICES]; 545 int offset[NUMVERTICES]; 546 double area; 547 double mean; 548 double values[3]; 549 550 /*First, get the area: */ 551 area=this->GetArea(); 552 553 /*Figure out the average for this element: */ 554 this->GetDofList1(&offset[0]); 555 mean=0; 556 for(i=0;i<NUMVERTICES;i++){ 557 partition[i]=(int)qmu_part[offset[i]]; 558 mean=mean+1.0/NUMVERTICES*vertex_response[offset[i]]; 559 } 560 561 /*Add contribution: */ 562 for(i=0;i<NUMVERTICES;i++){ 563 already=false; 564 for(j=0;j<i;j++){ 565 if (partition[i]==partition[j]){ 566 already=true; 567 break; 568 } 569 } 570 if(!already){ 571 VecSetValue(partition_contributions,partition[i],mean*area,ADD_VALUES); 572 VecSetValue(partition_areas,partition[i],area,ADD_VALUES); 573 }; 574 } 575 } 576 /*}}}*/ 577 /*FUNCTION Tria::ComputeBasalStress {{{1*/ 578 void Tria::ComputeBasalStress(Vec eps){ 579 ISSMERROR("Not Implemented yet"); 580 } 581 /*}}}*/ 582 /*FUNCTION Tria::ComputeStrainRate {{{1*/ 583 void Tria::ComputeStrainRate(Vec eps){ 584 ISSMERROR("Not Implemented yet"); 585 } 586 /*}}}*/ 587 /*FUNCTION Tria::SetCurrentConfiguration {{{1*/ 588 void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ 589 590 /*go into parameters and get the analysis_counter: */ 591 int analysis_counter; 592 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); 593 594 /*Get Element type*/ 595 this->element_type=this->element_type_list[analysis_counter]; 596 597 /*Pick up nodes*/ 598 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 599 else this->nodes=NULL; 600 601 } 602 /*}}}*/ 603 /*FUNCTION Tria::Configure {{{1*/ 604 void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ 605 606 /*go into parameters and get the analysis_counter: */ 607 int analysis_counter; 608 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); 609 610 /*Get Element type*/ 611 this->element_type=this->element_type_list[analysis_counter]; 612 613 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 614 * datasets, using internal ids and offsets hidden in hooks: */ 615 if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); 616 this->hmatice->configure(materialsin); 617 this->hmatpar->configure(materialsin); 618 619 /*Now, go pick up the objects inside the hooks: */ 620 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 621 else this->nodes=NULL; 622 this->matice=(Matice*)this->hmatice->delivers(); 623 this->matpar=(Matpar*)this->hmatpar->delivers(); 624 625 /*point parameters to real dataset: */ 626 this->parameters=parametersin; 627 628 } 629 /*}}}*/ 630 /*FUNCTION Tria::ControlInputGetGradient{{{1*/ 631 void Tria::ControlInputGetGradient(Vec gradient,int enum_type){ 632 633 int doflist1[NUMVERTICES]; 634 Input* input=NULL; 635 636 if(enum_type==RheologyBbarEnum){ 637 input=(Input*)matice->inputs->GetInput(enum_type); 638 } 639 else{ 640 input=inputs->GetInput(enum_type); 641 } 642 if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); 643 if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); 644 645 this->GetDofList1(&doflist1[0]); 646 ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); 647 648 }/*}}}*/ 649 /*FUNCTION Tria::ControlInputSetGradient{{{1*/ 650 void Tria::ControlInputSetGradient(double* gradient,int enum_type){ 651 652 int doflist1[NUMVERTICES]; 653 double grad_list[NUMVERTICES]; 654 Input* grad_input=NULL; 655 Input* input=NULL; 656 657 if(enum_type==RheologyBbarEnum){ 658 input=(Input*)matice->inputs->GetInput(enum_type); 659 } 660 else{ 661 input=inputs->GetInput(enum_type); 662 } 663 if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); 664 if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); 665 666 this->GetDofList1(&doflist1[0]); 667 for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]]; 668 grad_input=new TriaVertexInput(GradientEnum,grad_list); 669 670 ((ControlInput*)input)->SetGradient(grad_input); 671 672 }/*}}}*/ 673 /*FUNCTION Tria::RegularizeInversion {{{1*/ 674 double Tria::RegularizeInversion(){ 675 676 /*constants: */ 677 const int numdof=NDOF2*NUMVERTICES; 678 679 /* Intermediaries */ 680 int i,j,ig; 681 int num_controls; 682 int *control_type=NULL; 683 double Jelem = 0; 684 double cm_noisedmp; 685 double Jdet; 686 double xyz_list[NUMVERTICES][3]; 687 double dk[NDOF2],dB[NDOF2]; 688 GaussTria *gauss = NULL; 689 690 /*retrieve parameters and inputs*/ 691 this->parameters->FindParam(&num_controls,NumControlsEnum); 692 this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); 693 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); 694 695 /*If on water, return 0: */ 696 if(IsOnWater()) return 0; 697 698 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 699 700 for(i=0;i<num_controls;i++){ 701 /* Start looping on the number of gaussian points: */ 702 gauss=new GaussTria(2); 703 for (ig=gauss->begin();ig<gauss->end();ig++){ 704 705 gauss->GaussPoint(ig); 706 707 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 708 709 /*Add Tikhonov regularization term to misfit: 710 * 711 * WARNING: for now, the regularization is only taken into account by the gradient 712 * the misfit reflects the response only! 713 * 714 * */ 715 if (control_type[i]==DragCoefficientEnum){ 716 Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); 717 drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); 718 //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; 719 } 720 else if (control_type[i]==RheologyBbarEnum){ 721 //nothing 722 } 723 else if (control_type[i]==DhDtEnum){ 724 //nothing 725 } 726 else if (control_type[i]==VxEnum){ 727 //nothing 728 } 729 else if (control_type[i]==VyEnum){ 730 //nothing 731 } 732 else{ 733 ISSMERROR("unsupported control type: %s",EnumToString(control_type[i])); 734 } 735 } 736 737 delete gauss; 738 } 739 740 /*Clean up and return*/ 741 xfree((void**)&control_type); 742 return Jelem; 743 } 744 /*}}}*/ 745 /*FUNCTION Tria::CreateKMatrix {{{1*/ 746 void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ 747 748 /*retreive parameters: */ 749 ElementMatrix* Ke=NULL; 750 int analysis_type; 751 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 752 753 /*Checks in debugging mode{{{2*/ 754 ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); 755 /*}}}*/ 756 757 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ 758 switch(analysis_type){ 759 case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: 760 Ke=CreateKMatrixDiagnosticMacAyeal(); 761 break; 762 case DiagnosticHutterAnalysisEnum: 763 Ke=CreateKMatrixDiagnosticHutter(); 764 break; 765 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: 766 Ke=CreateKMatrixSlope(); 767 break; 768 case PrognosticAnalysisEnum: 769 Ke=CreateKMatrixPrognostic(); 770 break; 771 case BalancedthicknessAnalysisEnum: 772 Ke=CreateKMatrixBalancedthickness(); 773 break; 774 case AdjointBalancedthicknessAnalysisEnum: 775 Ke=CreateKMatrixAdjointBalancedthickness(); 776 break; 777 case BalancedvelocitiesAnalysisEnum: 778 Ke=CreateKMatrixBalancedvelocities(); 779 break; 780 default: 781 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); 782 } 783 784 /*Add to global matrix*/ 785 if(Ke){ 786 Ke->AddToGlobal(Kgg,Kff,Kfs); 787 delete Ke; 788 } 789 } 790 /*}}}*/ 791 /*FUNCTION Tria::CreatePVector {{{1*/ 792 void Tria::CreatePVector(Vec pg, Vec pf){ 793 794 /*retrive parameters: */ 795 ElementVector* pe=NULL; 796 int analysis_type; 797 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 798 799 /*asserts: {{{*/ 800 /*if debugging mode, check that all pointers exist*/ 801 ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); 802 /*}}}*/ 803 804 /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ 805 switch(analysis_type){ 806 case DiagnosticHorizAnalysisEnum: 807 pe=CreatePVectorDiagnosticMacAyeal(); 808 break; 809 case AdjointHorizAnalysisEnum: 810 pe=CreatePVectorAdjointHoriz(); 811 break; 812 case DiagnosticHutterAnalysisEnum: 813 pe=CreatePVectorDiagnosticHutter(); 814 break; 815 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: 816 pe=CreatePVectorSlope(); 817 break; 818 case PrognosticAnalysisEnum: 819 pe=CreatePVectorPrognostic(); 820 break; 821 case BalancedthicknessAnalysisEnum: 822 pe=CreatePVectorBalancedthickness(); 823 break; 824 case AdjointBalancedthicknessAnalysisEnum: 825 pe=CreatePVectorAdjointBalancedthickness(); 826 break; 827 case BalancedvelocitiesAnalysisEnum: 828 pe=CreatePVectorBalancedvelocities(); 829 break; 830 default: 831 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); 832 } 833 834 /*Add to global Vector*/ 835 if(pe){ 836 pe->AddToGlobal(pg,pf); 837 delete pe; 838 } 839 } 840 /*}}}*/ 841 /*FUNCTION Tria::DeleteResults {{{1*/ 842 void Tria::DeleteResults(void){ 843 844 /*Delete and reinitialize results*/ 845 delete this->results; 846 this->results=new Results(); 847 848 } 849 /*}}}*/ 850 /*FUNCTION Tria::GetNodeIndex {{{1*/ 851 int Tria::GetNodeIndex(Node* node){ 852 853 ISSMASSERT(nodes); 854 for(int i=0;i<NUMVERTICES;i++){ 855 if(node==nodes[i]) 856 return i; 857 } 858 ISSMERROR("Node provided not found among element nodes"); 859 } 860 /*}}}*/ 861 /*FUNCTION Tria::IsOnBed {{{1*/ 862 bool Tria::IsOnBed(){ 863 864 bool onbed; 865 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 866 return onbed; 867 } 868 /*}}}*/ 869 /*FUNCTION Tria::IsOnShelf {{{1*/ 870 bool Tria::IsOnShelf(){ 871 872 bool shelf; 873 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 874 return shelf; 875 } 876 /*}}}*/ 877 /*FUNCTION Tria::IsOnWater {{{1*/ 878 bool Tria::IsOnWater(){ 879 880 bool water; 881 inputs->GetParameterValue(&water,ElementOnWaterEnum); 882 return water; 883 } 884 /*}}}*/ 885 /*FUNCTION Tria::GetSolutionFromInputs{{{1*/ 886 void Tria::GetSolutionFromInputs(Vec solution){ 887 888 /*retrive parameters: */ 889 int analysis_type; 890 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 891 892 /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ 893 if (analysis_type==DiagnosticHorizAnalysisEnum) 894 GetSolutionFromInputsDiagnosticHoriz(solution); 895 else if (analysis_type==DiagnosticHutterAnalysisEnum) 896 GetSolutionFromInputsDiagnosticHutter(solution); 897 else 898 ISSMERROR("analysis: %s not supported yet",EnumToString(analysis_type)); 899 900 } 901 /*}}}*/ 902 /*FUNCTION Tria::GetVectorFromInputs{{{1*/ 903 void Tria::GetVectorFromInputs(Vec vector,int NameEnum){ 904 905 int doflist1[NUMVERTICES]; 906 907 /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ 908 for(int i=0;i<this->inputs->Size();i++){ 909 Input* input=(Input*)this->inputs->GetObjectByOffset(i); 910 if(input->EnumType()==NameEnum){ 911 /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ 912 this->GetDofList1(&doflist1[0]); 913 input->GetVectorFromInputs(vector,&doflist1[0]); 914 break; 915 } 916 } 917 } 918 /*}}}*/ 919 /*FUNCTION Tria::Gradj {{{1*/ 920 void Tria::Gradj(Vec gradient,int control_type){ 921 922 /*If on water, grad = 0: */ 923 if(IsOnWater())return; 924 925 switch(control_type){ 926 case DragCoefficientEnum: 927 GradjDrag(gradient); 928 break; 929 case RheologyBbarEnum: 930 GradjB(gradient); 931 break; 932 case DhDtEnum: 933 GradjDhDt(gradient); 934 break; 935 case VxEnum: 936 GradjVx(gradient); 937 break; 938 case VyEnum: 939 GradjVy(gradient); 940 break; 941 default: 942 ISSMERROR("%s%i","control type not supported yet: ",control_type); 943 } 944 } 945 /*}}}*/ 946 /*FUNCTION Tria::GradjB{{{1*/ 947 void Tria::GradjB(Vec gradient){ 948 949 /*Intermediaries*/ 950 int i,ig; 951 int doflist[NUMVERTICES]; 952 double vx,vy,lambda,mu,thickness,Jdet; 953 double cm_noisedmp,viscosity_complement; 954 double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; 955 double xyz_list[NUMVERTICES][3]; 956 double basis[3],epsilon[3]; 957 double dbasis[NDOF2][NUMVERTICES]; 958 double grad_g[NUMVERTICES]; 959 double grad[NUMVERTICES]={0.0}; 960 GaussTria *gauss = NULL; 961 962 /*retrieve some parameters: */ 963 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); 964 965 /* Get node coordinates and dof list: */ 966 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 967 GetDofList1(&doflist[0]); 968 969 /*Retrieve all inputs*/ 970 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 971 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 972 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 973 Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); 974 Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); 975 Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(rheologyb_input); 976 977 /* Start looping on the number of gaussian points: */ 978 gauss=new GaussTria(4); 979 for (ig=gauss->begin();ig<gauss->end();ig++){ 980 981 gauss->GaussPoint(ig); 982 983 thickness_input->GetParameterValue(&thickness,gauss); 984 rheologyb_input->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],gauss); 985 vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); 986 vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); 987 adjointx_input->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); 988 adjointy_input->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); 989 990 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 991 matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); 992 993 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 994 GetNodalFunctions(basis,gauss); 995 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); 996 997 /*standard gradient dJ/dki*/ 998 for (i=0;i<NUMVERTICES;i++){ 999 grad_g[i]=-viscosity_complement*thickness*( (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1])*Jdet*gauss->weight*basis[i]; 1000 } 1001 /*Add regularization term*/ 1002 for (i=0;i<NUMVERTICES;i++) grad_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); 1003 for(i=0;i<NUMVERTICES;i++) grad[i]+=grad_g[i]; 1004 } 1005 1006 VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES); 1007 1008 /*clean-up*/ 1009 delete gauss; 1010 } 1011 /*}}}*/ 1012 /*FUNCTION Tria::GradjDrag {{{1*/ 1013 void Tria::GradjDrag(Vec gradient){ 1014 1015 int i,ig; 1016 int drag_type,analysis_type; 1017 int doflist1[NUMVERTICES]; 1018 double vx,vy,lambda,mu,alpha_complement,Jdet; 1019 double bed,thickness,Neff,drag,cm_noisedmp; 1020 double xyz_list[NUMVERTICES][3]; 1021 double dh1dh3[NDOF2][NUMVERTICES]; 1022 double dk[NDOF2]; 1023 double grade_g[NUMVERTICES]={0.0}; 1024 double grade_g_gaussian[NUMVERTICES]; 1025 double l1l2l3[3]; 1026 double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ 1027 Friction* friction=NULL; 1028 GaussTria *gauss=NULL; 1029 1030 /*retrive parameters: */ 1031 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 1032 1033 /*retrieve some parameters ands return if iceshelf: */ 1034 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); 1035 if(IsOnShelf())return; 1036 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1037 GetDofList1(&doflist1[0]); 1038 1039 /*Build frictoin element, needed later: */ 1040 inputs->GetParameterValue(&drag_type,DragTypeEnum); 1041 friction=new Friction("2d",inputs,matpar,analysis_type); 1042 1043 /*Retrieve all inputs we will be needing: */ 1044 Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); 1045 Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); 1046 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 1047 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 1048 Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(dragcoefficient_input); 1049 1050 /* Start looping on the number of gaussian points: */ 1051 gauss=new GaussTria(4); 1052 for (ig=gauss->begin();ig<gauss->end();ig++){ 1053 1054 gauss->GaussPoint(ig); 1055 1056 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1057 GetNodalFunctions(l1l2l3, gauss); 1058 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); 1059 1060 /*Build alpha_complement_list: */ 1061 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); 1062 else alpha_complement=0; 1063 1064 dragcoefficient_input->GetParameterValue(&drag, gauss); 1065 adjointx_input->GetParameterValue(&lambda, gauss); 1066 adjointy_input->GetParameterValue(&mu, gauss); 1067 vx_input->GetParameterValue(&vx,gauss); 1068 vy_input->GetParameterValue(&vy,gauss); 1069 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); 1070 1071 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ 1072 for (i=0;i<NUMVERTICES;i++){ 1073 1074 //standard term dJ/dki 1075 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i]; 1076 1077 //noise dampening d/dki(1/2*(dk/dx)^2) 1078 grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); 1079 } 1080 1081 /*Add gradje_g_gaussian vector to gradje_g: */ 1082 for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i]; 1083 } 1084 1085 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 1086 1087 /*Clean up and return*/ 1088 delete gauss; 1089 delete friction; 1090 } 1091 /*}}}*/ 1092 /*FUNCTION Tria::GradjDhDt{{{1*/ 1093 void Tria::GradjDhDt(Vec gradient){ 1094 1095 /*Intermediaries*/ 1096 int doflist1[NUMVERTICES]; 1097 double lambda[NUMVERTICES]; 1098 double gradient_g[NUMVERTICES]; 1099 1100 GetDofList1(&doflist1[0]); 1101 1102 /*Compute Gradient*/ 1103 GetParameterListOnVertices(&lambda[0],AdjointEnum); 1104 for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i]; 1105 1106 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)gradient_g,INSERT_VALUES); 1107 } 1108 /*}}}*/ 1109 /*FUNCTION Tria::GradjVx{{{1*/ 1110 void Tria::GradjVx(Vec gradient){ 1111 1112 /*Intermediaries*/ 1113 int i,ig; 1114 int doflist1[NUMVERTICES]; 1115 double thickness,Jdet; 1116 double l1l2l3[3]; 1117 double Dlambda[2]; 1118 double xyz_list[NUMVERTICES][3]; 1119 double grade_g[NUMVERTICES] = {0.0}; 1120 GaussTria *gauss = NULL; 1121 1122 /* Get node coordinates and dof list: */ 1123 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1124 GetDofList1(&doflist1[0]); 1125 1126 /*Retrieve all inputs we will be needing: */ 1127 Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); 1128 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1129 1130 /* Start looping on the number of gaussian points: */ 1131 gauss=new GaussTria(2); 1132 for (ig=gauss->begin();ig<gauss->end();ig++){ 1133 1134 gauss->GaussPoint(ig); 1135 1136 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1137 GetNodalFunctions(l1l2l3, gauss); 1138 1139 adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); 1140 thickness_input->GetParameterValue(&thickness, gauss); 1141 1142 for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i]; 1143 } 1144 1145 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 1146 1147 /*Clean up and return*/ 1148 delete gauss; 1149 } 1150 /*}}}*/ 1151 /*FUNCTION Tria::GradjVy{{{1*/ 1152 void Tria::GradjVy(Vec gradient){ 1153 1154 /*Intermediaries*/ 1155 int i,ig; 1156 int doflist1[NUMVERTICES]; 1157 double thickness,Jdet; 1158 double l1l2l3[3]; 1159 double Dlambda[2]; 1160 double xyz_list[NUMVERTICES][3]; 1161 double grade_g[NUMVERTICES] = {0.0}; 1162 GaussTria *gauss = NULL; 1163 1164 /* Get node coordinates and dof list: */ 1165 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1166 GetDofList1(&doflist1[0]); 1167 1168 /*Retrieve all inputs we will be needing: */ 1169 Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); 1170 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1171 1172 /* Start looping on the number of gaussian points: */ 1173 gauss=new GaussTria(2); 1174 for (ig=gauss->begin();ig<gauss->end();ig++){ 1175 1176 gauss->GaussPoint(ig); 1177 1178 adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); 1179 thickness_input->GetParameterValue(&thickness, gauss); 1180 1181 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1182 GetNodalFunctions(l1l2l3, gauss); 1183 1184 for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[1]*Jdet*gauss->weight*l1l2l3[i]; 1185 } 1186 1187 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 1188 1189 /*Clean up and return*/ 1190 delete gauss; 1191 } 1192 /*}}}*/ 1193 /*FUNCTION Tria::InputControlUpdate{{{1*/ 1194 void Tria::InputControlUpdate(double scalar,bool save_parameter){ 1195 1196 /*Intermediary*/ 1197 int num_controls; 1198 int* control_type=NULL; 1199 Input* input=NULL; 1200 double cm_min,cm_max; 1201 1202 /*retrieve some parameters: */ 1203 this->parameters->FindParam(&cm_min,CmMinEnum); 1204 this->parameters->FindParam(&cm_max,CmMaxEnum); 1205 this->parameters->FindParam(&num_controls,NumControlsEnum); 1206 this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); 1207 1208 for(int i=0;i<num_controls;i++){ 1209 1210 if(control_type[i]==RheologyBbarEnum){ 1211 input=(Input*)matice->inputs->GetInput(control_type[i]); ISSMASSERT(input); 1212 } 1213 else{ 1214 input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input); 1215 } 1216 1217 if (input->Enum()!=ControlInputEnum){ 1218 ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i])); 1219 } 1220 1221 ((ControlInput*)input)->UpdateValue(scalar); 1222 input->Constrain(cm_min,cm_max); 1223 if (save_parameter) ((ControlInput*)input)->SaveValue(); 1224 1225 } 1226 1227 /*Clean up and return*/ 1228 xfree((void**)&control_type); 1229 } 1230 /*}}}*/ 1231 /*FUNCTION Tria::InputConvergence{{{1*/ 1232 bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ 1233 1234 bool converged=true; 1235 int i; 1236 Input** new_inputs=NULL; 1237 Input** old_inputs=NULL; 1238 1239 new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs 1240 old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs 1241 1242 for(i=0;i<num_enums/2;i++){ 1243 new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]); 1244 old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]); 1245 if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); 1246 if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); 1247 } 1248 1249 /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ 1250 for(i=0;i<num_criterionenums;i++){ 1251 IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]); 1252 if(eps[i]>criterionvalues[i]) converged=false; 1253 } 1254 1255 /*clean up and return*/ 1256 xfree((void**)&new_inputs); 1257 xfree((void**)&old_inputs); 1258 return converged; 1259 } 1260 /*}}}*/ 1261 /*FUNCTION Tria::InputDepthAverageAtBase {{{1*/ 1262 void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ 1263 1264 /*New input*/ 1265 Input* oldinput=NULL; 1266 Input* newinput=NULL; 1267 1268 /*copy input of enum_type*/ 1269 if (object_enum==ElementsEnum) 1270 oldinput=(Input*)this->inputs->GetInput(enum_type); 1271 else if (object_enum==MaterialsEnum) 1272 oldinput=(Input*)this->matice->inputs->GetInput(enum_type); 1273 else 1274 ISSMERROR("object %s not supported yet",EnumToString(object_enum)); 1275 if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumToString(enum_type)); 1276 newinput=(Input*)oldinput->copy(); 1277 1278 /*Assign new name (average)*/ 1279 newinput->ChangeEnum(average_enum_type); 1280 1281 /*Add new input to current element*/ 1282 if (object_enum==ElementsEnum) 1283 this->inputs->AddInput((Input*)newinput); 1284 else if (object_enum==MaterialsEnum) 1285 this->matice->inputs->AddInput((Input*)newinput); 1286 else 1287 ISSMERROR("object %s not supported yet",EnumToString(object_enum)); 1288 } 1289 /*}}}*/ 1290 /*FUNCTION Tria::InputDuplicate{{{1*/ 1291 void Tria::InputDuplicate(int original_enum,int new_enum){ 1292 1293 /*Call inputs method*/ 1294 if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); 1295 1296 } 1297 /*}}}*/ 1298 /*FUNCTION Tria::InputScale{{{1*/ 1299 void Tria::InputScale(int enum_type,double scale_factor){ 1300 1301 Input* input=NULL; 1302 1303 /*Make a copy of the original input: */ 1304 input=(Input*)this->inputs->GetInput(enum_type); 1305 if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); 1306 1307 /*Scale: */ 1308 input->Scale(scale_factor); 1309 } 1310 /*}}}*/ 1311 /*FUNCTION Tria::InputArtificialNoise{{{1*/ 1312 void Tria::InputArtificialNoise(int enum_type,double min,double max){ 1313 1314 Input* input=NULL; 1315 1316 /*Make a copy of the original input: */ 1317 input=(Input*)this->inputs->GetInput(enum_type); 1318 if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); 1319 1320 /*ArtificialNoise: */ 1321 input->ArtificialNoise(min,max); 1322 } 1323 /*}}}*/ 1324 /*FUNCTION Tria::InputToResult{{{1*/ 1325 void Tria::InputToResult(int enum_type,int step,double time){ 1326 1327 int i; 1328 Input *input = NULL; 1329 1330 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ 1331 if (enum_type==RheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); 1332 else input=this->inputs->GetInput(enum_type); 1333 if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumToString(enum_type)); 1334 1335 /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 1336 * object out of the input, with the additional step and time information: */ 1337 this->results->AddObject((Object*)input->SpawnResult(step,time)); 1338 if(input->Enum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); 1339 } 1340 /*}}}*/ 1341 /*FUNCTION Tria::MassFlux {{{1*/ 1342 double Tria::MassFlux( double* segment,bool process_units){ 1343 1344 const int numdofs=2; 1345 1346 int i; 1347 double mass_flux=0; 1348 double xyz_list[NUMVERTICES][3]; 1349 double normal[2]; 1350 double length,rho_ice; 1351 double x1,y1,x2,y2,h1,h2; 1352 double vx1,vx2,vy1,vy2; 1353 GaussTria* gauss_1=NULL; 1354 GaussTria* gauss_2=NULL; 1355 1356 /*Get material parameters :*/ 1357 rho_ice=matpar->GetRhoIce(); 1358 1359 /*First off, check that this segment belongs to this element: */ 1360 if ((int)*(segment+4)!=this->id)ISSMERROR("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id); 1361 1362 /*Recover segment node locations: */ 1363 x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); 1364 1365 /*Get xyz list: */ 1366 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1367 1368 /*get area coordinates of 0 and 1 locations: */ 1369 gauss_1=new GaussTria(); 1370 gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); 1371 gauss_2=new GaussTria(); 1372 gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); 1373 1374 normal[0]=cos(atan2(x1-x2,y2-y1)); 1375 normal[1]=sin(atan2(x1-x2,y2-y1)); 1376 1377 length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); 1378 1379 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1380 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 1381 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 1382 1383 thickness_input->GetParameterValue(&h1, gauss_1); 1384 thickness_input->GetParameterValue(&h2, gauss_2); 1385 vx_input->GetParameterValue(&vx1,gauss_1); 1386 vx_input->GetParameterValue(&vx2,gauss_2); 1387 vy_input->GetParameterValue(&vy1,gauss_1); 1388 vy_input->GetParameterValue(&vy2,gauss_2); 1389 1390 mass_flux= rho_ice*length*( 1391 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ 1392 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] 1393 ); 1394 1395 /*Process units: */ 1396 mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum,this->parameters); 1397 1398 /*clean up and return:*/ 1399 delete gauss_1; 1400 delete gauss_2; 1401 return mass_flux; 1402 } 1403 /*}}}*/ 1404 /*FUNCTION Tria::MaxAbsVx{{{1*/ 1405 void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ 1406 1407 /*Get maximum:*/ 1408 double maxabsvx=this->inputs->MaxAbs(VxEnum); 1409 1410 /*process units if requested: */ 1411 if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum,this->parameters); 1412 1413 /*Assign output pointers:*/ 1414 *pmaxabsvx=maxabsvx; 1415 } 1416 /*}}}*/ 1417 /*FUNCTION Tria::MaxAbsVy{{{1*/ 1418 void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ 1419 1420 /*Get maximum:*/ 1421 double maxabsvy=this->inputs->MaxAbs(VyEnum); 1422 1423 /*process units if requested: */ 1424 if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum,this->parameters); 1425 1426 /*Assign output pointers:*/ 1427 *pmaxabsvy=maxabsvy; 1428 } 1429 /*}}}*/ 1430 /*FUNCTION Tria::MaxAbsVz{{{1*/ 1431 void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ 1432 1433 /*Get maximum:*/ 1434 double maxabsvz=this->inputs->MaxAbs(VzEnum); 1435 1436 /*process units if requested: */ 1437 if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum,this->parameters); 1438 1439 /*Assign output pointers:*/ 1440 *pmaxabsvz=maxabsvz; 1441 } 1442 /*}}}*/ 1443 /*FUNCTION Tria::MaxVel{{{1*/ 1444 void Tria::MaxVel(double* pmaxvel, bool process_units){ 1445 1446 /*Get maximum:*/ 1447 double maxvel=this->inputs->Max(VelEnum); 1448 1449 /*process units if requested: */ 1450 if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum,this->parameters); 1451 1452 /*Assign output pointers:*/ 1453 *pmaxvel=maxvel; 1454 } 1455 /*}}}*/ 1456 /*FUNCTION Tria::MaxVx{{{1*/ 1457 void Tria::MaxVx(double* pmaxvx, bool process_units){ 1458 1459 /*Get maximum:*/ 1460 double maxvx=this->inputs->Max(VxEnum); 1461 1462 /*process units if requested: */ 1463 if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum,this->parameters); 1464 1465 /*Assign output pointers:*/ 1466 *pmaxvx=maxvx; 1467 } 1468 /*}}}*/ 1469 /*FUNCTION Tria::MaxVy{{{1*/ 1470 void Tria::MaxVy(double* pmaxvy, bool process_units){ 1471 1472 /*Get maximum:*/ 1473 double maxvy=this->inputs->Max(VyEnum); 1474 1475 /*process units if requested: */ 1476 if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum,this->parameters); 1477 1478 /*Assign output pointers:*/ 1479 *pmaxvy=maxvy; 1480 1481 } 1482 /*}}}*/ 1483 /*FUNCTION Tria::MaxVz{{{1*/ 1484 void Tria::MaxVz(double* pmaxvz, bool process_units){ 1485 1486 /*Get maximum:*/ 1487 double maxvz=this->inputs->Max(VzEnum); 1488 1489 /*process units if requested: */ 1490 if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum,this->parameters); 1491 1492 /*Assign output pointers:*/ 1493 *pmaxvz=maxvz; 1494 } 1495 /*}}}*/ 1496 /*FUNCTION Tria::MinVel{{{1*/ 1497 void Tria::MinVel(double* pminvel, bool process_units){ 1498 1499 /*Get minimum:*/ 1500 double minvel=this->inputs->Min(VelEnum); 1501 1502 /*process units if requested: */ 1503 if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum,this->parameters); 1504 1505 /*Assign output pointers:*/ 1506 *pminvel=minvel; 1507 } 1508 /*}}}*/ 1509 /*FUNCTION Tria::MinVx{{{1*/ 1510 void Tria::MinVx(double* pminvx, bool process_units){ 1511 1512 /*Get minimum:*/ 1513 double minvx=this->inputs->Min(VxEnum); 1514 1515 /*process units if requested: */ 1516 if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum,this->parameters); 1517 1518 /*Assign output pointers:*/ 1519 *pminvx=minvx; 1520 } 1521 /*}}}*/ 1522 /*FUNCTION Tria::MinVy{{{1*/ 1523 void Tria::MinVy(double* pminvy, bool process_units){ 1524 1525 /*Get minimum:*/ 1526 double minvy=this->inputs->Min(VyEnum); 1527 1528 /*process units if requested: */ 1529 if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum,this->parameters); 1530 1531 /*Assign output pointers:*/ 1532 *pminvy=minvy; 1533 } 1534 /*}}}*/ 1535 /*FUNCTION Tria::MinVz{{{1*/ 1536 void Tria::MinVz(double* pminvz, bool process_units){ 1537 1538 /*Get minimum:*/ 1539 double minvz=this->inputs->Min(VzEnum); 1540 1541 /*process units if requested: */ 1542 if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum,this->parameters); 1543 1544 /*Assign output pointers:*/ 1545 *pminvz=minvz; 1546 } 1547 /*}}}*/ 1548 /*FUNCTION Tria::TimeAdapt{{{1*/ 1549 double Tria::TimeAdapt(void){ 1550 1551 /*intermediary: */ 1552 int i; 1553 double C,dt; 1554 double dx,dy; 1555 double maxx,minx; 1556 double maxy,miny; 1557 double maxabsvx,maxabsvy; 1558 double xyz_list[NUMVERTICES][3]; 1559 1560 /*get CFL coefficient:*/ 1561 this->parameters->FindParam(&C,CflCoefficientEnum); 1562 1563 /*Get for Vx and Vy, the max of abs value: */ 1564 this->MaxAbsVx(&maxabsvx,false); 1565 this->MaxAbsVy(&maxabsvy,false); 1566 1567 /* Get node coordinates and dof list: */ 1568 GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES); 1569 1570 minx=xyz_list[0][0]; 1571 maxx=xyz_list[0][0]; 1572 miny=xyz_list[0][1]; 1573 maxy=xyz_list[0][1]; 1574 1575 for(i=1;i<NUMVERTICES;i++){ 1576 if (xyz_list[i][0]<minx)minx=xyz_list[i][0]; 1577 if (xyz_list[i][0]>maxx)maxx=xyz_list[i][0]; 1578 if (xyz_list[i][1]<miny)miny=xyz_list[i][1]; 1579 if (xyz_list[i][1]>maxy)maxy=xyz_list[i][1]; 1580 } 1581 dx=maxx-minx; 1582 dy=maxy-miny; 1583 1584 /*CFL criterion: */ 1585 dt=C/(maxabsvy/dx+maxabsvy/dy); 1586 1587 return dt; 1588 } 1589 /*}}}*/ 1590 /*FUNCTION Tria::ThicknessAbsMisfit {{{1*/ 1591 double Tria::ThicknessAbsMisfit(bool process_units){ 1592 1593 /* Constants */ 1594 const int numdof=1*NUMVERTICES; 1595 1596 /*Intermediaries*/ 1597 int i,ig; 1598 double thickness,thicknessobs,weight; 1599 double Jdet; 1600 double Jelem = 0; 1601 double xyz_list[NUMVERTICES][3]; 1602 GaussTria *gauss = NULL; 1603 1604 /*If on water, return 0: */ 1605 if(IsOnWater())return 0; 1606 1607 /* Get node coordinates and dof list: */ 1608 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1609 1610 /*Retrieve all inputs we will be needing: */ 1611 Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1612 Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); 1613 Input* weights_input =inputs->GetInput(WeightsEnum); ISSMASSERT(weights_input); 1614 1615 /* Start looping on the number of gaussian points: */ 1616 gauss=new GaussTria(2); 1617 for (ig=gauss->begin();ig<gauss->end();ig++){ 1618 1619 gauss->GaussPoint(ig); 1620 1621 /* Get Jacobian determinant: */ 1622 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1623 1624 /*Get parameters at gauss point*/ 1625 thickness_input->GetParameterValue(&thickness,gauss); 1626 thicknessobs_input->GetParameterValue(&thicknessobs,gauss); 1627 weights_input->GetParameterValue(&weight,gauss); 1628 1629 /*compute ThicknessAbsMisfit*/ 1630 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight; 1631 } 1632 1633 /* clean up and Return: */ 1634 delete gauss; 1635 return Jelem; 1636 } 1637 /*}}}*/ 1638 /*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ 1639 double Tria::SurfaceAbsVelMisfit(bool process_units){ 1640 1641 const int numdof=NDOF2*NUMVERTICES; 1642 1643 int i,ig; 1644 double Jelem=0; 1645 double velocity_mag,obs_velocity_mag; 1646 double meanvel, epsvel,misfit,Jdet; 1647 double xyz_list[NUMVERTICES][3]; 1648 double vx_list[NUMVERTICES]; 1649 double vy_list[NUMVERTICES]; 1650 double obs_vx_list[NUMVERTICES]; 1651 double obs_vy_list[NUMVERTICES]; 1652 double misfit_square_list[NUMVERTICES]; 1653 double misfit_list[NUMVERTICES]; 1654 double weights_list[NUMVERTICES]; 1655 GaussTria *gauss=NULL; 1656 1657 /*If on water, return 0: */ 1658 if(IsOnWater())return 0; 1659 1660 /* Get node coordinates and dof list: */ 1661 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1662 1663 /* Recover input data: */ 1664 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1665 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1666 GetParameterListOnVertices(&vx_list[0],VxEnum); 1667 GetParameterListOnVertices(&vy_list[0],VyEnum); 1668 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1669 1670 /*retrieve some parameters: */ 1671 this->parameters->FindParam(&meanvel,MeanVelEnum); 1672 this->parameters->FindParam(&epsvel,EpsVelEnum); 1673 1674 /* Compute SurfaceAbsVelMisfit at the 3 nodes 1675 * Here we integrate linearized functions: 1676 * 1677 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1678 * 1679 * where J_i are the misfits at the 3 nodes of the triangle 1680 * Phi_i is the nodal function (P1) with respect to 1681 * the vertex i 1682 */ 1683 1684 /*We are using an absolute misfit: 1685 * 1686 * 1 [ 2 2 ] 1687 * J = --- | (u - u ) + (v - v ) | 1688 * 2 [ obs obs ] 1689 * 1690 */ 1691 for (i=0;i<NUMVERTICES;i++){ 1692 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)); 1693 } 1694 /*Process units: */ 1695 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters); 1696 1697 /*Apply weights to misfits*/ 1698 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 1699 1700 /* Start looping on the number of gaussian points: */ 1701 gauss=new GaussTria(2); 1702 for (ig=gauss->begin();ig<gauss->end(); ig++){ 1703 1704 gauss->GaussPoint(ig); 1705 1706 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1707 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 1708 Jelem+=misfit*Jdet*gauss->weight; 1709 } 1710 1711 /*clean up and Return: */ 1712 delete gauss; 1713 return Jelem; 1714 } 1715 /*}}}*/ 1716 /*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ 1717 double Tria::SurfaceRelVelMisfit(bool process_units){ 1718 1719 const int numdof=2*NUMVERTICES; 1720 1721 int i,ig; 1722 double Jelem=0; 1723 double scalex=1,scaley=1; 1724 double meanvel, epsvel,misfit,Jdet; 1725 double velocity_mag,obs_velocity_mag; 1726 double xyz_list[NUMVERTICES][3]; 1727 double vx_list[NUMVERTICES]; 1728 double vy_list[NUMVERTICES]; 1729 double obs_vx_list[NUMVERTICES]; 1730 double obs_vy_list[NUMVERTICES]; 1731 double misfit_square_list[NUMVERTICES]; 1732 double misfit_list[NUMVERTICES]; 1733 double weights_list[NUMVERTICES]; 1734 GaussTria *gauss=NULL; 1735 1736 /*If on water, return 0: */ 1737 if(IsOnWater())return 0; 1738 1739 /* Get node coordinates and dof list: */ 1740 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1741 1742 /* Recover input data: */ 1743 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1744 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1745 GetParameterListOnVertices(&vx_list[0],VxEnum); 1746 GetParameterListOnVertices(&vy_list[0],VyEnum); 1747 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1748 1749 /*retrieve some parameters: */ 1750 this->parameters->FindParam(&meanvel,MeanVelEnum); 1751 this->parameters->FindParam(&epsvel,EpsVelEnum); 1752 1753 /* Compute SurfaceRelVelMisfit at the 3 nodes 1754 * Here we integrate linearized functions: 1755 * 1756 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1757 * 1758 * where J_i are the misfits at the 3 nodes of the triangle 1759 * Phi_i is the nodal function (P1) with respect to 1760 * the vertex i 1761 */ 1762 1763 /*We are using a relative misfit: 1764 * 1765 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1766 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1767 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1768 * obs obs 1769 */ 1770 for (i=0;i<NUMVERTICES;i++){ 1771 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2); 1772 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2); 1773 if(obs_vx_list[i]==0)scalex=0; 1774 if(obs_vy_list[i]==0)scaley=0; 1775 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)); 1776 } 1777 1778 /*Process units: */ 1779 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters); 1780 1781 /*Apply weights to misfits*/ 1782 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 1783 1784 /* Start looping on the number of gaussian points: */ 1785 gauss=new GaussTria(2); 1786 for (ig=gauss->begin();ig<gauss->end();ig++){ 1787 1788 gauss->GaussPoint(ig); 1789 1790 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1791 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 1792 Jelem+=misfit*Jdet*gauss->weight; 1793 } 1794 1795 /*clean up and Return: */ 1796 delete gauss; 1797 return Jelem; 1798 } 1799 /*}}}*/ 1800 /*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/ 1801 double Tria::SurfaceLogVelMisfit(bool process_units){ 1802 1803 const int numdof=NDOF2*NUMVERTICES; 1804 1805 int i,ig; 1806 double Jelem=0; 1807 double scalex=1,scaley=1; 1808 double meanvel, epsvel,misfit,Jdet; 1809 double velocity_mag,obs_velocity_mag; 1810 double xyz_list[NUMVERTICES][3]; 1811 double vx_list[NUMVERTICES]; 1812 double vy_list[NUMVERTICES]; 1813 double obs_vx_list[NUMVERTICES]; 1814 double obs_vy_list[NUMVERTICES]; 1815 double misfit_square_list[NUMVERTICES]; 1816 double misfit_list[NUMVERTICES]; 1817 double weights_list[NUMVERTICES]; 1818 GaussTria *gauss=NULL; 1819 1820 /*If on water, return 0: */ 1821 if(IsOnWater())return 0; 1822 1823 /* Get node coordinates and dof list: */ 1824 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1825 1826 /* Recover input data: */ 1827 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1828 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1829 GetParameterListOnVertices(&vx_list[0],VxEnum); 1830 GetParameterListOnVertices(&vy_list[0],VyEnum); 1831 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1832 1833 /*retrieve some parameters: */ 1834 this->parameters->FindParam(&meanvel,MeanVelEnum); 1835 this->parameters->FindParam(&epsvel,EpsVelEnum); 1836 1837 /* Compute SurfaceLogVelMisfit at the 3 nodes 1838 * Here we integrate linearized functions: 1839 * 1840 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1841 * 1842 * where J_i are the misfits at the 3 nodes of the triangle 1843 * Phi_i is the nodal function (P1) with respect to 1844 * the vertex i 1845 */ 1846 1847 /*We are using a logarithmic misfit: 1848 * 1849 * [ vel + eps ] 2 1850 * J = 4 \bar{v}^2 | log ( ----------- ) | 1851 * [ vel + eps ] 1852 * obs 1853 */ 1854 for (i=0;i<NUMVERTICES;i++){ 1855 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil. 1856 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. 1857 misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); 1858 } 1859 1860 /*Process units: */ 1861 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters); 1862 1863 /*Apply weights to misfits*/ 1864 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 1865 1866 /* Start looping on the number of gaussian points: */ 1867 gauss=new GaussTria(2); 1868 for (ig=gauss->begin();ig<gauss->end();ig++){ 1869 1870 gauss->GaussPoint(ig); 1871 1872 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1873 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 1874 Jelem+=misfit*Jdet*gauss->weight; 1875 } 1876 1877 /*clean-up and Return: */ 1878 delete gauss; 1879 return Jelem; 1880 } 1881 /*}}}*/ 1882 /*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/ 1883 double Tria::SurfaceLogVxVyMisfit(bool process_units){ 1884 1885 const int numdof=NDOF2*NUMVERTICES; 1886 1887 int i,ig; 1888 int fit=-1; 1889 double Jelem=0, S=0; 1890 double scalex=1,scaley=1; 1891 double meanvel, epsvel, misfit, Jdet; 1892 double velocity_mag,obs_velocity_mag; 1893 double xyz_list[NUMVERTICES][3]; 1894 double vx_list[NUMVERTICES]; 1895 double vy_list[NUMVERTICES]; 1896 double obs_vx_list[NUMVERTICES]; 1897 double obs_vy_list[NUMVERTICES]; 1898 double misfit_square_list[NUMVERTICES]; 1899 double misfit_list[NUMVERTICES]; 1900 double weights_list[NUMVERTICES]; 1901 GaussTria *gauss=NULL; 1902 1903 /*If on water, return 0: */ 1904 if(IsOnWater())return 0; 1905 1906 /* Get node coordinates and dof list: */ 1907 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1908 1909 /* Recover input data: */ 1910 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1911 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1912 GetParameterListOnVertices(&vx_list[0],VxEnum); 1913 GetParameterListOnVertices(&vy_list[0],VyEnum); 1914 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1915 1916 /*retrieve some parameters: */ 1917 this->parameters->FindParam(&meanvel,MeanVelEnum); 1918 this->parameters->FindParam(&epsvel,EpsVelEnum); 1919 1920 /* Compute SurfaceLogVxVyMisfit at the 3 nodes 1921 * Here we integrate linearized functions: 1922 * 1923 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1924 * 1925 * where J_i are the misfits at the 3 nodes of the triangle 1926 * Phi_i is the nodal function (P1) with respect to 1927 * the vertex i 1928 */ 1929 1930 /*We are using an logarithmic 2 misfit: 1931 * 1932 * 1 [ |u| + eps 2 |v| + eps 2 ] 1933 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1934 * 2 [ |u |+ eps |v |+ eps ] 1935 * obs obs 1936 */ 1937 for (i=0;i<NUMVERTICES;i++){ 1938 misfit_list[i]=0.5*pow(meanvel,(double)2)*( 1939 pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) + 1940 pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) ); 1941 } 1942 1943 /*Process units: */ 1944 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters); 1945 1946 /*Apply weights to misfits*/ 1947 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 1948 1949 /* Start looping on the number of gaussian points: */ 1950 gauss=new GaussTria(2); 1951 for (ig=gauss->begin();ig<gauss->end();ig++){ 1952 1953 gauss->GaussPoint(ig); 1954 1955 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1956 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 1957 Jelem+=misfit*Jdet*gauss->weight; 1958 } 1959 1960 /*clean-up and Return: */ 1961 delete gauss; 1962 return Jelem; 1963 } 1964 /*}}}*/ 1965 /*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/ 1966 double Tria::SurfaceAverageVelMisfit(bool process_units){ 1967 1968 const int numdof=2*NUMVERTICES; 1969 1970 int i,ig; 1971 int fit=-1; 1972 double Jelem=0,S=0; 1973 double scalex=1, scaley=1; 1974 double meanvel, epsvel,Jdet; 1975 double velocity_mag,obs_velocity_mag,misfit; 1976 double xyz_list[NUMVERTICES][3]; 1977 double vx_list[NUMVERTICES]; 1978 double vy_list[NUMVERTICES]; 1979 double obs_vx_list[NUMVERTICES]; 1980 double obs_vy_list[NUMVERTICES]; 1981 double misfit_square_list[NUMVERTICES]; 1982 double misfit_list[NUMVERTICES]; 1983 double weights_list[NUMVERTICES]; 1984 GaussTria *gauss=NULL; 1985 1986 /*If on water, return 0: */ 1987 if(IsOnWater())return 0; 1988 1989 /* Get node coordinates and dof list: */ 1990 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1991 1992 /* Recover input data: */ 1993 inputs->GetParameterValue(&S,SurfaceAreaEnum); 1994 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1995 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1996 GetParameterListOnVertices(&vx_list[0],VxEnum); 1997 GetParameterListOnVertices(&vy_list[0],VyEnum); 1998 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1999 2000 /*retrieve some parameters: */ 2001 this->parameters->FindParam(&meanvel,MeanVelEnum); 2002 this->parameters->FindParam(&epsvel,EpsVelEnum); 2003 2004 /* Compute SurfaceAverageVelMisfit at the 3 nodes 2005 * Here we integrate linearized functions: 2006 * 2007 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2008 * 2009 * where J_i are the misfits at the 3 nodes of the triangle 2010 * Phi_i is the nodal function (P1) with respect to 2011 * the vertex i 2012 */ 2013 2014 /*We are using a spacially average absolute misfit: 2015 * 2016 * 1 2 2 2017 * J = --- sqrt( (u - u ) + (v - v ) ) 2018 * S obs obs 2019 */ 2020 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); 2021 2022 /*Process units: */ 2023 if(process_units)UnitConversion(&misfit_square_list[0],NUMVERTICES,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters); 2024 2025 /*Take the square root, and scale by surface: */ 2026 for (i=0;i<NUMVERTICES;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S; 2027 2028 /*Apply weights to misfits*/ 2029 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 2030 2031 /* Start looping on the number of gaussian points: */ 2032 gauss=new GaussTria(2); 2033 for (ig=gauss->begin();ig<gauss->end();ig++){ 2034 2035 gauss->GaussPoint(ig); 2036 2037 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 2038 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 2039 Jelem+=misfit*Jdet*gauss->weight; 2040 } 2041 2042 /*clean-up and Return: */ 2043 delete gauss; 2044 return Jelem; 2045 } 2046 /*}}}*/ 2047 /*FUNCTION Tria::PatchFill{{{1*/ 2048 void Tria::PatchFill(int* prow, Patch* patch){ 2049 2050 int i,row; 2051 int vertices_ids[3]; 2052 2053 /*recover pointer: */ 2054 row=*prow; 2055 2056 for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. 2057 2058 for(i=0;i<this->results->Size();i++){ 2059 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2060 2061 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand 2062 *it to the result object, to fill the rest: */ 2063 patch->fillelementinfo(row,this->id,vertices_ids,3); 2064 elementresult->PatchFill(row,patch); 2065 2066 /*increment rower: */ 2067 row++; 2068 } 2069 2070 /*Assign output pointers:*/ 2071 *prow=row; 2072 } 2073 /*}}}*/ 2074 /*FUNCTION Tria::PatchSize{{{1*/ 2075 void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ 2076 2077 int i; 2078 int numrows = 0; 2079 int numnodes = 0; 2080 2081 /*Go through all the results objects, and update the counters: */ 2082 for (i=0;i<this->results->Size();i++){ 2083 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2084 /*first, we have one more result: */ 2085 numrows++; 2086 /*now, how many vertices and how many nodal values for this result? :*/ 2087 numnodes=elementresult->NumberOfNodalValues(); //ask result object. 2088 } 2089 2090 /*Assign output pointers:*/ 2091 *pnumrows=numrows; 2092 *pnumvertices=NUMVERTICES; 2093 *pnumnodes=numnodes; 2094 } 2095 /*}}}*/ 2096 /*FUNCTION Tria::ProcessResultsUnits{{{1*/ 2097 void Tria::ProcessResultsUnits(void){ 2098 2099 int i; 2100 2101 for(i=0;i<this->results->Size();i++){ 2102 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2103 elementresult->ProcessUnits(this->parameters); 2104 } 2105 } 2106 /*}}}*/ 2107 /*FUNCTION Tria::SurfaceArea {{{1*/ 2108 double Tria::SurfaceArea(void){ 2109 2110 int i; 2111 double S; 2112 double normal[3]; 2113 double v13[3],v23[3]; 2114 double xyz_list[NUMVERTICES][3]; 2115 2116 /*If on water, return 0: */ 2117 if(IsOnWater())return 0; 2118 2119 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 2120 2121 for (i=0;i<3;i++){ 2122 v13[i]=xyz_list[0][i]-xyz_list[2][i]; 2123 v23[i]=xyz_list[1][i]-xyz_list[2][i]; 2124 } 2125 2126 normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; 2127 normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; 2128 normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; 2129 2130 S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); 2131 2132 /*Return: */ 2133 return S; 2134 } 2135 /*}}}*/ 2136 /*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ 2137 void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index 2138 2139 /*Intermediaries*/ 2140 int i,j; 2141 int tria_node_ids[3]; 2142 int tria_vertex_ids[3]; 2143 int tria_type; 2144 double nodeinputs[3]; 2145 2146 /*Checks if debuging*/ 2147 /*{{{2*/ 2148 ISSMASSERT(iomodel->elements); 2149 /*}}}*/ 2150 2151 /*Recover element type*/ 2152 if ((analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum) && iomodel->prognostic_DG){ 2153 /*P1 Discontinuous Galerkin*/ 2154 tria_type=P1DGEnum; 2155 } 2156 else{ 2157 /*P1 Continuous Galerkin*/ 2158 tria_type=P1Enum; 2159 } 2160 this->SetElementType(tria_type,analysis_counter); 2161 2162 /*Recover vertices ids needed to initialize inputs*/ 2163 for(i=0;i<3;i++){ 2164 tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab 2165 } 2166 2167 /*Recover nodes ids needed to initialize the node hook.*/ 2168 if (tria_type==P1DGEnum){ 2169 /*Discontinuous Galerkin*/ 2170 tria_node_ids[0]=iomodel->nodecounter+3*index+1; 2171 tria_node_ids[1]=iomodel->nodecounter+3*index+2; 2172 tria_node_ids[2]=iomodel->nodecounter+3*index+3; 2173 } 2174 else{ 2175 /*Continuous Galerkin*/ 2176 for(i=0;i<3;i++){ 2177 tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab 2178 } 2179 } 2180 2181 /*hooks: */ 2182 this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type 2183 2184 /*Fill with IoModel*/ 2185 this->Update(index,iomodel); 2186 2187 /*Defaults if not provided in iomodel*/ 2188 switch(analysis_type){ 2189 2190 case DiagnosticHorizAnalysisEnum: 2191 2192 /*default vx,vy and vz: either observation or 0 */ 2193 if(!iomodel->vx){ 2194 if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; 2195 else for(i=0;i<3;i++)nodeinputs[i]=0; 2196 this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); 2197 this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); 2198 if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs)); 2199 } 2200 if(!iomodel->vy){ 2201 if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; 2202 else for(i=0;i<3;i++)nodeinputs[i]=0; 2203 this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); 2204 this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); 2205 if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs)); 2206 } 2207 if(!iomodel->vz){ 2208 if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; 2209 else for(i=0;i<3;i++)nodeinputs[i]=0; 2210 this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); 2211 this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); 2212 if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs)); 2213 } 2214 if(!iomodel->pressure){ 2215 for(i=0;i<3;i++)nodeinputs[i]=0; 2216 if(iomodel->qmu_analysis){ 2217 this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); 2218 this->inputs->AddInput(new TriaVertexInput(QmuPressureEnum,nodeinputs)); 2219 } 2220 } 2221 break; 2222 2223 default: 2224 /*No update for other solution types*/ 2225 break; 2226 2227 } 2228 2229 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. 2230 this->parameters=NULL; 2231 } 2232 /*}}}*/ 2233 /*FUNCTION Tria::Update(int index, IoModel* iomodel){{{1*/ 2234 void Tria::Update(int index, IoModel* iomodel){ //i is the element index 537 /*FUNCTION Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){{{1*/ 538 void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index 2235 539 2236 540 /*Intermediaries*/ … … 2375 679 } 2376 680 /*}}}*/ 681 682 /*Element virtual functions definitions: */ 683 /*FUNCTION Tria::AverageOntoPartition {{{1*/ 684 void Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){ 685 686 bool already=false; 687 int i,j; 688 int partition[NUMVERTICES]; 689 int offset[NUMVERTICES]; 690 double area; 691 double mean; 692 double values[3]; 693 694 /*First, get the area: */ 695 area=this->GetArea(); 696 697 /*Figure out the average for this element: */ 698 this->GetDofList1(&offset[0]); 699 mean=0; 700 for(i=0;i<NUMVERTICES;i++){ 701 partition[i]=(int)qmu_part[offset[i]]; 702 mean=mean+1.0/NUMVERTICES*vertex_response[offset[i]]; 703 } 704 705 /*Add contribution: */ 706 for(i=0;i<NUMVERTICES;i++){ 707 already=false; 708 for(j=0;j<i;j++){ 709 if (partition[i]==partition[j]){ 710 already=true; 711 break; 712 } 713 } 714 if(!already){ 715 VecSetValue(partition_contributions,partition[i],mean*area,ADD_VALUES); 716 VecSetValue(partition_areas,partition[i],area,ADD_VALUES); 717 }; 718 } 719 } 720 /*}}}*/ 721 /*FUNCTION Tria::ComputeBasalStress {{{1*/ 722 void Tria::ComputeBasalStress(Vec eps){ 723 ISSMERROR("Not Implemented yet"); 724 } 725 /*}}}*/ 726 /*FUNCTION Tria::ComputeStrainRate {{{1*/ 727 void Tria::ComputeStrainRate(Vec eps){ 728 ISSMERROR("Not Implemented yet"); 729 } 730 /*}}}*/ 731 /*FUNCTION Tria::SetCurrentConfiguration {{{1*/ 732 void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ 733 734 /*go into parameters and get the analysis_counter: */ 735 int analysis_counter; 736 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); 737 738 /*Get Element type*/ 739 this->element_type=this->element_type_list[analysis_counter]; 740 741 /*Pick up nodes*/ 742 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 743 else this->nodes=NULL; 744 745 } 746 /*}}}*/ 747 /*FUNCTION Tria::Configure {{{1*/ 748 void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ 749 750 /*go into parameters and get the analysis_counter: */ 751 int analysis_counter; 752 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); 753 754 /*Get Element type*/ 755 this->element_type=this->element_type_list[analysis_counter]; 756 757 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 758 * datasets, using internal ids and offsets hidden in hooks: */ 759 if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); 760 this->hmatice->configure(materialsin); 761 this->hmatpar->configure(materialsin); 762 763 /*Now, go pick up the objects inside the hooks: */ 764 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 765 else this->nodes=NULL; 766 this->matice=(Matice*)this->hmatice->delivers(); 767 this->matpar=(Matpar*)this->hmatpar->delivers(); 768 769 /*point parameters to real dataset: */ 770 this->parameters=parametersin; 771 772 } 773 /*}}}*/ 774 /*FUNCTION Tria::ControlInputGetGradient{{{1*/ 775 void Tria::ControlInputGetGradient(Vec gradient,int enum_type){ 776 777 int doflist1[NUMVERTICES]; 778 Input* input=NULL; 779 780 if(enum_type==RheologyBbarEnum){ 781 input=(Input*)matice->inputs->GetInput(enum_type); 782 } 783 else{ 784 input=inputs->GetInput(enum_type); 785 } 786 if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); 787 if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); 788 789 this->GetDofList1(&doflist1[0]); 790 ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); 791 792 }/*}}}*/ 793 /*FUNCTION Tria::ControlInputSetGradient{{{1*/ 794 void Tria::ControlInputSetGradient(double* gradient,int enum_type){ 795 796 int doflist1[NUMVERTICES]; 797 double grad_list[NUMVERTICES]; 798 Input* grad_input=NULL; 799 Input* input=NULL; 800 801 if(enum_type==RheologyBbarEnum){ 802 input=(Input*)matice->inputs->GetInput(enum_type); 803 } 804 else{ 805 input=inputs->GetInput(enum_type); 806 } 807 if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); 808 if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); 809 810 this->GetDofList1(&doflist1[0]); 811 for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]]; 812 grad_input=new TriaVertexInput(GradientEnum,grad_list); 813 814 ((ControlInput*)input)->SetGradient(grad_input); 815 816 }/*}}}*/ 817 /*FUNCTION Tria::RegularizeInversion {{{1*/ 818 double Tria::RegularizeInversion(){ 819 820 /*constants: */ 821 const int numdof=NDOF2*NUMVERTICES; 822 823 /* Intermediaries */ 824 int i,j,ig; 825 int num_controls; 826 int *control_type=NULL; 827 double Jelem = 0; 828 double cm_noisedmp; 829 double Jdet; 830 double xyz_list[NUMVERTICES][3]; 831 double dk[NDOF2],dB[NDOF2]; 832 GaussTria *gauss = NULL; 833 834 /*retrieve parameters and inputs*/ 835 this->parameters->FindParam(&num_controls,NumControlsEnum); 836 this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); 837 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); 838 839 /*If on water, return 0: */ 840 if(IsOnWater()) return 0; 841 842 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 843 844 for(i=0;i<num_controls;i++){ 845 /* Start looping on the number of gaussian points: */ 846 gauss=new GaussTria(2); 847 for (ig=gauss->begin();ig<gauss->end();ig++){ 848 849 gauss->GaussPoint(ig); 850 851 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 852 853 /*Add Tikhonov regularization term to misfit: 854 * 855 * WARNING: for now, the regularization is only taken into account by the gradient 856 * the misfit reflects the response only! 857 * 858 * */ 859 if (control_type[i]==DragCoefficientEnum){ 860 Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); 861 drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); 862 //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; 863 } 864 else if (control_type[i]==RheologyBbarEnum){ 865 //nothing 866 } 867 else if (control_type[i]==DhDtEnum){ 868 //nothing 869 } 870 else if (control_type[i]==VxEnum){ 871 //nothing 872 } 873 else if (control_type[i]==VyEnum){ 874 //nothing 875 } 876 else{ 877 ISSMERROR("unsupported control type: %s",EnumToString(control_type[i])); 878 } 879 } 880 881 delete gauss; 882 } 883 884 /*Clean up and return*/ 885 xfree((void**)&control_type); 886 return Jelem; 887 } 888 /*}}}*/ 889 /*FUNCTION Tria::CreateKMatrix {{{1*/ 890 void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ 891 892 /*retreive parameters: */ 893 ElementMatrix* Ke=NULL; 894 int analysis_type; 895 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 896 897 /*Checks in debugging mode{{{2*/ 898 ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); 899 /*}}}*/ 900 901 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ 902 switch(analysis_type){ 903 case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: 904 Ke=CreateKMatrixDiagnosticMacAyeal(); 905 break; 906 case DiagnosticHutterAnalysisEnum: 907 Ke=CreateKMatrixDiagnosticHutter(); 908 break; 909 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: 910 Ke=CreateKMatrixSlope(); 911 break; 912 case PrognosticAnalysisEnum: 913 Ke=CreateKMatrixPrognostic(); 914 break; 915 case BalancedthicknessAnalysisEnum: 916 Ke=CreateKMatrixBalancedthickness(); 917 break; 918 case AdjointBalancedthicknessAnalysisEnum: 919 Ke=CreateKMatrixAdjointBalancedthickness(); 920 break; 921 case BalancedvelocitiesAnalysisEnum: 922 Ke=CreateKMatrixBalancedvelocities(); 923 break; 924 default: 925 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); 926 } 927 928 /*Add to global matrix*/ 929 if(Ke){ 930 Ke->AddToGlobal(Kgg,Kff,Kfs); 931 delete Ke; 932 } 933 } 934 /*}}}*/ 935 /*FUNCTION Tria::CreatePVector {{{1*/ 936 void Tria::CreatePVector(Vec pg, Vec pf){ 937 938 /*retrive parameters: */ 939 ElementVector* pe=NULL; 940 int analysis_type; 941 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 942 943 /*asserts: {{{*/ 944 /*if debugging mode, check that all pointers exist*/ 945 ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); 946 /*}}}*/ 947 948 /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ 949 switch(analysis_type){ 950 case DiagnosticHorizAnalysisEnum: 951 pe=CreatePVectorDiagnosticMacAyeal(); 952 break; 953 case AdjointHorizAnalysisEnum: 954 pe=CreatePVectorAdjointHoriz(); 955 break; 956 case DiagnosticHutterAnalysisEnum: 957 pe=CreatePVectorDiagnosticHutter(); 958 break; 959 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: 960 pe=CreatePVectorSlope(); 961 break; 962 case PrognosticAnalysisEnum: 963 pe=CreatePVectorPrognostic(); 964 break; 965 case BalancedthicknessAnalysisEnum: 966 pe=CreatePVectorBalancedthickness(); 967 break; 968 case AdjointBalancedthicknessAnalysisEnum: 969 pe=CreatePVectorAdjointBalancedthickness(); 970 break; 971 case BalancedvelocitiesAnalysisEnum: 972 pe=CreatePVectorBalancedvelocities(); 973 break; 974 default: 975 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); 976 } 977 978 /*Add to global Vector*/ 979 if(pe){ 980 pe->AddToGlobal(pg,pf); 981 delete pe; 982 } 983 } 984 /*}}}*/ 985 /*FUNCTION Tria::DeleteResults {{{1*/ 986 void Tria::DeleteResults(void){ 987 988 /*Delete and reinitialize results*/ 989 delete this->results; 990 this->results=new Results(); 991 992 } 993 /*}}}*/ 994 /*FUNCTION Tria::GetNodeIndex {{{1*/ 995 int Tria::GetNodeIndex(Node* node){ 996 997 ISSMASSERT(nodes); 998 for(int i=0;i<NUMVERTICES;i++){ 999 if(node==nodes[i]) 1000 return i; 1001 } 1002 ISSMERROR("Node provided not found among element nodes"); 1003 } 1004 /*}}}*/ 1005 /*FUNCTION Tria::IsOnBed {{{1*/ 1006 bool Tria::IsOnBed(){ 1007 1008 bool onbed; 1009 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1010 return onbed; 1011 } 1012 /*}}}*/ 1013 /*FUNCTION Tria::IsOnShelf {{{1*/ 1014 bool Tria::IsOnShelf(){ 1015 1016 bool shelf; 1017 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 1018 return shelf; 1019 } 1020 /*}}}*/ 1021 /*FUNCTION Tria::IsOnWater {{{1*/ 1022 bool Tria::IsOnWater(){ 1023 1024 bool water; 1025 inputs->GetParameterValue(&water,ElementOnWaterEnum); 1026 return water; 1027 } 1028 /*}}}*/ 1029 /*FUNCTION Tria::GetSolutionFromInputs{{{1*/ 1030 void Tria::GetSolutionFromInputs(Vec solution){ 1031 1032 /*retrive parameters: */ 1033 int analysis_type; 1034 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 1035 1036 /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ 1037 if (analysis_type==DiagnosticHorizAnalysisEnum) 1038 GetSolutionFromInputsDiagnosticHoriz(solution); 1039 else if (analysis_type==DiagnosticHutterAnalysisEnum) 1040 GetSolutionFromInputsDiagnosticHutter(solution); 1041 else 1042 ISSMERROR("analysis: %s not supported yet",EnumToString(analysis_type)); 1043 1044 } 1045 /*}}}*/ 1046 /*FUNCTION Tria::GetVectorFromInputs{{{1*/ 1047 void Tria::GetVectorFromInputs(Vec vector,int NameEnum){ 1048 1049 int doflist1[NUMVERTICES]; 1050 1051 /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ 1052 for(int i=0;i<this->inputs->Size();i++){ 1053 Input* input=(Input*)this->inputs->GetObjectByOffset(i); 1054 if(input->EnumType()==NameEnum){ 1055 /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ 1056 this->GetDofList1(&doflist1[0]); 1057 input->GetVectorFromInputs(vector,&doflist1[0]); 1058 break; 1059 } 1060 } 1061 } 1062 /*}}}*/ 1063 /*FUNCTION Tria::Gradj {{{1*/ 1064 void Tria::Gradj(Vec gradient,int control_type){ 1065 1066 /*If on water, grad = 0: */ 1067 if(IsOnWater())return; 1068 1069 switch(control_type){ 1070 case DragCoefficientEnum: 1071 GradjDrag(gradient); 1072 break; 1073 case RheologyBbarEnum: 1074 GradjB(gradient); 1075 break; 1076 case DhDtEnum: 1077 GradjDhDt(gradient); 1078 break; 1079 case VxEnum: 1080 GradjVx(gradient); 1081 break; 1082 case VyEnum: 1083 GradjVy(gradient); 1084 break; 1085 default: 1086 ISSMERROR("%s%i","control type not supported yet: ",control_type); 1087 } 1088 } 1089 /*}}}*/ 1090 /*FUNCTION Tria::GradjB{{{1*/ 1091 void Tria::GradjB(Vec gradient){ 1092 1093 /*Intermediaries*/ 1094 int i,ig; 1095 int doflist[NUMVERTICES]; 1096 double vx,vy,lambda,mu,thickness,Jdet; 1097 double cm_noisedmp,viscosity_complement; 1098 double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; 1099 double xyz_list[NUMVERTICES][3]; 1100 double basis[3],epsilon[3]; 1101 double dbasis[NDOF2][NUMVERTICES]; 1102 double grad_g[NUMVERTICES]; 1103 double grad[NUMVERTICES]={0.0}; 1104 GaussTria *gauss = NULL; 1105 1106 /*retrieve some parameters: */ 1107 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); 1108 1109 /* Get node coordinates and dof list: */ 1110 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1111 GetDofList1(&doflist[0]); 1112 1113 /*Retrieve all inputs*/ 1114 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1115 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 1116 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 1117 Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); 1118 Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); 1119 Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(rheologyb_input); 1120 1121 /* Start looping on the number of gaussian points: */ 1122 gauss=new GaussTria(4); 1123 for (ig=gauss->begin();ig<gauss->end();ig++){ 1124 1125 gauss->GaussPoint(ig); 1126 1127 thickness_input->GetParameterValue(&thickness,gauss); 1128 rheologyb_input->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],gauss); 1129 vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); 1130 vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); 1131 adjointx_input->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); 1132 adjointy_input->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); 1133 1134 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 1135 matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); 1136 1137 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1138 GetNodalFunctions(basis,gauss); 1139 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); 1140 1141 /*standard gradient dJ/dki*/ 1142 for (i=0;i<NUMVERTICES;i++){ 1143 grad_g[i]=-viscosity_complement*thickness*( (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1])*Jdet*gauss->weight*basis[i]; 1144 } 1145 /*Add regularization term*/ 1146 for (i=0;i<NUMVERTICES;i++) grad_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); 1147 for(i=0;i<NUMVERTICES;i++) grad[i]+=grad_g[i]; 1148 } 1149 1150 VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES); 1151 1152 /*clean-up*/ 1153 delete gauss; 1154 } 1155 /*}}}*/ 1156 /*FUNCTION Tria::GradjDrag {{{1*/ 1157 void Tria::GradjDrag(Vec gradient){ 1158 1159 int i,ig; 1160 int drag_type,analysis_type; 1161 int doflist1[NUMVERTICES]; 1162 double vx,vy,lambda,mu,alpha_complement,Jdet; 1163 double bed,thickness,Neff,drag,cm_noisedmp; 1164 double xyz_list[NUMVERTICES][3]; 1165 double dh1dh3[NDOF2][NUMVERTICES]; 1166 double dk[NDOF2]; 1167 double grade_g[NUMVERTICES]={0.0}; 1168 double grade_g_gaussian[NUMVERTICES]; 1169 double l1l2l3[3]; 1170 double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ 1171 Friction* friction=NULL; 1172 GaussTria *gauss=NULL; 1173 1174 /*retrive parameters: */ 1175 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 1176 1177 /*retrieve some parameters ands return if iceshelf: */ 1178 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); 1179 if(IsOnShelf())return; 1180 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1181 GetDofList1(&doflist1[0]); 1182 1183 /*Build frictoin element, needed later: */ 1184 inputs->GetParameterValue(&drag_type,DragTypeEnum); 1185 friction=new Friction("2d",inputs,matpar,analysis_type); 1186 1187 /*Retrieve all inputs we will be needing: */ 1188 Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); 1189 Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); 1190 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 1191 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 1192 Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(dragcoefficient_input); 1193 1194 /* Start looping on the number of gaussian points: */ 1195 gauss=new GaussTria(4); 1196 for (ig=gauss->begin();ig<gauss->end();ig++){ 1197 1198 gauss->GaussPoint(ig); 1199 1200 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1201 GetNodalFunctions(l1l2l3, gauss); 1202 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); 1203 1204 /*Build alpha_complement_list: */ 1205 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); 1206 else alpha_complement=0; 1207 1208 dragcoefficient_input->GetParameterValue(&drag, gauss); 1209 adjointx_input->GetParameterValue(&lambda, gauss); 1210 adjointy_input->GetParameterValue(&mu, gauss); 1211 vx_input->GetParameterValue(&vx,gauss); 1212 vy_input->GetParameterValue(&vy,gauss); 1213 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); 1214 1215 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ 1216 for (i=0;i<NUMVERTICES;i++){ 1217 1218 //standard term dJ/dki 1219 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i]; 1220 1221 //noise dampening d/dki(1/2*(dk/dx)^2) 1222 grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); 1223 } 1224 1225 /*Add gradje_g_gaussian vector to gradje_g: */ 1226 for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i]; 1227 } 1228 1229 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 1230 1231 /*Clean up and return*/ 1232 delete gauss; 1233 delete friction; 1234 } 1235 /*}}}*/ 1236 /*FUNCTION Tria::GradjDhDt{{{1*/ 1237 void Tria::GradjDhDt(Vec gradient){ 1238 1239 /*Intermediaries*/ 1240 int doflist1[NUMVERTICES]; 1241 double lambda[NUMVERTICES]; 1242 double gradient_g[NUMVERTICES]; 1243 1244 GetDofList1(&doflist1[0]); 1245 1246 /*Compute Gradient*/ 1247 GetParameterListOnVertices(&lambda[0],AdjointEnum); 1248 for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i]; 1249 1250 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)gradient_g,INSERT_VALUES); 1251 } 1252 /*}}}*/ 1253 /*FUNCTION Tria::GradjVx{{{1*/ 1254 void Tria::GradjVx(Vec gradient){ 1255 1256 /*Intermediaries*/ 1257 int i,ig; 1258 int doflist1[NUMVERTICES]; 1259 double thickness,Jdet; 1260 double l1l2l3[3]; 1261 double Dlambda[2]; 1262 double xyz_list[NUMVERTICES][3]; 1263 double grade_g[NUMVERTICES] = {0.0}; 1264 GaussTria *gauss = NULL; 1265 1266 /* Get node coordinates and dof list: */ 1267 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1268 GetDofList1(&doflist1[0]); 1269 1270 /*Retrieve all inputs we will be needing: */ 1271 Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); 1272 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1273 1274 /* Start looping on the number of gaussian points: */ 1275 gauss=new GaussTria(2); 1276 for (ig=gauss->begin();ig<gauss->end();ig++){ 1277 1278 gauss->GaussPoint(ig); 1279 1280 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1281 GetNodalFunctions(l1l2l3, gauss); 1282 1283 adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); 1284 thickness_input->GetParameterValue(&thickness, gauss); 1285 1286 for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i]; 1287 } 1288 1289 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 1290 1291 /*Clean up and return*/ 1292 delete gauss; 1293 } 1294 /*}}}*/ 1295 /*FUNCTION Tria::GradjVy{{{1*/ 1296 void Tria::GradjVy(Vec gradient){ 1297 1298 /*Intermediaries*/ 1299 int i,ig; 1300 int doflist1[NUMVERTICES]; 1301 double thickness,Jdet; 1302 double l1l2l3[3]; 1303 double Dlambda[2]; 1304 double xyz_list[NUMVERTICES][3]; 1305 double grade_g[NUMVERTICES] = {0.0}; 1306 GaussTria *gauss = NULL; 1307 1308 /* Get node coordinates and dof list: */ 1309 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1310 GetDofList1(&doflist1[0]); 1311 1312 /*Retrieve all inputs we will be needing: */ 1313 Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); 1314 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1315 1316 /* Start looping on the number of gaussian points: */ 1317 gauss=new GaussTria(2); 1318 for (ig=gauss->begin();ig<gauss->end();ig++){ 1319 1320 gauss->GaussPoint(ig); 1321 1322 adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); 1323 thickness_input->GetParameterValue(&thickness, gauss); 1324 1325 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1326 GetNodalFunctions(l1l2l3, gauss); 1327 1328 for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[1]*Jdet*gauss->weight*l1l2l3[i]; 1329 } 1330 1331 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 1332 1333 /*Clean up and return*/ 1334 delete gauss; 1335 } 1336 /*}}}*/ 1337 /*FUNCTION Tria::InputControlUpdate{{{1*/ 1338 void Tria::InputControlUpdate(double scalar,bool save_parameter){ 1339 1340 /*Intermediary*/ 1341 int num_controls; 1342 int* control_type=NULL; 1343 Input* input=NULL; 1344 double cm_min,cm_max; 1345 1346 /*retrieve some parameters: */ 1347 this->parameters->FindParam(&cm_min,CmMinEnum); 1348 this->parameters->FindParam(&cm_max,CmMaxEnum); 1349 this->parameters->FindParam(&num_controls,NumControlsEnum); 1350 this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); 1351 1352 for(int i=0;i<num_controls;i++){ 1353 1354 if(control_type[i]==RheologyBbarEnum){ 1355 input=(Input*)matice->inputs->GetInput(control_type[i]); ISSMASSERT(input); 1356 } 1357 else{ 1358 input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input); 1359 } 1360 1361 if (input->Enum()!=ControlInputEnum){ 1362 ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i])); 1363 } 1364 1365 ((ControlInput*)input)->UpdateValue(scalar); 1366 input->Constrain(cm_min,cm_max); 1367 if (save_parameter) ((ControlInput*)input)->SaveValue(); 1368 1369 } 1370 1371 /*Clean up and return*/ 1372 xfree((void**)&control_type); 1373 } 1374 /*}}}*/ 1375 /*FUNCTION Tria::InputConvergence{{{1*/ 1376 bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ 1377 1378 bool converged=true; 1379 int i; 1380 Input** new_inputs=NULL; 1381 Input** old_inputs=NULL; 1382 1383 new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs 1384 old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs 1385 1386 for(i=0;i<num_enums/2;i++){ 1387 new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]); 1388 old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]); 1389 if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); 1390 if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); 1391 } 1392 1393 /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ 1394 for(i=0;i<num_criterionenums;i++){ 1395 IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]); 1396 if(eps[i]>criterionvalues[i]) converged=false; 1397 } 1398 1399 /*clean up and return*/ 1400 xfree((void**)&new_inputs); 1401 xfree((void**)&old_inputs); 1402 return converged; 1403 } 1404 /*}}}*/ 1405 /*FUNCTION Tria::InputDepthAverageAtBase {{{1*/ 1406 void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ 1407 1408 /*New input*/ 1409 Input* oldinput=NULL; 1410 Input* newinput=NULL; 1411 1412 /*copy input of enum_type*/ 1413 if (object_enum==ElementsEnum) 1414 oldinput=(Input*)this->inputs->GetInput(enum_type); 1415 else if (object_enum==MaterialsEnum) 1416 oldinput=(Input*)this->matice->inputs->GetInput(enum_type); 1417 else 1418 ISSMERROR("object %s not supported yet",EnumToString(object_enum)); 1419 if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumToString(enum_type)); 1420 newinput=(Input*)oldinput->copy(); 1421 1422 /*Assign new name (average)*/ 1423 newinput->ChangeEnum(average_enum_type); 1424 1425 /*Add new input to current element*/ 1426 if (object_enum==ElementsEnum) 1427 this->inputs->AddInput((Input*)newinput); 1428 else if (object_enum==MaterialsEnum) 1429 this->matice->inputs->AddInput((Input*)newinput); 1430 else 1431 ISSMERROR("object %s not supported yet",EnumToString(object_enum)); 1432 } 1433 /*}}}*/ 1434 /*FUNCTION Tria::InputDuplicate{{{1*/ 1435 void Tria::InputDuplicate(int original_enum,int new_enum){ 1436 1437 /*Call inputs method*/ 1438 if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); 1439 1440 } 1441 /*}}}*/ 1442 /*FUNCTION Tria::InputScale{{{1*/ 1443 void Tria::InputScale(int enum_type,double scale_factor){ 1444 1445 Input* input=NULL; 1446 1447 /*Make a copy of the original input: */ 1448 input=(Input*)this->inputs->GetInput(enum_type); 1449 if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); 1450 1451 /*Scale: */ 1452 input->Scale(scale_factor); 1453 } 1454 /*}}}*/ 1455 /*FUNCTION Tria::InputArtificialNoise{{{1*/ 1456 void Tria::InputArtificialNoise(int enum_type,double min,double max){ 1457 1458 Input* input=NULL; 1459 1460 /*Make a copy of the original input: */ 1461 input=(Input*)this->inputs->GetInput(enum_type); 1462 if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); 1463 1464 /*ArtificialNoise: */ 1465 input->ArtificialNoise(min,max); 1466 } 1467 /*}}}*/ 1468 /*FUNCTION Tria::InputToResult{{{1*/ 1469 void Tria::InputToResult(int enum_type,int step,double time){ 1470 1471 int i; 1472 Input *input = NULL; 1473 1474 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ 1475 if (enum_type==RheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); 1476 else input=this->inputs->GetInput(enum_type); 1477 if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumToString(enum_type)); 1478 1479 /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 1480 * object out of the input, with the additional step and time information: */ 1481 this->results->AddObject((Object*)input->SpawnResult(step,time)); 1482 if(input->Enum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); 1483 } 1484 /*}}}*/ 1485 /*FUNCTION Tria::MassFlux {{{1*/ 1486 double Tria::MassFlux( double* segment,bool process_units){ 1487 1488 const int numdofs=2; 1489 1490 int i; 1491 double mass_flux=0; 1492 double xyz_list[NUMVERTICES][3]; 1493 double normal[2]; 1494 double length,rho_ice; 1495 double x1,y1,x2,y2,h1,h2; 1496 double vx1,vx2,vy1,vy2; 1497 GaussTria* gauss_1=NULL; 1498 GaussTria* gauss_2=NULL; 1499 1500 /*Get material parameters :*/ 1501 rho_ice=matpar->GetRhoIce(); 1502 1503 /*First off, check that this segment belongs to this element: */ 1504 if ((int)*(segment+4)!=this->id)ISSMERROR("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id); 1505 1506 /*Recover segment node locations: */ 1507 x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); 1508 1509 /*Get xyz list: */ 1510 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1511 1512 /*get area coordinates of 0 and 1 locations: */ 1513 gauss_1=new GaussTria(); 1514 gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); 1515 gauss_2=new GaussTria(); 1516 gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); 1517 1518 normal[0]=cos(atan2(x1-x2,y2-y1)); 1519 normal[1]=sin(atan2(x1-x2,y2-y1)); 1520 1521 length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); 1522 1523 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1524 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 1525 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 1526 1527 thickness_input->GetParameterValue(&h1, gauss_1); 1528 thickness_input->GetParameterValue(&h2, gauss_2); 1529 vx_input->GetParameterValue(&vx1,gauss_1); 1530 vx_input->GetParameterValue(&vx2,gauss_2); 1531 vy_input->GetParameterValue(&vy1,gauss_1); 1532 vy_input->GetParameterValue(&vy2,gauss_2); 1533 1534 mass_flux= rho_ice*length*( 1535 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ 1536 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] 1537 ); 1538 1539 /*Process units: */ 1540 mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum,this->parameters); 1541 1542 /*clean up and return:*/ 1543 delete gauss_1; 1544 delete gauss_2; 1545 return mass_flux; 1546 } 1547 /*}}}*/ 1548 /*FUNCTION Tria::MaxAbsVx{{{1*/ 1549 void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ 1550 1551 /*Get maximum:*/ 1552 double maxabsvx=this->inputs->MaxAbs(VxEnum); 1553 1554 /*process units if requested: */ 1555 if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum,this->parameters); 1556 1557 /*Assign output pointers:*/ 1558 *pmaxabsvx=maxabsvx; 1559 } 1560 /*}}}*/ 1561 /*FUNCTION Tria::MaxAbsVy{{{1*/ 1562 void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ 1563 1564 /*Get maximum:*/ 1565 double maxabsvy=this->inputs->MaxAbs(VyEnum); 1566 1567 /*process units if requested: */ 1568 if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum,this->parameters); 1569 1570 /*Assign output pointers:*/ 1571 *pmaxabsvy=maxabsvy; 1572 } 1573 /*}}}*/ 1574 /*FUNCTION Tria::MaxAbsVz{{{1*/ 1575 void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ 1576 1577 /*Get maximum:*/ 1578 double maxabsvz=this->inputs->MaxAbs(VzEnum); 1579 1580 /*process units if requested: */ 1581 if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum,this->parameters); 1582 1583 /*Assign output pointers:*/ 1584 *pmaxabsvz=maxabsvz; 1585 } 1586 /*}}}*/ 1587 /*FUNCTION Tria::MaxVel{{{1*/ 1588 void Tria::MaxVel(double* pmaxvel, bool process_units){ 1589 1590 /*Get maximum:*/ 1591 double maxvel=this->inputs->Max(VelEnum); 1592 1593 /*process units if requested: */ 1594 if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum,this->parameters); 1595 1596 /*Assign output pointers:*/ 1597 *pmaxvel=maxvel; 1598 } 1599 /*}}}*/ 1600 /*FUNCTION Tria::MaxVx{{{1*/ 1601 void Tria::MaxVx(double* pmaxvx, bool process_units){ 1602 1603 /*Get maximum:*/ 1604 double maxvx=this->inputs->Max(VxEnum); 1605 1606 /*process units if requested: */ 1607 if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum,this->parameters); 1608 1609 /*Assign output pointers:*/ 1610 *pmaxvx=maxvx; 1611 } 1612 /*}}}*/ 1613 /*FUNCTION Tria::MaxVy{{{1*/ 1614 void Tria::MaxVy(double* pmaxvy, bool process_units){ 1615 1616 /*Get maximum:*/ 1617 double maxvy=this->inputs->Max(VyEnum); 1618 1619 /*process units if requested: */ 1620 if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum,this->parameters); 1621 1622 /*Assign output pointers:*/ 1623 *pmaxvy=maxvy; 1624 1625 } 1626 /*}}}*/ 1627 /*FUNCTION Tria::MaxVz{{{1*/ 1628 void Tria::MaxVz(double* pmaxvz, bool process_units){ 1629 1630 /*Get maximum:*/ 1631 double maxvz=this->inputs->Max(VzEnum); 1632 1633 /*process units if requested: */ 1634 if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum,this->parameters); 1635 1636 /*Assign output pointers:*/ 1637 *pmaxvz=maxvz; 1638 } 1639 /*}}}*/ 1640 /*FUNCTION Tria::MinVel{{{1*/ 1641 void Tria::MinVel(double* pminvel, bool process_units){ 1642 1643 /*Get minimum:*/ 1644 double minvel=this->inputs->Min(VelEnum); 1645 1646 /*process units if requested: */ 1647 if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum,this->parameters); 1648 1649 /*Assign output pointers:*/ 1650 *pminvel=minvel; 1651 } 1652 /*}}}*/ 1653 /*FUNCTION Tria::MinVx{{{1*/ 1654 void Tria::MinVx(double* pminvx, bool process_units){ 1655 1656 /*Get minimum:*/ 1657 double minvx=this->inputs->Min(VxEnum); 1658 1659 /*process units if requested: */ 1660 if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum,this->parameters); 1661 1662 /*Assign output pointers:*/ 1663 *pminvx=minvx; 1664 } 1665 /*}}}*/ 1666 /*FUNCTION Tria::MinVy{{{1*/ 1667 void Tria::MinVy(double* pminvy, bool process_units){ 1668 1669 /*Get minimum:*/ 1670 double minvy=this->inputs->Min(VyEnum); 1671 1672 /*process units if requested: */ 1673 if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum,this->parameters); 1674 1675 /*Assign output pointers:*/ 1676 *pminvy=minvy; 1677 } 1678 /*}}}*/ 1679 /*FUNCTION Tria::MinVz{{{1*/ 1680 void Tria::MinVz(double* pminvz, bool process_units){ 1681 1682 /*Get minimum:*/ 1683 double minvz=this->inputs->Min(VzEnum); 1684 1685 /*process units if requested: */ 1686 if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum,this->parameters); 1687 1688 /*Assign output pointers:*/ 1689 *pminvz=minvz; 1690 } 1691 /*}}}*/ 1692 /*FUNCTION Tria::TimeAdapt{{{1*/ 1693 double Tria::TimeAdapt(void){ 1694 1695 /*intermediary: */ 1696 int i; 1697 double C,dt; 1698 double dx,dy; 1699 double maxx,minx; 1700 double maxy,miny; 1701 double maxabsvx,maxabsvy; 1702 double xyz_list[NUMVERTICES][3]; 1703 1704 /*get CFL coefficient:*/ 1705 this->parameters->FindParam(&C,CflCoefficientEnum); 1706 1707 /*Get for Vx and Vy, the max of abs value: */ 1708 this->MaxAbsVx(&maxabsvx,false); 1709 this->MaxAbsVy(&maxabsvy,false); 1710 1711 /* Get node coordinates and dof list: */ 1712 GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES); 1713 1714 minx=xyz_list[0][0]; 1715 maxx=xyz_list[0][0]; 1716 miny=xyz_list[0][1]; 1717 maxy=xyz_list[0][1]; 1718 1719 for(i=1;i<NUMVERTICES;i++){ 1720 if (xyz_list[i][0]<minx)minx=xyz_list[i][0]; 1721 if (xyz_list[i][0]>maxx)maxx=xyz_list[i][0]; 1722 if (xyz_list[i][1]<miny)miny=xyz_list[i][1]; 1723 if (xyz_list[i][1]>maxy)maxy=xyz_list[i][1]; 1724 } 1725 dx=maxx-minx; 1726 dy=maxy-miny; 1727 1728 /*CFL criterion: */ 1729 dt=C/(maxabsvy/dx+maxabsvy/dy); 1730 1731 return dt; 1732 } 1733 /*}}}*/ 1734 /*FUNCTION Tria::ThicknessAbsMisfit {{{1*/ 1735 double Tria::ThicknessAbsMisfit(bool process_units){ 1736 1737 /* Constants */ 1738 const int numdof=1*NUMVERTICES; 1739 1740 /*Intermediaries*/ 1741 int i,ig; 1742 double thickness,thicknessobs,weight; 1743 double Jdet; 1744 double Jelem = 0; 1745 double xyz_list[NUMVERTICES][3]; 1746 GaussTria *gauss = NULL; 1747 1748 /*If on water, return 0: */ 1749 if(IsOnWater())return 0; 1750 1751 /* Get node coordinates and dof list: */ 1752 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1753 1754 /*Retrieve all inputs we will be needing: */ 1755 Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 1756 Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); 1757 Input* weights_input =inputs->GetInput(WeightsEnum); ISSMASSERT(weights_input); 1758 1759 /* Start looping on the number of gaussian points: */ 1760 gauss=new GaussTria(2); 1761 for (ig=gauss->begin();ig<gauss->end();ig++){ 1762 1763 gauss->GaussPoint(ig); 1764 1765 /* Get Jacobian determinant: */ 1766 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1767 1768 /*Get parameters at gauss point*/ 1769 thickness_input->GetParameterValue(&thickness,gauss); 1770 thicknessobs_input->GetParameterValue(&thicknessobs,gauss); 1771 weights_input->GetParameterValue(&weight,gauss); 1772 1773 /*compute ThicknessAbsMisfit*/ 1774 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight; 1775 } 1776 1777 /* clean up and Return: */ 1778 delete gauss; 1779 return Jelem; 1780 } 1781 /*}}}*/ 1782 /*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ 1783 double Tria::SurfaceAbsVelMisfit(bool process_units){ 1784 1785 const int numdof=NDOF2*NUMVERTICES; 1786 1787 int i,ig; 1788 double Jelem=0; 1789 double velocity_mag,obs_velocity_mag; 1790 double meanvel, epsvel,misfit,Jdet; 1791 double xyz_list[NUMVERTICES][3]; 1792 double vx_list[NUMVERTICES]; 1793 double vy_list[NUMVERTICES]; 1794 double obs_vx_list[NUMVERTICES]; 1795 double obs_vy_list[NUMVERTICES]; 1796 double misfit_square_list[NUMVERTICES]; 1797 double misfit_list[NUMVERTICES]; 1798 double weights_list[NUMVERTICES]; 1799 GaussTria *gauss=NULL; 1800 1801 /*If on water, return 0: */ 1802 if(IsOnWater())return 0; 1803 1804 /* Get node coordinates and dof list: */ 1805 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1806 1807 /* Recover input data: */ 1808 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1809 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1810 GetParameterListOnVertices(&vx_list[0],VxEnum); 1811 GetParameterListOnVertices(&vy_list[0],VyEnum); 1812 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1813 1814 /*retrieve some parameters: */ 1815 this->parameters->FindParam(&meanvel,MeanVelEnum); 1816 this->parameters->FindParam(&epsvel,EpsVelEnum); 1817 1818 /* Compute SurfaceAbsVelMisfit at the 3 nodes 1819 * Here we integrate linearized functions: 1820 * 1821 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1822 * 1823 * where J_i are the misfits at the 3 nodes of the triangle 1824 * Phi_i is the nodal function (P1) with respect to 1825 * the vertex i 1826 */ 1827 1828 /*We are using an absolute misfit: 1829 * 1830 * 1 [ 2 2 ] 1831 * J = --- | (u - u ) + (v - v ) | 1832 * 2 [ obs obs ] 1833 * 1834 */ 1835 for (i=0;i<NUMVERTICES;i++){ 1836 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)); 1837 } 1838 /*Process units: */ 1839 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters); 1840 1841 /*Apply weights to misfits*/ 1842 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 1843 1844 /* Start looping on the number of gaussian points: */ 1845 gauss=new GaussTria(2); 1846 for (ig=gauss->begin();ig<gauss->end(); ig++){ 1847 1848 gauss->GaussPoint(ig); 1849 1850 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1851 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 1852 Jelem+=misfit*Jdet*gauss->weight; 1853 } 1854 1855 /*clean up and Return: */ 1856 delete gauss; 1857 return Jelem; 1858 } 1859 /*}}}*/ 1860 /*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ 1861 double Tria::SurfaceRelVelMisfit(bool process_units){ 1862 1863 const int numdof=2*NUMVERTICES; 1864 1865 int i,ig; 1866 double Jelem=0; 1867 double scalex=1,scaley=1; 1868 double meanvel, epsvel,misfit,Jdet; 1869 double velocity_mag,obs_velocity_mag; 1870 double xyz_list[NUMVERTICES][3]; 1871 double vx_list[NUMVERTICES]; 1872 double vy_list[NUMVERTICES]; 1873 double obs_vx_list[NUMVERTICES]; 1874 double obs_vy_list[NUMVERTICES]; 1875 double misfit_square_list[NUMVERTICES]; 1876 double misfit_list[NUMVERTICES]; 1877 double weights_list[NUMVERTICES]; 1878 GaussTria *gauss=NULL; 1879 1880 /*If on water, return 0: */ 1881 if(IsOnWater())return 0; 1882 1883 /* Get node coordinates and dof list: */ 1884 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1885 1886 /* Recover input data: */ 1887 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1888 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1889 GetParameterListOnVertices(&vx_list[0],VxEnum); 1890 GetParameterListOnVertices(&vy_list[0],VyEnum); 1891 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1892 1893 /*retrieve some parameters: */ 1894 this->parameters->FindParam(&meanvel,MeanVelEnum); 1895 this->parameters->FindParam(&epsvel,EpsVelEnum); 1896 1897 /* Compute SurfaceRelVelMisfit at the 3 nodes 1898 * Here we integrate linearized functions: 1899 * 1900 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1901 * 1902 * where J_i are the misfits at the 3 nodes of the triangle 1903 * Phi_i is the nodal function (P1) with respect to 1904 * the vertex i 1905 */ 1906 1907 /*We are using a relative misfit: 1908 * 1909 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1910 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1911 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1912 * obs obs 1913 */ 1914 for (i=0;i<NUMVERTICES;i++){ 1915 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2); 1916 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2); 1917 if(obs_vx_list[i]==0)scalex=0; 1918 if(obs_vy_list[i]==0)scaley=0; 1919 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)); 1920 } 1921 1922 /*Process units: */ 1923 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters); 1924 1925 /*Apply weights to misfits*/ 1926 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 1927 1928 /* Start looping on the number of gaussian points: */ 1929 gauss=new GaussTria(2); 1930 for (ig=gauss->begin();ig<gauss->end();ig++){ 1931 1932 gauss->GaussPoint(ig); 1933 1934 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1935 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 1936 Jelem+=misfit*Jdet*gauss->weight; 1937 } 1938 1939 /*clean up and Return: */ 1940 delete gauss; 1941 return Jelem; 1942 } 1943 /*}}}*/ 1944 /*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/ 1945 double Tria::SurfaceLogVelMisfit(bool process_units){ 1946 1947 const int numdof=NDOF2*NUMVERTICES; 1948 1949 int i,ig; 1950 double Jelem=0; 1951 double scalex=1,scaley=1; 1952 double meanvel, epsvel,misfit,Jdet; 1953 double velocity_mag,obs_velocity_mag; 1954 double xyz_list[NUMVERTICES][3]; 1955 double vx_list[NUMVERTICES]; 1956 double vy_list[NUMVERTICES]; 1957 double obs_vx_list[NUMVERTICES]; 1958 double obs_vy_list[NUMVERTICES]; 1959 double misfit_square_list[NUMVERTICES]; 1960 double misfit_list[NUMVERTICES]; 1961 double weights_list[NUMVERTICES]; 1962 GaussTria *gauss=NULL; 1963 1964 /*If on water, return 0: */ 1965 if(IsOnWater())return 0; 1966 1967 /* Get node coordinates and dof list: */ 1968 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 1969 1970 /* Recover input data: */ 1971 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 1972 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 1973 GetParameterListOnVertices(&vx_list[0],VxEnum); 1974 GetParameterListOnVertices(&vy_list[0],VyEnum); 1975 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 1976 1977 /*retrieve some parameters: */ 1978 this->parameters->FindParam(&meanvel,MeanVelEnum); 1979 this->parameters->FindParam(&epsvel,EpsVelEnum); 1980 1981 /* Compute SurfaceLogVelMisfit at the 3 nodes 1982 * Here we integrate linearized functions: 1983 * 1984 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 1985 * 1986 * where J_i are the misfits at the 3 nodes of the triangle 1987 * Phi_i is the nodal function (P1) with respect to 1988 * the vertex i 1989 */ 1990 1991 /*We are using a logarithmic misfit: 1992 * 1993 * [ vel + eps ] 2 1994 * J = 4 \bar{v}^2 | log ( ----------- ) | 1995 * [ vel + eps ] 1996 * obs 1997 */ 1998 for (i=0;i<NUMVERTICES;i++){ 1999 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil. 2000 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. 2001 misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); 2002 } 2003 2004 /*Process units: */ 2005 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters); 2006 2007 /*Apply weights to misfits*/ 2008 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 2009 2010 /* Start looping on the number of gaussian points: */ 2011 gauss=new GaussTria(2); 2012 for (ig=gauss->begin();ig<gauss->end();ig++){ 2013 2014 gauss->GaussPoint(ig); 2015 2016 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 2017 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 2018 Jelem+=misfit*Jdet*gauss->weight; 2019 } 2020 2021 /*clean-up and Return: */ 2022 delete gauss; 2023 return Jelem; 2024 } 2025 /*}}}*/ 2026 /*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/ 2027 double Tria::SurfaceLogVxVyMisfit(bool process_units){ 2028 2029 const int numdof=NDOF2*NUMVERTICES; 2030 2031 int i,ig; 2032 int fit=-1; 2033 double Jelem=0, S=0; 2034 double scalex=1,scaley=1; 2035 double meanvel, epsvel, misfit, Jdet; 2036 double velocity_mag,obs_velocity_mag; 2037 double xyz_list[NUMVERTICES][3]; 2038 double vx_list[NUMVERTICES]; 2039 double vy_list[NUMVERTICES]; 2040 double obs_vx_list[NUMVERTICES]; 2041 double obs_vy_list[NUMVERTICES]; 2042 double misfit_square_list[NUMVERTICES]; 2043 double misfit_list[NUMVERTICES]; 2044 double weights_list[NUMVERTICES]; 2045 GaussTria *gauss=NULL; 2046 2047 /*If on water, return 0: */ 2048 if(IsOnWater())return 0; 2049 2050 /* Get node coordinates and dof list: */ 2051 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 2052 2053 /* Recover input data: */ 2054 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 2055 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 2056 GetParameterListOnVertices(&vx_list[0],VxEnum); 2057 GetParameterListOnVertices(&vy_list[0],VyEnum); 2058 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 2059 2060 /*retrieve some parameters: */ 2061 this->parameters->FindParam(&meanvel,MeanVelEnum); 2062 this->parameters->FindParam(&epsvel,EpsVelEnum); 2063 2064 /* Compute SurfaceLogVxVyMisfit at the 3 nodes 2065 * Here we integrate linearized functions: 2066 * 2067 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2068 * 2069 * where J_i are the misfits at the 3 nodes of the triangle 2070 * Phi_i is the nodal function (P1) with respect to 2071 * the vertex i 2072 */ 2073 2074 /*We are using an logarithmic 2 misfit: 2075 * 2076 * 1 [ |u| + eps 2 |v| + eps 2 ] 2077 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 2078 * 2 [ |u |+ eps |v |+ eps ] 2079 * obs obs 2080 */ 2081 for (i=0;i<NUMVERTICES;i++){ 2082 misfit_list[i]=0.5*pow(meanvel,(double)2)*( 2083 pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) + 2084 pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) ); 2085 } 2086 2087 /*Process units: */ 2088 if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters); 2089 2090 /*Apply weights to misfits*/ 2091 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 2092 2093 /* Start looping on the number of gaussian points: */ 2094 gauss=new GaussTria(2); 2095 for (ig=gauss->begin();ig<gauss->end();ig++){ 2096 2097 gauss->GaussPoint(ig); 2098 2099 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 2100 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 2101 Jelem+=misfit*Jdet*gauss->weight; 2102 } 2103 2104 /*clean-up and Return: */ 2105 delete gauss; 2106 return Jelem; 2107 } 2108 /*}}}*/ 2109 /*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/ 2110 double Tria::SurfaceAverageVelMisfit(bool process_units){ 2111 2112 const int numdof=2*NUMVERTICES; 2113 2114 int i,ig; 2115 int fit=-1; 2116 double Jelem=0,S=0; 2117 double scalex=1, scaley=1; 2118 double meanvel, epsvel,Jdet; 2119 double velocity_mag,obs_velocity_mag,misfit; 2120 double xyz_list[NUMVERTICES][3]; 2121 double vx_list[NUMVERTICES]; 2122 double vy_list[NUMVERTICES]; 2123 double obs_vx_list[NUMVERTICES]; 2124 double obs_vy_list[NUMVERTICES]; 2125 double misfit_square_list[NUMVERTICES]; 2126 double misfit_list[NUMVERTICES]; 2127 double weights_list[NUMVERTICES]; 2128 GaussTria *gauss=NULL; 2129 2130 /*If on water, return 0: */ 2131 if(IsOnWater())return 0; 2132 2133 /* Get node coordinates and dof list: */ 2134 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 2135 2136 /* Recover input data: */ 2137 inputs->GetParameterValue(&S,SurfaceAreaEnum); 2138 GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); 2139 GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); 2140 GetParameterListOnVertices(&vx_list[0],VxEnum); 2141 GetParameterListOnVertices(&vy_list[0],VyEnum); 2142 GetParameterListOnVertices(&weights_list[0],WeightsEnum); 2143 2144 /*retrieve some parameters: */ 2145 this->parameters->FindParam(&meanvel,MeanVelEnum); 2146 this->parameters->FindParam(&epsvel,EpsVelEnum); 2147 2148 /* Compute SurfaceAverageVelMisfit at the 3 nodes 2149 * Here we integrate linearized functions: 2150 * 2151 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2152 * 2153 * where J_i are the misfits at the 3 nodes of the triangle 2154 * Phi_i is the nodal function (P1) with respect to 2155 * the vertex i 2156 */ 2157 2158 /*We are using a spacially average absolute misfit: 2159 * 2160 * 1 2 2 2161 * J = --- sqrt( (u - u ) + (v - v ) ) 2162 * S obs obs 2163 */ 2164 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); 2165 2166 /*Process units: */ 2167 if(process_units)UnitConversion(&misfit_square_list[0],NUMVERTICES,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters); 2168 2169 /*Take the square root, and scale by surface: */ 2170 for (i=0;i<NUMVERTICES;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S; 2171 2172 /*Apply weights to misfits*/ 2173 for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i]; 2174 2175 /* Start looping on the number of gaussian points: */ 2176 gauss=new GaussTria(2); 2177 for (ig=gauss->begin();ig<gauss->end();ig++){ 2178 2179 gauss->GaussPoint(ig); 2180 2181 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 2182 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); 2183 Jelem+=misfit*Jdet*gauss->weight; 2184 } 2185 2186 /*clean-up and Return: */ 2187 delete gauss; 2188 return Jelem; 2189 } 2190 /*}}}*/ 2191 /*FUNCTION Tria::PatchFill{{{1*/ 2192 void Tria::PatchFill(int* prow, Patch* patch){ 2193 2194 int i,row; 2195 int vertices_ids[3]; 2196 2197 /*recover pointer: */ 2198 row=*prow; 2199 2200 for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. 2201 2202 for(i=0;i<this->results->Size();i++){ 2203 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2204 2205 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand 2206 *it to the result object, to fill the rest: */ 2207 patch->fillelementinfo(row,this->id,vertices_ids,3); 2208 elementresult->PatchFill(row,patch); 2209 2210 /*increment rower: */ 2211 row++; 2212 } 2213 2214 /*Assign output pointers:*/ 2215 *prow=row; 2216 } 2217 /*}}}*/ 2218 /*FUNCTION Tria::PatchSize{{{1*/ 2219 void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ 2220 2221 int i; 2222 int numrows = 0; 2223 int numnodes = 0; 2224 2225 /*Go through all the results objects, and update the counters: */ 2226 for (i=0;i<this->results->Size();i++){ 2227 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2228 /*first, we have one more result: */ 2229 numrows++; 2230 /*now, how many vertices and how many nodal values for this result? :*/ 2231 numnodes=elementresult->NumberOfNodalValues(); //ask result object. 2232 } 2233 2234 /*Assign output pointers:*/ 2235 *pnumrows=numrows; 2236 *pnumvertices=NUMVERTICES; 2237 *pnumnodes=numnodes; 2238 } 2239 /*}}}*/ 2240 /*FUNCTION Tria::ProcessResultsUnits{{{1*/ 2241 void Tria::ProcessResultsUnits(void){ 2242 2243 int i; 2244 2245 for(i=0;i<this->results->Size();i++){ 2246 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2247 elementresult->ProcessUnits(this->parameters); 2248 } 2249 } 2250 /*}}}*/ 2251 /*FUNCTION Tria::SurfaceArea {{{1*/ 2252 double Tria::SurfaceArea(void){ 2253 2254 int i; 2255 double S; 2256 double normal[3]; 2257 double v13[3],v23[3]; 2258 double xyz_list[NUMVERTICES][3]; 2259 2260 /*If on water, return 0: */ 2261 if(IsOnWater())return 0; 2262 2263 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 2264 2265 for (i=0;i<3;i++){ 2266 v13[i]=xyz_list[0][i]-xyz_list[2][i]; 2267 v23[i]=xyz_list[1][i]-xyz_list[2][i]; 2268 } 2269 2270 normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; 2271 normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; 2272 normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; 2273 2274 S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); 2275 2276 /*Return: */ 2277 return S; 2278 } 2279 /*}}}*/ 2280 /*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ 2281 void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index 2282 2283 /*Intermediaries*/ 2284 int i,j; 2285 int tria_node_ids[3]; 2286 int tria_vertex_ids[3]; 2287 int tria_type; 2288 double nodeinputs[3]; 2289 2290 /*Checks if debuging*/ 2291 /*{{{2*/ 2292 ISSMASSERT(iomodel->elements); 2293 /*}}}*/ 2294 2295 /*Recover element type*/ 2296 if ((analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum) && iomodel->prognostic_DG){ 2297 /*P1 Discontinuous Galerkin*/ 2298 tria_type=P1DGEnum; 2299 } 2300 else{ 2301 /*P1 Continuous Galerkin*/ 2302 tria_type=P1Enum; 2303 } 2304 this->SetElementType(tria_type,analysis_counter); 2305 2306 /*Recover vertices ids needed to initialize inputs*/ 2307 for(i=0;i<3;i++){ 2308 tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab 2309 } 2310 2311 /*Recover nodes ids needed to initialize the node hook.*/ 2312 if (tria_type==P1DGEnum){ 2313 /*Discontinuous Galerkin*/ 2314 tria_node_ids[0]=iomodel->nodecounter+3*index+1; 2315 tria_node_ids[1]=iomodel->nodecounter+3*index+2; 2316 tria_node_ids[2]=iomodel->nodecounter+3*index+3; 2317 } 2318 else{ 2319 /*Continuous Galerkin*/ 2320 for(i=0;i<3;i++){ 2321 tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab 2322 } 2323 } 2324 2325 /*hooks: */ 2326 this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type 2327 2328 /*Fill with IoModel*/ 2329 this->InputUpdateFromIoModel(index,iomodel); 2330 2331 /*Defaults if not provided in iomodel*/ 2332 switch(analysis_type){ 2333 2334 case DiagnosticHorizAnalysisEnum: 2335 2336 /*default vx,vy and vz: either observation or 0 */ 2337 if(!iomodel->vx){ 2338 if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; 2339 else for(i=0;i<3;i++)nodeinputs[i]=0; 2340 this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); 2341 this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); 2342 if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs)); 2343 } 2344 if(!iomodel->vy){ 2345 if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; 2346 else for(i=0;i<3;i++)nodeinputs[i]=0; 2347 this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); 2348 this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); 2349 if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs)); 2350 } 2351 if(!iomodel->vz){ 2352 if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; 2353 else for(i=0;i<3;i++)nodeinputs[i]=0; 2354 this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); 2355 this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); 2356 if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs)); 2357 } 2358 if(!iomodel->pressure){ 2359 for(i=0;i<3;i++)nodeinputs[i]=0; 2360 if(iomodel->qmu_analysis){ 2361 this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); 2362 this->inputs->AddInput(new TriaVertexInput(QmuPressureEnum,nodeinputs)); 2363 } 2364 } 2365 break; 2366 2367 default: 2368 /*No update for other solution types*/ 2369 break; 2370 2371 } 2372 2373 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. 2374 this->parameters=NULL; 2375 } 2376 /*}}}*/ 2377 2377 /*FUNCTION Tria::UpdateGeometry{{{1*/ 2378 2378 void Tria::UpdateGeometry(void){ -
issm/trunk/src/c/objects/Elements/Tria.h
r6213 r6216 66 66 void InputUpdateFromConstant(int constant, int name); 67 67 void InputUpdateFromConstant(bool constant, int name); 68 void InputUpdateFromIoModel(int index, IoModel* iomodel); 68 69 /*}}}*/ 69 70 /*Element virtual functions definitions: {{{1*/ … … 122 123 double SurfaceArea(void); 123 124 void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); 124 void Update(int index, IoModel* iomodel);125 125 void UpdateGeometry(void); 126 126 double TimeAdapt(); -
issm/trunk/src/c/objects/Loads/Icefront.h
r6007 r6216 66 66 void InputUpdateFromConstant(bool constant, int name); 67 67 void InputUpdateFromSolution(double* solution); 68 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; 68 69 /*}}}*/ 69 70 /*Load virtual functions definitions: {{{1*/ -
issm/trunk/src/c/objects/Loads/Numericalflux.h
r6029 r6216 62 62 void InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} 63 63 void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} 64 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; 64 65 /*}}}*/ 65 66 /*Load virtual functions definitions: {{{1*/ -
issm/trunk/src/c/objects/Loads/Pengrid.h
r5937 r6216 67 67 void InputUpdateFromConstant(bool constant, int name); 68 68 void InputUpdateFromSolution(double* solution); 69 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; 69 70 /*}}}*/ 70 71 /*Load virtual functions definitions: {{{1*/ -
issm/trunk/src/c/objects/Loads/Penpair.h
r5940 r6216 54 54 void InputUpdateFromConstant(bool constant, int name); 55 55 void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} 56 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; 56 57 /*}}}*/ 57 58 /*Load virtual functions definitions: {{{1*/ -
issm/trunk/src/c/objects/Loads/Riftfront.h
r5941 r6216 73 73 void InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} 74 74 void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} 75 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; 75 76 /*}}}*/ 76 77 /*Load virtual functions definitions: {{{1*/ -
issm/trunk/src/c/objects/Materials/Material.h
r6213 r6216 22 22 virtual void InputDuplicate(int original_enum,int new_enum)=0; 23 23 virtual void Configure(Elements* elements)=0; 24 virtual void Update(int index, IoModel* iomodel)=0;25 24 26 25 }; -
issm/trunk/src/c/objects/Materials/Matice.cpp
r6213 r6216 248 248 249 249 return matice; 250 }251 /*}}}*/252 /*FUNCTION Matice::Update{{{1*/253 void Matice::Update(int index, IoModel* iomodel){254 255 int i,j;256 257 /*if 2d*/258 if(iomodel->dim==2){259 260 /*Intermediaries*/261 const int num_vertices = 3; //Tria has 3 vertices262 double nodeinputs[num_vertices];263 264 /*Control Inputs*/265 if (iomodel->control_analysis && iomodel->control_type){266 for(i=0;i<iomodel->num_control_type;i++){267 switch((int)iomodel->control_type[i]){268 case RheologyBbarEnum:269 if (iomodel->rheology_B){270 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+j]-1)];271 this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs));272 }273 break;274 }275 }276 }277 }278 279 /*if 3d*/280 else if(iomodel->dim==3){281 282 /*Intermediaries*/283 const int num_vertices = 6; //Penta has 6 vertices284 double nodeinputs[num_vertices];285 286 /*Control Inputs*/287 if (iomodel->control_analysis && iomodel->control_type){288 for(i=0;i<iomodel->num_control_type;i++){289 switch((int)iomodel->control_type[i]){290 case RheologyBbarEnum:291 if (iomodel->rheology_B){292 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+j]-1)];293 this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs));294 }295 break;296 }297 }298 }299 }300 else{301 ISSMERROR(" Mesh type not supported yet!");302 }303 304 return;305 250 } 306 251 /*}}}*/ … … 721 666 } 722 667 /*}}}*/ 668 /*FUNCTION Matice::InputUpdateFromIoModel{{{1*/ 669 void Matice::InputUpdateFromIoModel(int index, IoModel* iomodel){ 670 671 int i,j; 672 673 /*if 2d*/ 674 if(iomodel->dim==2){ 675 676 /*Intermediaries*/ 677 const int num_vertices = 3; //Tria has 3 vertices 678 double nodeinputs[num_vertices]; 679 680 /*Control Inputs*/ 681 if (iomodel->control_analysis && iomodel->control_type){ 682 for(i=0;i<iomodel->num_control_type;i++){ 683 switch((int)iomodel->control_type[i]){ 684 case RheologyBbarEnum: 685 if (iomodel->rheology_B){ 686 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; 687 this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs)); 688 } 689 break; 690 } 691 } 692 } 693 } 694 695 /*if 3d*/ 696 else if(iomodel->dim==3){ 697 698 /*Intermediaries*/ 699 const int num_vertices = 6; //Penta has 6 vertices 700 double nodeinputs[num_vertices]; 701 702 /*Control Inputs*/ 703 if (iomodel->control_analysis && iomodel->control_type){ 704 for(i=0;i<iomodel->num_control_type;i++){ 705 switch((int)iomodel->control_type[i]){ 706 case RheologyBbarEnum: 707 if (iomodel->rheology_B){ 708 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; 709 this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs)); 710 } 711 break; 712 } 713 } 714 } 715 } 716 else{ 717 ISSMERROR(" Mesh type not supported yet!"); 718 } 719 720 return; 721 } 722 /*}}}*/ 723 723 /*FUNCTION Matice::IsInput{{{1*/ 724 724 bool Matice::IsInput(int name){ -
issm/trunk/src/c/objects/Materials/Matice.h
r6213 r6216 52 52 void InputUpdateFromConstant(bool constant, int name); 53 53 void InputUpdateFromSolution(double* solution); 54 void Update(int index, IoModel* iomodel);54 void InputUpdateFromIoModel(int index, IoModel* iomodel); 55 55 /*}}}*/ 56 56 /*Material virtual functions resolution: {{{1*/ -
issm/trunk/src/c/objects/Materials/Matpar.h
r6213 r6216 56 56 void InputUpdateFromConstant(bool constant, int name); 57 57 void InputUpdateFromSolution(double* solution); 58 void Update(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};58 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; 59 59 /*}}}*/ 60 60 /*Material virtual functions resolution: {{{1*/ -
issm/trunk/src/c/objects/Node.h
r5843 r6216 62 62 void InputUpdateFromConstant(bool constant, int name); 63 63 void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} 64 void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("Not implemented yet!");} 64 65 /*}}}*/ 65 66 /*Node numerical routines {{{1*/ -
issm/trunk/src/c/objects/Update.h
r5311 r6216 25 25 virtual void InputUpdateFromConstant(bool constant, int name)=0; 26 26 virtual void InputUpdateFromSolution(double* solution)=0; 27 virtual void InputUpdateFromIoModel(int index, IoModel* iomodel)=0; 27 28 28 29 };
Note:
See TracChangeset
for help on using the changeset viewer.