Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 6237) +++ /issm/trunk/src/c/Makefile.am (revision 6238) @@ -551,4 +551,6 @@ ./modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp\ ./modules/ControlInputSetGradientx/ControlInputSetGradientx.h\ + ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\ + ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\ ./modules/SystemMatricesx/SystemMatricesx.cpp\ ./modules/SystemMatricesx/SystemMatricesx.h\ @@ -1115,4 +1117,6 @@ ./modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp\ ./modules/ControlInputSetGradientx/ControlInputSetGradientx.h\ + ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\ + ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\ ./modules/SystemMatricesx/SystemMatricesx.cpp\ ./modules/SystemMatricesx/SystemMatricesx.h\ Index: /issm/trunk/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp =================================================================== --- /issm/trunk/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp (revision 6238) +++ /issm/trunk/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp (revision 6238) @@ -0,0 +1,18 @@ +/*!\file ControlInputScaleGradientx + * \brief retrieve gradient from inputs in elements + */ + +#include "./ControlInputScaleGradientx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int NameEnum,double scale){ + + for(int i=0;iSize();i++){ + Element* element=(Element*)elements->GetObjectByOffset(i); + element->ControlInputScaleGradient(NameEnum,scale); + } + +} Index: /issm/trunk/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h =================================================================== --- /issm/trunk/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h (revision 6238) +++ /issm/trunk/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h (revision 6238) @@ -0,0 +1,12 @@ +/*!\file: ControlInputScaleGradientx.h + */ + +#ifndef _CONTROLINPUTSSCALEGRADIENTX_H +#define _CONTROLINPUTSSCALEGRADIENTX_H + +#include "../../objects/objects.h" +#include "../../Container/Container.h" + +void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int NameEnum,double scale); + +#endif Index: /issm/trunk/src/c/modules/Orthx/Orthx.cpp =================================================================== --- /issm/trunk/src/c/modules/Orthx/Orthx.cpp (revision 6237) +++ /issm/trunk/src/c/modules/Orthx/Orthx.cpp (revision 6238) @@ -25,14 +25,4 @@ } - /*scale to 1: gradient=gradient/max(abs(gradient))*/ - VecNorm(newgradj,NORM_INFINITY,&norm_new); - if (norm_new<=0){ - ISSMERROR("||∂J/∂α||∞ = 0 gradient is zero"); - } - if(isnan(norm_new)){ - ISSMERROR("input vector norm is NaN\n"); - } - VecScale(newgradj,1.0/norm_new); - /*Assign correct pointer*/ *pnewgradj=newgradj; Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 6237) +++ /issm/trunk/src/c/modules/modules.h (revision 6238) @@ -21,4 +21,5 @@ #include "./ControlInputGetGradientx/ControlInputGetGradientx.h" #include "./ControlInputSetGradientx/ControlInputSetGradientx.h" +#include "./ControlInputScaleGradientx/ControlInputScaleGradientx.h" #include "./CostFunctionx/CostFunctionx.h" #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h" Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 6237) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 6238) @@ -60,4 +60,5 @@ virtual void ControlInputGetGradient(Vec gradient,int enum_type)=0; virtual void ControlInputSetGradient(double* gradient,int enum_type)=0; + virtual void ControlInputScaleGradient(int enum_type, double scale)=0; virtual void ProcessResultsUnits(void)=0; virtual void MinVel(double* pminvel, bool process_units)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6237) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6238) @@ -816,4 +816,20 @@ ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); +}/*}}}*/ +/*FUNCTION Penta::ControlInputScaleGradient{{{1*/ +void Penta::ControlInputScaleGradient(int enum_type,double scale){ + + Input* input=NULL; + + if(enum_type==RheologyBbarEnum){ + input=(Input*)matice->inputs->GetInput(enum_type); + } + else{ + input=inputs->GetInput(enum_type); + } + if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); + + ((ControlInput*)input)->ScaleGradient(scale); }/*}}}*/ /*FUNCTION Penta::ControlInputSetGradient{{{1*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 6237) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 6238) @@ -95,4 +95,5 @@ void InputScale(int enum_type,double scale_factor); void ControlInputGetGradient(Vec gradient,int enum_type); + void ControlInputScaleGradient(int enum_type,double scale); void ControlInputSetGradient(double* gradient,int enum_type); void InputToResult(int enum_type,int step,double time); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6237) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6238) @@ -793,4 +793,20 @@ }/*}}}*/ +/*FUNCTION Tria::ControlInputScaleGradient{{{1*/ +void Tria::ControlInputScaleGradient(int enum_type,double scale){ + + Input* input=NULL; + + if(enum_type==RheologyBbarEnum){ + input=(Input*)matice->inputs->GetInput(enum_type); + } + else{ + input=inputs->GetInput(enum_type); + } + if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); + + ((ControlInput*)input)->ScaleGradient(scale); +}/*}}}*/ /*FUNCTION Tria::ControlInputSetGradient{{{1*/ void Tria::ControlInputSetGradient(double* gradient,int enum_type){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 6237) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 6238) @@ -96,4 +96,5 @@ void InputScale(int enum_type,double scale_factor); void ControlInputGetGradient(Vec gradient,int enum_type); + void ControlInputScaleGradient(int enum_type,double scale); void ControlInputSetGradient(double* gradient,int enum_type); void InputToResult(int enum_type,int step,double time); Index: /issm/trunk/src/c/objects/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/ControlInput.cpp (revision 6237) +++ /issm/trunk/src/c/objects/Inputs/ControlInput.cpp (revision 6238) @@ -266,4 +266,9 @@ if(gradient) gradient->GetVectorFromInputs(gradient_vec,doflist); }/*}}}*/ +/*FUNCTION ControlInput::ScaleGradient{{{1*/ +void ControlInput::ScaleGradient(double scaling_factor){ + if(!gradient) ISSMERROR("Gradient of ControlInput %s not found",EnumToString(enum_type)); + gradient->Scale(scaling_factor); +}/*}}}*/ /*FUNCTION ControlInput::SetGradient{{{1*/ void ControlInput::SetGradient(Input* gradient_in){ Index: /issm/trunk/src/c/objects/Inputs/ControlInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/ControlInput.h (revision 6237) +++ /issm/trunk/src/c/objects/Inputs/ControlInput.h (revision 6238) @@ -78,4 +78,5 @@ ElementResult* SpawnGradient(int step, double time); void GetGradient(Vec gradient_vec,int* doflist); + void ScaleGradient(double scale); void SetGradient(Input* gradient_in); void UpdateValue(double scalar); Index: /issm/trunk/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradient_core.cpp (revision 6237) +++ /issm/trunk/src/c/solutions/gradient_core.cpp (revision 6238) @@ -21,4 +21,5 @@ int *control_type = NULL; double *optscal_list = NULL; + double optscal,norm_grad; /*Intermediaries*/ @@ -34,4 +35,5 @@ femmodel->parameters->FindParam(&verbose,VerboseEnum); + /*Compute and norm gradient of all controls*/ for (int i=0;ielements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]); - Orthx(&new_gradient,gradient,old_gradient); - VecFree(&old_gradient); + Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient); } else{ _printf_("%s"," normalizing directions...\n"); - Orthx(&new_gradient,gradient,NULL); + Orthx(&new_gradient,gradient,NULL); VecFree(&gradient); } - VecFree(&gradient); - /*Scale gradient for current step and current parameter*/ - VecScale(new_gradient,optscal_list[num_controls*step+i]); + /*Get scaling factor of current control:*/ + VecNorm(new_gradient,NORM_INFINITY,&norm_grad); + if(norm_grad<=0) ISSMERROR("||∂J/∂α||∞ = 0 gradient norm of J with respect to %s is zero",EnumToString(control_type[i])); + if(isnan(norm_grad))ISSMERROR("||∂J/∂α||∞ = NaN gradient norm of J with respect to %s is NaN" ,EnumToString(control_type[i])); + if(i==0 || (optscal_list[num_controls*step+i]/norm_grad) elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type[i],new_gradient); - - /*Free ressources and return:*/ VecFree(&new_gradient); } + printf("norm grad = %g\n",norm_grad); + printf("optscal = %g\n",optscal); + + /*Scale Gradients*/ + for (int i=0;i elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type[i],optscal); /*Clean up and return*/ Index: /issm/trunk/src/m/solutions/gradient_core.m =================================================================== --- /issm/trunk/src/m/solutions/gradient_core.m (revision 6237) +++ /issm/trunk/src/m/solutions/gradient_core.m (revision 6238) @@ -25,5 +25,5 @@ control_type=femmodel.parameters.ControlType; control_steady=femmodel.parameters.ControlSteady; - optscal=femmodel.parameters.OptScal; + optscal_list=femmodel.parameters.OptScal; for i=1:num_controls, @@ -32,6 +32,5 @@ grad=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters); - - if control_steady; + if control_steady, femmodel=diagnostic_core(femmodel); end @@ -46,8 +45,16 @@ end - %Scale Gradient - new_gradient=optscal(step,i)*new_gradient; + %Get scaling factor of current control: + norm_grad=norm(new_gradient,inf); + if(norm_grad<=0), error(['||∂J/∂α||∞ = 0 gradient norm of J with respect to ' EnumToString(control_type(i)) ' is zero']); end + if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i)) ' is NaN' ]); end + if(i==1 | (optscal_list(step,i)/norm_grad) Configure(elements,loads, nodes,vertices, materials,parameters); +nodes-> Configure(elements,loads, nodes,vertices, materials,parameters); +loads-> Configure(elements, loads, nodes,vertices, materials,parameters); + +/*call "x" code layer*/ +ControlInputScaleGradientx(elements,nodes,vertices,loads, materials,parameters,control_type,scaling_factor); + +/*write output datasets: */ +WriteData(ELEMENTS,elements); +WriteData(NODES,nodes); +WriteData(VERTICES,vertices); +WriteData(LOADS,loads); +WriteData(MATERIALS,materials); +WriteParams(PARAMETERS,parameters); + +/*Free ressources: */ +delete elements; +delete nodes; +delete vertices; +delete loads; +delete materials; +delete parameters; + +/*end module: */ +MODULEEND(); +} + +void ControlInputScaleGradientUsage(void) +{ + _printf_("\n"); + _printf_(" usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,control_type,scaling_factor);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h =================================================================== --- /issm/trunk/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h (revision 6238) +++ /issm/trunk/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h (revision 6238) @@ -0,0 +1,43 @@ +/* + ControlInputScaleGradient.h +*/ + +#ifndef _CONTROLINPUTSCALEGRADIENT_H +#define _CONTROLINPUTSCALEGRADIENT_H + +/* local prototypes: */ +void ControlInputScaleGradientUsage(void); + +#include "../../c/modules/modules.h" +#include "../../c/Container/Container.h" +#include "../../c/shared/shared.h" +#include "../../c/EnumDefinitions/EnumDefinitions.h" + +#undef __FUNCT__ +#define __FUNCT__ "ControlInputScaleGradient" + +/* serial input macros: */ +#define ELEMENTSIN (mxArray*)prhs[0] +#define NODESIN (mxArray*)prhs[1] +#define VERTICESIN (mxArray*)prhs[2] +#define LOADSIN (mxArray*)prhs[3] +#define MATERIALSIN (mxArray*)prhs[4] +#define PARAMETERSIN (mxArray*)prhs[5] +#define CONTROLTYPE (mxArray*)prhs[6] +#define SCALE (mxArray*)prhs[7] + +/* serial output macros: */ +#define ELEMENTS (mxArray**)&plhs[0] +#define NODES (mxArray**)&plhs[1] +#define VERTICES (mxArray**)&plhs[2] +#define LOADS (mxArray**)&plhs[3] +#define MATERIALS (mxArray**)&plhs[4] +#define PARAMETERS (mxArray**)&plhs[5] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 6 +#undef NRHS +#define NRHS 8 + +#endif Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 6237) +++ /issm/trunk/src/mex/Makefile.am (revision 6238) @@ -15,4 +15,5 @@ ControlOptimization\ ControlInputGetGradient\ + ControlInputScaleGradient\ ControlInputSetGradient\ ContourToMesh \ @@ -133,4 +134,7 @@ ControlInputGetGradient/ControlInputGetGradient.h +ControlInputScaleGradient_SOURCES = ControlInputScaleGradient/ControlInputScaleGradient.cpp\ + ControlInputScaleGradient/ControlInputScaleGradient.h + ControlInputSetGradient_SOURCES = ControlInputSetGradient/ControlInputSetGradient.cpp\ ControlInputSetGradient/ControlInputSetGradient.h