Changeset 6213


Ignore:
Timestamp:
10/08/10 16:49:12 (14 years ago)
Author:
Mathieu Morlighem
Message:

Prepared ISSM for multivariable CM, Added modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

Location:
issm/trunk/src/c
Files:
1 added
28 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Parameters.cpp

    r4873 r6213  
    175175}
    176176/*}}}*/
     177/*FUNCTION Parameters::FindParam(int** pintarray,int* pM,int enum_type){{{1*/
     178int   Parameters::FindParam(int** pintarray,int* pM, int enum_type){
     179
     180        /*Go through a dataset, and find a Param* object
     181         *which parameter name is "name" : */
     182
     183        vector<Object*>::iterator object;
     184        Param* param=NULL;
     185
     186        int found=0;
     187
     188        for ( object=objects.begin() ; object < objects.end(); object++ ){
     189
     190                /*Ok, this object is a parameter, recover it and ask which name it has: */
     191                param=(Param*)(*object);
     192
     193                if(param->EnumType()==enum_type){
     194                        /*Ok, this is the one! Recover the value of this parameter: */
     195                        param->GetParameterValue(pintarray,pM);
     196                        found=1;
     197                        break;
     198                }
     199        }
     200        return found;
     201
     202}
     203/*}}}*/
    177204/*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM,int enum_type){{{1*/
    178205int   Parameters::FindParam(double** pdoublearray,int* pM, int enum_type){

  • issm/trunk/src/c/Container/Parameters.h

    r4873 r6213  
    3131                int   FindParam(char** pstring,int enum_type);
    3232                int   FindParam(char*** pstringarray,int* pM,int enum_type);
     33                int   FindParam(int** pintarray,int* pM,int enum_type);
    3334                int   FindParam(double** pdoublearray,int* pM,int enum_type);
    3435                int   FindParam(double** pdoublearray,int* pM,int* pN,int enum_type);

  • issm/trunk/src/c/Makefile.am

    r6200 r6213  
    383383                                        ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\
    384384                                        ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\
     385                                        ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\
    385386                                        ./modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp\
    386387                                        ./modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp\
    … …  
    952953                                        ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\
    953954                                        ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\
     955                                        ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\
    954956                                        ./modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp\
    955957                                        ./modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp\

  • issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp

    r6201 r6213  
    3737        }
    3838        if(iomodel->control_analysis){
     39                IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type");
    3940                IoModelFetchData(&iomodel->thickness_obs,NULL,NULL,iomodel_handle,"thickness_obs");
    4041                IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
    … …  
    6667        xfree((void**)&iomodel->thickness_obs);
    6768        xfree((void**)&iomodel->weights);
     69        xfree((void**)&iomodel->control_type);
    6870}

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

    r6201 r6213  
    2727
    2828                /*What control type?*/
    29                 if (iomodel->control_type!=DragCoefficientEnum &&
    30                                         iomodel->control_type!=RheologyBbarEnum &&
    31                                         iomodel->control_type!=DhDtEnum &&
    32                                         iomodel->control_type!=VxEnum &&
    33                                         iomodel->control_type!=VyEnum
    34                                         ) ISSMERROR("control_type %s not supported yet!",EnumToString(iomodel->control_type));
    35 
    36                 parameters->AddObject(new IntParam(ControlTypeEnum,iomodel->control_type));
     29                IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type");
     30                parameters->AddObject(new IntVecParam(ControlTypeEnum,iomodel->control_type,iomodel->num_control_type));
     31                xfree((void**)&iomodel->control_type);
    3732
    3833                /*What solution type?*/

  • issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r5579 r6213  
    1818void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){
    1919
    20         bool continuous=true;
    21         Elements* elements=NULL;
     20        bool       continuous = true;
     21        Elements  *elements   = NULL;
     22        Materials *materials  = NULL;
    2223                       
    2324        /*Create elements, vertices and materials, independent of analysis_type: */
    2425        CreateElementsVerticesAndMaterials(pelements, pvertices, pmaterials, iomodel,iomodel_handle,nummodels);
    2526
    26         /*Recover elements, for future update: */
     27        /*Recover elements and materials, for future update: */
    2728        elements=*pelements;
     29        materials=*pmaterials;
    2830
    2931        /*Now, branch onto analysis dependent model generation: */
    … …  
    102104        }
    103105
     106        /*Update Elements and Materials For Control methods*/
     107        UpdateElementsAndMaterialsControl(elements,materials,iomodel,iomodel_handle);
     108
    104109        /*Generate objects that are not dependent on any analysis_type: */
    105110        CreateParameters(pparameters,iomodel,iomodel_handle,solution_type,analysis_type,analysis_counter);

  • issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r6201 r6213  
    4040        IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
    4141        IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
     42        IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type");
    4243       
    4344        /*Create elements and materials: */
    … …  
    6162        xfree((void**)&iomodel->rheology_B);
    6263        xfree((void**)&iomodel->rheology_n);
     64        xfree((void**)&iomodel->control_type);
    6365
    6466        /*Add new constrant material property tgo materials, at the end: */
    … …  
    9597        *pvertices=vertices;
    9698        *pmaterials=materials;
    97 
    9899}

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r6201 r6213  
    4747        }
    4848        if(iomodel->control_analysis){
     49                IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type");
    4950                IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs");
    5051                IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
    … …  
    8788        xfree((void**)&iomodel->vy_obs);
    8889        xfree((void**)&iomodel->weights);
     90        xfree((void**)&iomodel->control_type);
    8991}

  • issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r5579 r6213  
    2121void  CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int solution_type,int analysis_type);
    2222void  CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int solution_type,int analysis_type);
     23void  UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    2324
    2425/*Creation of fem datasets: specialised drivers: */

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

    r6200 r6213  
    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;
    5859                virtual void   UpdateGeometry(void)=0;
    5960                virtual void   InputToResult(int enum_type,int step,double time)=0;

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

    r6212 r6213  
    641641        }
    642642        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));
    643644
    644645        this->GetDofList1(&doflist1[0]);
    … …  
    661662        }
    662663        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));
    663665
    664666        this->GetDofList1(&doflist1[0]);
    … …  
    959961
    960962        /*Intermediary*/
    961         int    control_type;
     963        int    num_controls;
     964        int*   control_type=NULL;
    962965        Input* input=NULL;
    963966        double cm_min,cm_max;
    … …  
    966969        this->parameters->FindParam(&cm_min,CmMinEnum);
    967970        this->parameters->FindParam(&cm_max,CmMaxEnum);
    968         this->parameters->FindParam(&control_type,ControlTypeEnum);
    969 
    970         if(control_type==RheologyBbarEnum){
    971                 if (!IsOnBed()) return;
    972                 input=(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input);
    973         }
    974         else{
    975                 input=(Input*)this->inputs->GetInput(control_type);   ISSMASSERT(input);
    976         }
    977 
    978         if (input->Enum()!=ControlInputEnum) ISSMERROR("input %s is not a ControlInput",EnumToString(control_type));
    979 
    980         ((ControlInput*)input)->UpdateValue(scalar);
    981         input->Constrain(cm_min,cm_max);
    982         if (save_parameter) ((ControlInput*)input)->SaveValue();
    983 
    984         if(control_type==RheologyBbarEnum){
    985                 this->InputExtrude(RheologyBEnum,MaterialsEnum);
    986         }
     971        this->parameters->FindParam(&num_controls,NumControlsEnum);
     972        this->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
     973
     974        for(int i=0;i<num_controls;i++){
     975
     976                if(control_type[i]==RheologyBbarEnum){
     977                        if (!IsOnBed()) return;
     978                        input=(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input);
     979                }
     980                else{
     981                        input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input);
     982                }
     983
     984                if (input->Enum()!=ControlInputEnum) ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i]));
     985
     986                ((ControlInput*)input)->UpdateValue(scalar);
     987                input->Constrain(cm_min,cm_max);
     988                if (save_parameter) ((ControlInput*)input)->SaveValue();
     989
     990                if(control_type[i]==RheologyBbarEnum){
     991                        this->InputExtrude(RheologyBEnum,MaterialsEnum);
     992                }
     993        }
     994
     995        /*Clean up and return*/
     996        xfree((void**)&control_type);
    987997}
    988998/*}}}*/
    … …  
    16881698}
    16891699/*}}}*/
    1690 /*FUNCTION Penta::Update {{{1*/
     1700/*FUNCTION Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type) {{{1*/
    16911701void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type){
    16921702
    16931703        /*Intermediaries*/
    1694         IssmInt i;
     1704        IssmInt i,j;
    16951705        int     penta_type;
    16961706        int     penta_node_ids[6];
    … …  
    17291739        this->SetHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
    17301740
    1731         //add as many inputs per element as requested:
    1732         if (iomodel->thickness) {
    1733                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1];
    1734                 this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));
    1735         }
    1736         if (iomodel->surface) {
    1737                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1];
    1738                 this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));
    1739         }
    1740         if (iomodel->bed) {
    1741                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1];
    1742                 this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs));
    1743         }
    1744         if (iomodel->drag_coefficient) {
    1745                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1];
    1746                 this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs));
    1747 
    1748                 if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
    1749                 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
    1750                 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
    1751 
    1752         }
    1753         if (iomodel->melting_rate) {
    1754                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts;
    1755                 this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));
    1756         }
    1757         if (iomodel->accumulation_rate) {
    1758                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts;
    1759                 this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));
    1760         }
    1761         if (iomodel->geothermalflux) {
    1762                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1];
    1763                 this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs));
    1764         }       
    1765         if (iomodel->pressure) {
    1766                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1];
    1767                 this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));
    1768         }
    1769         if (iomodel->temperature) {
    1770                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1];
    1771                 this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs));
    1772         }
    1773         if (iomodel->dhdt) {
    1774                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts;
    1775                 this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs));
    1776         }
    1777         /*vx,vy and vz: */
    1778         if (iomodel->vx) {
    1779                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts;
    1780                 this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
    1781                 this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));
    1782                 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs));
    1783         }
    1784         if (iomodel->vy) {
    1785                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts;
    1786                 this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
    1787                 this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));
    1788                 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs));
    1789         }
    1790         if (iomodel->vz) {
    1791                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts;
    1792                 this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
    1793                 this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
    1794                 if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs));
    1795         }
    1796         if (iomodel->vx_obs) {
    1797                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
    1798                 this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs));
    1799         }
    1800         if (iomodel->vy_obs) {
    1801                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
    1802                 this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs));
    1803         }
    1804         if (iomodel->vz_obs) {
    1805                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
    1806                 this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs));
    1807         }
    1808         if (iomodel->weights) {
    1809                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];
    1810                 this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs));
    1811         }
    1812         if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
    1813         if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
    1814         if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
    1815         if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
     1741        /*Fill with IoModel*/
     1742        this->Update(index,iomodel);
    18161743
    18171744        /*Defaults if not provided in iomodel*/
    … …  
    18881815        }
    18891816
     1817}
     1818/*}}}*/
     1819/*FUNCTION Penta::Update(int index,IoModel* iomodel) {{{1*/
     1820void 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
    18901923        /*Control Inputs*/
    1891         if (iomodel->control_analysis){
    1892                 switch(iomodel->control_type){
    1893                         case DhDtEnum:
    1894                                 if (iomodel->dhdt){
    1895                                         for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts;
    1896                                         this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs));
    1897                                 }
    1898                                 break;
    1899                         case VxEnum:
    1900                                 if (iomodel->vx){
    1901                                         for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts;
    1902                                         this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs));
    1903                                 }
    1904                                 break;
    1905                         case VyEnum:
    1906                                 if (iomodel->vy){
    1907                                         for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts;
    1908                                         this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs));
    1909                                 }
    1910                                 break;
    1911                         case DragCoefficientEnum:
    1912                                 if (iomodel->drag_coefficient){
    1913                                         for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1];
    1914                                         this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs));
    1915                                 }
    1916                                 break;
    1917                         case RheologyBbarEnum:
    1918                                 /*Matice will take care of it*/ break;
    1919                         default:
    1920                                 ISSMERROR("Control %s not implemented yet",EnumToString(iomodel->control_type));
     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                        }
    19211956                }
    19221957        }

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

    r6212 r6213  
    119119                double SurfaceArea(void);
    120120                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
     121                void   Update(int index, IoModel* iomodel);
    121122                void   UpdateGeometry(void);
    122123                double TimeAdapt();

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

    r6212 r6213  
    641641        }
    642642        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));
    643644
    644645        this->GetDofList1(&doflist1[0]);
    … …  
    661662        }
    662663        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));
    663665
    664666        this->GetDofList1(&doflist1[0]);
    … …  
    677679        /* Intermediaries */
    678680        int        i,j,ig;
    679         int        control_type;
     681        int        num_controls;
     682        int       *control_type=NULL;
    680683        double     Jelem = 0;
    681684        double     cm_noisedmp;
    … …  
    686689
    687690        /*retrieve parameters and inputs*/
    688         this->parameters->FindParam(&control_type,ControlTypeEnum);
     691        this->parameters->FindParam(&num_controls,NumControlsEnum);
     692        this->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    689693        this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    690694
    … …  
    694698        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    695699
    696         /* Start looping on the number of gaussian points: */
    697         gauss=new GaussTria(2);
    698         for (ig=gauss->begin();ig<gauss->end();ig++){
    699 
    700                 gauss->GaussPoint(ig);
    701 
    702                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    703 
    704                 /*Add Tikhonov regularization term to misfit:
    705                  *
    706                  * WARNING: for now, the regularization is only taken into account by the gradient
    707                  * the misfit reflects the response only!
    708                  *
    709                  * */
    710                 if (control_type==DragCoefficientEnum){
    711                         Input* drag_input=inputs->GetInput(DragCoefficientEnum);      ISSMASSERT(drag_input);
    712                         drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    713                         //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
     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                        }
    714735                }
    715                 else if (control_type==RheologyBbarEnum){
    716                         //nothing
    717                 }
    718                 else if (control_type==DhDtEnum){
    719                         //nothing
    720                 }
    721                 else if (control_type==VxEnum){
    722                         //nothing
    723                 }
    724                 else if (control_type==VyEnum){
    725                         //nothing
    726                 }
    727                 else{
    728                         ISSMERROR("unsupported control type: %s",EnumToString(control_type));
    729                 }
    730         }
    731 
    732         /*clean-up and return: */
    733         delete gauss;
     736
     737                delete gauss;
     738        }
     739
     740        /*Clean up and return*/
     741        xfree((void**)&control_type);
    734742        return Jelem;
    735743}
    … …  
    11871195
    11881196        /*Intermediary*/
    1189         int    control_type;
     1197        int    num_controls;
     1198        int*   control_type=NULL;
    11901199        Input* input=NULL;
    11911200        double cm_min,cm_max;
    … …  
    11941203        this->parameters->FindParam(&cm_min,CmMinEnum);
    11951204        this->parameters->FindParam(&cm_max,CmMaxEnum);
    1196         this->parameters->FindParam(&control_type,ControlTypeEnum);
    1197 
    1198         if(control_type==RheologyBbarEnum){
    1199                 input=(Input*)matice->inputs->GetInput(control_type); ISSMASSERT(input);
    1200         }
    1201         else{
    1202                 input=(Input*)this->inputs->GetInput(control_type);   ISSMASSERT(input);
    1203         }
    1204 
    1205         if (input->Enum()!=ControlInputEnum){
    1206                 ISSMERROR("input %s is not a ControlInput",EnumToString(control_type));
    1207         }
    1208 
    1209         ((ControlInput*)input)->UpdateValue(scalar);
    1210         input->Constrain(cm_min,cm_max);
    1211         if (save_parameter) ((ControlInput*)input)->SaveValue();
     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);
    12121229}
    12131230/*}}}*/
    … …  
    21172134}
    21182135/*}}}*/
    2119 /*FUNCTION Tria::Update{{{1*/
     2136/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/
    21202137void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index
    21212138
    21222139        /*Intermediaries*/
    2123         int    i;
     2140        int    i,j;
    21242141        int    tria_node_ids[3];
    21252142        int    tria_vertex_ids[3];
    … …  
    21642181        /*hooks: */
    21652182        this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
    2166        
    2167         /*add as many inputs per element as requested:*/
    2168         if (iomodel->thickness) {
    2169                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1];
    2170                 this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs));
    2171         }
    2172         if (iomodel->surface) {
    2173                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1];
    2174                 this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs));
    2175         }
    2176         if (iomodel->bed) {
    2177                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1];
    2178                 this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs));
    2179         }
    2180         if (iomodel->drag_coefficient) {
    2181                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1];
    2182                 this->inputs->AddInput(new TriaVertexInput(DragCoefficientEnum,nodeinputs));
    2183 
    2184                 if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
    2185                 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
    2186                 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
    2187         }
    2188         if (iomodel->thickness_obs) {
    2189                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness_obs[tria_vertex_ids[i]-1];
    2190                 this->inputs->AddInput(new TriaVertexInput(ThicknessObsEnum,nodeinputs));
    2191         }
    2192         if (iomodel->melting_rate) {
    2193                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts;
    2194                 this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs));
    2195         }
    2196         if (iomodel->accumulation_rate) {
    2197                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
    2198                 this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs));
    2199         }
    2200         if (iomodel->geothermalflux) {
    2201                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1];
    2202                 this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs));
    2203         }
    2204         if (iomodel->dhdt){
    2205                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]/iomodel->yts;
    2206                 this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs));
    2207         }
    2208         if (iomodel->pressure){
    2209                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1];
    2210                 this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs));
    2211         }
    2212         if (iomodel->temperature) {
    2213                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1];
    2214                 this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs));
    2215         }
    2216         /*vx,vy and vz: */
    2217         if (iomodel->vx) {
    2218                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts;
    2219                 this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
    2220                 this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs));
    2221                 if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs));
    2222         }
    2223         if (iomodel->vy) {
    2224                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts;
    2225                 this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
    2226                 this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs));
    2227                 if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs));
    2228         }
    2229         if (iomodel->vz) {
    2230                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts;
    2231                 this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs));
    2232                 this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs));
    2233                 if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs));
    2234         }
    2235         if (iomodel->vx_obs) {
    2236                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts;
    2237                 this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs));
    2238         }
    2239         if (iomodel->vy_obs) {
    2240                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts;
    2241                 this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs));
    2242         }
    2243         if (iomodel->vz_obs) {
    2244                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts;
    2245                 this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs));
    2246         }
    2247         if (iomodel->weights) {
    2248                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1];
    2249                 this->inputs->AddInput(new TriaVertexInput(WeightsEnum,nodeinputs));
    2250         }
    2251         if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
    2252         if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
    2253         if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
    2254         if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
     2183
     2184        /*Fill with IoModel*/
     2185        this->Update(index,iomodel);
    22552186
    22562187        /*Defaults if not provided in iomodel*/
    … …  
    22962227        }
    22972228
    2298         /*Control Inputs*/
    2299         if (iomodel->control_analysis){
    2300                 switch(iomodel->control_type){
    2301                         case DhDtEnum:
    2302                                 if (iomodel->dhdt){
    2303                                         for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]/iomodel->yts;
    2304                                         this->inputs->AddInput(new ControlInput(DhDtEnum,TriaVertexInputEnum,nodeinputs));
    2305                                 }
    2306                                 break;
    2307                         case VxEnum:
    2308                                 if (iomodel->vx){
    2309                                         for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts;
    2310                                         this->inputs->AddInput(new ControlInput(VxEnum,TriaVertexInputEnum,nodeinputs));
    2311                                 }
    2312                                 break;
    2313                         case VyEnum:
    2314                                 if (iomodel->vy){
    2315                                         for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts;
    2316                                         this->inputs->AddInput(new ControlInput(VyEnum,TriaVertexInputEnum,nodeinputs));
    2317                                 }
    2318                                 break;
    2319                         case DragCoefficientEnum:
    2320                                 if (iomodel->drag_coefficient){
    2321                                         for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1];
    2322                                         this->inputs->AddInput(new ControlInput(DragCoefficientEnum,TriaVertexInputEnum,nodeinputs));
    2323                                 }
    2324                                 break;
    2325                         case RheologyBbarEnum:
    2326                                 /*Matice will take care of it*/ break;
    2327                         default:
    2328                                 ISSMERROR("Control %s not implemented yet",EnumToString(iomodel->control_type));
    2329                 }
    2330 
    2331         }
    2332 
    23332229        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    23342230        this->parameters=NULL;
    2335 
     2231}
     2232/*}}}*/
     2233/*FUNCTION Tria::Update(int index, IoModel* iomodel){{{1*/
     2234void Tria::Update(int index, IoModel* iomodel){ //i is the element index
     2235
     2236        /*Intermediaries*/
     2237        int    i,j;
     2238        int    tria_vertex_ids[3];
     2239        double nodeinputs[3];
     2240
     2241        /*Checks if debuging*/
     2242        /*{{{2*/
     2243        ISSMASSERT(iomodel->elements);
     2244        /*}}}*/
     2245
     2246        /*Recover vertices ids needed to initialize inputs*/
     2247        for(i=0;i<3;i++){
     2248                tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab
     2249        }
     2250
     2251        /*add as many inputs per element as requested:*/
     2252        if (iomodel->thickness) {
     2253                for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1];
     2254                this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs));
     2255        }
     2256        if (iomodel->surface) {
     2257                for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1];
     2258                this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs));
     2259        }
     2260        if (iomodel->bed) {
     2261                for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1];
     2262                this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs));
     2263        }
     2264        if (iomodel->drag_coefficient) {
     2265                for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1];
     2266                this->inputs->AddInput(new TriaVertexInput(DragCoefficientEnum,nodeinputs));
     2267
     2268                if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
     2269                if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
     2270                this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
     2271        }
     2272        if (iomodel->thickness_obs) {
     2273                for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness_obs[tria_vertex_ids[i]-1];
     2274                this->inputs->AddInput(new TriaVertexInput(ThicknessObsEnum,nodeinputs));
     2275        }
     2276        if (iomodel->melting_rate) {
     2277                for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts;
     2278                this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs));
     2279        }
     2280        if (iomodel->accumulation_rate) {
     2281                for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
     2282                this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs));
     2283        }
     2284        if (iomodel->geothermalflux) {
     2285                for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1];
     2286                this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs));
     2287        }
     2288        if (iomodel->dhdt){
     2289                for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]/iomodel->yts;
     2290                this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs));
     2291        }
     2292        if (iomodel->pressure){
     2293                for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1];
     2294                this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs));
     2295        }
     2296        if (iomodel->temperature) {
     2297                for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1];
     2298                this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs));
     2299        }
     2300        /*vx,vy and vz: */
     2301        if (iomodel->vx) {
     2302                for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts;
     2303                this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
     2304                this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs));
     2305                if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs));
     2306        }
     2307        if (iomodel->vy) {
     2308                for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts;
     2309                this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
     2310                this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs));
     2311                if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs));
     2312        }
     2313        if (iomodel->vz) {
     2314                for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts;
     2315                this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs));
     2316                this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs));
     2317                if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs));
     2318        }
     2319        if (iomodel->vx_obs) {
     2320                for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts;
     2321                this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs));
     2322        }
     2323        if (iomodel->vy_obs) {
     2324                for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts;
     2325                this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs));
     2326        }
     2327        if (iomodel->vz_obs) {
     2328                for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts;
     2329                this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs));
     2330        }
     2331        if (iomodel->weights) {
     2332                for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1];
     2333                this->inputs->AddInput(new TriaVertexInput(WeightsEnum,nodeinputs));
     2334        }
     2335        if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
     2336        if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
     2337        if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
     2338        if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
     2339
     2340        /*Control Inputs*/
     2341        if (iomodel->control_analysis && iomodel->control_type){
     2342                for(i=0;i<iomodel->num_control_type;i++){
     2343                        switch((int)iomodel->control_type[i]){
     2344                                case DhDtEnum:
     2345                                        if (iomodel->dhdt){
     2346                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->dhdt[tria_vertex_ids[j]-1]/iomodel->yts;
     2347                                                this->inputs->AddInput(new ControlInput(DhDtEnum,TriaVertexInputEnum,nodeinputs));
     2348                                        }
     2349                                        break;
     2350                                case VxEnum:
     2351                                        if (iomodel->vx){
     2352                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->vx[tria_vertex_ids[j]-1]/iomodel->yts;
     2353                                                this->inputs->AddInput(new ControlInput(VxEnum,TriaVertexInputEnum,nodeinputs));
     2354                                        }
     2355                                        break;
     2356                                case VyEnum:
     2357                                        if (iomodel->vy){
     2358                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->vy[tria_vertex_ids[j]-1]/iomodel->yts;
     2359                                                this->inputs->AddInput(new ControlInput(VyEnum,TriaVertexInputEnum,nodeinputs));
     2360                                        }
     2361                                        break;
     2362                                case DragCoefficientEnum:
     2363                                        if (iomodel->drag_coefficient){
     2364                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->drag_coefficient[tria_vertex_ids[j]-1];
     2365                                                this->inputs->AddInput(new ControlInput(DragCoefficientEnum,TriaVertexInputEnum,nodeinputs));
     2366                                        }
     2367                                        break;
     2368                                case RheologyBbarEnum:
     2369                                        /*Matice will take care of it*/ break;
     2370                                default:
     2371                                        ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i]));
     2372                        }
     2373                }
     2374        }
    23362375}
    23372376/*}}}*/

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

    r6200 r6213  
    122122                double SurfaceArea(void);
    123123                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
     124                void   Update(int index, IoModel* iomodel);
    124125                void   UpdateGeometry(void);
    125126                double TimeAdapt();

  • issm/trunk/src/c/objects/IoModel.cpp

    r6201 r6213  
    8484        xfree((void**)&this->rheology_B);
    8585        xfree((void**)&this->rheology_n);
     86        xfree((void**)&this->control_type);
    8687        xfree((void**)&this->cm_responses);
    8788        xfree((void**)&this->weights);
    … …  
    152153        /*Get control parameters: */
    153154        IoModelFetchData(&this->num_control_type,iomodel_handle,"num_control_type");
    154         IoModelFetchData(&this->control_type,iomodel_handle,"control_type");
    155155
    156156        /*!Get solution parameters: */
    … …  
    295295
    296296        /*!solution parameters: */
     297        this->control_type=NULL;
    297298        this->cm_responses=NULL;
    298299        this->weights=NULL;

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

    r6201 r6213  
    122122                /*control methods: */
    123123                int      num_control_type;
    124                 int      control_type;
     124                double*  control_type;
    125125
    126126                /*solution parameters: */

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

    r5320 r6213  
    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;
    2425
    2526};

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

    r6200 r6213  
    5757
    5858                /*Control Inputs*/
    59                 if (iomodel->control_analysis){
    60                         switch(iomodel->control_type){
    61                                 case RheologyBbarEnum:
    62                                         if (iomodel->rheology_B){
    63                                                 for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+i]-1)];
    64                                                 this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs));
    65                                         }
    66                                         break;
    67                                 default:
    68                                         /*Nothing*/;
     59                if (iomodel->control_analysis && iomodel->control_type){
     60                        for(i=0;i<iomodel->num_control_type;i++){
     61                                switch((int)iomodel->control_type[i]){
     62                                        case RheologyBbarEnum:
     63                                                if (iomodel->rheology_B){
     64                                                        for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+i]-1)];
     65                                                        this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs));
     66                                                }
     67                                                break;
     68                                }
    6969                        }
    7070                }
    … …  
    9191
    9292                /*Control Inputs*/
    93                 if (iomodel->control_analysis){
    94                         switch(iomodel->control_type){
    95                                 case RheologyBbarEnum:
    96                                         if (iomodel->rheology_B){
    97                                                 for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+i]-1)];
    98                                                 this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs));
    99                                         }
    100                                         break;
     93                if (iomodel->control_analysis && iomodel->control_type){
     94                        for(i=0;i<iomodel->num_control_type;i++){
     95                                switch((int)iomodel->control_type[i]){
     96                                        case RheologyBbarEnum:
     97                                                if (iomodel->rheology_B){
     98                                                        for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+i]-1)];
     99                                                        this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs));
     100                                                }
     101                                                break;
     102                                }
    101103                        }
    102104                }
    … …  
    246248
    247249        return matice;
     250}
     251/*}}}*/
     252/*FUNCTION Matice::Update{{{1*/
     253void 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;
    248305}
    249306/*}}}*/

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

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

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

    r5320 r6213  
    5656                void   InputUpdateFromConstant(bool constant, int name);
    5757                void   InputUpdateFromSolution(double* solution);
     58                void  Update(int index, IoModel* iomodel){ISSMERROR("not implemented yet");};
    5859                /*}}}*/
    5960                /*Material virtual functions resolution: {{{1*/

  • issm/trunk/src/c/objects/Params/DoubleParam.cpp

    r5103 r6213  
    129129}
    130130/*}}}*/
    131 /*FUNCTION DoubleParam::GetParameterValue{{{1*/
     131/*FUNCTION DoubleParam::GetParameterValue(int* pinteger){{{1*/
    132132void DoubleParam::GetParameterValue(int* pinteger){
    133133#ifdef _SERIAL_
    … …  
    138138}
    139139/*}}}*/
    140 /*FUNCTION DoubleParam::GetParameterValue{{{1*/
     140/*FUNCTION DoubleParam::GetParameterValue(bool* pbool){{{1*/
    141141void DoubleParam::GetParameterValue(bool* pbool){
    142142#ifdef _SERIAL_
    … …  
    148148#else
    149149        ISSMERROR("Double param of enum %i (%s) cannot return an bool",enum_type,EnumToString(enum_type));
     150#endif
     151}
     152/*}}}*/
     153/*FUNCTION DoubleParam::GetParameterValue(int** pintarray,int* pM){{{1*/
     154void DoubleParam::GetParameterValue(int** pintarray,int* pM){
     155#ifdef _SERIAL_
     156        int* output=NULL;
     157
     158        output=(int*)xmalloc(1*sizeof(int));
     159        *output=(int)value;
     160
     161        /*Assign output pointers:*/
     162        if(pM) *pM=1;
     163        *pintarray=output;
     164#else
     165        ISSMERROR("Double param of enum %i (%s) cannot return an array of integers",enum_type,EnumToString(enum_type));
    150166#endif
    151167}

  • issm/trunk/src/c/objects/Params/DoubleParam.h

    r6165 r6213  
    5252                void  GetParameterValue(bool* pbool);
    5353                void  GetParameterValue(int* pinteger);
    54                 void  GetParameterValue(int** pintarray,int* pM){ISSMERROR("Double param of enum %i (%s) cannot return an array of integers",enum_type,EnumToString(enum_type));}
     54                void  GetParameterValue(int** pintarray,int* pM);
    5555                void  GetParameterValue(double* pdouble){*pdouble=value;}
    5656                void  GetParameterValue(char** pstring){ISSMERROR("Double param of enum %i (%s) cannot return a string",enum_type,EnumToString(enum_type));}

  • issm/trunk/src/c/objects/Params/IntVecParam.cpp

    r6165 r6213  
    3636}
    3737/*}}}*/
     38/*FUNCTION IntVecParam::IntVecParam(int enum_type,double* values,int M){{{1*/
     39IntVecParam::IntVecParam(int in_enum_type,double* in_values, int in_M){
     40
     41        enum_type=in_enum_type;
     42        M=in_M;
     43
     44        values=(int*)xmalloc(M*sizeof(int));
     45        for(int i=0;i<in_M;i++) values[i]=(int)in_values[i];
     46}
     47/*}}}*/
    3848/*FUNCTION IntVecParam::~IntVecParam(){{{1*/
    3949IntVecParam::~IntVecParam(){
    … …  
    101111int   IntVecParam::MarshallSize(){
    102112       
    103         return sizeof(M)
     113        return sizeof(M)+
    104114                +M*sizeof(int)
    105115                +sizeof(enum_type)+

  • issm/trunk/src/c/objects/Params/IntVecParam.h

    r6169 r6213  
    3535                IntVecParam();
    3636                IntVecParam(int enum_type,int* values,int M);
     37                IntVecParam(int enum_type,double* values,int M);
    3738                ~IntVecParam();
    3839                /*}}}*/
    … …  
    6970                void  SetValue(char* string){ISSMERROR("IntVec param of enum %i (%s) cannot hold a string",enum_type,EnumToString(enum_type));}
    7071                void  SetValue(char** stringarray,int M){ISSMERROR("IntVec param of enum %i (%s) cannot hold a string array",enum_type,EnumToString(enum_type));}
    71                 void  SetValue(double* doublearray,int M);
     72                void  SetValue(double* doublearray,int M){ISSMERROR("IntVec param of enum %i (%s) cannot hold a double mat array",enum_type,EnumToString(enum_type));}
    7273                void  SetValue(double* pdoublearray,int M,int N){ISSMERROR("IntVec param of enum %i (%s) cannot hold a double mat array",enum_type,EnumToString(enum_type));}
    7374                void  SetValue(Vec vec){ISSMERROR("IntVec param of enum %i (%s) cannot hold a Vec",enum_type,EnumToString(enum_type));}

  • issm/trunk/src/c/solutions/control_core.cpp

    r6200 r6213  
    1818        /*parameters: */
    1919        int     verbose=0;
    20         int     control_type;
     20        int     num_controls;
    2121        int     nsteps;
    2222        double  eps_cm;
    … …  
    2828        bool    qmu_analysis=false;
    2929
     30        int*    control_type = NULL;
    3031        double* responses=NULL;
    3132        double* optscal=NULL;
    … …  
    4849        /*Recover parameters used throughout the solution:{{{1*/
    4950        femmodel->parameters->FindParam(&nsteps,NStepsEnum);
    50         femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
     51        femmodel->parameters->FindParam(&num_controls,NumControlsEnum);
     52        femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    5153        femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    5254        femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
    … …  
    122124        /*some results not computed by steadystate_core or diagnostic_core: */
    123125        if(!qmu_analysis){ //do not save this if we are running the control core from a qmu run!
    124                 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
     126                for(i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);
    125127                femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
    126                 femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0));
     128                //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0));
    127129        }
    128130
    129131        cleanup_and_return:
    130132        /*Free ressources: */
     133        xfree((void**)&control_type);
    131134        xfree((void**)&responses);
    132135        xfree((void**)&optscal);

  • issm/trunk/src/c/solutions/controlrestart.cpp

    r5530 r6213  
    99void controlrestart(FemModel* femmodel,double* J){
    1010
    11         int      control_type;
     11        int      num_controls;
     12        int*     control_type = NULL;
    1213        int      nsteps;
    1314        bool     qmu_analysis=true;
    1415
    1516        /*retrieve output file name: */
    16         femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
     17        femmodel->parameters->FindParam(&num_controls,NumControlsEnum);
     18        femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    1719        femmodel->parameters->FindParam(&nsteps,NStepsEnum);
    1820        femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
    … …  
    2224        if(!qmu_analysis){
    2325                /*we essentially want J and the parameter: */
    24                 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
     26                for(int i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);
    2527                femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
    26                 femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0));
     28                //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0));
    2729
    2830                /*write to disk: */
    … …  
    3032        }
    3133
     34        /*Clean up and return*/
     35        xfree((void**)&control_type);
    3236}

  • issm/trunk/src/c/solutions/gradient_core.cpp

    r6200 r6213  
    1717        /*parameters: */
    1818        bool control_steady;
    19         int  control_type;
     19        int  num_controls;
     20        int* control_type=NULL;
    2021        int  verbose;
    2122
    … …  
    2728        /*retrieve parameters:*/
    2829        femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
    29         femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
     30        femmodel->parameters->FindParam(&num_controls,NumControlsEnum);
     31        femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    3032        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    3133
    32         if(verbose)_printf_("%s\n","      compute gradient:");
    33         Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type);
     34        for (int i=0;i<num_controls;i++){
    3435
    35         if(control_steady)diagnostic_core(femmodel);
    36        
    37         if (step>0 && search_scalar==0){
    38                 _printf_("%s","      orthogonalization...\n");
    39                 if(verbose)_printf_("%s\n","      retrieve old gradient:");
    40                 ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type);
    41                 Orthx(&new_gradient,gradient,old_gradient);
    42                 VecFree(&old_gradient);
     36                if(verbose)_printf_("      compute gradient of J with respect to %s\n",EnumToString(control_type[i]));
     37                Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type[i]);
     38
     39                if(control_steady)diagnostic_core(femmodel);
     40
     41                if (step>0 && search_scalar==0){
     42                        _printf_("%s","      orthogonalization...\n");
     43                        ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]);
     44                        Orthx(&new_gradient,gradient,old_gradient);
     45                        VecFree(&old_gradient);
     46                }
     47                else{
     48                        _printf_("%s","      normalizing directions...\n");
     49                        Orthx(&new_gradient,gradient,NULL);
     50                }
     51                VecFree(&gradient);
     52
     53                /*plug back into inputs: */
     54                ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type[i],new_gradient);
     55
     56                /*Free ressources and return:*/
     57                VecFree(&new_gradient);
    4358        }
    44         else{
    45                 _printf_("%s","      normalizing directions...\n");
    46                 Orthx(&new_gradient,gradient,NULL);
    47         }
    48         VecFree(&gradient);
    4959
    50         /*plug back into inputs: */
    51         ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type,new_gradient);
    52 
    53         /*Free ressources and return:*/
    54         VecFree(&new_gradient);
     60        /*Clean up and return*/
     61        xfree((void**)&control_type);
    5562}

  • issm/trunk/src/c/solutions/objectivefunctionC.cpp

    r6200 r6213  
    3131        int        n;
    3232        double    *optscal        = NULL;
    33         double    *responses      =NULL;
    34         int        control_type;
     33        double    *responses      = NULL;
    3534        int        solution_type;
    3635        int        analysis_type;
    … …  
    4746        femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    4847        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    49         femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
    5048        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    5149        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
Note: See TracChangeset for help on using the changeset viewer.