Changeset 6216


Ignore:
Timestamp:
10/08/10 17:03:53 (14 years ago)
Author:
Mathieu Morlighem
Message:

Added method InputUpdateFromIoModel, cleaner that Update

Location:
issm/trunk/src/c
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r6213 r6216  
    5858                if(iomodel->my_elements[i]){
    5959                        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.
    6161
    6262                        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.
    6464                        counter++;
    6565                }

  • issm/trunk/src/c/objects/Elements/Element.h

    r6213 r6216  
    5656                virtual void   DeleteResults(void)=0;
    5757                virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0;
    58                 virtual void   Update(int index, IoModel* iomodel)=0;
    5958                virtual void   UpdateGeometry(void)=0;
    6059                virtual void   InputToResult(int enum_type,int step,double time)=0;

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

    r6213 r6216  
    470470void  Penta::InputUpdateFromVectorDakota(bool* vector, int name, int type){
    471471        ISSMERROR(" not supported yet!");
     472}
     473/*}}}*/
     474/*FUNCTION Penta::InputUpdateFromIoModel(int index,IoModel* iomodel) {{{1*/
     475void 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
    472642}
    473643/*}}}*/
    … …  
    17401910
    17411911        /*Fill with IoModel*/
    1742         this->Update(index,iomodel);
     1912        this->InputUpdateFromIoModel(index,iomodel);
    17431913
    17441914        /*Defaults if not provided in iomodel*/
    … …  
    18171987}
    18181988/*}}}*/
    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 Matlab
    1835         }
    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 element
    1960         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 /*}}}*/
    19891989/*FUNCTION Penta::UpdateGeometry{{{1*/
    19901990void  Penta::UpdateGeometry(void){

  • issm/trunk/src/c/objects/Elements/Penta.h

    r6213 r6216  
    6969                void  InputUpdateFromVectorDakota(double* vector, int name, int type);
    7070                void  InputUpdateFromVectorDakota(int* vector, int name, int type);
     71                void  InputUpdateFromIoModel(int index, IoModel* iomodel);
    7172                /*}}}*/
    7273                /*Element virtual functions definitions: {{{1*/
    … …  
    119120                double SurfaceArea(void);
    120121                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    121                 void   Update(int index, IoModel* iomodel);
    122122                void   UpdateGeometry(void);
    123123                double TimeAdapt();

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

    r6213 r6216  
    535535}
    536536/*}}}*/
    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*/
     538void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index
    2235539
    2236540        /*Intermediaries*/
    … …  
    2375679}
    2376680/*}}}*/
     681
     682/*Element virtual functions definitions: */
     683/*FUNCTION Tria::AverageOntoPartition {{{1*/
     684void  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*/
     722void  Tria::ComputeBasalStress(Vec eps){
     723        ISSMERROR("Not Implemented yet");
     724}
     725/*}}}*/
     726/*FUNCTION Tria::ComputeStrainRate {{{1*/
     727void  Tria::ComputeStrainRate(Vec eps){
     728        ISSMERROR("Not Implemented yet");
     729}
     730/*}}}*/
     731/*FUNCTION Tria::SetCurrentConfiguration {{{1*/
     732void  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*/
     748void  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*/
     775void 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*/
     794void 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*/
     818double 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*/
     890void  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*/
     936void  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*/
     986void  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*/
     995int 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*/
     1006bool Tria::IsOnBed(){
     1007       
     1008        bool onbed;
     1009        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1010        return onbed;
     1011}
     1012/*}}}*/
     1013/*FUNCTION Tria::IsOnShelf {{{1*/
     1014bool   Tria::IsOnShelf(){
     1015
     1016        bool shelf;
     1017        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     1018        return shelf;
     1019}
     1020/*}}}*/
     1021/*FUNCTION Tria::IsOnWater {{{1*/
     1022bool   Tria::IsOnWater(){
     1023
     1024        bool water;
     1025        inputs->GetParameterValue(&water,ElementOnWaterEnum);
     1026        return water;
     1027}
     1028/*}}}*/
     1029/*FUNCTION Tria::GetSolutionFromInputs{{{1*/
     1030void  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*/
     1047void  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*/
     1064void  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*/
     1091void  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*/
     1157void  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*/
     1237void  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*/
     1254void  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*/
     1296void  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*/
     1338void  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*/
     1376bool 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*/
     1406void  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*/
     1435void  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*/
     1443void  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*/
     1456void  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*/
     1469void  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*/
     1486double 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*/
     1549void  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*/
     1562void  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*/
     1575void  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*/
     1588void  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*/
     1601void  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*/
     1614void  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*/
     1628void  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*/
     1641void  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*/
     1654void  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*/
     1667void  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*/
     1680void  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*/
     1693double  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*/
     1735double 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*/
     1783double 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*/
     1861double 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*/
     1945double 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*/
     2027double 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*/
     2110double 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*/
     2192void  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*/
     2219void  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*/
     2241void  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*/
     2252double 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*/
     2281void 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/*}}}*/
    23772377/*FUNCTION Tria::UpdateGeometry{{{1*/
    23782378void  Tria::UpdateGeometry(void){

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

    r6213 r6216  
    6666                void  InputUpdateFromConstant(int constant, int name);
    6767                void  InputUpdateFromConstant(bool constant, int name);
     68                void  InputUpdateFromIoModel(int index, IoModel* iomodel);
    6869                /*}}}*/
    6970                /*Element virtual functions definitions: {{{1*/
    … …  
    122123                double SurfaceArea(void);
    123124                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    124                 void   Update(int index, IoModel* iomodel);
    125125                void   UpdateGeometry(void);
    126126                double TimeAdapt();

  • issm/trunk/src/c/objects/Loads/Icefront.h

    r6007 r6216  
    6666                void  InputUpdateFromConstant(bool constant, int name);
    6767                void  InputUpdateFromSolution(double* solution);
     68                void  InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};
    6869                /*}}}*/
    6970                /*Load virtual functions definitions: {{{1*/

  • issm/trunk/src/c/objects/Loads/Numericalflux.h

    r6029 r6216  
    6262                void    InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
    6363                void    InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");}
     64                void  InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};
    6465                /*}}}*/
    6566                /*Load virtual functions definitions: {{{1*/

  • issm/trunk/src/c/objects/Loads/Pengrid.h

    r5937 r6216  
    6767                void  InputUpdateFromConstant(bool constant, int name);
    6868                void  InputUpdateFromSolution(double* solution);
     69                void  InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};
    6970                /*}}}*/
    7071                /*Load virtual functions definitions: {{{1*/

  • issm/trunk/src/c/objects/Loads/Penpair.h

    r5940 r6216  
    5454                void  InputUpdateFromConstant(bool constant, int name);
    5555                void  InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");}
     56                void  InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};
    5657                /*}}}*/
    5758                        /*Load virtual functions definitions: {{{1*/

  • issm/trunk/src/c/objects/Loads/Riftfront.h

    r5941 r6216  
    7373                void    InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
    7474                void    InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");}
     75                void  InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};
    7576                /*}}}*/
    7677                /*Load virtual functions definitions: {{{1*/

  • issm/trunk/src/c/objects/Materials/Material.h

    r6213 r6216  
    2222                virtual void   InputDuplicate(int original_enum,int new_enum)=0;
    2323                virtual void   Configure(Elements* elements)=0;
    24                 virtual void   Update(int index, IoModel* iomodel)=0;
    2524
    2625};

  • issm/trunk/src/c/objects/Materials/Matice.cpp

    r6213 r6216  
    248248
    249249        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 vertices
    262                 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 vertices
    284                 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;
    305250}
    306251/*}}}*/
    … …  
    721666}
    722667/*}}}*/
     668/*FUNCTION Matice::InputUpdateFromIoModel{{{1*/
     669void 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/*}}}*/
    723723/*FUNCTION Matice::IsInput{{{1*/
    724724bool Matice::IsInput(int name){

  • issm/trunk/src/c/objects/Materials/Matice.h

    r6213 r6216  
    5252                void  InputUpdateFromConstant(bool constant, int name);
    5353                void  InputUpdateFromSolution(double* solution);
    54                 void  Update(int index, IoModel* iomodel);
     54                void  InputUpdateFromIoModel(int index, IoModel* iomodel);
    5555                /*}}}*/
    5656                /*Material virtual functions resolution: {{{1*/

  • issm/trunk/src/c/objects/Materials/Matpar.h

    r6213 r6216  
    5656                void   InputUpdateFromConstant(bool constant, int name);
    5757                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");};
    5959                /*}}}*/
    6060                /*Material virtual functions resolution: {{{1*/

  • issm/trunk/src/c/objects/Node.h

    r5843 r6216  
    6262                void  InputUpdateFromConstant(bool constant, int name);
    6363                void  InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");}
     64                void  InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("Not implemented yet!");}
    6465                /*}}}*/
    6566                /*Node numerical routines {{{1*/

  • issm/trunk/src/c/objects/Update.h

    r5311 r6216  
    2525                virtual void  InputUpdateFromConstant(bool constant, int name)=0;
    2626                virtual void  InputUpdateFromSolution(double* solution)=0;
     27                virtual void  InputUpdateFromIoModel(int index, IoModel* iomodel)=0;
    2728
    2829};
Note: See TracChangeset for help on using the changeset viewer.