Changeset 4967
- Timestamp:
- 08/04/10 11:01:57 (15 years ago)
- Location:
- issm/trunk/src
- Files:
-
- 6 added
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Container/Inputs.cpp
r4943 r4967 426 426 } 427 427 /*}}}*/ 428 /*FUNCTION Inputs::AXPY{{{1*/ 429 void Inputs::AXPY(int YEnum, double scalar, int XEnum){ 430 431 Input* xinput=NULL; 432 Input* yinput=NULL; 433 434 /*Find x and y inputs: */ 435 xinput=(Input*)this->GetInput(XEnum); 436 yinput=(Input*)this->GetInput(YEnum); 437 438 /*some checks: */ 439 if(!xinput) ISSMERROR(" input %s could not be found!",EnumAsString(XEnum)); 440 if(!yinput) ISSMERROR(" input %s could not be found!",EnumAsString(YEnum)); 441 442 /*Apply AXPY: */ 443 yinput->AXPY(xinput,scalar); 444 } 445 /*}}}*/ -
issm/trunk/src/c/Container/Inputs.h
r4943 r4967 34 34 Input* GetInput(int enum_name); 35 35 Inputs* SpawnTriaInputs(int* indices); 36 void AXPY(int YEnum, double scalar, int XEnum); 36 37 37 38 void GetParameterAverage(double* pvalue, int enum_type); -
issm/trunk/src/c/Makefile.am
r4939 r4967 435 435 ./modules/InputControlConstrainx/InputControlConstrainx.h\ 436 436 ./modules/InputControlConstrainx/InputControlConstrainx.cpp\ 437 ./modules/InputControlUpdatex/InputControlUpdatex.h\ 438 ./modules/InputControlUpdatex/InputControlUpdatex.cpp\ 437 439 ./modules/SurfaceAreax/SurfaceAreax.h\ 438 440 ./modules/SurfaceAreax/SurfaceAreax.cpp\ … … 978 980 ./modules/InputControlConstrainx/InputControlConstrainx.h\ 979 981 ./modules/InputControlConstrainx/InputControlConstrainx.cpp\ 982 ./modules/InputControlUpdatex/InputControlUpdatex.h\ 983 ./modules/InputControlUpdatex/InputControlUpdatex.cpp\ 980 984 ./modules/SurfaceAreax/SurfaceAreax.h\ 981 985 ./modules/SurfaceAreax/SurfaceAreax.cpp\ -
issm/trunk/src/c/modules/Gradjx/Gradjx.cpp
r4592 r4967 15 15 int dim; 16 16 int numberofvertices; 17 bool extrude;18 17 Vec gradient = NULL; 19 18 20 19 /*retrieve some parameters: */ 21 20 parameters->FindParam(&dim,DimEnum); 22 parameters->FindParam(&extrude,ExtrudeParamEnum);23 21 numberofvertices=vertices->NumberOfVertices(); 24 22 … … 36 34 VecAssemblyEnd(gradient); 37 35 38 /*Extrude if needed: */39 if(dim==3 && extrude) VecExtrudex(gradient, elements,nodes, vertices, loads, materials, parameters,0);40 41 36 /*Assign output pointers: */ 42 37 *pgradient=gradient; -
issm/trunk/src/c/modules/InputDuplicatex/InputDuplicatex.cpp
r4572 r4967 12 12 13 13 /*Go through elemnets, and ask to reinitialie the input: */ 14 15 14 int i; 16 15 for(i=0;i<elements->Size();i++){ … … 18 17 element->InputDuplicate(original_enum,new_enum); 19 18 } 19 for(i=0;i<materials->Size();i++){ 20 Material* material=(Material*)materials->GetObjectByOffset(i); 21 //material->InputDuplicate(original_enum,new_enum); 22 } 20 23 21 24 } -
issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
r4919 r4967 40 40 IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); 41 41 IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); 42 if(iomodel->control_analysis && strcmp(iomodel->control_type,"rheology_B")==0){ 43 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel 44 } 42 45 43 46 /*Create elements and materials: */ … … 61 64 xfree((void**)&iomodel->rheology_B); 62 65 xfree((void**)&iomodel->rheology_n); 66 xfree((void**)&iomodel->control_parameter); 63 67 64 68 /*Add new constrant material property tgo materials, at the end: */ -
issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
r4578 r4967 47 47 IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); 48 48 IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); 49 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel 49 if(strcmp(iomodel->control_type,"drag_coefficient")==0){ 50 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel 51 } 50 52 } 51 53 -
issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/UpdateElementsDiagnosticStokes.cpp
r4578 r4967 48 48 IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); 49 49 IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); 50 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel 50 if(strcmp(iomodel->control_type,"drag_coefficient")==0){ 51 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel 52 } 51 53 } 52 54 -
issm/trunk/src/c/modules/modules.h
r4939 r4967 28 28 #include "./InputAXPYx/InputAXPYx.h" 29 29 #include "./InputControlConstrainx/InputControlConstrainx.h" 30 #include "./InputControlUpdatex/InputControlUpdatex.h" 30 31 #include "./InputConvergencex/InputConvergencex.h" 31 32 #include "./InputDuplicatex/InputDuplicatex.h" -
issm/trunk/src/c/objects/Elements/Element.h
r4942 r4967 72 72 virtual void GetVectorFromInputs(Vec vector,int NameEnum)=0; 73 73 virtual void InputAXPY(int YEnum, double scalar, int XEnum)=0; 74 virtual void InputControlUpdate(double scalar,bool save_parameter)=0; 74 75 virtual void InputControlConstrain(int control_type,double cm_min, double cm_max)=0; 75 76 virtual bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r4945 r4967 1026 1026 } 1027 1027 /*}}}*/ 1028 /*FUNCTION Penta::InputControlUpdate{{{1*/ 1029 void Penta::InputControlUpdate(double scalar,bool save_parameter){ 1030 1031 /*Intermediary*/ 1032 Input* input=NULL; 1033 double cm_min,cm_max; 1034 int control_type; 1035 1036 /*retrieve some parameters: */ 1037 this->parameters->FindParam(&cm_min,CmMinEnum); 1038 this->parameters->FindParam(&cm_max,CmMaxEnum); 1039 this->parameters->FindParam(&control_type,ControlTypeEnum); 1040 1041 1042 /*Rheology*/ 1043 if(control_type==RheologyB2dEnum){ 1044 1045 /*First, get revert to previous parameter value (kept in ControlParameter input)*/ 1046 matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBEnum); 1047 1048 /*Now, update using the gradient new = old + scalar * gradient*/ 1049 //matice->inputs->AXPY(RheologyB2dEnum,scalar,GradientEnum); 1050 // For now: Gradient is in element (TO BE CHANGED) and parameter in matice 1051 Input* input_B =(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input_B); 1052 Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad); 1053 input_B->AXPY(input_Grad,scalar); 1054 1055 /*Constrain constrain input*/ 1056 input=(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input); 1057 input->Constrain(cm_min,cm_max); 1058 1059 /*Finally: save input if requested*/ 1060 if (save_parameter) matice->inputs->DuplicateInput(RheologyBEnum,ControlParameterEnum); 1061 1062 } 1063 else if(control_type==DragCoefficientEnum){ 1064 1065 /*First, get revert to previous parameter value (kept in ControlParameter input)*/ 1066 this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum); 1067 1068 /*Now, update using the gradient new = old + scalar * gradient*/ 1069 this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum); 1070 1071 /*Constrain input*/ 1072 input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input); 1073 input->Constrain(cm_min,cm_max); 1074 1075 /*Finally: save input if requested*/ 1076 if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum); 1077 1078 } 1079 else{ 1080 ISSMERROR("control type %s not implemented yet",EnumAsString(control_type)); 1081 } 1082 } 1083 /*}}}*/ 1028 1084 /*FUNCTION Penta::InputDepthAverageAtBase{{{1*/ 1029 1085 void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ … … 1156 1212 1157 1213 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ 1158 for (i=0;i<this->inputs->Size();i++){1159 input= (Input*)this->inputs->GetObjectByOffset(i);1160 if (input->EnumType()==enum_type){1161 found=true; break;1162 }1163 } 1164 if (! found) ISSMERROR("Input %s not found in penta->inputs",EnumAsString(enum_type));1214 if (enum_type==RheologyB2dEnum){ 1215 input=this->matice->inputs->GetInput(enum_type); 1216 } 1217 else{ 1218 input=this->inputs->GetInput(enum_type); 1219 } 1220 if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type)); 1165 1221 1166 1222 /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result … … 1748 1804 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); 1749 1805 1750 }1751 if (iomodel->rheology_B) {1752 for(i=0;i<6;i++)nodeinputs[i]=iomodel->rheology_B[penta_vertex_ids[i]-1];1753 this->inputs->AddInput(new PentaVertexInput(RheologyBEnum,nodeinputs));1754 1806 } 1755 1807 if (iomodel->control_parameter) { … … 5392 5444 name==VzEnum || 5393 5445 name==TemperatureEnum || 5394 name==RheologyBEnum || 5395 name==RheologyNEnum || 5446 name==ControlParameterEnum || 5396 5447 name==FitEnum || 5397 5448 name==DragCoefficientEnum || -
issm/trunk/src/c/objects/Elements/Penta.h
r4942 r4967 87 87 void GradjDrag(Vec gradient); 88 88 void InputAXPY(int YEnum, double scalar, int XEnum); 89 void InputControlUpdate(double scalar,bool save_parameter); 89 90 void InputControlConstrain(int control_type,double cm_min, double cm_max); 90 91 bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r4943 r4967 416 416 417 417 /*update input*/ 418 this->inputs->AddInput(new TriaVertexInput(name,values)); 418 if (name==RheologyB2dEnum || name==RheologyBEnum){ 419 matice->inputs->AddInput(new TriaVertexInput(name,values)); 420 } 421 else{ 422 this->inputs->AddInput(new TriaVertexInput(name,values)); 423 } 419 424 return; 420 425 … … 628 633 } 629 634 else if (control_type==RheologyB2dEnum){ 630 inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);635 matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum); 631 636 Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; 632 637 } … … 965 970 adjointx_input=inputs->GetInput(AdjointxEnum); 966 971 adjointy_input=inputs->GetInput(AdjointyEnum); 967 rheologyb_input= inputs->GetInput(RheologyB2dEnum);972 rheologyb_input=matice->inputs->GetInput(RheologyB2dEnum); 968 973 969 974 /* Start looping on the number of gaussian points: */ … … 1229 1234 yinput->AXPY(xinput,scalar); 1230 1235 1231 /*Move input to Material if required (needed if control method) TO BE IMPROVED*/ 1232 if (YEnum==RheologyBEnum || YEnum==RheologyB2dEnum){ 1233 this->matice->inputs->AddInput((Input*)yinput->copy()); 1236 } 1237 /*}}}*/ 1238 /*FUNCTION Tria::InputControlUpdate{{{1*/ 1239 void Tria::InputControlUpdate(double scalar,bool save_parameter){ 1240 1241 /*Intermediary*/ 1242 Input* input=NULL; 1243 double cm_min,cm_max; 1244 int control_type; 1245 1246 /*retrieve some parameters: */ 1247 this->parameters->FindParam(&cm_min,CmMinEnum); 1248 this->parameters->FindParam(&cm_max,CmMaxEnum); 1249 this->parameters->FindParam(&control_type,ControlTypeEnum); 1250 1251 /*Rheology*/ 1252 if(control_type==RheologyB2dEnum){ 1253 1254 /*First, get revert to previous parameter value (kept in ControlParameter input)*/ 1255 matice->inputs->DuplicateInput(ControlParameterEnum,RheologyB2dEnum); 1256 1257 /*Now, update using the gradient new = old + scalar * gradient*/ 1258 //matice->inputs->AXPY(RheologyB2dEnum,scalar,GradientEnum); 1259 // For now: Gradient is in element (TO BE CHANGED) and parameter in matice 1260 Input* input_B =(Input*)matice->inputs->GetInput(RheologyB2dEnum); ISSMASSERT(input_B); 1261 Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad); 1262 input_B->AXPY(input_Grad,scalar); 1263 1264 /*Constrain constrain input*/ 1265 input=(Input*)matice->inputs->GetInput(RheologyB2dEnum); ISSMASSERT(input); 1266 input->Constrain(cm_min,cm_max); 1267 1268 /*Finally: save input if requested*/ 1269 if (save_parameter) matice->inputs->DuplicateInput(RheologyB2dEnum,ControlParameterEnum); 1270 1271 } 1272 else if(control_type==DragCoefficientEnum){ 1273 1274 /*First, get revert to previous parameter value (kept in ControlParameter input)*/ 1275 this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum); 1276 1277 /*Now, update using the gradient new = old + scalar * gradient*/ 1278 this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum); 1279 1280 /*Constrain input*/ 1281 input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input); 1282 input->Constrain(cm_min,cm_max); 1283 1284 /*Finally: save input if requested*/ 1285 if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum); 1286 1287 } 1288 else{ 1289 ISSMERROR("control type %s not implemented yet",EnumAsString(control_type)); 1234 1290 } 1235 1291 … … 1319 1375 1320 1376 /*Call inputs method*/ 1321 i nputs->DuplicateInput(original_enum,new_enum);1377 if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); 1322 1378 1323 1379 } … … 1344 1400 1345 1401 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ 1346 for (i=0;i<this->inputs->Size();i++){ 1347 input=(Input*)this->inputs->GetObjectByOffset(i); 1348 if (input->EnumType()==enum_type){ 1349 found=true; 1350 break; 1351 } 1352 } 1353 if (!found) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type)); 1402 if (enum_type==RheologyB2dEnum){ 1403 input=this->matice->inputs->GetInput(enum_type); 1404 } 1405 else{ 1406 input=this->inputs->GetInput(enum_type); 1407 } 1408 if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type)); 1354 1409 1355 1410 /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result … … 2122 2177 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); 2123 2178 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); 2124 }2125 if (iomodel->rheology_B) {2126 for(i=0;i<3;i++)nodeinputs[i]=iomodel->rheology_B[tria_vertex_ids[i]-1];2127 this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs));2128 2179 } 2129 2180 if (iomodel->control_parameter) { … … 6044 6095 name==MeltingRateEnum || 6045 6096 name==AccumulationRateEnum || 6097 name==ControlParameterEnum || 6046 6098 name==VxEnum || 6047 6099 name==VyEnum || 6048 name==RheologyBEnum ||6049 name==RheologyB2dEnum ||6050 name==RheologyNEnum ||6051 6100 name==FitEnum || 6052 6101 name==DragCoefficientEnum || -
issm/trunk/src/c/objects/Elements/Tria.h
r4942 r4967 83 83 void GradjDrag(Vec gradient); 84 84 void InputAXPY(int YEnum, double scalar, int XEnum); 85 void InputControlUpdate(double scalar,bool save_parameter); 85 86 void InputControlConstrain(int control_type,double cm_min, double cm_max); 86 87 bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); -
issm/trunk/src/c/objects/Materials/Material.h
r4248 r4967 18 18 public: 19 19 virtual ~Material(){}; 20 21 /*Numerics*/ 22 virtual void InputDuplicate(int original_enum,int new_enum)=0; 23 20 24 }; 21 25 #endif -
issm/trunk/src/c/objects/Materials/Matice.cpp
r4931 r4967 55 55 this->inputs->AddInput(new TriaVertexInput(RheologyNEnum,nodeinputs)); 56 56 } 57 58 /*Get control_parameter*/ 59 if (iomodel->control_parameter) { 60 for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->control_parameter[int(iomodel->elements[num_vertices*index+i]-1)]; 61 this->inputs->AddInput(new TriaVertexInput(ControlParameterEnum,nodeinputs)); 62 } 57 63 } 58 64 … … 74 80 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->rheology_n[index]; 75 81 this->inputs->AddInput(new PentaVertexInput(RheologyNEnum,nodeinputs)); 82 } 83 84 /*Get control_parameter*/ 85 if (iomodel->control_parameter) { 86 for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->control_parameter[int(iomodel->elements[num_vertices*index+i]-1)]; 87 this->inputs->AddInput(new PentaVertexInput(ControlParameterEnum,nodeinputs)); 76 88 } 77 89 } … … 531 543 } 532 544 /*}}}*/ 545 /*FUNCTION Matice::InputDuplicate{{{1*/ 546 void Matice::InputDuplicate(int original_enum,int new_enum){ 547 548 /*Call inputs method*/ 549 printf("%s\n",EnumAsString(original_enum)); 550 if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); 551 552 } 553 /*}}}*/ 533 554 /*FUNCTION Matice::InputUpdateFromVector(double* vector, int name, int type) {{{1*/ 534 555 void Matice::InputUpdateFromVector(double* vector, int name, int type){ 535 /*Nothing updated yet*/ 556 557 /*Intermediaries*/ 558 Element *element = NULL; 559 560 /*Recover element*/ 561 element=(Element*)helement->delivers(); 562 563 /*Check that name is an element input*/ 564 if (!IsInput(name)) return; 565 566 switch(type){ 567 568 case VertexEnum: 569 570 switch(element->Enum()){ 571 572 case TriaEnum: 573 double values[3]; 574 for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexDof()]; 575 this->inputs->AddInput(new TriaVertexInput(name,values)); 576 return; 577 578 default: ISSMERROR("element %s not implemented yet",EnumAsString(element->Enum())); 579 } 580 default: ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); 581 } 536 582 } 537 583 /*}}}*/ … … 566 612 } 567 613 /*}}}*/ 614 /*FUNCTION Matice::IsInput{{{1*/ 615 bool Matice::IsInput(int name){ 616 if ( 617 name==RheologyBEnum || 618 name==RheologyB2dEnum || 619 name==RheologyNEnum || 620 name==ControlParameterEnum 621 ){ 622 return true; 623 } 624 else return false; 625 } 626 /*}}}*/ -
issm/trunk/src/c/objects/Materials/Matice.h
r4931 r4967 50 50 void InputUpdateFromSolution(double* solution); 51 51 /*}}}*/ 52 /*Material virtual functions resolution: {{{1*/ 53 void InputDuplicate(int original_enum,int new_enum); 54 /*}}}*/ 52 55 /*Matice Numerics: {{{1*/ 53 56 void Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); … … 60 63 double GetB2d(); 61 64 double GetN(); 65 bool IsInput(int name); 62 66 /*}}}*/ 63 67 }; -
issm/trunk/src/c/objects/Materials/Matpar.h
r4966 r4967 54 54 void InputUpdateFromSolution(double* solution); 55 55 /*}}}*/ 56 /*Material virtual functions resolution: {{{1*/ 57 void InputDuplicate(int original_enum,int new_enum){ISSMERROR("not implemented yet");}; 58 /*}}}*/ 56 59 /*Numerics: {{{1*/ 57 60 double GetG(); -
issm/trunk/src/c/solutions/control_core.cpp
r4894 r4967 23 23 double eps_cm; 24 24 double tolx; 25 double cm_min;26 double cm_max;27 25 bool cm_gradient; 28 26 int dim; … … 54 52 femmodel->parameters->FindParam(&eps_cm,EpsCmEnum); 55 53 femmodel->parameters->FindParam(&tolx,TolXEnum); 56 femmodel->parameters->FindParam(&cm_min,CmMinEnum);57 femmodel->parameters->FindParam(&cm_max,CmMaxEnum);58 54 femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum); 59 55 femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); … … 90 86 BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs); 91 87 92 _printf_("%s\n"," updating parameter using optimized search scalar..."); 93 InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type); 94 InputAXPYx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],GradientEnum); 95 96 _printf_("%s\n"," constraining the new distribution..."); 97 InputControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max); 98 99 _printf_("%s\n"," save new parameter..."); 100 InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum); 101 102 _printf_("%s%i%s%g\n"," value of misfit J after optimization #",n+1,": ",J[n]); 88 _printf_("%s\n"," updating parameter using optimized search scalar..."); //true means update parameter and copy it onto ControlParameter input 89 InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true); 103 90 104 91 if(controlconvergence(J,fit,eps_cm,n)) break; -
issm/trunk/src/c/solutions/objectivefunctionC.cpp
r4839 r4967 32 32 double* optscal=NULL; 33 33 double* fit=NULL; 34 double cm_min;35 double cm_max;36 34 int control_type; 37 35 int analysis_type; … … 48 46 femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); 49 47 femmodel->parameters->FindParam(&fit,NULL,FitEnum); 50 femmodel->parameters->FindParam(&cm_min,CmMinEnum);51 femmodel->parameters->FindParam(&cm_max,CmMaxEnum);52 48 femmodel->parameters->FindParam(&isstokes,IsStokesEnum); 53 49 femmodel->parameters->FindParam(&control_type,ControlTypeEnum); … … 59 55 else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); 60 56 61 /*Use ControlParameterEnum input to reinitialize our input parameter: */ 62 InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type); 63 64 /*Use search scalar to shoot parameter in the gradient direction: */ 65 InputAXPYx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],GradientEnum); 66 67 /*Constrain:*/ 68 InputControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max); 57 /*update parameter according to scalar: */ //false means: do not copy updated parameter onto ControlParameter input 58 InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],false); 69 59 70 60 /*Run diagnostic with updated inputs: */ -
issm/trunk/src/m/solutions/control_core.m
r4895 r4967 19 19 eps_cm=femmodel.parameters.EpsCm; 20 20 tolx=femmodel.parameters.TolX; 21 cm_min=femmodel.parameters.CmMin;22 cm_max=femmodel.parameters.CmMax;23 21 cm_gradient=femmodel.parameters.CmGradient; 24 22 control_steady=femmodel.parameters.ControlSteady; … … 55 53 56 54 displaystring('\n%s',[' updating parameter using optimized search scalar:']); 57 scalar=search_scalar*optscal(n); 58 femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type); 59 [femmodel.elements femmodel.materials]=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,GradientEnum); 60 61 displaystring('\n%s',[' constraning the new distribution...']); 62 [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max); 63 64 displaystring('\n%s',[' save new parameter...']); 65 femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,ControlParameterEnum); 55 [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal(n),1); 66 56 67 57 disp([' value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]); -
issm/trunk/src/m/solutions/objectivefunctionC.m
r4804 r4967 10 10 analysis_type=femmodel.parameters.AnalysisType; 11 11 isstokes=femmodel.parameters.IsStokes; 12 cm_min=femmodel.parameters.CmMin;13 cm_max=femmodel.parameters.CmMax;14 12 15 13 %set current configuration … … 20 18 end 21 19 22 %Use ControlParameterEnum input to reinitialize our input parameter:23 femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);24 25 20 %Use search scalar to shoot parameter in the gradient direction: 26 scalar=search_scalar*optscal; 27 [femmodel.elements femmodel.materials]=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,GradientEnum); 28 29 %Constrain: 30 [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max); 21 [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal,0); 31 22 32 23 %Run diagnostic with updated inputs: -
issm/trunk/src/mex/Makefile.am
r4938 r4967 22 22 ElementConnectivity\ 23 23 InputControlConstrain \ 24 InputControlUpdate \ 24 25 InputConvergence\ 25 26 GetSolutionFromInputs\ … … 134 135 InputControlConstrain/InputControlConstrain.h 135 136 137 InputControlUpdate_SOURCES = InputControlUpdate/InputControlUpdate.cpp\ 138 InputControlUpdate/InputControlUpdate.h 139 136 140 InputConvergence_SOURCES = InputConvergence/InputConvergence.cpp\ 137 141 InputConvergence/InputConvergence.h
Note:
See TracChangeset
for help on using the changeset viewer.