Changeset 11317
- Timestamp:
- 02/03/12 13:50:06 (13 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 6 added
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Makefile.am
r11295 r11317 423 423 ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\ 424 424 ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\ 425 ./modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp\ 426 ./modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.h\ 427 ./modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp\ 428 ./modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h\ 425 429 ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\ 426 430 ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\ -
issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
r8387 r11317 10 10 11 11 void GetVectorFromInputsx( Vec* pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name, int type){ 12 13 12 14 13 int i; … … 41 40 /*Assign output pointers:*/ 42 41 *pvector=vector; 43 44 42 } 45 43 -
issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
r11313 r11317 27 27 /*Allocate gradient_list */ 28 28 gradient_list = (Vec*)xmalloc(num_controls*sizeof(Vec)); 29 norm_list 29 norm_list = (double*)xmalloc(num_controls*sizeof(double)); 30 30 for(i=0;i<num_controls;i++){ 31 31 gradient_list[i]=NewVec(num_controls*numberofvertices); … … 59 59 60 60 /*Clean-up and assign output pointer*/ 61 *pnorm_list=norm_list; 62 *pgradient=gradient; 61 if(pnorm_list){ 62 *pnorm_list=norm_list; 63 } 64 else{ 65 xfree((void**)&norm_list); 66 } 67 if(pgradient) *pgradient=gradient; 63 68 xfree((void**)&gradient_list); 64 69 xfree((void**)&control_type); -
issm/trunk-jpl/src/c/modules/modules.h
r11295 r11317 31 31 #include "./GetSolutionFromInputsx/GetSolutionFromInputsx.h" 32 32 #include "./GetVectorFromInputsx/GetVectorFromInputsx.h" 33 #include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h" 34 #include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h" 33 35 #include "./Gradjx/Gradjx.h" 34 36 #include "./GroundinglineMigrationx/GroundinglineMigrationx.h" -
issm/trunk-jpl/src/c/objects/Elements/Element.h
r11313 r11317 103 103 virtual void ControlInputSetGradient(double* gradient,int enum_type,int control_index)=0; 104 104 virtual void ControlInputScaleGradient(int enum_type, double scale)=0; 105 virtual void GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index)=0; 106 virtual void SetControlInputsFromVector(double* vector,int control_enum,int control_index)=0; 105 107 virtual void InputControlUpdate(double scalar,bool save_parameter)=0; 106 108 #endif -
issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
r11314 r11317 1607 1607 this->results->AddObject((Object*)input->SpawnResult(step,time)); 1608 1608 #ifdef _HAVE_CONTROL_ 1609 if(input->ObjectEnum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); 1609 if(input->ObjectEnum()==ControlInputEnum){ 1610 if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); 1611 } 1610 1612 #endif 1611 1613 … … 5129 5131 } 5130 5132 /*}}}*/ 5133 /*FUNCTION Penta::GetVectorFromControlInputs{{{1*/ 5134 void Penta::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index){ 5135 5136 int doflist1[NUMVERTICES]; 5137 5138 /*Get out if this is not an element input*/ 5139 if(!IsInput(control_enum)) return; 5140 5141 /*Prepare index list*/ 5142 GradientIndexing(&doflist1[0],control_index); 5143 5144 /*Get input (either in element or material)*/ 5145 Input* input=inputs->GetInput(control_enum); 5146 if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum)); 5147 5148 /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ 5149 input->GetVectorFromInputs(vector,&doflist1[0]); 5150 } 5151 /*}}}*/ 5152 /*FUNCTION Penta::SetControlInputsFromVector{{{1*/ 5153 void Penta::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ 5154 5155 double values[NUMVERTICES]; 5156 int doflist1[NUMVERTICES]; 5157 Input *input = NULL; 5158 Input *new_input = NULL; 5159 5160 /*Get out if this is not an element input*/ 5161 if(!IsInput(control_enum)) return; 5162 5163 /*Prepare index list*/ 5164 GradientIndexing(&doflist1[0],control_index); 5165 5166 /*Get values on vertices*/ 5167 for (int i=0;i<NUMVERTICES;i++){ 5168 values[i]=vector[doflist1[i]]; 5169 } 5170 new_input = new PentaP1Input(control_enum,values); 5171 5172 5173 if(control_enum==MaterialsRheologyBbarEnum){ 5174 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); 5175 } 5176 else{ 5177 input=(Input*)this->inputs->GetInput(control_enum); _assert_(input); 5178 } 5179 5180 if (input->ObjectEnum()!=ControlInputEnum){ 5181 _error_("input %s is not a ControlInput",EnumToStringx(control_enum)); 5182 } 5183 5184 ((ControlInput*)input)->SetInput(new_input); 5185 } 5186 /*}}}*/ 5131 5187 #endif 5132 5188 -
issm/trunk-jpl/src/c/objects/Elements/Penta.h
r11313 r11317 149 149 void GradjBbarPattyn(Vec gradient,int control_index); 150 150 void GradjBbarStokes(Vec gradient,int control_index); 151 void GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index); 152 void SetControlInputsFromVector(double* vector,int control_enum,int control_index); 151 153 void ControlInputGetGradient(Vec gradient,int enum_type,int control_index); 152 154 void ControlInputScaleGradient(int enum_type,double scale); -
issm/trunk-jpl/src/c/objects/Elements/Tria.cpp
r11313 r11317 1450 1450 1451 1451 #ifdef _HAVE_CONTROL_ 1452 if(input->ObjectEnum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); 1452 if(input->ObjectEnum()==ControlInputEnum){ 1453 if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); 1454 } 1453 1455 #endif 1454 1456 } … … 4728 4730 } 4729 4731 /*}}}*/ 4732 /*FUNCTION Tria::GetVectorFromControlInputs{{{1*/ 4733 void Tria::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index){ 4734 4735 int doflist1[NUMVERTICES]; 4736 4737 /*Get out if this is not an element input*/ 4738 if(!IsInput(control_enum)) return; 4739 4740 /*Prepare index list*/ 4741 GradientIndexing(&doflist1[0],control_index); 4742 4743 /*Get input (either in element or material)*/ 4744 Input* input=inputs->GetInput(control_enum); 4745 if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum)); 4746 4747 /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ 4748 input->GetVectorFromInputs(vector,&doflist1[0]); 4749 } 4750 /*}}}*/ 4751 /*FUNCTION Tria::SetControlInputsFromVector{{{1*/ 4752 void Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){ 4753 4754 double values[NUMVERTICES]; 4755 int doflist1[NUMVERTICES]; 4756 Input *input = NULL; 4757 Input *new_input = NULL; 4758 4759 /*Get out if this is not an element input*/ 4760 if(!IsInput(control_enum)) return; 4761 4762 /*Prepare index list*/ 4763 GradientIndexing(&doflist1[0],control_index); 4764 4765 /*Get values on vertices*/ 4766 for (int i=0;i<NUMVERTICES;i++){ 4767 values[i]=vector[doflist1[i]]; 4768 } 4769 new_input = new TriaP1Input(control_enum,values); 4770 4771 4772 if(control_enum==MaterialsRheologyBbarEnum){ 4773 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); 4774 } 4775 else{ 4776 input=(Input*)this->inputs->GetInput(control_enum); _assert_(input); 4777 } 4778 4779 if (input->ObjectEnum()!=ControlInputEnum){ 4780 _error_("input %s is not a ControlInput",EnumToStringx(control_enum)); 4781 } 4782 4783 ((ControlInput*)input)->SetInput(new_input); 4784 } 4785 /*}}}*/ 4730 4786 #endif 4731 4787 -
issm/trunk-jpl/src/c/objects/Elements/Tria.h
r11313 r11317 150 150 void GradjVxBalancedthickness(Vec gradient,int control_index); 151 151 void GradjVyBalancedthickness(Vec gradient,int control_index); 152 void GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index); 153 void SetControlInputsFromVector(double* vector,int control_enum,int control_index); 152 154 void ControlInputGetGradient(Vec gradient,int enum_type,int control_index); 153 155 void ControlInputScaleGradient(int enum_type,double scale); -
issm/trunk-jpl/src/c/objects/Inputs/ControlInput.cpp
r11291 r11317 386 386 387 387 }/*}}}*/ 388 /*FUNCTION ControlInput::SetInput{{{1*/ 389 void ControlInput::SetInput(Input* in_input){ 390 391 delete values; this->values=in_input; 392 this->SaveValue(); //because this is what SpawnResult saves FIXME 393 394 }/*}}}*/ 388 395 /*FUNCTION ControlInput::SpawnResult{{{1*/ 389 396 ElementResult* ControlInput::SpawnResult(int step, double time){ … … 396 403 /*FUNCTION ControlInput::SpawnGradient{{{1*/ 397 404 ElementResult* ControlInput::SpawnGradient(int step, double time){ 405 _assert_(gradient); 398 406 return gradient->SpawnResult(step,time); 399 407 }/*}}}*/ -
issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h
r10135 r11317 54 54 /*}}}*/ 55 55 /*numerics: {{{1*/ 56 void SetInput(Input* in_input); 56 57 void GetInputValue(bool* pvalue); 57 58 void GetInputValue(int* pvalue); -
issm/trunk-jpl/src/c/solutions/control_core.cpp
r11306 r11317 70 70 /*Initialize responses: */ 71 71 J=(double*)xmalloc(nsteps*sizeof(double)); 72 step_responses=(int*)xmalloc(num_responses*sizeof( double));72 step_responses=(int*)xmalloc(num_responses*sizeof(int)); 73 73 74 74 /*Initialize some of the BrentSearch arguments: */ -
issm/trunk-jpl/src/c/solutions/controltao_core.cpp
r11313 r11317 25 25 26 26 /*TAO*/ 27 int ierr,numberofvertices; 28 AppCtx user; 29 TaoSolver tao; 30 Vec X = NULL; 31 Vec XL = NULL; 32 Vec XU = NULL; 27 int ierr,numberofvertices; 28 int num_controls; 29 AppCtx user; 30 TaoSolver tao; 31 int *control_list = NULL; 32 Vec X = NULL; 33 Vec XL = NULL; 34 Vec XU = NULL; 33 35 34 36 /*Initialize TAO*/ … … 37 39 ierr = TaoInitialize(&argc,&args,(char*)0,""); 38 40 if(ierr) _error_("Could not initialize Tao"); 41 42 /*Recover some parameters*/ 43 femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); 44 femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum); 39 45 40 46 /*Initialize TAO*/ … … 52 58 XL=NewVec(numberofvertices); VecSet(XL,1.); 53 59 XU=NewVec(numberofvertices); VecSet(XU,200.); 54 TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);60 //TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU); 55 61 56 GetVectorFrom Inputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum,VertexEnum);62 GetVectorFromControlInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); 57 63 TaoSetInitialVector(tao,X); 58 64 … … 64 70 TaoView(tao,PETSC_VIEWER_STDOUT_WORLD); 65 71 TaoGetSolutionVector(tao,&X); 66 InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum); 67 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum); 72 SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); 73 for(int i=0;i<num_controls;i++){ 74 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]); 75 } 68 76 69 77 /*Clean up and return*/ 78 xfree((void**)&control_list); 70 79 VecFree(&X); 71 80 TaoDestroy(&tao); … … 75 84 76 85 /*Retreive arguments*/ 77 AppCtx *user = (AppCtx*)userCtx; 78 FemModel *femmodel = user->femmodel; 79 Vec gradient = NULL; 86 int solution_type,num_cost_functions; 87 AppCtx *user = (AppCtx *)userCtx; 88 FemModel *femmodel = user->femmodel; 89 int *cost_functions = NULL; 90 double *cost_functionsd= NULL; 91 Vec gradient = NULL; 80 92 81 int costfunction=SurfaceAbsVelMisfitEnum; 82 femmodel->parameters->SetParam(&costfunction,1,1,StepResponsesEnum); 83 //InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum); 84 InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum); 85 adjointdiagnostic_core(user->femmodel); 86 //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum); 93 /*Set new variable*/ 94 SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); 95 96 /*Recover some parameters*/ 97 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); 98 femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum); 99 femmodel->parameters->FindParam(&cost_functionsd,NULL,NULL,InversionCostFunctionsEnum); 100 101 /*Prepare objective function*/ 102 cost_functions=(int*)xmalloc(num_cost_functions*sizeof(int)); 103 for(int i=0;i<num_cost_functions;i++) cost_functions[i]=(int)cost_functionsd[i]; //FIXME 104 femmodel->parameters->SetParam(cost_functions,1,num_cost_functions,StepResponsesEnum); 105 106 /*Compute solution and adjoint*/ 107 void (*solutioncore)(FemModel*)=NULL; 108 void (*adjointcore)(FemModel*)=NULL; 109 CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); 110 AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); 111 solutioncore(femmodel); 112 adjointcore(femmodel); 113 114 /*Compute objective function*/ 115 CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); 116 117 /*Compute gradient*/ 87 118 Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); 88 VecScale(gradient,-1.); 89 VecCopy(gradient,G); 90 VecFree(&gradient); 91 CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); 119 VecCopy(gradient,G); VecFree(&gradient); 120 VecScale(G,-1.); 121 122 /*Clean-up and return*/ 123 xfree((void**)&cost_functions); 124 xfree((void**)&cost_functionsd); 92 125 return 0; 93 126 }
Note:
See TracChangeset
for help on using the changeset viewer.