Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 11316) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 11317) @@ -423,4 +423,8 @@ ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\ ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\ + ./modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp\ + ./modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.h\ + ./modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp\ + ./modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h\ ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\ ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\ Index: /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 11317) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 11317) @@ -0,0 +1,55 @@ +/*!\file GetVectorFromControlInputsx + * \brief retrieve vector from inputs in elements + */ + +#include "./GetVectorFromControlInputsx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void GetVectorFromControlInputsx(Vec* pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){ + + int num_controls; + int *control_type = NULL; + Vec vector=NULL; + + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + /*Allocate and populate gradient*/ + vector=NewVec(num_controls*vertices->NumberOfVertices()); + + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->GetVectorFromControlInputs(vector,control_type[i],i); + } + } + + VecAssemblyBegin(vector); + VecAssemblyEnd(vector); + + /*Assign output pointers:*/ + xfree((void**)&control_type); + *pvector=vector; +} + +void GetVectorFromControlInputsx( double** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){ + + /*output: */ + double* vector=NULL; + + /*intermediary: */ + Vec vec_vector=NULL; + + GetVectorFromControlInputsx( &vec_vector, elements,nodes, vertices, loads, materials, parameters); + VecToMPISerial(&vector,vec_vector); + + /*Free ressources:*/ + VecFree(&vec_vector); + + /*Assign output pointers:*/ + *pvector=vector; +} Index: /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.h (revision 11317) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.h (revision 11317) @@ -0,0 +1,15 @@ +/*!\file: GetVectorFromControlInputsx.h + */ + +#ifndef _GETVECTORFROMCONTROLINPUTSXX_H +#define _GETVECTORFROMCONTROLINPUTSXX_H + +#include "../../objects/objects.h" +#include "../../Container/Container.h" + +/* local prototypes: */ +void GetVectorFromControlInputsx( Vec* pvector, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters); +void GetVectorFromControlInputsx( double** pvector, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters); + +#endif /* _GETVECTORFROMCONTROLINPUTSXX_H */ + Index: /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 11317) @@ -10,5 +10,4 @@ void GetVectorFromInputsx( Vec* pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name, int type){ - int i; @@ -41,5 +40,4 @@ /*Assign output pointers:*/ *pvector=vector; - } Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 11317) @@ -27,5 +27,5 @@ /*Allocate gradient_list */ gradient_list = (Vec*)xmalloc(num_controls*sizeof(Vec)); - norm_list = (double*)xmalloc(num_controls*sizeof(double)); + norm_list = (double*)xmalloc(num_controls*sizeof(double)); for(i=0;iFindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->SetControlInputsFromVector(vector,control_type[i],i); + } + } + + xfree((void**)&control_type); +} + +void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vec vector){ + + double* serial_vector=NULL; + + VecToMPISerial(&serial_vector,vector); + + SetControlInputsFromVectorx(elements,nodes, vertices, loads, materials, parameters,serial_vector); + + /*Free ressources:*/ + xfree((void**)&serial_vector); +} Index: /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h (revision 11317) +++ /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h (revision 11317) @@ -0,0 +1,14 @@ +/*!\file: SetControlInputsFromVectorx.h + */ + +#ifndef _SETCONTROLINPUTSXFROMVECTOR_H +#define _SETCONTROLINPUTSXFROMVECTOR_H + +#include "../../objects/objects.h" +#include "../../Container/Container.h" + +/* local prototypes: */ +void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vec vector); +void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,double* vector); + +#endif Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 11316) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 11317) @@ -31,4 +31,6 @@ #include "./GetSolutionFromInputsx/GetSolutionFromInputsx.h" #include "./GetVectorFromInputsx/GetVectorFromInputsx.h" +#include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h" +#include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h" #include "./Gradjx/Gradjx.h" #include "./GroundinglineMigrationx/GroundinglineMigrationx.h" Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 11317) @@ -103,4 +103,6 @@ virtual void ControlInputSetGradient(double* gradient,int enum_type,int control_index)=0; virtual void ControlInputScaleGradient(int enum_type, double scale)=0; + virtual void GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index)=0; + virtual void SetControlInputsFromVector(double* vector,int control_enum,int control_index)=0; virtual void InputControlUpdate(double scalar,bool save_parameter)=0; #endif Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11317) @@ -1607,5 +1607,7 @@ this->results->AddObject((Object*)input->SpawnResult(step,time)); #ifdef _HAVE_CONTROL_ - if(input->ObjectEnum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); + if(input->ObjectEnum()==ControlInputEnum){ + if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); + } #endif @@ -5129,4 +5131,58 @@ } /*}}}*/ +/*FUNCTION Penta::GetVectorFromControlInputs{{{1*/ +void Penta::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index){ + + int doflist1[NUMVERTICES]; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&doflist1[0],control_index); + + /*Get input (either in element or material)*/ + Input* input=inputs->GetInput(control_enum); + if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum)); + + /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ + input->GetVectorFromInputs(vector,&doflist1[0]); +} +/*}}}*/ +/*FUNCTION Penta::SetControlInputsFromVector{{{1*/ +void Penta::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ + + double values[NUMVERTICES]; + int doflist1[NUMVERTICES]; + Input *input = NULL; + Input *new_input = NULL; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&doflist1[0],control_index); + + /*Get values on vertices*/ + for (int i=0;iinputs->GetInput(control_enum); _assert_(input); + } + else{ + input=(Input*)this->inputs->GetInput(control_enum); _assert_(input); + } + + if (input->ObjectEnum()!=ControlInputEnum){ + _error_("input %s is not a ControlInput",EnumToStringx(control_enum)); + } + + ((ControlInput*)input)->SetInput(new_input); +} +/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11317) @@ -149,4 +149,6 @@ void GradjBbarPattyn(Vec gradient,int control_index); void GradjBbarStokes(Vec gradient,int control_index); + void GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index); + void SetControlInputsFromVector(double* vector,int control_enum,int control_index); void ControlInputGetGradient(Vec gradient,int enum_type,int control_index); void ControlInputScaleGradient(int enum_type,double scale); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11317) @@ -1450,5 +1450,7 @@ #ifdef _HAVE_CONTROL_ - if(input->ObjectEnum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); + if(input->ObjectEnum()==ControlInputEnum){ + if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); + } #endif } @@ -4728,4 +4730,58 @@ } /*}}}*/ +/*FUNCTION Tria::GetVectorFromControlInputs{{{1*/ +void Tria::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index){ + + int doflist1[NUMVERTICES]; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&doflist1[0],control_index); + + /*Get input (either in element or material)*/ + Input* input=inputs->GetInput(control_enum); + if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum)); + + /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ + input->GetVectorFromInputs(vector,&doflist1[0]); +} +/*}}}*/ +/*FUNCTION Tria::SetControlInputsFromVector{{{1*/ +void Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ + + double values[NUMVERTICES]; + int doflist1[NUMVERTICES]; + Input *input = NULL; + Input *new_input = NULL; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&doflist1[0],control_index); + + /*Get values on vertices*/ + for (int i=0;iinputs->GetInput(control_enum); _assert_(input); + } + else{ + input=(Input*)this->inputs->GetInput(control_enum); _assert_(input); + } + + if (input->ObjectEnum()!=ControlInputEnum){ + _error_("input %s is not a ControlInput",EnumToStringx(control_enum)); + } + + ((ControlInput*)input)->SetInput(new_input); +} +/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11317) @@ -150,4 +150,6 @@ void GradjVxBalancedthickness(Vec gradient,int control_index); void GradjVyBalancedthickness(Vec gradient,int control_index); + void GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index); + void SetControlInputsFromVector(double* vector,int control_enum,int control_index); void ControlInputGetGradient(Vec gradient,int enum_type,int control_index); void ControlInputScaleGradient(int enum_type,double scale); Index: /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.cpp (revision 11317) @@ -386,4 +386,11 @@ }/*}}}*/ +/*FUNCTION ControlInput::SetInput{{{1*/ +void ControlInput::SetInput(Input* in_input){ + + delete values; this->values=in_input; + this->SaveValue(); //because this is what SpawnResult saves FIXME + +}/*}}}*/ /*FUNCTION ControlInput::SpawnResult{{{1*/ ElementResult* ControlInput::SpawnResult(int step, double time){ @@ -396,4 +403,5 @@ /*FUNCTION ControlInput::SpawnGradient{{{1*/ ElementResult* ControlInput::SpawnGradient(int step, double time){ + _assert_(gradient); return gradient->SpawnResult(step,time); }/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h (revision 11316) +++ /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h (revision 11317) @@ -54,4 +54,5 @@ /*}}}*/ /*numerics: {{{1*/ + void SetInput(Input* in_input); void GetInputValue(bool* pvalue); void GetInputValue(int* pvalue); Index: /issm/trunk-jpl/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/control_core.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/solutions/control_core.cpp (revision 11317) @@ -70,5 +70,5 @@ /*Initialize responses: */ J=(double*)xmalloc(nsteps*sizeof(double)); - step_responses=(int*)xmalloc(num_responses*sizeof(double)); + step_responses=(int*)xmalloc(num_responses*sizeof(int)); /*Initialize some of the BrentSearch arguments: */ Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 11316) +++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 11317) @@ -25,10 +25,12 @@ /*TAO*/ - int ierr,numberofvertices; - AppCtx user; - TaoSolver tao; - Vec X = NULL; - Vec XL = NULL; - Vec XU = NULL; + int ierr,numberofvertices; + int num_controls; + AppCtx user; + TaoSolver tao; + int *control_list = NULL; + Vec X = NULL; + Vec XL = NULL; + Vec XU = NULL; /*Initialize TAO*/ @@ -37,4 +39,8 @@ ierr = TaoInitialize(&argc,&args,(char*)0,""); if(ierr) _error_("Could not initialize Tao"); + + /*Recover some parameters*/ + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum); /*Initialize TAO*/ @@ -52,7 +58,7 @@ XL=NewVec(numberofvertices); VecSet(XL,1.); XU=NewVec(numberofvertices); VecSet(XU,200.); - TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU); + //TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU); - GetVectorFromInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum,VertexEnum); + GetVectorFromControlInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); TaoSetInitialVector(tao,X); @@ -64,8 +70,11 @@ TaoView(tao,PETSC_VIEWER_STDOUT_WORLD); TaoGetSolutionVector(tao,&X); - InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum); + SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); + for(int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]); + } /*Clean up and return*/ + xfree((void**)&control_list); VecFree(&X); TaoDestroy(&tao); @@ -75,19 +84,43 @@ /*Retreive arguments*/ - AppCtx *user = (AppCtx*)userCtx; - FemModel *femmodel = user->femmodel; - Vec gradient = NULL; + int solution_type,num_cost_functions; + AppCtx *user = (AppCtx *)userCtx; + FemModel *femmodel = user->femmodel; + int *cost_functions = NULL; + double *cost_functionsd= NULL; + Vec gradient = NULL; - int costfunction=SurfaceAbsVelMisfitEnum; - femmodel->parameters->SetParam(&costfunction,1,1,StepResponsesEnum); - //InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum); - InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum); - adjointdiagnostic_core(user->femmodel); - //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum); + /*Set new variable*/ + SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); + + /*Recover some parameters*/ + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum); + femmodel->parameters->FindParam(&cost_functionsd,NULL,NULL,InversionCostFunctionsEnum); + + /*Prepare objective function*/ + cost_functions=(int*)xmalloc(num_cost_functions*sizeof(int)); + for(int i=0;iparameters->SetParam(cost_functions,1,num_cost_functions,StepResponsesEnum); + + /*Compute solution and adjoint*/ + void (*solutioncore)(FemModel*)=NULL; + void (*adjointcore)(FemModel*)=NULL; + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); + solutioncore(femmodel); + adjointcore(femmodel); + + /*Compute objective function*/ + CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + /*Compute gradient*/ Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - VecScale(gradient,-1.); - VecCopy(gradient,G); - VecFree(&gradient); - CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + VecCopy(gradient,G); VecFree(&gradient); + VecScale(G,-1.); + + /*Clean-up and return*/ + xfree((void**)&cost_functions); + xfree((void**)&cost_functionsd); return 0; }