Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22739) @@ -1387,21 +1387,40 @@ /*Get number of controls*/ int num_controls; + bool isautodiff; parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + /*Get number of vertices*/ int numvertices = this->GetNumberOfVertices(); - + if(isautodiff){ + int* N; + int* M; + int start = 0; + parameters->FindParam(&N,NULL,ControlInputSizeNEnum); + parameters->FindParam(&M,NULL,ControlInputSizeMEnum); + if(control_index>0) for(int n=0;nvertices[i]->Sid() + start + n*M[control_index]; + } + } + } + else{ + int M; + parameters->FindParam(&M,ControlInputSizeMEnum); /*get gradient indices*/ - if(onsid){ - for(int i=0;ivertices[i]->Sid() + control_index; - } - } - else{ - for(int i=0;ivertices[i]->Pid() + control_index; - } - } - + if(onsid){ + for(int i=0;ivertices[i]->Sid() + (control_index)*M; + } + } + else{ + for(int i=0;ivertices[i]->Pid() + (control_index)*M; + } + } + } } /*}}}*/ @@ -1619,6 +1638,6 @@ /*For now we only support nodal vectors*/ - if(M!=iomodel->numberofvertices) _error_("not supported"); - if(N!=1) _error_("not supported"); + //if(M!=iomodel->numberofvertices) _error_("not supported"); + //if(N!=1) _error_("not supported"); /*Recover vertices ids needed to initialize inputs*/ @@ -1637,4 +1656,44 @@ this->AddControlInput(input_enum,values,min_vector,max_vector,P1Enum,id); } + + else if(M==iomodel->numberofvertices+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(int t=0;tObjectEnum()){ + case TriaEnum: + values_input->AddTimeInput(new TriaInput(ControlInputValuesEnum,values,P1Enum)); + mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum)); + maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum)); + break; + case PentaEnum: + printf("-------------- file: Element.cpp line: %i\n",__LINE__); + values_input->AddTimeInput(new PentaInput(ControlInputValuesEnum,values,P1Enum)); + mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum)); + maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum)); + break; + case TetraEnum: + values_input->AddTimeInput(new TetraInput(ControlInputValuesEnum,values,P1Enum)); + mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum)); + maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum)); + break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id)); + xDelete(times); + } + else _error_("not currently supported type of M and N attempted"); /*clean up*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22739) @@ -202,4 +202,5 @@ virtual void ComputeEsaStrainAndVorticity(void)=0; virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; + virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M)=0; virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; virtual void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum)=0; @@ -224,5 +225,6 @@ virtual int GetNumberOfNodes(int enum_type)=0; virtual int GetNumberOfVertices(void)=0; - virtual void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid)=0; + virtual void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,int offset, bool onsid)=0; + virtual void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data, bool onsid)=0; virtual void GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0; virtual void GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0; @@ -285,4 +287,5 @@ virtual void ResetHooks()=0; virtual void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");}; + virtual void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M)=0; virtual void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22739) @@ -561,4 +561,30 @@ } /*}}}*/ +void Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){/*{{{*/ + + int vertexpidlist[NUMVERTICES]; + IssmDouble grad_list[NUMVERTICES]; + Input* grad_input=NULL; + Input* input=NULL; + + if(enum_type==MaterialsRheologyBbarEnum){ + input=(Input*)inputs->GetInput(MaterialsRheologyBEnum); + } + else if(enum_type==DamageDbarEnum){ + input=(Input*)inputs->GetInput(DamageDEnum); + } + else{ + input=inputs->GetInput(enum_type); + } + if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found"); + if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput"); + + GradientIndexing(&vertexpidlist[0],control_index); + for(int i=0;iSetGradient(grad_input); + +} +/*}}}*/ void Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/ @@ -1149,4 +1175,83 @@ ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data); } +/*}}}*/ +void Penta::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data,int offset,bool onsid){/*{{{*/ + + int* idlist = NULL; + IssmDouble* values = NULL; + int* M; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) _error_("Enum "<inputs->GetInput(control_enum); _assert_(input); + + + /*Cast to Controlinput*/ + if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); + ControlInput* controlinput = xDynamicCast(input); + + if(strcmp(data,"value")==0){ + input = controlinput->values; + } + else if (strcmp(data,"lowerbound")==0){ + input = controlinput->minvalues; + } + else if (strcmp(data,"upperbound")==0){ + input = controlinput->maxvalues; + } + else if (strcmp(data,"gradient")==0){ + input = controlinput->gradient; + } + else{ + _error_("Data " << data << " not supported yet"); + } + /*Check what input we are dealing with*/ + + switch(input->ObjectEnum()){ + case PentaInputEnum: + { + PentaInput* pentainput = xDynamicCast(input); + if(pentainput->interpolation_type!=P1Enum) _error_("not supported yet"); + + /*Create list of indices and values for global vector*/ + idlist = xNew(NUMVERTICES); + values = xNew(NUMVERTICES); + GradientIndexing(&idlist[0],control_index,true); + for(int i=0;ivalues[i]; + //if(this->Id()<=10)_printf_("index "<SetValues(NUMVERTICES,idlist,values,INS_VAL); + break; + } + + case TransientInputEnum: + { + parameters->FindParam(&M,NULL,ControlInputSizeMEnum); + TransientInput* transientinput = xDynamicCast(input); + int N = transientinput->numtimesteps; + idlist = xNew(NUMVERTICES*N); + values = xNew(NUMVERTICES*N); + for(int t=0;tnumtimesteps;t++) { + IssmDouble time = transientinput->GetTimeByOffset(t); + input = transientinput->GetTimeInput(time); + TriaInput* timeinput = xDynamicCast(input); + if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet"); + /*Create list of indices and values for global vector*/ + for(int i=0;ivertices[i]->Sid()+t*M[control_index]; + values[N*i+t] = timeinput->values[i]; + } + } + + vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL); + break; + } + default: _error_("input "<ObjectEnum()<<" not supported yet"); + + } + xDelete(idlist); + xDelete(values); + } /*}}}*/ void Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/ @@ -2252,4 +2357,46 @@ } +/*}}}*/ +void Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N, int M){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + int vertexpidlist[NUMVERTICES],control_init; + + /*Specific case for depth averaged quantities*/ + control_init=control_enum; + if(control_enum==MaterialsRheologyBbarEnum){ + control_enum=MaterialsRheologyBEnum; + if(!IsOnBase()) return; + } + if(control_enum==DamageDbarEnum){ + control_enum=DamageDEnum; + if(!IsOnBase()) return; + } + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&vertexpidlist[0],control_index); + + /*Get values on vertices*/ + for(int i=0;iinputs->GetInput(control_enum); _assert_(input); + if(input->ObjectEnum()!=ControlInputEnum){ + _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); + } + + ((ControlInput*)input)->SetInput(new_input); + + if(control_init==MaterialsRheologyBbarEnum){ + this->InputExtrude(control_enum,-1); + } + if(control_init==DamageDbarEnum){ + this->InputExtrude(control_enum,-1); + } +} /*}}}*/ void Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22739) @@ -58,4 +58,5 @@ void ComputeStressTensor(); void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M); void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); @@ -86,4 +87,5 @@ Penta* GetUpperPenta(void); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid); + void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); @@ -148,4 +150,5 @@ void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); void SetClone(int* minranks); + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22739) @@ -52,4 +52,5 @@ void ComputeStressTensor(){_error_("not implemented yet");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");}; void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; @@ -72,4 +73,5 @@ int GetNumberOfNodes(int enum_type){_error_("not implemented yet");}; int GetNumberOfVertices(void); + void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid){_error_("not implemented yet");}; void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid){_error_("not implemented yet");}; void GetVerticesCoordinates(IssmDouble** pxyz_list); @@ -135,4 +137,5 @@ void ResetFSBasalBoundaryCondition(void){_error_("not implemented yet");}; void ResetHooks(){_error_("not implemented yet");}; + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M){_error_("not implemented yet");}; void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");}; void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22739) @@ -52,4 +52,5 @@ void ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");}; void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; @@ -78,4 +79,5 @@ int GetNumberOfNodes(int enum_type){_error_("not implemented yet");}; int GetNumberOfVertices(void); + void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid){_error_("not implemented yet");}; void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid){_error_("not implemented yet");}; void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); @@ -143,4 +145,5 @@ void ResetHooks(); void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){_error_("not implemented yet");}; void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");}; void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22739) @@ -742,4 +742,37 @@ } /*}}}*/ +void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M){/*{{{*/ + + int idlist[NUMVERTICES]; + int gradidlist[NUMVERTICES]; + IssmDouble grad_list[NUMVERTICES]; + Input* grad_input=NULL; + + Input* input=inputs->GetInput(enum_type); + if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found"); + if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput"); + + GradientIndexing(&gradidlist[0],control_index,true); + + for(int n=0;nvertices[i]->Sid()+n*M; + grad_list[i]=gradient[idlist[i]]; + } + + ControlInput* controlinput = xDynamicCast(input); + if(controlinput->layout_enum!=TransientInputEnum){ + grad_input=new TriaInput(GradientEnum,grad_list,P1Enum); + controlinput->SetGradient(grad_input); + } + else{ + grad_input = new TriaInput(GradientEnum,grad_list,P1Enum); + controlinput->SetGradient(grad_input,n); + controlinput->Configure(parameters); + } + } + + +}/*}}}*/ void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/ @@ -1676,4 +1709,80 @@ ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data); +} +/*}}}*/ +void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data,int offset, bool onsid){/*{{{*/ + + int* idlist = NULL; + IssmDouble* values = NULL; + int* M; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) _error_("Enum "<inputs->GetInput(control_enum); _assert_(input); + + parameters->FindParam(&M,NULL,ControlInputSizeMEnum); + + /*Cast to Controlinput*/ + if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); + ControlInput* controlinput = xDynamicCast(input); + + if(strcmp(data,"value")==0){ + input = controlinput->values; + } + else if (strcmp(data,"lowerbound")==0){ + input = controlinput->minvalues; + } + else if (strcmp(data,"upperbound")==0){ + input = controlinput->maxvalues; + } + else if (strcmp(data,"gradient")==0){ + input = controlinput->gradient; + } + else{ + _error_("Data " << data << " not supported yet"); + } + /*Check what input we are dealing with*/ + + switch(input->ObjectEnum()){ + case TriaInputEnum: + { + TriaInput* triainput = xDynamicCast(input); + if(triainput->interpolation_type!=P1Enum) _error_("not supported yet"); + + /*Create list of indices and values for global vector*/ + idlist = xNew(NUMVERTICES); + values = xNew(NUMVERTICES); + GradientIndexing(&idlist[0],control_index,true); + for(int i=0;ivalues[i]; + } + vector->SetValues(NUMVERTICES,idlist,values,INS_VAL); + break; + } + + case TransientInputEnum: + { + TransientInput* transientinput = xDynamicCast(input); + int N = transientinput->numtimesteps; + idlist = xNew(NUMVERTICES*N); + values = xNew(NUMVERTICES*N); + for(int t=0;tnumtimesteps;t++) { + IssmDouble time = transientinput->GetTimeByOffset(t); + input = transientinput->GetTimeInput(time); + TriaInput* timeinput = xDynamicCast(input); + if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet"); + /*Create list of indices and values for global vector*/ + for(int i=0;ivertices[i]->Sid()+t*M[control_index]; + values[N*i+t] = timeinput->values[i]; + } + } + vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL); + break; + } + default: _error_("input "<ObjectEnum()<<" not supported yet"); + } + xDelete(idlist); + xDelete(values); } /*}}}*/ @@ -1952,5 +2061,5 @@ /*Get parameters: */ - iomodel->FindConstant(&yts,"md.constants.yts"); + iomodel->FindConstant(&yts,"md.constants.yts"); iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol"); iomodel->FindConstant(&ad_analysis, "md.autodiff.isautodiff"); @@ -1963,5 +2072,5 @@ /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<3;i++){ + for(i=0;i<3;i++){ tria_vertex_ids[i]=reCast(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab } @@ -2063,9 +2172,9 @@ /*DatasetInputs*/ if (control_analysis && iomodel->Data("md.inversion.cost_functions_coefficients")){ - + /*Generate cost functions associated with the iomodel*/ - char** cost_functions = NULL; - int* cost_functions_enums = NULL; - int num_cost_functions; + char** cost_functions = NULL; + int* cost_functions_enums = NULL; + int num_cost_functions; iomodel->FindConstant(&num_cost_functions,"md.inversion.num_cost_functions"); @@ -2091,5 +2200,4 @@ } } -/*}}}*/ void Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){/*{{{*/ @@ -3149,4 +3257,58 @@ _error_("not implemented yet"); +} +/*}}}*/ +void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + int idlist[NUMVERTICES],control_init; + + + /*Get Domain type*/ + int domaintype; + parameters->FindParam(&domaintype,DomainTypeEnum); + + /*Specific case for depth averaged quantities*/ + control_init=control_enum; + if(domaintype==Domain2DverticalEnum){ + if(control_enum==MaterialsRheologyBbarEnum){ + control_enum=MaterialsRheologyBEnum; + if(!IsOnBase()) return; + } + if(control_enum==DamageDbarEnum){ + control_enum=DamageDEnum; + if(!IsOnBase()) return; + } + } + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + Input* input = (Input*)this->inputs->GetInput(control_enum); _assert_(input); + if(input->ObjectEnum()!=ControlInputEnum){ + _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); + } + + ControlInput* controlinput = xDynamicCast(input); + input = controlinput->values; + + /*Get values on vertices*/ + for(int n=0;nvertices[i]->Sid()+n*M; + values[i]=vector[idlist[i]]; + //_printf_("index: "<ObjectEnum()==TriaInputEnum){ + Input* new_input = new TriaInput(control_enum,values,P1Enum); + controlinput->SetInput(new_input); + } + else if(input->ObjectEnum()==TransientInputEnum){ + Input* new_input = new TriaInput(control_enum,values,P1Enum); + controlinput->SetInput(new_input,n); + controlinput->Configure(parameters); + } + else _error_("Type not supported"); + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22739) @@ -63,4 +63,5 @@ void ComputeSurfaceNormalVelocity(); void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M); void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); @@ -85,4 +86,5 @@ int GetNumberOfNodes(int enum_type); int GetNumberOfVertices(void); + void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); @@ -118,4 +120,5 @@ void ResetHooks(); void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 22739) @@ -26,4 +26,5 @@ control_id=id; enum_type=in_enum_type; + layout_enum = input_layout_enum; _assert_(interp==P1Enum); @@ -48,10 +49,26 @@ } /*}}}*/ +ControlInput::ControlInput(int in_enum_type,int input_layout_enum,Input* input_pvalues,Input* input_pmin,Input* input_pmax,Input* input_pgrad,int interp,int id){/*{{{*/ + + this->control_id=id; + this->enum_type=in_enum_type; + this->layout_enum = input_layout_enum; + + _assert_(interp==P1Enum); + if(input_layout_enum!=TransientInputEnum) _error_("Wrong type of layout_enum, needs to be a TransientInputEnum"); + + this->values = input_pvalues; + this->savedvalues= NULL; + this->minvalues = input_pmin; + this->maxvalues = input_pmax; + this->gradient = input_pgrad; + } +/*}}}*/ ControlInput::~ControlInput(){/*{{{*/ - delete values; - delete savedvalues; - delete minvalues; - delete maxvalues; - delete gradient; + if(values) delete values; + if(savedvalues) delete savedvalues; + if(minvalues) delete minvalues; + if(maxvalues) delete maxvalues; + if(gradient) delete gradient; } /*}}}*/ @@ -65,4 +82,5 @@ output->enum_type=this->enum_type; output->control_id=this->control_id; + output->layout_enum = this->control_id; if(values) output->values = xDynamicCast(this->values->copy()); @@ -79,4 +97,5 @@ _printf_("ControlInput:\n"); _printf_(setw(15)<<" ControlInput "<enum_type)<<"\n"); + _printf_(setw(15)<<" ControlInput "<layout_enum)<<"\n"); _printf_("---values: \n"); if (values) values->Echo(); _printf_("---savedvalues: \n");if (savedvalues) savedvalues->Echo(); @@ -98,4 +117,5 @@ MARSHALLING(enum_type); MARSHALLING(control_id); + MARSHALLING(layout_enum); if (marshall_direction == MARSHALLING_BACKWARD){ @@ -115,4 +135,11 @@ gradient =new PentaInput(); break; + case TransientInputEnum: + values =new TransientInput(); + savedvalues=new TransientInput(); + minvalues =new TransientInput(); + maxvalues =new TransientInput(); + gradient =new TransientInput(); + break; default: _error_("Input of Enum " << EnumToStringx(enum_type) << " not supported yet by ControlInput"); @@ -146,4 +173,10 @@ }/*}}}*/ void ControlInput::Configure(Parameters* parameters){/*{{{*/ + if(this->values->ObjectEnum()==TransientInputEnum){ + this->values->Configure(parameters); + this->minvalues->Configure(parameters); + this->maxvalues->Configure(parameters); + this->gradient->Configure(parameters); + } /*do nothing: */ } @@ -194,5 +227,5 @@ }/*}}}*/ void ControlInput::GetVectorFromInputs(Vector* vector,int* doflist,const char* data){/*{{{*/ - if(strcmp(data,"value")==0){ + if(strcmp(data,"value")==0){ _assert_(values); values->GetVectorFromInputs(vector,doflist); @@ -226,5 +259,6 @@ this->savedvalues=xDynamicCast(this->values->copy()); }/*}}}*/ -void ControlInput::SetGradient(Input* gradient_in){/*{{{*/ +void ControlInput::SetGradient(Input* gradient_in,int timestep){/*{{{*/ +if(this->values->ObjectEnum()!=TransientInputEnum)_error_("you are in the wrong place, go home"); /*Get enum for current gradient*/ @@ -243,4 +277,40 @@ } + TransientInput* transient_input = xDynamicCast(this->gradient); + TransientInput* values_input = xDynamicCast(this->values); + if(values_input->numtimesteps==transient_input->numtimesteps){ + TransientInput* new_trans_input = new TransientInput(ControlInputGradEnum); + IssmDouble time = transient_input->GetTimeByOffset(timestep); + for(int i=0;inumtimesteps;i++){ + if(transient_input->timesteps[i]==time) new_trans_input->AddTimeInput(xDynamicCast(gradient_in),time); + else { + Input* input = transient_input->GetTimeInput(transient_input->timesteps[i]); + new_trans_input->AddTimeInput(xDynamicCast(input),transient_input->timesteps[i]); + } + } + this->gradient=new_trans_input; + } + else{ + IssmDouble time = values_input->GetTimeByOffset(timestep); + transient_input->AddTimeInput(gradient_in,time); + } +}/*}}}*/ +void ControlInput::SetGradient(Input* gradient_in){/*{{{*/ + + /*Get enum for current gradient*/ + switch(this->control_id){ + case 1: + gradient_in->ChangeEnum(Gradient1Enum); + break; + case 2: + gradient_in->ChangeEnum(Gradient2Enum); + break; + case 3: + gradient_in->ChangeEnum(Gradient3Enum); + break; + default: + _error_("more than 3 controls not implemented yet (Gradient " << this->control_id << " was requested). EnumDefinitions.h needs to be updated."); + } + /*Delete old gradient and assign new gradient*/ if(gradient) delete gradient; @@ -250,6 +320,27 @@ void ControlInput::SetInput(Input* in_input){/*{{{*/ + if(layout_enum==TransientInputEnum)_error_("need two arguments in SetInput for TransientInput Controls"); delete values; this->values=in_input; this->SaveValue(); //because this is what SpawnResult saves FIXME + +}/*}}}*/ +void ControlInput::SetInput(Input* in_input,int timeoffset){/*{{{*/ + Input* input = this->values; + if(input->ObjectEnum()!=TransientInputEnum)_error_("cannot have timeoffset argument if not TransientInput Control"); + TransientInput* transient_input = xDynamicCast(input); + IssmDouble time = transient_input->GetTimeByOffset(timeoffset); + TransientInput* new_trans_input = new TransientInput(ControlInputValuesEnum); + for(int i=0;inumtimesteps;i++){ + if(transient_input->timesteps[i]==time) new_trans_input->AddTimeInput(xDynamicCast(in_input),time); + else { + input = transient_input->GetTimeInput(transient_input->timesteps[i]); + new_trans_input->AddTimeInput(xDynamicCast(input),transient_input->timesteps[i]); + } + } + this->values=new_trans_input; + + // this->values->Echo(); + //this->values->Echo(); + //new_trans_input->Echo(); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 22739) @@ -18,4 +18,5 @@ int control_id; int enum_type; + int layout_enum; Input* gradient; Input* maxvalues; @@ -27,4 +28,5 @@ ControlInput(); ControlInput(int enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int interp,int id); + ControlInput(int enum_type,int enum_input,Input* pvalues,Input* pmin,Input* pmax,Input* pgrad,int interp,int id); ~ControlInput(); /*}}}*/ @@ -75,6 +77,8 @@ void Scale(IssmDouble scale_factor){_error_("not implemented yet");}; void Set(IssmDouble setvalue){_error_("Set not implemented yet");}; + void SetGradient(Input* gradient_in,int timestep); void SetGradient(Input* gradient_in); void SetInput(Input* in_input); + void SetInput(Input* in_input, int timeoffset); void UpdateValue(IssmDouble scalar); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 22739) @@ -83,4 +83,5 @@ int i; + printf("-------------- file: TransientInput.cpp line: %i\n",__LINE__); _printf_("TransientInput:\n"); @@ -232,5 +233,5 @@ /*First, recover current time from parameters: */ - this->parameters->FindParam(&time,TimeEnum); + parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ @@ -268,8 +269,6 @@ /*}}}*/ IssmDouble TransientInput::GetTimeByOffset(int offset){/*{{{*/ - if (offset < 0) offset=0; - _assert_(offsetnumtimesteps); - + _assert_(offset<(this->numtimesteps)); return this->timesteps[offset]; } @@ -469,5 +468,5 @@ int found; int offset; - bool interp; + bool interp=false; /*First, recover interp bool: */ Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22739) @@ -69,8 +69,10 @@ int num_responses,num_controls,numberofvertices,solution_type; femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - + int* N = NULL; + int N_add = 0; + if (solution_type == TransientSolutionEnum){ femmodel = input_struct->femmodel->copy(); - } + } IssmPDouble *Jlist = input_struct->Jlist; @@ -85,4 +87,5 @@ femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); + femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum); numberofvertices=femmodel->vertices->NumberOfVertices(); @@ -92,7 +95,9 @@ GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - for(int i=0;i(scaling_factors[c]); if(X[index]>XU[index]) X[index]=XU[index]; @@ -103,10 +108,9 @@ /*Start Tracing*/ simul_starttrace(femmodel); - /*Set X as our new control input and as INDEPENDENT*/ #ifdef _HAVE_AD_ - IssmDouble* aX=xNew(num_controls*numberofvertices,"t"); + IssmDouble* aX=xNew(intn,"t"); #else - IssmDouble* aX=xNew(num_controls*numberofvertices); + IssmDouble* aX=xNew(intn); #endif if(my_rank==0){ @@ -115,4 +119,5 @@ } } + ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); SetControlInputsFromVectorx(femmodel,aX); @@ -127,5 +132,4 @@ if(solution_type==TransientSolutionEnum){ IssmDouble restart_time; - femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum); femmodel->parameters->SetParam(restart_time,TimeEnum); @@ -139,5 +143,4 @@ DataSet* dependent_objects=NULL; IssmDouble J=0.; - femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); @@ -262,5 +265,7 @@ /*Add to totalgradient: */ - if(my_rank==0) for(int i=0;i=XU[index]) G[index]=0.; if(X[index]<=XL[index]) G[index]=0.; @@ -339,4 +345,7 @@ double *X = NULL; double *G = NULL; + int* N = NULL; + int offset = 0; + int N_add; /*Recover some parameters*/ @@ -350,4 +359,5 @@ femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); femmodel->parameters->SetParam(false,SaveResultsEnum); + femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum); numberofvertices=femmodel->vertices->NumberOfVertices(); @@ -377,5 +387,5 @@ Xpetsc->GetSize(&intn); delete Xpetsc; - _assert_(intn==numberofvertices*num_controls); + //_assert_(intn==numberofvertices*num_controls); /*Get problem dimension and initialize gradient and initial guess*/ @@ -384,11 +394,13 @@ /*Scale control for M1QN3*/ - for(int i=0;i(scaling_factors[c]); } - } - + N_add+=N[c]; + } + /*Allocate m1qn3 working arrays (see documentation)*/ long m = 100; @@ -437,16 +449,20 @@ GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - for(int i=0;i(scaling_factors[c]); if(X[index]>XU[index]) X[index]=XU[index]; if(X[index](intn); - for(int i=0;i(X[i]); + for(int i=0;i(X[i]); + } SetControlInputsFromVectorx(femmodel,aX); xDelete(aX); @@ -454,9 +470,8 @@ /*Set final gradient in inputs*/ IssmDouble* aG=xNew(intn); - for(int i=0;i(G[i]); + for(int i=0;i(G[i]); + } ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG); - xDelete(aG); - - /*Add last cost function to results*/ if (solution_type == TransientSolutionEnum){ @@ -465,4 +480,6 @@ femmodel->OutputControlsx(&femmodel->results); femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,1,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,Gradient1Enum,G,numberofvertices,intn/numberofvertices,1,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,InversionControlParametersEnum,X,numberofvertices,intn/numberofvertices,1,0)); } else{ @@ -471,4 +488,9 @@ femmodel->OutputControlsx(&femmodel->results); } + + xDelete(aG); + + /*Add last cost function to results*/ + /*Finalize*/ if(VerboseControl()) _printf0_(" preparing final solution\n"); Index: /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp (revision 22739) @@ -9,21 +9,50 @@ void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* gradient){ - /*Intermediaries*/ - int num_controls; - int *control_type = NULL; + bool isautodiff; + parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + if(isautodiff){ + /*Intermediaries*/ + int num_controls; + int *control_type = NULL; + int* M_all; + int* N_all; - /*Retrieve some parameters*/ - parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + parameters->FindParam(&M_all,NULL,ControlInputSizeMEnum); + parameters->FindParam(&N_all,NULL,ControlInputSizeNEnum); - for(int i=0;iSize();j++){ - Element* element=(Element*)elements->GetObjectByOffset(j); - element->ControlInputSetGradient(gradient,control_type[i],i); + int offset = 0; + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ControlInputSetGradient(gradient,control_type[i],i,offset,N_all[i],M_all[i]); + } + offset+=M_all[i]*N_all[i]; } + + /*Clean up and return*/ + xDelete(control_type); } + else{ + int num_controls; + int *control_type = NULL; - /*Clean up and return*/ - xDelete(control_type); + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + int offset = 0; + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ControlInputSetGradient(gradient,control_type[i],i); + } + } + + /*Clean up and return*/ + xDelete(control_type); + } } Index: /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 22739) @@ -9,32 +9,76 @@ void GetVectorFromControlInputsx(Vector** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,const char* data,bool onsid){/*{{{*/ - int num_controls; - int *control_type = NULL; - Vector* vector=NULL; + bool isautodiff; + parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + if(isautodiff){ + int* N = NULL; + int* M = NULL; + int num_controls; + int* control_type = NULL; + Vector* vector=NULL; - /*Retrieve some parameters*/ - parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + parameters->FindParam(&N,NULL,ControlInputSizeNEnum); + parameters->FindParam(&M,NULL,ControlInputSizeMEnum); - /*Allocate and populate gradient*/ - vector=new Vector(num_controls*vertices->NumberOfVertices()); + /*1. Get vector size*/ + int size = 0; + for(int i=0;iSize();j++){ - Element* element=(Element*)elements->GetObjectByOffset(j); - element->GetVectorFromControlInputs(vector,control_type[i],i,data,onsid); + + /*2. Allocate vector*/ + vector=new Vector(size); + + /*3. Populate vector*/ + int offset = 0; + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->GetVectorFromControlInputs(vector,control_type[i],i,data,offset,onsid); + } + offset += M[i]*N[i]; } + + vector->Assemble(); + + /*Assign output pointers:*/ + xDelete(control_type); + xDelete(M); + xDelete(N); + + *pvector=vector; + } + else{ + int num_controls; + int* control_type = NULL; + Vector* vector=NULL; + + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + + /*2. Allocate vector*/ + vector=new Vector(num_controls*vertices->NumberOfVertices()); + + /*3. Populate vector*/ + int offset = 0; + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->GetVectorFromControlInputs(vector,control_type[i],i,data,onsid); + } + } + vector->Assemble(); + + /*Assign output pointers:*/ + xDelete(control_type); + *pvector=vector; } - vector->Assemble(); - - /*Assign output pointers:*/ - xDelete(control_type); - *pvector=vector; }/*}}}*/ void GetVectorFromControlInputsx( IssmDouble** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, const char* data,bool onsid){/*{{{*/ - - /*output: */ - IssmDouble* vector=NULL; /*intermediary: */ @@ -42,5 +86,5 @@ GetVectorFromControlInputsx( &vec_vector, elements,nodes, vertices, loads, materials, parameters,data,onsid); - vector=vec_vector->ToMPISerial(); + IssmDouble* vector=vec_vector->ToMPISerial(); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 22739) @@ -10,7 +10,9 @@ -#if !defined(_HAVE_ADOLC_) -void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){ - +//#if !defined(_HAVE_ADOLC_) +void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){ + bool isautodiff; + iomodel->FindConstant(&isautodiff,"md.autodiff.isautodiff"); + if(!isautodiff){ /*Intermediary*/ bool control_analysis; @@ -65,4 +67,5 @@ } } + parameters->AddObject(new IntParam(ControlInputSizeMEnum,iomodel->numberofvertices)); for(int i=0;i(controls); } -#else -void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){ - +//} +//#else +//void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){ +else{ /*Intermediaries*/ - int num_independent_objects,M,N; + int num_independent_objects,M,N,M_par,N_par; char** names = NULL; int* types = NULL; + int* control_sizes = NULL; + int* M_all = NULL; IssmDouble* independent = NULL; - IssmDouble** independents_min = NULL; - IssmDouble** independents_max = NULL; + IssmDouble* independents_fullmin = NULL; + IssmDouble* independents_fullmax = NULL; bool control_analysis =false; iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol"); - + /*Now, return if no control*/ if(!control_analysis) return; @@ -154,12 +160,20 @@ - /*TODO: fetch min and max*/ - independents_min = xNew(num_independent_objects); - independents_max = xNew(num_independent_objects); - for(int i=0;i(num_independent_objects); /*create independent objects, and at the same time, fetch the corresponding independent variables, *and declare them as such in ADOLC: */ + iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters"); + iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters"); + iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes"); + + int* start_point = NULL; + start_point = xNew(num_independent_objects); + int counter = 0; + for(int i=0;iFetchData(&independent,&M,&N,iofieldname); - iomodel->FetchData(&independents_min[i],&M,&N,"md.autodiff.independent_min_parameters"); - iomodel->FetchData(&independents_max[i],&M,&N,"md.autodiff.independent_max_parameters"); - - _assert_(independent); + _assert_(independent); + _assert_(N==control_sizes[i]); + _printf_("N control size: "<(M*N); + independents_max = NULL; independents_max = xNew(M*N); + for(int m=0;mSize();j++){ Element* element=xDynamicCast(elements->GetObjectByOffset(j)); - element->ControlInputCreate(independent,independents_min[i],independents_max[i],iomodel,M,N,input_enum,i+1); + element->ControlInputCreate(independent,independents_min,independents_max,iomodel,M,N,input_enum,i+1); } xDelete(independent); + xDelete(independents_min); + xDelete(independents_max); + } else{ @@ -188,19 +218,21 @@ } } + parameters->AddObject(new IntVecParam(ControlInputSizeNEnum,control_sizes,num_independent_objects)); + parameters->AddObject(new IntVecParam(ControlInputSizeMEnum,M_all,num_independent_objects)); /*cleanup*/ for(int i=0;i(names[i]); - xDelete(independents_min[i]); - xDelete(independents_max[i]); - } - xDelete(independents_min); - xDelete(independents_max); + } xDelete(names); xDelete(types); - + xDelete(M_all); + xDelete(independents_fullmin); + xDelete(independents_fullmax); + xDelete(start_point); /*Step2: create cost functions (dependents)*/ return; } -#endif +} +//#endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 22739) @@ -69,5 +69,5 @@ /*Solution specific updates*/ if(VerboseMProcessor()) _printf0_(" updating elements and materials for control parameters" << "\n"); - UpdateElementsAndMaterialsControl(elements,materials,iomodel); + UpdateElementsAndMaterialsControl(elements,parameters,materials,iomodel); #ifdef _HAVE_DAKOTA_ if(VerboseMProcessor()) _printf0_(" updating elements and materials for uncertainty quantification" << "\n"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 22738) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 22739) @@ -18,5 +18,5 @@ void CreateParametersDakota(Parameters* parameters,IoModel* iomodel,char* rootpath); void CreateOutputDefinitions(Elements* elements, Parameters* parameters,IoModel* iomodel); -void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel); +void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel); void UpdateElementsAndMaterialsDakota(Elements* elements,Materials* materials, IoModel* iomodel); void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_type); Index: /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 22738) +++ /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 22739) @@ -9,19 +9,45 @@ void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){ - int num_controls; - int *control_type = NULL; + bool isautodiff; + femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + if(isautodiff){ + int num_controls; + int* control_type = NULL; + int* M = NULL; + int* N = NULL; - /*Retrieve some parameters*/ - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + /*Retrieve some parameters*/ + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum); + femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum); - for(int i=0;ielements->Size();j++){ - Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); - element->SetControlInputsFromVector(vector,control_type[i],i); + int offset = 0; + for(int i=0;ielements->Size();j++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + element->SetControlInputsFromVector(vector,control_type[i],i,offset,N[i],M[i]); + } + offset += M[i]*N[i]; } + + + xDelete(control_type); } + else{ - xDelete(control_type); + int num_controls; + int* control_type = NULL; + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + int offset = 0; + for(int i=0;ielements->Size();j++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + element->SetControlInputsFromVector(vector,control_type[i],i); + } + } + xDelete(control_type); + } }