Changeset 25406


Ignore:
Timestamp:
08/16/20 15:34:50 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on transient calibration. Not working yet

Location:
issm/trunk-jpl/src/c
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp

    r25379 r25406  
    101101IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/
    102102
    103          /*diverse: */
    104          IssmDouble time;
    105 
    106          /*recover time parameters: */
    107          femmodel->parameters->FindParam(&time,TimeEnum);
    108          if(time < last_time) timepassedflag = false;
    109          last_time = time;
    110 
    111                  int i;
    112                  IssmDouble J=0.;
    113                  IssmDouble J_sum=0.;
    114 
    115          if(datatime<=time && !timepassedflag){
    116                  for(i=0;i<femmodel->elements->Size();i++){
    117                          Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
    118                          J+=this->Cfsurfacelogvel_Calculation(element,definitionenum);
    119                  }
    120 
    121                  ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    122                  ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    123                  J=J_sum;
    124 
    125                  timepassedflag = true;
    126                  return J;
    127                 }
    128          else return J;
    129  }
    130         /*}}}*/
     103
     104        /*recover time parameters: */
     105        IssmDouble time;
     106        femmodel->parameters->FindParam(&time,TimeEnum);
     107        if(time < last_time) this->timepassedflag = false; //FIXME: do we really need this? looks like an old piece of code before copy
     108        this->last_time = time;
     109
     110        IssmDouble J=0.;
     111        IssmDouble J_sum=0.;
     112        if(this->datatime<=time && !timepassedflag){
     113                for(int i=0;i<femmodel->elements->Size();i++){
     114                        Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
     115                        J+=this->Cfsurfacelogvel_Calculation(element,definitionenum);
     116                }
     117
     118                ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     119                ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     120                J=J_sum;
     121
     122                this->timepassedflag = true;
     123        }
     124        return J;
     125}/*}}}*/
    131126IssmDouble Cfsurfacelogvel::Cfsurfacelogvel_Calculation(Element* element, int definitionenum){/*{{{*/
    132127
    … …  
    165160        Input *vy_input = NULL;
    166161        if(numcomponents==2){
    167               vy_input = topelement->GetInput(VyEnum); _assert_(vy_input);
     162                vy_input = topelement->GetInput(VyEnum); _assert_(vy_input);
    168163        }
    169164
    … …  
    192187                 *        *                 [        vel + eps     ] 2
    193188                 *               * J = 4 \bar{v}^2 | log ( -----------  ) |
    194                  *                      *                 [       vel   + eps    ]
    195                  *                             *                            obs
    196                  *                                    */
     189                 *                          [       vel   + eps    ]
     190                 *                                      obs
     191                 */
    197192                if(numcomponents==1){
    198193                        velocity_mag    =fabs(vx)+epsvel;
    … …  
    205200
    206201                misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
     202                _printf_(velocity_mag<<" "<<obs_velocity_mag<<"\n");
     203                _error_("S");
    207204
    208205                /*Add to cost function*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r25396 r25406  
    2020#include "../Inputs/PentaInput.h"
    2121#include "../Inputs/DatasetInput.h"
     22#include "../Inputs/ControlInput.h"
    2223#include "../Inputs/ArrayInput.h"
    2324/*}}}*/
    … …  
    17471748        /*Intermediaries*/
    17481749        const int numvertices = this->GetNumberOfVertices();
     1750        IssmDouble  value,value_min,value_max;
    17491751
    17501752        int        *vertexids   = xNew<int>(numvertices);
    … …  
    17791781                this->AddControlInput(input_enum,inputs,iomodel,values,values_min,values_max,P0Enum,id);
    17801782        }
    1781 
     1783        else if(M==iomodel->numberofelements+1 && N>1){
     1784
     1785                /*create transient input: */
     1786                IssmDouble* times = xNew<IssmDouble>(N);
     1787                for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1788                inputs->SetTransientControlInput(input_enum,id,times,N);
     1789
     1790                /*Get input*/
     1791                ControlInput* control_input = inputs->GetControlInput(input_enum); _assert_(control_input);
     1792                TransientInput* transientinput_val = control_input->GetTransientInput("value");
     1793                TransientInput* transientinput_min = control_input->GetTransientInput("lowerbound");
     1794                TransientInput* transientinput_max = control_input->GetTransientInput("upperbound");
     1795                for(int t=0;t<N;t++){
     1796                        value     = vector[N*this->Sid()+t];
     1797                        value_min = scale*min_vector[N*this->Sid()+t];
     1798                        value_max = scale*max_vector[N*this->Sid()+t];
     1799                        switch(this->ObjectEnum()){
     1800                                case TriaEnum:
     1801                                        transientinput_val->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum);
     1802                                        transientinput_min->AddTriaTimeInput( t,1,&(this->lid),&value_min,P0Enum);
     1803                                        transientinput_max->AddTriaTimeInput( t,1,&(this->lid),&value_max,P0Enum);
     1804                                        break;
     1805                                case PentaEnum:
     1806                                        //transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum);
     1807                                        _error_("to be implemented");
     1808                                        break;
     1809                                default: _error_("Not implemented yet");
     1810                        }
     1811                }
     1812                xDelete<IssmDouble>(times);
     1813        }
    17821814        else if(M==iomodel->numberofvertices+1 && N>1){
    17831815                _error_("not supported tet");

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r25384 r25406  
    259259                virtual int        GetNumberOfNodes(int enum_type)=0;
    260260                virtual int        GetNumberOfVertices(void)=0;
    261                 virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset)=0;
    262                 virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data)=0;
     261                virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,int N,const char* data,int offset)=0;
    263262                virtual void       GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
    264263                virtual void       GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r25379 r25406  
    215215        }
    216216
     217        /*Create Control Input*/
     218        inputs->SetControlInput(input_enum,PentaInputEnum,interpolation_enum,id);
     219        ControlInput* control_input = inputs->GetControlInput(input_enum); _assert_(input_enum);
     220
    217221        /*Call inputs method*/
    218222        switch(interpolation_enum){
    219223                case P1Enum:
    220                         inputs->SetPentaControlInput(input_enum,PentaInputEnum,interpolation_enum,id,NUMVERTICES,vertexlids,values,values_min,values_max);
     224                        control_input->SetControl(interpolation_enum,NUMVERTICES,&vertexlids[0],values,values_min,values_max);
    221225                        break;
    222226                default:
    … …  
    18301834}
    18311835/*}}}*/
    1832 void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    1833 
    1834         _error_("NOT NEEDED ANYMORE");
     1836void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,int N,const char* data,int offset){/*{{{*/
     1837
    18351838        /*Get input*/
    18361839        if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
    … …  
    18781881                                        for(int t=0;t<transientinput->numtimesteps;t++) {
    18791882                                                IssmDouble time = transientinput->GetTimeByOffset(t);
    1880                                                 _error_("not implemented");
    1881                                                 //PentaInput* timeinput = xDynamicCast<PentaInput*>(transientinput->GetTimeInput(time));
    1882                                                 //if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
    1883                                                 //input->Serve(NUMVERTICES,&lidlist[0]);
    1884                                                 ///*Create list of indices and values for global vector*/
    1885                                                 //for(int i=0;i<NUMVERTICES;i++){
    1886                                                 //              idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
    1887                                                 //              values[N*i+t] = timeinput->values[i];
    1888                                                 //}
    1889                                         }
    1890                                         vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
    1891                                         xDelete<int>(M);
    1892                                         xDelete<int>(idlist);
    1893                                         xDelete<IssmDouble>(values);
    1894                                         break;
    1895                                 }
    1896                 default: _error_("input "<<EnumToStringx(input->ObjectEnum())<<" not supported yet");
    1897         }
    1898 }
    1899 /*}}}*/
    1900 void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset){/*{{{*/
    1901 
    1902         /*Get input*/
    1903         if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
    1904         ElementInput* input=this->inputs->GetControlInputData(control_enum,data);   _assert_(input);
    1905 
    1906         /*Lid list once for all*/
    1907         int lidlist[NUMVERTICES];
    1908         for(int i=0;i<NUMVERTICES;i++) lidlist[i] = vertices[i]->lid;
    1909 
    1910         /*Check what input we are dealing with*/
    1911         switch(input->ObjectEnum()){
    1912                 case PentaInputEnum:
    1913                           {
    1914                                 IssmDouble values[NUMVERTICES];
    1915                                 int        idlist[NUMVERTICES];
    1916 
    1917                                 PentaInput* triainput = xDynamicCast<PentaInput*>(input);
    1918 
    1919                                 /*Create list of indices and values for global vector*/
    1920                                 GradientIndexing(&idlist[0],control_index);
    1921 
    1922                                 if(triainput->GetInputInterpolationType()==P1Enum){
    1923                                         input->Serve(NUMVERTICES,&lidlist[0]);
    1924                                         for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
    1925                                         vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
    1926                                 }
    1927                                 else if(triainput->GetInputInterpolationType()==P0Enum){
    1928                                         input->Serve(1,&this->lid);
    1929                                         vector->SetValue(idlist[0],triainput->element_values[0],INS_VAL);
    1930                                 }
    1931                                 else{
    1932                                         _error_("not supported yet");
    1933                                 }
    1934                                 break;
    1935                           }
    1936 
    1937                 case TransientInputEnum:
    1938                                 {
    1939                                         TransientInput* transientinput = xDynamicCast<TransientInput*>(input);
    1940                                         int  N = transientinput->numtimesteps;
    1941                                         int* M = NULL;
    1942                                         parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    1943                                         int* idlist = xNew<int>(NUMVERTICES*N);
    1944                                         IssmDouble* values = xNew<IssmDouble>(NUMVERTICES*N);
    1945                                         for(int t=0;t<transientinput->numtimesteps;t++) {
    1946                                                 IssmDouble time = transientinput->GetTimeByOffset(t);
    1947                                                 _error_("not implemented");
     1883                                                _error_("not implemented, SEE TRIA!");
    19481884                                                //PentaInput* timeinput = xDynamicCast<PentaInput*>(transientinput->GetTimeInput(time));
    19491885                                                //if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r25379 r25406  
    9999                Penta*         GetLowerPenta(void);
    100100                Penta*         GetUpperPenta(void);
    101                 void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
    102                 void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset);
     101                void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,int N,const char* data,int offset);
    103102                int            GetVertexIndex(Vertex* vertex);
    104103                void           GetVerticesCoordinatesBase(IssmDouble** pxyz_list);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r25379 r25406  
    7373                int         GetNumberOfNodes(int enum_type){_error_("not implemented yet");};
    7474                int         GetNumberOfVertices(void);
    75                 void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset){_error_("not implemented yet");};
    76                 void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data){_error_("not implemented yet");};
     75                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,int N,const char* data,int offset){_error_("not implemented yet");};
    7776                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
    7877                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r25379 r25406  
    7777                int         GetNumberOfNodes(int enum_type){_error_("not implemented yet");};
    7878                int         GetNumberOfVertices(void);
    79                 void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset){_error_("not implemented yet");};
    80                 void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data){_error_("not implemented yet");};
     79                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,int N,const char* data,int offset){_error_("not implemented yet");};
    8180                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    8281                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r25379 r25406  
    220220        }
    221221
     222        /*Create Control Input*/
     223        inputs->SetControlInput(input_enum,TriaInputEnum,interpolation_enum,id);
     224        ControlInput* control_input = inputs->GetControlInput(input_enum); _assert_(input_enum);
     225
    222226        /*Call inputs method*/
    223227        switch(interpolation_enum){
    224228                case P0Enum:
    225                         inputs->SetTriaControlInput(input_enum,TriaInputEnum,interpolation_enum,id,1,&this->lid,values,values_min,values_max);
     229                        control_input->SetControl(interpolation_enum,1,&this->lid,values,values_min,values_max);
    226230                        break;
    227231                case P1Enum:
    228                         inputs->SetTriaControlInput(input_enum,TriaInputEnum,interpolation_enum,id,NUMVERTICES,vertexlids,values,values_min,values_max);
     232                        control_input->SetControl(interpolation_enum,NUMVERTICES,&vertexlids[0],values,values_min,values_max);
    229233                        break;
    230234                default:
    … …  
    11041108        IssmDouble  values[NUMVERTICES];
    11051109        int         lidlist[NUMVERTICES];
    1106         int         idlist[NUMVERTICES];
    1107 
    1108         ElementInput* input=this->inputs->GetControlInputData(control_enum,"gradient");   _assert_(input);
     1110
     1111        /*Get list of ids for this element and this control*/
     1112        int* idlist = xNew<int>(NUMVERTICES*N);
     1113        GradientIndexing(&idlist[0],control_index);
     1114
     1115        ControlInput* control_input=this->inputs->GetControlInput(control_enum); _assert_(control_input);
    11091116        this->GetVerticesLidList(&lidlist[0]);
    1110         GradientIndexing(&idlist[0],control_index);
    11111117
    11121118        /*Get values on vertices*/
    1113         if(input->ObjectEnum()==TriaInputEnum && input->GetInputInterpolationType()==P1Enum){
    1114                 _assert_(N==1);
    1115                 for(int i=0;i<NUMVERTICES;i++){
    1116                         values[i] = gradient[idlist[i]];
    1117                 }
    1118                 input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
    1119         }
    1120         else if(input->ObjectEnum()==TriaInputEnum && input->GetInputInterpolationType()==P0Enum){
    1121                 _assert_(N==1);
    1122                 input->SetInput(P0Enum,this->lid,gradient[idlist[0]]);
    1123         }
    1124         else if(input->ObjectEnum()==TransientInputEnum){
     1119        if(control_input->layout_enum==TriaInputEnum){
     1120                ElementInput* gradient_input = control_input->GetInput("gradient"); _assert_(gradient_input);
     1121                if(gradient_input->GetInputInterpolationType()==P1Enum){
     1122                        _assert_(N==1);
     1123                        for(int i=0;i<NUMVERTICES;i++){
     1124                                values[i] = gradient[idlist[i]];
     1125                        }
     1126                        gradient_input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
     1127                }
     1128                else if(gradient_input->GetInputInterpolationType()==P0Enum){
     1129                        _assert_(N==1);
     1130                        gradient_input->SetInput(P0Enum,this->lid,gradient[idlist[0]]);
     1131                }
     1132                else{
     1133                        _error_("not implemented yet");
     1134                }
     1135        }
     1136        else if(control_input->layout_enum==TransientInputEnum){
     1137                _assert_(N>1);
     1138                TransientInput* gradient_input = control_input->GetTransientInput("gradient"); _assert_(gradient_input);
     1139
    11251140                for(int n=0;n<N;n++){
    1126                         _error_("not implemented");
    1127                         //Input* new_input = new TriaInput(control_enum,gradient,P1Enum);
    1128                         //controlinput->SetInput(new_input,n);
    1129                         //controlinput->Configure(parameters);
     1141                        TriaInput* input_n = gradient_input->GetTriaInput(n); _assert_(input_n);
     1142                        if(input_n->GetInputInterpolationType()==P1Enum){
     1143                                _error_("not implemented");
     1144                        }
     1145                        else if(input_n->GetInputInterpolationType()==P0Enum){
     1146                                input_n->SetInput(P0Enum,this->lid,gradient[idlist[n]]);
     1147                        }
     1148                        else{
     1149                                _error_("not implemented yet");
     1150                        }
    11301151                }
    11311152        }
    11321153        else _error_("Type not supported");
     1154
     1155        /*Clean up*/
     1156        xDelete<int>(idlist);
    11331157
    11341158}/*}}}*/
    … …  
    24072431}
    24082432/*}}}*/
    2409 void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    2410 
    2411         _error_("DO WE NEED THIS FUNCTION?? REALLY?? see below");
    2412 
    2413         /*Get out if this is not an element input*/
    2414         if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
    2415 
    2416         /*Prepare index list*/
    2417         int idlist[NUMVERTICES];
    2418         GradientIndexing(&idlist[0],control_index);
    2419 
    2420         /*Get input (either in element or material)*/
    2421         ElementInput* input=this->inputs->GetControlInputData(control_enum,data);   _assert_(input);
    2422 
    2423         /*Intermediaries*/
    2424         int numindices;
    2425         int indices[NUMVERTICES];
    2426 
    2427         /*Check interpolation*/
    2428         int interpolation = input->GetInterpolation();
    2429         switch(interpolation){
    2430                 case P0Enum:
     2433void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,int N,const char* data,int offset){/*{{{*/
     2434
     2435        /*Get list of ids for this element and this control*/
     2436        int* idlist = xNew<int>(NUMVERTICES*N);
     2437        int  lidlist[NUMVERTICES];
     2438        this->GradientIndexing(&idlist[0],control_index);
     2439        this->GetVerticesLidList(&lidlist[0]);
     2440
     2441        /*Get Control*/
     2442        ControlInput* control_input=this->inputs->GetControlInput(control_enum); _assert_(control_input);
     2443
     2444        /*Get values on vertices*/
     2445        if(control_input->layout_enum==TriaInputEnum){
     2446                _assert_(N==1);
     2447                TriaInput* input = xDynamicCast<TriaInput*>(control_input->GetInput("value")); _assert_(input);
     2448
     2449                if(input->GetInputInterpolationType()==P1Enum){
     2450                        IssmDouble values[NUMVERTICES];
     2451                        input->Serve(NUMVERTICES,&lidlist[0]);
     2452                        for(int i=0;i<NUMVERTICES;i++) values[i] = input->element_values[i];
     2453                        vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     2454                }
     2455                else if(input->GetInputInterpolationType()==P0Enum){
    24312456                        input->Serve(1,&this->lid);
    2432                         break;
    2433                 case P1Enum:
    2434                         numindices = NUMVERTICES;
    2435                         for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
    2436                         input->Serve(numindices,&indices[0]);
    2437                         break;
    2438                 default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
    2439         }
    2440 
    2441         /*Flag as collapsed for later use*/
    2442         if(this->iscollapsed){
    2443                 xDynamicCast<PentaInput*>(input)->SetServeCollapsed(true);
    2444         }
    2445 
    2446         /* Start looping on the number of vertices: */
    2447         IssmDouble values[NUMVERTICES];
    2448         Gauss*gauss=this->NewGauss();
    2449         for(int iv=0;iv<NUMVERTICES;iv++){
    2450                 gauss->GaussVertex(iv);
    2451                 input->GetInputValue(&values[iv],gauss);
    2452         }
    2453         delete gauss;
    2454 
    2455         vector->SetValues(NUMVERTICES,idlist,&values[0],INS_VAL);
    2456 }
    2457 /*}}}*/
    2458 void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset){/*{{{*/
    2459 
    2460         /*Get input*/
    2461         ElementInput* input=this->inputs->GetControlInputData(control_enum,data);   _assert_(input);
    2462 
    2463         /*Lid list once for all*/
    2464         int lidlist[NUMVERTICES];
    2465         for(int i=0;i<NUMVERTICES;i++) lidlist[i] = vertices[i]->lid;
    2466 
    2467         /*Check what input we are dealing with*/
    2468         switch(input->ObjectEnum()){
    2469                 case TriaInputEnum:
    2470                           {
    2471                                 IssmDouble values[NUMVERTICES];
    2472                                 int        idlist[NUMVERTICES];
    2473 
    2474                                 TriaInput* triainput = xDynamicCast<TriaInput*>(input);
    2475 
    2476                                 /*Create list of indices and values for global vector*/
    2477                                 GradientIndexing(&idlist[0],control_index);
    2478 
    2479                                 if(triainput->GetInputInterpolationType()==P1Enum){
    2480                                         input->Serve(NUMVERTICES,&lidlist[0]);
    2481                                         for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
    2482                                         vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
    2483                                 }
    2484                                 else if(triainput->GetInputInterpolationType()==P0Enum){
    2485                                         input->Serve(1,&this->lid);
    2486                                         vector->SetValue(idlist[0],triainput->element_values[0],INS_VAL);
    2487                                 }
    2488                                 else{
    2489                                         _error_("not supported yet");
    2490                                 }
    2491                                 break;
    2492                           }
    2493 
    2494                 case TransientInputEnum:
    2495                                 {
    2496                                         TransientInput* transientinput = xDynamicCast<TransientInput*>(input);
    2497                                         int  N = transientinput->numtimesteps;
    2498                                         int* M = NULL;
    2499                                         parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    2500                                         int* idlist = xNew<int>(NUMVERTICES*N);
    2501                                         IssmDouble* values = xNew<IssmDouble>(NUMVERTICES*N);
    2502                                         for(int t=0;t<transientinput->numtimesteps;t++) {
    2503                                                 IssmDouble time = transientinput->GetTimeByOffset(t);
    2504                                                 _error_("not implemented");
    2505                                                 //TriaInput* timeinput = xDynamicCast<TriaInput*>(transientinput->GetTimeInput(time));
    2506                                                 //if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
    2507                                                 //input->Serve(NUMVERTICES,&lidlist[0]);
    2508                                                 ///*Create list of indices and values for global vector*/
    2509                                                 //for(int i=0;i<NUMVERTICES;i++){
    2510                                                 //              idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
    2511                                                 //              values[N*i+t] = timeinput->values[i];
    2512                                                 //}
    2513                                         }
    2514                                         vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
    2515                                         xDelete<int>(M);
    2516                                         xDelete<int>(idlist);
    2517                                         xDelete<IssmDouble>(values);
    2518                                         break;
    2519                                 }
    2520                 default: _error_("input "<<EnumToStringx(input->ObjectEnum())<<" not supported yet");
    2521         }
     2457                        vector->SetValue(idlist[0],input->element_values[0],INS_VAL);
     2458                }
     2459                else{
     2460                        _error_("not supported yet");
     2461                }
     2462        }
     2463        else if(control_input->layout_enum==TransientInputEnum){
     2464                TransientInput* input = control_input->GetTransientInput("value"); _assert_(input);
     2465                _assert_(N==input->numtimesteps);
     2466                IssmDouble* values = xNew<IssmDouble>(NUMVERTICES*N);
     2467
     2468                int count = 0;
     2469                for(int n=0;n<N;n++){
     2470                        TriaInput* input_n = input->GetTriaInput(n); _assert_(input_n);
     2471                        if(input_n->GetInputInterpolationType()==P1Enum){
     2472                                _error_("not implemented");
     2473                        }
     2474                        else if(input_n->GetInputInterpolationType()==P0Enum){
     2475                                input_n->Serve(1,&this->lid);
     2476                                values[n] = input_n->element_values[0];
     2477                                count++;
     2478                        }
     2479                        else{
     2480                                _error_("not implemented yet");
     2481                        }
     2482                }
     2483                vector->SetValues(count,idlist,values,INS_VAL);
     2484                xDelete<IssmDouble>(values);
     2485        }
     2486        else _error_("Type not supported");
     2487
     2488        /*Clean up*/
     2489        xDelete<int>(idlist);
    25222490}
    25232491/*}}}*/
    … …  
    39443912        IssmDouble  values[NUMVERTICES];
    39453913        int         lidlist[NUMVERTICES];
    3946         int         idlist[NUMVERTICES];
    39473914
    39483915        /*Get Domain type*/
    … …  
    39663933        if(!IsInputEnum(control_enum)) return;
    39673934
    3968         ElementInput* input=this->inputs->GetControlInputData(control_enum,"value");   _assert_(input);
     3935        /*Get list of ids for this element and this control*/
     3936        int* idlist = xNew<int>(NUMVERTICES*N);
     3937        GradientIndexing(&idlist[0],control_index);
     3938
     3939        ControlInput* control_input=this->inputs->GetControlInput(control_enum); _assert_(control_input);
    39693940        this->GetVerticesLidList(&lidlist[0]);
    3970         GradientIndexing(&idlist[0],control_index);
    3971 
    3972         /*Get values on vertices*/
    3973         if(input->ObjectEnum()==TriaInputEnum && input->GetInputInterpolationType()==P1Enum){
    3974                 _assert_(N==1);
    3975                 for(int i=0;i<NUMVERTICES;i++){
    3976                         values[i] = vector[idlist[i]];
    3977                 }
    3978                 input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
    3979         }
    3980         else if(input->ObjectEnum()==TriaInputEnum && input->GetInputInterpolationType()==P0Enum){
    3981                 _assert_(N==1);
    3982                 input->SetInput(P0Enum,this->lid,vector[idlist[0]]);
    3983         }
    3984         else if(input->ObjectEnum()==TransientInputEnum){
     3941
     3942        if(control_input->layout_enum==TriaInputEnum){
     3943                ElementInput* input = control_input->GetInput("value"); _assert_(input);
     3944                if(input->GetInputInterpolationType()==P1Enum){
     3945                        _assert_(N==1);
     3946                        for(int i=0;i<NUMVERTICES;i++){
     3947                                values[i] = vector[idlist[i]];
     3948                        }
     3949                        input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
     3950                }
     3951                else if(input->GetInputInterpolationType()==P0Enum){
     3952                        _assert_(N==1);
     3953                        input->SetInput(P0Enum,this->lid,vector[idlist[0]]);
     3954                }
     3955                else{
     3956                        _error_("not implemented yet");
     3957                }
     3958        }
     3959        else if(control_input->layout_enum==TransientInputEnum){
     3960                _assert_(N>1);
     3961                TransientInput* input = control_input->GetTransientInput("value"); _assert_(input);
     3962
    39853963                for(int n=0;n<N;n++){
    3986                         _error_("not implemented");
    3987                         //Input* new_input = new TriaInput(control_enum,values,P1Enum);
    3988                         //controlinput->SetInput(new_input,n);
    3989                         //controlinput->Configure(parameters);
     3964                        TriaInput* input_n = input->GetTriaInput(n); _assert_(input_n);
     3965                        if(input_n->GetInputInterpolationType()==P1Enum){
     3966                                _error_("not implemented");
     3967                        }
     3968                        else if(input_n->GetInputInterpolationType()==P0Enum){
     3969                                input_n->SetInput(P0Enum,this->lid,vector[idlist[n]]);
     3970                        }
     3971                        else{
     3972                                _error_("not implemented yet");
     3973                        }
    39903974                }
    39913975        }
    39923976        else _error_("Type not supported");
     3977
     3978        /*Clean up*/
     3979        xDelete<int>(idlist);
    39933980}
    39943981/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r25379 r25406  
    9191                int         GetNumberOfNodes(int enum_type);
    9292                int         GetNumberOfVertices(void);
    93                 void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset);
    94                 void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
     93                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,int N,const char* data,int offset);
    9594                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    9695                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp

    r25379 r25406  
    1313#include "./TriaInput.h"
    1414#include "./PentaInput.h"
     15#include "./TransientInput.h"
    1516//#include "../../toolkits/objects/Vector.h"
    1617
    … …  
    5051}
    5152/*}}}*/
     53ControlInput::ControlInput(int enum_in,int nbe, int nbv,int id,IssmDouble* times, int numtimes){/*{{{*/
     54
     55        this->enum_type   = enum_in;
     56        this->control_id  = id;
     57        this->layout_enum = TransientInputEnum; /*Tria or Penta?*/
     58
     59        this->values     =new TransientInput(enum_in,nbe,nbv,times,numtimes);
     60        this->savedvalues=new TransientInput(enum_in,nbe,nbv,times,numtimes);
     61        this->minvalues  =new TransientInput(enum_in,nbe,nbv,times,numtimes);
     62        this->maxvalues  =new TransientInput(enum_in,nbe,nbv,times,numtimes);
     63        this->gradient   =new TransientInput(enum_in,nbe,nbv,times,numtimes);
     64}
     65/*}}}*/
    5266ControlInput::~ControlInput(){/*{{{*/
    5367        if(values)      delete values;
    … …  
    6781        output->enum_type=this->enum_type;
    6882        output->control_id=this->control_id;
    69         output->layout_enum = this->control_id;
    70 
    71         if(values)      output->values      = xDynamicCast<ElementInput*>(this->values->copy());
    72         if(savedvalues) output->savedvalues = xDynamicCast<ElementInput*>(this->savedvalues->copy());
    73         if(minvalues)   output->minvalues   = xDynamicCast<ElementInput*>(this->minvalues->copy());
    74         if(maxvalues)   output->maxvalues   = xDynamicCast<ElementInput*>(this->maxvalues->copy());
    75         if(gradient)    output->gradient    = xDynamicCast<ElementInput*>(this->gradient->copy());
     83        output->layout_enum = this->layout_enum;
     84
     85        if(values)      output->values      = this->values->copy();
     86        if(savedvalues) output->savedvalues = this->savedvalues->copy();
     87        if(minvalues)   output->minvalues   = this->minvalues->copy();
     88        if(maxvalues)   output->maxvalues   = this->maxvalues->copy();
     89        if(gradient)    output->gradient    = this->gradient->copy();
    7690
    7791        return output;
     92}
     93/*}}}*/
     94void ControlInput::Configure(Parameters* params){/*{{{*/
     95        this->values->Configure(params);
     96        this->savedvalues->Configure(params);
     97        this->minvalues->Configure(params);
     98        this->maxvalues->Configure(params);
     99        this->gradient->Configure(params);
    78100}
    79101/*}}}*/
    … …  
    82104        _printf_("ControlInput:\n");
    83105        _printf_(setw(15)<<"   ControlInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<"\n");
    84         _printf_(setw(15)<<"   ControlInput "<<setw(25)<<left<<EnumToStringx(this->layout_enum)<<"\n");
     106        _printf_(setw(15)<<"   Layout      "<<setw(25)<<left<<EnumToStringx(this->layout_enum)<<"\n");
    85107        _printf_("---values: \n");     if (values)      values->Echo();
    86108        _printf_("---savedvalues: \n");if (savedvalues) savedvalues->Echo();
    … …  
    114136
    115137        /*Set input*/
    116         //TriaInput* input = xDynamicCast<TriaInput*>(this->inputs[id]);
    117         this->values->SetInput(interp,numindices,indices,values_in);
    118         this->minvalues->SetInput(interp,numindices,indices,values_min);
    119         this->maxvalues->SetInput(interp,numindices,indices,values_max);
     138        if(this->values->ObjectEnum()==TriaInputEnum || this->values->ObjectEnum()==PentaInputEnum){
     139                xDynamicCast<ElementInput*>(this->values)->SetInput(interp,numindices,indices,values_in);
     140                xDynamicCast<ElementInput*>(this->minvalues)->SetInput(interp,numindices,indices,values_min);
     141                xDynamicCast<ElementInput*>(this->maxvalues)->SetInput(interp,numindices,indices,values_max);
     142        }
     143        else{
     144                _error_("not supported");
     145        }
    120146}
    121147/*}}}*/
    … …  
    124150        _assert_(this);
    125151        _assert_(this->gradient);
    126         this->gradient->SetInput(interp,numindices,indices,values_in);
     152        if(this->gradient->ObjectEnum()==TriaInputEnum || this->gradient->ObjectEnum()==PentaInputEnum){
     153                xDynamicCast<ElementInput*>(this->gradient)->SetInput(interp,numindices,indices,values_in);
     154        }
     155        else{
     156                _error_("not supported");
     157        }
    127158}
    128159/*}}}*/
    … …  
    161192
    162193        /*Cast and return*/
    163         if(this->values->ObjectEnum()!=TriaInputEnum){
     194        if(this->values->ObjectEnum()==TriaInputEnum){
     195                return xDynamicCast<TriaInput*>(this->values);
     196        }
     197        else if(this->values->ObjectEnum()==TransientInputEnum){
     198                return xDynamicCast<TransientInput*>(this->values)->GetTriaInput();
     199        }
     200        else{
    164201                _error_("Cannot return a TriaInput");
    165202        }
    166         return xDynamicCast<TriaInput*>(this->values);
    167203
    168204}
    … …  
    180216ElementInput* ControlInput::GetInput(const char* data){/*{{{*/
    181217
    182         if(strcmp(data,"value")==0){
    183                 _assert_(values);
    184                 return values;
    185         }
    186         else if(strcmp(data,"savedvalues")==0){
    187                 _assert_(savedvalues);
    188                 return values;
    189         }
    190         else if (strcmp(data,"lowerbound")==0){
    191                 _assert_(minvalues);
    192                 return minvalues;
    193         }
    194         else if (strcmp(data,"upperbound")==0){
    195                 _assert_(maxvalues);
    196                 return maxvalues;
    197         }
    198         else if (strcmp(data,"gradient")==0){
    199                 _assert_(gradient);
    200                 return gradient;
    201         }
    202         else{
    203                 _error_("Data " << data << " not supported yet");
    204         }
    205 
    206 }
    207 /*}}}*/
     218        if(this->gradient->ObjectEnum()==TriaInputEnum || this->gradient->ObjectEnum()==PentaInputEnum){
     219                if(strcmp(data,"value")==0){
     220                        _assert_(values);
     221                        return xDynamicCast<ElementInput*>(values);
     222                }
     223                else if(strcmp(data,"savedvalues")==0){
     224                        _assert_(savedvalues);
     225                        return xDynamicCast<ElementInput*>(values);
     226                }
     227                else if (strcmp(data,"lowerbound")==0){
     228                        _assert_(minvalues);
     229                        return xDynamicCast<ElementInput*>(minvalues);
     230                }
     231                else if (strcmp(data,"upperbound")==0){
     232                        _assert_(maxvalues);
     233                        return xDynamicCast<ElementInput*>(maxvalues);
     234                }
     235                else if (strcmp(data,"gradient")==0){
     236                        _assert_(gradient);
     237                        return xDynamicCast<ElementInput*>(gradient);
     238                }
     239                else{
     240                        _error_("Data " << data << " not supported yet");
     241                }
     242        }
     243        else{
     244                int* temp = xNew<int>(3);
     245                _error_("not supported");
     246        }
     247
     248}
     249/*}}}*/
     250TransientInput* ControlInput::GetTransientInput(const char* data){/*{{{*/
     251
     252        if(this->values->ObjectEnum()==TransientInputEnum){
     253                if(strcmp(data,"value")==0){
     254                        _assert_(values);
     255                        return xDynamicCast<TransientInput*>(values);
     256                }
     257                else if(strcmp(data,"savedvalues")==0){
     258                        _assert_(savedvalues);
     259                        return xDynamicCast<TransientInput*>(values);
     260                }
     261                else if (strcmp(data,"lowerbound")==0){
     262                        _assert_(minvalues);
     263                        return xDynamicCast<TransientInput*>(minvalues);
     264                }
     265                else if (strcmp(data,"upperbound")==0){
     266                        _assert_(maxvalues);
     267                        return xDynamicCast<TransientInput*>(maxvalues);
     268                }
     269                else if (strcmp(data,"gradient")==0){
     270                        _assert_(gradient);
     271                        return xDynamicCast<TransientInput*>(gradient);
     272                }
     273                else{
     274                        _error_("Data " << data << " not supported yet");
     275                }
     276        }
     277        else{
     278                _error_("not supported");
     279        }
     280
     281}
     282/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r25379 r25406  
    1010class Gauss;
    1111class ElementInput;
     12class TransientInput;
    1213
    1314class ControlInput: public Input{
    1415
    1516        public:
    16                 int            control_id;
    17                 int            enum_type;
    18                 int            layout_enum;
    19                 ElementInput *gradient;
    20                 ElementInput *maxvalues;
    21                 ElementInput *minvalues;
    22                 ElementInput *savedvalues;
    23                 ElementInput *values;
     17                int    control_id;
     18                int    enum_type;
     19                int    layout_enum;
     20                Input *gradient;
     21                Input *maxvalues;
     22                Input *minvalues;
     23                Input *savedvalues;
     24                Input *values;
    2425
    2526                /*ControlInput constructors, destructors: {{{*/
    2627                ControlInput();
    2728                ControlInput(int nbe, int nbv,int input_layout_enum,int interp,int id);
     29                ControlInput(int enum_in,int nbe, int nbv,int id,IssmDouble* times, int numtimes);
    2830                ~ControlInput();
    2931                /*}}}*/
    3032                /*Object virtual functions definitions:{{{ */
    3133                Input* copy();
    32                 void    DeepEcho();
    33                 void    Echo();
    34                 int     Id();
    35                 void    Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
    36                 int     ObjectEnum();
     34                void   Configure(Parameters* params);
     35                void   DeepEcho();
     36                void   Echo();
     37                int    Id();
     38                void   Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
     39                int    ObjectEnum();
    3740                /*}}}*/
    3841                void SetInput(Input* in_input){_error_("not impelemented");};
    3942                void SetInput(Input* in_input,int timeoffset){_error_("not impelemented");};
    4043                ElementInput* GetInput(const char* data);
     44                TransientInput* GetTransientInput(const char* data);
    4145                void SetControl(int interp,int numindices,int* indices,IssmDouble* values_in,IssmDouble* values_min,IssmDouble* values_max);
    4246                void SetGradient(int interp,int numindices,int* indices,IssmDouble* values_in);

  • issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp

    r25379 r25406  
    690690        }
    691691}/*}}}*/
    692 void Inputs::SetTriaControlInput(int enum_in,int layout,int interpolation,int control_id,int numindices,int* indices,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max){/*{{{*/
     692void Inputs::SetControlInput(int enum_in,int layout,int interpolation,int control_id){/*{{{*/
    693693
    694694        bool recreate = false;
     695
    695696        /*Get input id*/
    696697        int id = EnumToIndex(enum_in);
    … …  
    711712        }
    712713
    713         /*Set input*/
    714         ControlInput* input = xDynamicCast<ControlInput*>(this->inputs[id]);
    715         input->SetControl(interpolation,numindices,indices,values,values_min,values_max);
     714}/*}}}*/
     715void Inputs::SetTransientControlInput(int enum_in,int control_id,IssmDouble* times,int numtimes){/*{{{*/
     716
     717        bool recreate = false;
     718
     719        /*Get input id*/
     720        int id = EnumToIndex(enum_in);
     721
     722        /*Create it if necessary*/
     723        if(this->inputs[id]){
     724                if(this->inputs[id]->ObjectEnum()!=ControlInputEnum){
     725                        delete this->inputs[id];
     726                        recreate = true;
     727                }
     728        }
     729        else{
     730                recreate = true;
     731        }
     732
     733        if(recreate){
     734                this->inputs[id] = new ControlInput(enum_in,this->numberofelements_local,this->numberofvertices_local,control_id,times,numtimes);
     735        }
     736
    716737}/*}}}*/
    717738void Inputs::SetTriaControlInputGradient(int enum_in,int interpolation,int numindices,int* indices,IssmDouble* values){/*{{{*/
    … …  
    822843        input->SetInput(interpolation,row,numindices,values);
    823844}/*}}}*/
    824 void Inputs::SetPentaControlInput(int enum_in,int layout,int interpolation,int control_id,int numindices,int* indices,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max){/*{{{*/
    825 
    826         bool recreate = false;
    827         /*Get input id*/
    828         int id = EnumToIndex(enum_in);
    829 
    830         /*Create it if necessary*/
    831         if(this->inputs[id]){
    832                 if(this->inputs[id]->ObjectEnum()!=ControlInputEnum){
    833                         delete this->inputs[id];
    834                         recreate = true;
    835                 }
    836         }
    837         else{
    838                 recreate = true;
    839         }
    840 
    841         if(recreate){
    842                 this->inputs[id] = new ControlInput(this->numberofelements_local,this->numberofvertices_local,layout,interpolation,control_id);
    843         }
    844 
    845         /*Set input*/
    846         ControlInput* input = xDynamicCast<ControlInput*>(this->inputs[id]);
    847         input->SetControl(interpolation,numindices,indices,values,values_min,values_max);
    848 }/*}}}*/
    849845void Inputs::SetPentaControlInputGradient(int enum_in,int interpolation,int numindices,int* indices,IssmDouble* values){/*{{{*/
    850846

  • issm/trunk-jpl/src/c/classes/Inputs/Inputs.h

    r25379 r25406  
    7373                void  SetDoubleInput(int enum_in,int index,IssmDouble value);
    7474                void  SetTransientInput(int enum_in,IssmDouble* times,int numtimes);
     75                void  SetControlInput(int enum_in,int layout,int interpolation,int control_id);
     76                void  SetTransientControlInput(int enum_in,int control_id,IssmDouble* times,int numtimes);
    7577                TransientInput* SetDatasetTransientInput(int enum_in,int id,IssmDouble* times,int numtimes);
    7678                void  SetArrayInput(int enum_in,int row,IssmDouble* layers,int numlayers);
    77                 void  SetTriaControlInput(int enum_in,int layout,int interpolation,int id,int numindices,int* indices,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max);
    7879                void  SetTriaControlInputGradient(int enum_in,int interpolation,int numindices,int* indices,IssmDouble* values);
    7980                void  SetTriaControlInputGradient(int enum_in,int interpolation,int numindices,int* indices,IssmDouble* values,int n);
    … …  
    8283                void  SetTriaInput(int enum_in,int interpolation,int numindices,int* indices,IssmDouble* values);
    8384                void  SetTriaInput(int enum_in,int interpolation,int row,int numindices,IssmDouble* values);
    84                 void  SetPentaControlInput(int enum_in,int layout,int interpolation,int id,int numindices,int* indices,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max);
    8585                void  SetPentaControlInputGradient(int enum_in,int interpolation,int numindices,int* indices,IssmDouble* values);
    8686                void  SetPentaDatasetInput(int enum_in,int id,int interpolation,int numindices,int* indices,IssmDouble* values);

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

    r25379 r25406  
    8080                output->inputs = xNew<Input*>(this->numtimesteps);
    8181                for(int i=0;i<this->numtimesteps;i++){
    82                         output->inputs[i] = this->inputs[i]->copy();
     82                        if(this->inputs[i]){
     83                                output->inputs[i] = this->inputs[i]->copy();
     84                        }
     85                        else{
     86                                output->inputs[i] = NULL;
     87                        }
    8388                }
    8489        }

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

    r25379 r25406  
    3838                /*}}}*/
    3939                /*Object virtual functions definitions:{{{*/
    40                 Input* copy();
     40                Input*  copy();
    4141                void    Configure(Parameters* params);
    4242                void    DeepEcho();

  • issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp

    r25317 r25406  
    3232                for(int j=0;j<elements->Size();j++){
    3333                        Element* element=(Element*)elements->GetObjectByOffset(j);
    34                         element->GetVectorFromControlInputs(vector,control_type[i],i,data,offset);
     34                        element->GetVectorFromControlInputs(vector,control_type[i],i,N[i],data,offset);
    3535                }
    3636                offset += M[i]*N[i];
    … …  
    4747
    4848        /*intermediary: */
    49         int                 N;
     49        int N;
    5050        Vector<IssmDouble>* vec_vector=NULL;
    5151
Note: See TracChangeset for help on using the changeset viewer.