Index: /issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.cpp (revision 11313) @@ -9,12 +9,23 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name){ +void ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){ - Vec gradient=NULL; - gradient=NewVec(vertices->NumberOfVertices()); + /*Intermediaries*/ + int num_controls; + int *control_type = NULL; + Vec gradient=NULL; - for(int i=0;iSize();i++){ - Element* element=(Element*)elements->GetObjectByOffset(i); - element->ControlInputGetGradient(gradient,name); + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + /*Allocate and populate gradient*/ + gradient=NewVec(num_controls*vertices->NumberOfVertices()); + + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ControlInputGetGradient(gradient,control_type[i],i); + } } @@ -22,5 +33,6 @@ VecAssemblyEnd(gradient); - /*Assign output pointers:*/ + /*Clean up and return*/ + xfree((void**)&control_type); *pgradient=gradient; } Index: /issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.h (revision 11312) +++ /issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.h (revision 11313) @@ -8,5 +8,5 @@ #include "../../Container/Container.h" -void ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int name); +void ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters); #endif Index: /issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp (revision 11313) @@ -9,10 +9,33 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name,double scale){ +void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* norm_list,int step){ - for(int i=0;iSize();i++){ - Element* element=(Element*)elements->GetObjectByOffset(i); - element->ControlInputScaleGradient(name,scale); + /*Intermediaries*/ + int i,j,num_controls; + int *control_type = NULL; + double *scalar_list = NULL; + double scalar; + + + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum); + + /*Compute scaling factor*/ + scalar=scalar_list[0]/norm_list[0]; + for(i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ControlInputScaleGradient(control_type[i],scalar); + } + } + + /*Clean up and return*/ + xfree((void**)&control_type); + xfree((void**)&scalar_list); } Index: /issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h (revision 11312) +++ /issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h (revision 11313) @@ -8,5 +8,5 @@ #include "../../Container/Container.h" -void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int name,double scale); +void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,double* norm_list,int step); #endif Index: /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp (revision 11313) @@ -9,13 +9,26 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name,double* gradient){ +void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* gradient){ - for(int i=0;iSize();i++){ - Element* element=(Element*)elements->GetObjectByOffset(i); - element->ControlInputSetGradient(gradient,name); + /*Intermediaries*/ + int num_controls; + int *control_type = NULL; + + /*Retrieve some parameters*/ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + for(int i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ControlInputSetGradient(gradient,control_type[i],i); + } } + /*Clean up and return*/ + xfree((void**)&control_type); + } -void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name,Vec gradient){ +void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vec gradient){ /*Serialize gradient*/ @@ -23,5 +36,5 @@ VecToMPISerial(&serial_gradient,gradient); - ControlInputSetGradientx(elements,nodes,vertices, loads, materials, parameters, name,serial_gradient); + ControlInputSetGradientx(elements,nodes,vertices, loads, materials, parameters,serial_gradient); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.h (revision 11312) +++ /issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.h (revision 11313) @@ -8,6 +8,6 @@ #include "../../Container/Container.h" -void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int name,double* gradient); -void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int name,Vec gradient); +void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,double* gradient); +void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vec gradient); #endif Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 11313) @@ -10,35 +10,56 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void Gradjx( Vec* pgradient,double* pnorm_grad, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int control_index){ +void Gradjx(Vec* pgradient,double** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){ - int i; - int numberofvertices; - double norm_grad; - int *control_type = NULL; - Vec gradient = NULL; + int i,j,numberofvertices; + int num_controls; + double norm_inf; + double *norm_list = NULL; + int *control_type = NULL; + Vec gradient = NULL; + Vec *gradient_list = NULL; /*retrieve some parameters: */ + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); _assert_(num_controls); parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); numberofvertices=vertices->NumberOfVertices(); - /*Allocate gradient: */ - gradient=NewVec(numberofvertices); + /*Allocate gradient_list */ + gradient_list = (Vec*)xmalloc(num_controls*sizeof(Vec)); + norm_list = (double*)xmalloc(num_controls*sizeof(double)); + for(i=0;iSize();i++){ - Element* element=(Element*)elements->GetObjectByOffset(i); - element->Gradj(gradient,control_type[control_index]); + /*Compute all gradient_list*/ + for(i=0;iSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->Gradj(gradient_list[i],control_type[i],i); + } + + VecAssemblyBegin(gradient_list[i]); + VecAssemblyEnd(gradient_list[i]); + + VecNorm(gradient_list[i],NORM_INFINITY,&norm_list[i]); } - /*Assemble vector: */ - VecAssemblyBegin(gradient); - VecAssemblyEnd(gradient); + /*Add all gradient_list together*/ + for(i=0;iObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type)); - this->GetDofList1(&doflist1[0]); + GradientIndexing(&doflist1[0],control_index); ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); @@ -4350,5 +4350,5 @@ }/*}}}*/ /*FUNCTION Penta::ControlInputSetGradient{{{1*/ -void Penta::ControlInputSetGradient(double* gradient,int enum_type){ +void Penta::ControlInputSetGradient(double* gradient,int enum_type,int control_index){ int doflist1[NUMVERTICES]; @@ -4366,5 +4366,5 @@ if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type)); - this->GetDofList1(&doflist1[0]); + GradientIndexing(&doflist1[0],control_index); for(int i=0;iFindParam(&num_controls,InversionNumControlParametersEnum); + + /*get gradient indices*/ + for(int i=0;inodes[i]->GetVertexDof() + control_index; + } + +} +/*}}}*/ /*FUNCTION Penta::Gradj {{{1*/ -void Penta::Gradj(Vec gradient,int control_type){ +void Penta::Gradj(Vec gradient,int control_type,int control_index){ /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ @@ -4451,11 +4465,11 @@ switch(approximation){ case MacAyealApproximationEnum: - GradjDragMacAyeal(gradient); + GradjDragMacAyeal(gradient,control_index); break; case PattynApproximationEnum: - GradjDragPattyn(gradient); + GradjDragPattyn(gradient,control_index); break; case StokesApproximationEnum: - GradjDragStokes(gradient); + GradjDragStokes(gradient,control_index); break; case NoneApproximationEnum: @@ -4471,11 +4485,11 @@ switch(approximation){ case MacAyealApproximationEnum: - GradjBbarMacAyeal(gradient); + GradjBbarMacAyeal(gradient,control_index); break; case PattynApproximationEnum: - GradjBbarPattyn(gradient); + GradjBbarPattyn(gradient,control_index); break; case StokesApproximationEnum: - GradjBbarStokes(gradient); + GradjBbarStokes(gradient,control_index); break; case NoneApproximationEnum: @@ -4510,5 +4524,5 @@ if (!IsOnBed()) return; tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjDragGradient(gradient,resp); + tria->GradjDragGradient(gradient,resp,control_index); delete tria->matice; delete tria; break; @@ -4516,5 +4530,5 @@ if (!IsOnBed()) return; tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjBGradient(gradient,resp); + tria->GradjBGradient(gradient,resp,control_index); delete tria->matice; delete tria; break; @@ -4526,5 +4540,5 @@ /*}}}*/ /*FUNCTION Penta::GradjDragMacAyeal {{{1*/ -void Penta::GradjDragMacAyeal(Vec gradient){ +void Penta::GradjDragMacAyeal(Vec gradient,int control_index){ /*Gradient is 0 if on shelf or not on bed*/ @@ -4533,10 +4547,10 @@ /*Spawn tria*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjDragMacAyeal(gradient); + tria->GradjDragMacAyeal(gradient,control_index); delete tria->matice; delete tria; } /*}}}*/ /*FUNCTION Penta::GradjDragPattyn {{{1*/ -void Penta::GradjDragPattyn(Vec gradient){ +void Penta::GradjDragPattyn(Vec gradient,int control_index){ int i,j,ig; @@ -4559,7 +4573,7 @@ /*Retrieve all inputs and parameters*/ parameters->FindParam(&analysis_type,AnalysisTypeEnum); + GradientIndexing(&doflist1[0],control_index); GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); for(i=0;iGetInput(AdjointxEnum); _assert_(adjointx_input); Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); @@ -4609,5 +4623,5 @@ /*}}}*/ /*FUNCTION Penta::GradjDragStokes {{{1*/ -void Penta::GradjDragStokes(Vec gradient){ +void Penta::GradjDragStokes(Vec gradient,int control_index){ int i,j,ig; @@ -4634,12 +4648,12 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); for(i=0;iGetInput(FrictionCoefficientEnum); _assert_(drag_input); - Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input =inputs->GetInput(VzEnum); _assert_(vz_input); - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* adjointz_input=inputs->GetInput(AdjointzEnum); _assert_(adjointz_input); + Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input =inputs->GetInput(VzEnum); _assert_(vz_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* adjointz_input=inputs->GetInput(AdjointzEnum); _assert_(adjointz_input); /*Build frictoin element, needed later: */ @@ -4701,5 +4715,5 @@ /*}}}*/ /*FUNCTION Penta::GradjBbarMacAyeal {{{1*/ -void Penta::GradjBbarMacAyeal(Vec gradient){ +void Penta::GradjBbarMacAyeal(Vec gradient,int control_index){ /*This element should be collapsed into a tria element at its base*/ @@ -4711,5 +4725,5 @@ /*Collapse element to the base*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). - tria->GradjBMacAyeal(gradient); + tria->GradjBMacAyeal(gradient,control_index); delete tria->matice; delete tria; @@ -4719,5 +4733,5 @@ } /*}}}*/ /*FUNCTION Penta::GradjBbarPattyn {{{1*/ -void Penta::GradjBbarPattyn(Vec gradient){ +void Penta::GradjBbarPattyn(Vec gradient,int control_index){ /*Gradient is computed on bed only (Bbar)*/ @@ -4729,5 +4743,5 @@ /*Collapse element to the base*/ Tria* tria=(Tria*)SpawnTria(0,1,2); - tria->GradjBMacAyeal(gradient); //We use MacAyeal as an estimate for now + tria->GradjBMacAyeal(gradient,control_index); //We use MacAyeal as an estimate for now delete tria->matice; delete tria; @@ -4736,5 +4750,5 @@ } /*}}}*/ /*FUNCTION Penta::GradjBbarStokes {{{1*/ -void Penta::GradjBbarStokes(Vec gradient){ +void Penta::GradjBbarStokes(Vec gradient,int control_index){ /*Gradient is computed on bed only (Bbar)*/ @@ -4746,5 +4760,5 @@ /*Collapse element to the base*/ Tria* tria=(Tria*)SpawnTria(0,1,2); - tria->GradjBMacAyeal(gradient); //We use MacAyeal as an estimate for now + tria->GradjBMacAyeal(gradient,control_index); //We use MacAyeal as an estimate for now delete tria->matice; delete tria; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11312) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11313) @@ -141,14 +141,15 @@ #ifdef _HAVE_CONTROL_ double DragCoefficientAbsGradient(bool process_units,int weight_index); - void Gradj(Vec gradient,int control_type); - void GradjDragMacAyeal(Vec gradient); - void GradjDragPattyn(Vec gradient); - void GradjDragStokes(Vec gradient); - void GradjBbarMacAyeal(Vec gradient); - void GradjBbarPattyn(Vec gradient); - void GradjBbarStokes(Vec gradient); - void ControlInputGetGradient(Vec gradient,int enum_type); + void GradientIndexing(int* indexing,int control_index); + void Gradj(Vec gradient,int control_type,int control_index); + void GradjDragMacAyeal(Vec gradient,int control_index); + void GradjDragPattyn(Vec gradient,int control_index); + void GradjDragStokes(Vec gradient,int control_index); + void GradjBbarMacAyeal(Vec gradient,int control_index); + void GradjBbarPattyn(Vec gradient,int control_index); + void GradjBbarStokes(Vec gradient,int control_index); + void ControlInputGetGradient(Vec gradient,int enum_type,int control_index); void ControlInputScaleGradient(int enum_type,double scale); - void ControlInputSetGradient(double* gradient,int enum_type); + void ControlInputSetGradient(double* gradient,int enum_type,int control_index); double RheologyBbarAbsGradient(bool process_units,int weight_index); double ThicknessAbsMisfit( bool process_units,int weight_index); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11313) @@ -3266,5 +3266,5 @@ /*}}}*/ /*FUNCTION Tria::ControlInputGetGradient{{{1*/ -void Tria::ControlInputGetGradient(Vec gradient,int enum_type){ +void Tria::ControlInputGetGradient(Vec gradient,int enum_type,int control_index){ int doflist1[NUMVERTICES]; @@ -3280,5 +3280,5 @@ if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type)); - this->GetDofList1(&doflist1[0]); + GradientIndexing(&doflist1[0],control_index); ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); @@ -3301,5 +3301,5 @@ }/*}}}*/ /*FUNCTION Tria::ControlInputSetGradient{{{1*/ -void Tria::ControlInputSetGradient(double* gradient,int enum_type){ +void Tria::ControlInputSetGradient(double* gradient,int enum_type,int control_index){ int doflist1[NUMVERTICES]; @@ -3317,5 +3317,5 @@ if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type)); - this->GetDofList1(&doflist1[0]); + GradientIndexing(&doflist1[0],control_index); for(int i=0;iinputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); @@ -3424,5 +3424,5 @@ /*}}}*/ /*FUNCTION Tria::GradjBMacAyeal{{{1*/ -void Tria::GradjBMacAyeal(Vec gradient){ +void Tria::GradjBMacAyeal(Vec gradient,int control_index){ /*Intermediaries*/ @@ -3439,12 +3439,12 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist[0]); + GradientIndexing(&doflist[0],control_index); /*Retrieve all inputs*/ - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); @@ -3481,5 +3481,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDragMacAyeal {{{1*/ -void Tria::GradjDragMacAyeal(Vec gradient){ +void Tria::GradjDragMacAyeal(Vec gradient,int control_index){ int i,ig; @@ -3502,5 +3502,5 @@ parameters->FindParam(&analysis_type,AnalysisTypeEnum); GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist1[0]); + GradientIndexing(&doflist1[0],control_index); /*Build frictoin element, needed later: */ @@ -3508,8 +3508,8 @@ /*Retrieve all inputs we will be needing: */ - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); @@ -3568,5 +3568,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDragGradient{{{1*/ -void Tria::GradjDragGradient(Vec gradient, int weight_index){ +void Tria::GradjDragGradient(Vec gradient, int weight_index,int control_index){ int i,ig; @@ -3582,5 +3582,5 @@ if(IsFloating())return; GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist1[0]); + GradientIndexing(&doflist1[0],control_index); Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); @@ -3612,5 +3612,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDhDtBalancedthickness{{{1*/ -void Tria::GradjDhDtBalancedthickness(Vec gradient){ +void Tria::GradjDhDtBalancedthickness(Vec gradient,int control_index){ /*Intermediaries*/ @@ -3619,7 +3619,6 @@ double gradient_g[NUMVERTICES]; - GetDofList1(&doflist1[0]); - /*Compute Gradient*/ + GradientIndexing(&doflist1[0],control_index); GetInputListOnVertices(&lambda[0],AdjointEnum); for(int i=0;iFindParam(&num_controls,InversionNumControlParametersEnum); + + /*get gradient indices*/ + for(int i=0;inodes[i]->GetVertexDof() + control_index; + } + } /*}}}*/ @@ -4225,4 +4238,5 @@ /*Clean up and return*/ delete gauss; + xfree((void**)&responses); return pe; } Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11312) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11313) @@ -140,16 +140,17 @@ #ifdef _HAVE_CONTROL_ double DragCoefficientAbsGradient(bool process_units,int weight_index); - void Gradj(Vec gradient,int control_type); - void GradjBGradient(Vec gradient,int weight_index); - void GradjBMacAyeal(Vec gradient); - void GradjDragMacAyeal(Vec gradient); - void GradjDragStokes(Vec gradient); - void GradjDragGradient(Vec gradient,int weight_index); - void GradjDhDtBalancedthickness(Vec gradient); - void GradjVxBalancedthickness(Vec gradient); - void GradjVyBalancedthickness(Vec gradient); - void ControlInputGetGradient(Vec gradient,int enum_type); + void GradientIndexing(int* indexing,int control_index); + void Gradj(Vec gradient,int control_type,int control_index); + void GradjBGradient(Vec gradient,int weight_index,int control_index); + void GradjBMacAyeal(Vec gradient,int control_index); + void GradjDragMacAyeal(Vec gradient,int control_index); + void GradjDragStokes(Vec gradient,int control_index); + void GradjDragGradient(Vec gradient,int weight_index,int control_index); + void GradjDhDtBalancedthickness(Vec gradient,int control_index); + void GradjVxBalancedthickness(Vec gradient,int control_index); + void GradjVyBalancedthickness(Vec gradient,int control_index); + void ControlInputGetGradient(Vec gradient,int enum_type,int control_index); void ControlInputScaleGradient(int enum_type,double scale); - void ControlInputSetGradient(double* gradient,int enum_type); + void ControlInputSetGradient(double* gradient,int enum_type,int control_index); double RheologyBbarAbsGradient(bool process_units,int weight_index); double ThicknessAbsMisfit( bool process_units,int weight_index); Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 11313) @@ -85,5 +85,5 @@ adjointdiagnostic_core(user->femmodel); //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum); - Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,0); + Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); VecScale(gradient,-1.); VecCopy(gradient,G); Index: /issm/trunk-jpl/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/gradient_core.cpp (revision 11312) +++ /issm/trunk-jpl/src/c/solutions/gradient_core.cpp (revision 11313) @@ -14,54 +14,38 @@ void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ - - /*parameters: */ - bool control_steady; - int num_controls; - int *control_type = NULL; - double *optscal_list = NULL; - double optscal=1.e+100,norm_grad; /*Intermediaries*/ - Vec new_gradient=NULL; - Vec gradient=NULL; - Vec old_gradient=NULL; + double norm_inf; + double *norm_list = NULL; + Vec new_gradient = NULL; + Vec gradient = NULL; + Vec old_gradient = NULL; - /*retrieve parameters:*/ - femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); - femmodel->parameters->FindParam(&optscal_list,NULL,NULL,InversionGradientScalingEnum); + /*Compute gradient*/ + _printf_(VerboseControl()," compute cost function gradient\n"); + Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters); - /*Compute and norm gradient of all controls*/ - for (int i=0;ielements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,i); - - if (orthogonalize){ - _printf_(VerboseControl()," orthogonalization\n"); - ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]); - Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient); - } - else{ - new_gradient=gradient; - } - - /*Get scaling factor of current control:*/ - VecNorm(new_gradient,NORM_INFINITY,&norm_grad); - if(norm_grad<=0) _error_("||∂J/∂α||∞ = 0 gradient norm of J with respect to %s is zero",EnumToStringx(control_type[i])); - if(isnan(norm_grad))_error_("||∂J/∂α||∞ = NaN gradient norm of J with respect to %s is NaN" ,EnumToStringx(control_type[i])); - optscal=min(optscal,optscal_list[num_controls*step+i]/norm_grad); - - /*plug back into inputs: */ - ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],new_gradient); - VecFree(&new_gradient); + if (orthogonalize){ + _printf_(VerboseControl()," orthogonalization\n"); + ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters); + Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient); + } + else{ + new_gradient=gradient; } + /*Check that gradient is clean*/ + VecNorm(new_gradient,NORM_INFINITY,&norm_inf); + if(norm_inf<=0) _error_("||∂J/∂α||∞ = 0 gradient norm is zero"); + if(isnan(norm_inf))_error_("||∂J/∂α||∞ = NaN gradient norm is NaN"); + + /*plug back into inputs: */ + ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,new_gradient); + VecFree(&new_gradient); + /*Scale Gradients*/ - for (int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],optscal); + ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,norm_list,step); /*Clean up and return*/ - xfree((void**)&control_type); - xfree((void**)&optscal_list); + xfree((void**)&norm_list); } Index: /issm/trunk-jpl/src/m/solutions/control_core.m =================================================================== --- /issm/trunk-jpl/src/m/solutions/control_core.m (revision 11312) +++ /issm/trunk-jpl/src/m/solutions/control_core.m (revision 11313) @@ -50,5 +50,5 @@ eval(['femmodel=' adjointcore '(femmodel);']); - femmodel=gradient_core(femmodel,n,search_scalar==0); + femmodel=gradient_core(femmodel,n-1,search_scalar==0); %Return gradient if asked Index: /issm/trunk-jpl/src/m/solutions/gradient_core.m =================================================================== --- /issm/trunk-jpl/src/m/solutions/gradient_core.m (revision 11312) +++ /issm/trunk-jpl/src/m/solutions/gradient_core.m (revision 11313) @@ -20,37 +20,23 @@ end -%recover parameters common to all solutions -num_controls=femmodel.parameters.InversionNumControlParameters; -control_type=femmodel.parameters.InversionControlParameters; -control_steady=femmodel.parameters.ControlSteady; -gradient_scaling_list=femmodel.parameters.InversionGradientScaling; -gradient_scaling=10^100; +issmprintf(VerboseControl,[' compute cost function gradient']); +[grad norm_list]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters); -for i=1:num_controls, - - issmprintf(VerboseControl,[' compute gradient of J with respect to %s'],EnumToString(control_type(i))); - [grad norm_grad]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,i-1); - - if orthogonalize, - issmprintf(VerboseControl,'%s',[' orthogonalization']); - old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i)); - new_gradient=Orth(grad,old_gradient); - else - new_gradient=grad; - end - - %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 - gradient_scaling=min(gradient_scaling,gradient_scaling_list(step,i)/norm_grad); - - %plug back into inputs: - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputSetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,control_type(i),new_gradient); - +if orthogonalize, + issmprintf(VerboseControl,'%s',[' orthogonalization']); + old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters); + new_gradient=Orth(grad,old_gradient); +else + new_gradient=grad; end +%Check that gradient is clean +norm_grad=norm(new_gradient,inf); +if(norm_grad<=0), error(['||∂J/∂α||∞ = 0 gradient norm is zero']); end +if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm is NaN' ]); end + +%plug back into inputs: +[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputSetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,new_gradient); + %Scale all gradients -for i=1:num_controls, - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputScaleGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,control_type(i),gradient_scaling); -end +[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputScaleGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,norm_list,step); Index: /issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.cpp =================================================================== --- /issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.cpp (revision 11312) +++ /issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.cpp (revision 11313) @@ -14,5 +14,4 @@ Materials *materials = NULL; Parameters *parameters = NULL; - int control_type; Vec gradient = NULL; @@ -32,5 +31,4 @@ FetchMatlabData((DataSet**)&materials,MATERIALS); FetchMatlabData(¶meters,PARAMETERS); - FetchMatlabData(&control_type,CONTROLTYPE); /*configure: */ @@ -41,5 +39,5 @@ /*!core code:*/ - ControlInputGetGradientx(&gradient,elements, nodes,vertices,loads, materials,parameters,control_type); + ControlInputGetGradientx(&gradient,elements, nodes,vertices,loads, materials,parameters); /*write output datasets: */ @@ -61,5 +59,5 @@ { _printf_(true,"\n"); - _printf_(true," usage: [gradient] = %s(elements,nodes,vertices,loads, materials,parameters,control_type);\n",__FUNCT__); + _printf_(true," usage: [gradient] = %s(elements,nodes,vertices,loads, materials,parameters);\n",__FUNCT__); _printf_(true,"\n"); } Index: /issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.h =================================================================== --- /issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.h (revision 11312) +++ /issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.h (revision 11313) @@ -23,5 +23,4 @@ #define MATERIALS (mxArray*)prhs[4] #define PARAMETERS (mxArray*)prhs[5] -#define CONTROLTYPE (mxArray*)prhs[6] /* serial output macros: */ @@ -32,5 +31,5 @@ #define NLHS 1 #undef NRHS -#define NRHS 7 +#define NRHS 6 #endif Index: /issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.cpp =================================================================== --- /issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.cpp (revision 11312) +++ /issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.cpp (revision 11313) @@ -8,4 +8,5 @@ /*input datasets: */ +int step; Elements *elements = NULL; Nodes *nodes = NULL; @@ -14,6 +15,5 @@ Materials *materials = NULL; Parameters *parameters = NULL; -int control_type; -double scaling_factor; +double *norm_list = NULL; /*Boot module: */ @@ -30,6 +30,6 @@ FetchMatlabData((DataSet**)&materials,MATERIALSIN); FetchMatlabData(¶meters,PARAMETERSIN); -FetchMatlabData(&control_type,CONTROLTYPE); -FetchMatlabData(&scaling_factor,SCALE); +FetchMatlabData(&norm_list,NULL,NULL,NORMLIST); +FetchMatlabData(&step,STEP); /*configure: */ @@ -39,5 +39,5 @@ /*call "x" code layer*/ -ControlInputScaleGradientx(elements,nodes,vertices,loads, materials,parameters,control_type,scaling_factor); +ControlInputScaleGradientx(elements,nodes,vertices,loads, materials,parameters,norm_list,step); /*write output datasets: */ @@ -64,5 +64,5 @@ { _printf_(true,"\n"); - _printf_(true," usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,control_type,scaling_factor);\n",__FUNCT__); + _printf_(true," usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,norm_list,step);\n",__FUNCT__); _printf_(true,"\n"); } Index: /issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h =================================================================== --- /issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h (revision 11312) +++ /issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h (revision 11313) @@ -24,6 +24,6 @@ #define MATERIALSIN (mxArray*)prhs[4] #define PARAMETERSIN (mxArray*)prhs[5] -#define CONTROLTYPE (mxArray*)prhs[6] -#define SCALE (mxArray*)prhs[7] +#define NORMLIST (mxArray*)prhs[6] +#define STEP (mxArray*)prhs[7] /* serial output macros: */ Index: /issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.cpp =================================================================== --- /issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.cpp (revision 11312) +++ /issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.cpp (revision 11313) @@ -14,5 +14,4 @@ Materials *materials = NULL; Parameters *parameters = NULL; -int control_type; double *gradient = NULL; @@ -30,5 +29,4 @@ FetchMatlabData((DataSet**)&materials,MATERIALSIN); FetchMatlabData(¶meters,PARAMETERSIN); -FetchMatlabData(&control_type,CONTROLTYPE); FetchMatlabData(&gradient,NULL,GRADIENT); @@ -39,5 +37,5 @@ /*call "x" code layer*/ -ControlInputSetGradientx(elements,nodes,vertices,loads, materials,parameters,control_type,gradient); +ControlInputSetGradientx(elements,nodes,vertices,loads, materials,parameters,gradient); /*write output datasets: */ @@ -64,5 +62,5 @@ { _printf_(true,"\n"); - _printf_(true," usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,control_type,gradient);\n",__FUNCT__); + _printf_(true," usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,gradient);\n",__FUNCT__); _printf_(true,"\n"); } Index: /issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.h =================================================================== --- /issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.h (revision 11312) +++ /issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.h (revision 11313) @@ -24,6 +24,5 @@ #define MATERIALSIN (mxArray*)prhs[4] #define PARAMETERSIN (mxArray*)prhs[5] -#define CONTROLTYPE (mxArray*)prhs[6] -#define GRADIENT (mxArray*)prhs[7] +#define GRADIENT (mxArray*)prhs[6] /* serial output macros: */ @@ -39,5 +38,5 @@ #define NLHS 6 #undef NRHS -#define NRHS 8 +#define NRHS 7 #endif Index: /issm/trunk-jpl/src/mex/Gradj/Gradj.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Gradj/Gradj.cpp (revision 11312) +++ /issm/trunk-jpl/src/mex/Gradj/Gradj.cpp (revision 11313) @@ -8,6 +8,6 @@ /*input datasets: */ - int control_index; - double grad_norm; + int control_index,num_controls; + double *norm_list = NULL; Elements *elements = NULL; Nodes *nodes = NULL; @@ -33,5 +33,5 @@ FetchMatlabData((DataSet**)&materials,MATERIALS); FetchMatlabData(¶meters,PARAMETERS); - FetchMatlabData(&control_index,CONTROLINDEX); + parameters->FindParam(&num_controls,InversionNumControlParametersEnum); /*configure: */ @@ -41,8 +41,8 @@ /*!Call core code: */ - Gradjx(&gradient,&grad_norm,elements,nodes, vertices,loads, materials,parameters, control_index); + Gradjx(&gradient,&norm_list,elements,nodes, vertices,loads, materials,parameters); /*write output : */ - WriteMatlabData(GRADNORM,grad_norm); + WriteMatlabData(NORMLIST,norm_list,num_controls,1); WriteMatlabData(GRADG,gradient); Index: /issm/trunk-jpl/src/mex/Gradj/Gradj.h =================================================================== --- /issm/trunk-jpl/src/mex/Gradj/Gradj.h (revision 11312) +++ /issm/trunk-jpl/src/mex/Gradj/Gradj.h (revision 11313) @@ -24,9 +24,8 @@ #define MATERIALS (mxArray*)prhs[4] #define PARAMETERS (mxArray*)prhs[5] -#define CONTROLINDEX (mxArray*)prhs[6] /* serial output macros: */ #define GRADG (mxArray**)&plhs[0] -#define GRADNORM (mxArray**)&plhs[1] +#define NORMLIST (mxArray**)&plhs[1] /* serial arg counts: */ @@ -34,5 +33,5 @@ #define NLHS 2 #undef NRHS -#define NRHS 7 +#define NRHS 6 #endif /* _GRADJ_H */