Changeset 11313


Ignore:
Timestamp:
02/03/12 07:55:52 (13 years ago)
Author:
Mathieu Morlighem
Message:

Objective function gradient is now ONE unique vector

Location:
issm/trunk-jpl/src
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.cpp

    r7177 r11313  
    99#include "../../EnumDefinitions/EnumDefinitions.h"
    1010
    11 void ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name){
     11void ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
    1212
    13         Vec gradient=NULL;
    14         gradient=NewVec(vertices->NumberOfVertices());
     13        /*Intermediaries*/
     14        int  num_controls;
     15        int *control_type = NULL;
     16        Vec  gradient=NULL;
    1517
    16         for(int i=0;i<elements->Size();i++){
    17                 Element* element=(Element*)elements->GetObjectByOffset(i);
    18                 element->ControlInputGetGradient(gradient,name);
     18        /*Retrieve some parameters*/
     19        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     20        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     21
     22        /*Allocate and populate gradient*/
     23        gradient=NewVec(num_controls*vertices->NumberOfVertices());
     24
     25        for(int i=0;i<num_controls;i++){
     26                for(int j=0;j<elements->Size();j++){
     27                        Element* element=(Element*)elements->GetObjectByOffset(j);
     28                        element->ControlInputGetGradient(gradient,control_type[i],i);
     29                }
    1930        }
    2031
    … …  
    2233        VecAssemblyEnd(gradient);
    2334
    24         /*Assign output pointers:*/
     35        /*Clean up and return*/
     36        xfree((void**)&control_type);
    2537        *pgradient=gradient;
    2638}

  • issm/trunk-jpl/src/c/modules/ControlInputGetGradientx/ControlInputGetGradientx.h

    r7177 r11313  
    88#include "../../Container/Container.h"
    99
    10 void    ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,int name);
     10void    ControlInputGetGradientx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters);
    1111
    1212#endif

  • issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp

    r7177 r11313  
    99#include "../../EnumDefinitions/EnumDefinitions.h"
    1010
    11 void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name,double scale){
     11void ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* norm_list,int step){
    1212
    13         for(int i=0;i<elements->Size();i++){
    14                 Element* element=(Element*)elements->GetObjectByOffset(i);
    15                 element->ControlInputScaleGradient(name,scale);
     13        /*Intermediaries*/
     14        int     i,j,num_controls;
     15        int    *control_type = NULL;
     16        double *scalar_list = NULL;
     17        double  scalar;
     18
     19
     20        /*Retrieve some parameters*/
     21        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     22        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     23        parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum);
     24
     25        /*Compute scaling factor*/
     26        scalar=scalar_list[0]/norm_list[0];
     27        for(i=0;i<num_controls;i++){
     28                scalar=min(scalar,scalar_list[num_controls*step+i]/norm_list[i]);
    1629        }
    1730
     31        for(i=0;i<num_controls;i++){
     32                for(j=0;j<elements->Size();j++){
     33                        Element* element=(Element*)elements->GetObjectByOffset(j);
     34                        element->ControlInputScaleGradient(control_type[i],scalar);
     35                }
     36        }
     37
     38        /*Clean up and return*/
     39        xfree((void**)&control_type);
     40        xfree((void**)&scalar_list);
    1841}

  • issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h

    r7177 r11313  
    88#include "../../Container/Container.h"
    99
    10 void    ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,int name,double scale);
     10void    ControlInputScaleGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,double* norm_list,int step);
    1111
    1212#endif

  • issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp

    r7177 r11313  
    99#include "../../EnumDefinitions/EnumDefinitions.h"
    1010
    11 void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name,double* gradient){
     11void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* gradient){
    1212
    13         for(int i=0;i<elements->Size();i++){
    14                 Element* element=(Element*)elements->GetObjectByOffset(i);
    15                 element->ControlInputSetGradient(gradient,name);
     13        /*Intermediaries*/
     14        int  num_controls;
     15        int *control_type = NULL;
     16
     17        /*Retrieve some parameters*/
     18        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     19        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     20
     21        for(int i=0;i<num_controls;i++){
     22                for(int j=0;j<elements->Size();j++){
     23                        Element* element=(Element*)elements->GetObjectByOffset(j);
     24                        element->ControlInputSetGradient(gradient,control_type[i],i);
     25                }
    1626        }
    1727
     28        /*Clean up and return*/
     29        xfree((void**)&control_type);
     30
    1831}
    19 void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name,Vec gradient){
     32void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vec gradient){
    2033
    2134        /*Serialize gradient*/
    … …  
    2336        VecToMPISerial(&serial_gradient,gradient);
    2437
    25         ControlInputSetGradientx(elements,nodes,vertices, loads, materials, parameters, name,serial_gradient);
     38        ControlInputSetGradientx(elements,nodes,vertices, loads, materials, parameters,serial_gradient);
    2639
    2740        /*Clean up and return*/

  • issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.h

    r7177 r11313  
    88#include "../../Container/Container.h"
    99
    10 void    ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,int name,double* gradient);
    11 void    ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,int name,Vec gradient);
     10void    ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,double* gradient);
     11void    ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,Vec gradient);
    1212
    1313#endif

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp

    r11310 r11313  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void Gradjx( Vec* pgradient,double* pnorm_grad, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int control_index){
     12void Gradjx(Vec* pgradient,double** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
    1313
    14         int    i;
    15         int    numberofvertices;
    16         double norm_grad;
    17         int   *control_type   = NULL;
    18         Vec    gradient = NULL;
     14        int     i,j,numberofvertices;
     15        int     num_controls;
     16        double  norm_inf;
     17        double *norm_list       = NULL;
     18        int    *control_type    = NULL;
     19        Vec     gradient        = NULL;
     20        Vec    *gradient_list   = NULL;
    1921       
    2022        /*retrieve some parameters: */
     23        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
    2124        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    2225        numberofvertices=vertices->NumberOfVertices();
    2326
    24         /*Allocate gradient: */
    25         gradient=NewVec(numberofvertices);
     27        /*Allocate gradient_list */
     28        gradient_list = (Vec*)xmalloc(num_controls*sizeof(Vec));
     29        norm_list     = (double*)xmalloc(num_controls*sizeof(double));
     30        for(i=0;i<num_controls;i++){
     31                gradient_list[i]=NewVec(num_controls*numberofvertices);
     32        }
     33        gradient=NewVec(num_controls*numberofvertices);
    2634
    27         /*Compute gradients: */
    28         for (i=0;i<elements->Size();i++){
    29                 Element* element=(Element*)elements->GetObjectByOffset(i);
    30                 element->Gradj(gradient,control_type[control_index]);
     35        /*Compute all gradient_list*/
     36        for(i=0;i<num_controls;i++){
     37
     38                for(j=0;j<elements->Size();j++){
     39                        Element* element=(Element*)elements->GetObjectByOffset(j);
     40                        element->Gradj(gradient_list[i],control_type[i],i);
     41                }
     42
     43                VecAssemblyBegin(gradient_list[i]);
     44                VecAssemblyEnd(gradient_list[i]);
     45
     46                VecNorm(gradient_list[i],NORM_INFINITY,&norm_list[i]);
    3147        }
    3248
    33         /*Assemble vector: */
    34         VecAssemblyBegin(gradient);
    35         VecAssemblyEnd(gradient);
     49        /*Add all gradient_list together*/
     50        for(i=0;i<num_controls;i++){
     51                VecAXPY(gradient,1.,gradient_list[i]);
     52                VecFree(&gradient_list[i]);
     53        }
    3654
    37         /*Clean up and return*/
     55        /*Check that gradient is clean*/
     56        VecNorm(gradient,NORM_INFINITY,&norm_inf);
     57        if(norm_inf<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm is zero");
     58        if(isnan(norm_inf))_error_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
     59
     60        /*Clean-up and assign output pointer*/
     61        *pnorm_list=norm_list;
     62        *pgradient=gradient;
     63        xfree((void**)&gradient_list);
    3864        xfree((void**)&control_type);
    39         if(pnorm_grad){
    40                 VecNorm(gradient,NORM_INFINITY,&norm_grad);
    41                 *pnorm_grad=norm_grad;
    42         }
    43         *pgradient=gradient;
    4465}

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h

    r11310 r11313  
    1010
    1111/* local prototypes: */
    12 void Gradjx(Vec* pgrad_g,double* pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,  Parameters* parameters, int control_type);
     12void Gradjx(Vec* pgrad_g,double** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters);
    1313
    1414#endif  /* _GRADJX_H */
    15 

  • issm/trunk-jpl/src/c/objects/Elements/Element.h

    r11260 r11313  
    9090
    9191                #ifdef _HAVE_CONTROL_
    92                 virtual void   Gradj(Vec gradient,int control_type)=0;
     92                virtual void   Gradj(Vec gradient,int control_type,int control_index)=0;
    9393                virtual double ThicknessAbsMisfit(bool process_units  ,int weight_index)=0;
    9494                virtual double SurfaceAbsVelMisfit(bool process_units ,int weight_index)=0;
    … …  
    100100                virtual double RheologyBbarAbsGradient(bool process_units,int weight_index)=0;
    101101                virtual double DragCoefficientAbsGradient(bool process_units,int weight_index)=0;
    102                 virtual void   ControlInputGetGradient(Vec gradient,int enum_type)=0;
    103                 virtual void   ControlInputSetGradient(double* gradient,int enum_type)=0;
     102                virtual void   ControlInputGetGradient(Vec gradient,int enum_type,int control_index)=0;
     103                virtual void   ControlInputSetGradient(double* gradient,int enum_type,int control_index)=0;
    104104                virtual void   ControlInputScaleGradient(int enum_type, double scale)=0;
    105105                virtual void   InputControlUpdate(double scalar,bool save_parameter)=0;

  • issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

    r11299 r11313  
    43144314#ifdef _HAVE_CONTROL_
    43154315/*FUNCTION Penta::ControlInputGetGradient{{{1*/
    4316 void Penta::ControlInputGetGradient(Vec gradient,int enum_type){
     4316void Penta::ControlInputGetGradient(Vec gradient,int enum_type,int control_index){
    43174317
    43184318        int doflist1[NUMVERTICES];
    … …  
    43294329        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
    43304330
    4331         this->GetDofList1(&doflist1[0]);
     4331        GradientIndexing(&doflist1[0],control_index);
    43324332        ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]);
    43334333
    … …  
    43504350}/*}}}*/
    43514351/*FUNCTION Penta::ControlInputSetGradient{{{1*/
    4352 void Penta::ControlInputSetGradient(double* gradient,int enum_type){
     4352void Penta::ControlInputSetGradient(double* gradient,int enum_type,int control_index){
    43534353
    43544354        int    doflist1[NUMVERTICES];
    … …  
    43664366        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
    43674367
    4368         this->GetDofList1(&doflist1[0]);
     4368        GradientIndexing(&doflist1[0],control_index);
    43694369        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]];
    43704370        grad_input=new PentaP1Input(GradientEnum,grad_list);
    … …  
    44344434}
    44354435/*}}}*/
     4436/*FUNCTION Penta::GradientIndexing{{{1*/
     4437void Penta::GradientIndexing(int* indexing,int control_index){
     4438
     4439        /*Get some parameters*/
     4440        int num_controls;
     4441        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     4442
     4443        /*get gradient indices*/
     4444        for(int i=0;i<NUMVERTICES;i++){
     4445                indexing[i]=num_controls*this->nodes[i]->GetVertexDof() + control_index;
     4446        }
     4447
     4448}
     4449/*}}}*/
    44364450/*FUNCTION Penta::Gradj {{{1*/
    4437 void  Penta::Gradj(Vec gradient,int control_type){
     4451void  Penta::Gradj(Vec gradient,int control_type,int control_index){
    44384452        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    44394453
    … …  
    44514465                        switch(approximation){
    44524466                                case MacAyealApproximationEnum:
    4453                                         GradjDragMacAyeal(gradient);
     4467                                        GradjDragMacAyeal(gradient,control_index);
    44544468                                        break;
    44554469                                case PattynApproximationEnum:
    4456                                         GradjDragPattyn(gradient);
     4470                                        GradjDragPattyn(gradient,control_index);
    44574471                                        break;
    44584472                                case StokesApproximationEnum:
    4459                                         GradjDragStokes(gradient);
     4473                                        GradjDragStokes(gradient,control_index);
    44604474                                        break;
    44614475                                case NoneApproximationEnum:
    … …  
    44714485                        switch(approximation){
    44724486                                case MacAyealApproximationEnum:
    4473                                         GradjBbarMacAyeal(gradient);
     4487                                        GradjBbarMacAyeal(gradient,control_index);
    44744488                                        break;
    44754489                                case PattynApproximationEnum:
    4476                                         GradjBbarPattyn(gradient);
     4490                                        GradjBbarPattyn(gradient,control_index);
    44774491                                        break;
    44784492                                case StokesApproximationEnum:
    4479                                         GradjBbarStokes(gradient);
     4493                                        GradjBbarStokes(gradient,control_index);
    44804494                                        break;
    44814495                                case NoneApproximationEnum:
    … …  
    45104524                        if (!IsOnBed()) return;
    45114525                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    4512                         tria->GradjDragGradient(gradient,resp);
     4526                        tria->GradjDragGradient(gradient,resp,control_index);
    45134527                        delete tria->matice; delete tria;
    45144528                        break;
    … …  
    45164530                        if (!IsOnBed()) return;
    45174531                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    4518                         tria->GradjBGradient(gradient,resp);
     4532                        tria->GradjBGradient(gradient,resp,control_index);
    45194533                        delete tria->matice; delete tria;
    45204534                        break;
    … …  
    45264540/*}}}*/
    45274541/*FUNCTION Penta::GradjDragMacAyeal {{{1*/
    4528 void  Penta::GradjDragMacAyeal(Vec gradient){
     4542void  Penta::GradjDragMacAyeal(Vec gradient,int control_index){
    45294543
    45304544        /*Gradient is 0 if on shelf or not on bed*/
    … …  
    45334547        /*Spawn tria*/
    45344548        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    4535         tria->GradjDragMacAyeal(gradient);
     4549        tria->GradjDragMacAyeal(gradient,control_index);
    45364550        delete tria->matice; delete tria;
    45374551
    45384552} /*}}}*/
    45394553/*FUNCTION Penta::GradjDragPattyn {{{1*/
    4540 void  Penta::GradjDragPattyn(Vec gradient){
     4554void  Penta::GradjDragPattyn(Vec gradient,int control_index){
    45414555
    45424556        int        i,j,ig;
    … …  
    45594573        /*Retrieve all inputs and parameters*/
    45604574        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     4575        GradientIndexing(&doflist1[0],control_index);
    45614576        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    45624577        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4563         GetDofList1(&doflist1[0]);
    45644578        Input* adjointx_input=inputs->GetInput(AdjointxEnum);               _assert_(adjointx_input);
    45654579        Input* adjointy_input=inputs->GetInput(AdjointyEnum);               _assert_(adjointy_input);
    … …  
    46094623/*}}}*/
    46104624/*FUNCTION Penta::GradjDragStokes {{{1*/
    4611 void  Penta::GradjDragStokes(Vec gradient){
     4625void  Penta::GradjDragStokes(Vec gradient,int control_index){
    46124626
    46134627        int        i,j,ig;
    … …  
    46344648        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    46354649        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4636         GetDofList1(&doflist1[0]);
     4650        GradientIndexing(&doflist1[0],control_index);
    46374651        Input* drag_input    =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);
    4638         Input* vx_input      =inputs->GetInput(VxEnum);              _assert_(vx_input);
    4639         Input* vy_input      =inputs->GetInput(VyEnum);              _assert_(vy_input);
    4640         Input* vz_input      =inputs->GetInput(VzEnum);              _assert_(vz_input);
    4641         Input* adjointx_input=inputs->GetInput(AdjointxEnum);        _assert_(adjointx_input);
    4642         Input* adjointy_input=inputs->GetInput(AdjointyEnum);        _assert_(adjointy_input);
    4643         Input* adjointz_input=inputs->GetInput(AdjointzEnum);        _assert_(adjointz_input);
     4652        Input* vx_input      =inputs->GetInput(VxEnum);                  _assert_(vx_input);
     4653        Input* vy_input      =inputs->GetInput(VyEnum);                  _assert_(vy_input);
     4654        Input* vz_input      =inputs->GetInput(VzEnum);                  _assert_(vz_input);
     4655        Input* adjointx_input=inputs->GetInput(AdjointxEnum);            _assert_(adjointx_input);
     4656        Input* adjointy_input=inputs->GetInput(AdjointyEnum);            _assert_(adjointy_input);
     4657        Input* adjointz_input=inputs->GetInput(AdjointzEnum);            _assert_(adjointz_input);
    46444658
    46454659        /*Build frictoin element, needed later: */
    … …  
    47014715/*}}}*/
    47024716/*FUNCTION Penta::GradjBbarMacAyeal {{{1*/
    4703 void  Penta::GradjBbarMacAyeal(Vec gradient){
     4717void  Penta::GradjBbarMacAyeal(Vec gradient,int control_index){
    47044718
    47054719        /*This element should be collapsed into a tria element at its base*/
    … …  
    47114725        /*Collapse element to the base*/
    47124726        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    4713         tria->GradjBMacAyeal(gradient);
     4727        tria->GradjBMacAyeal(gradient,control_index);
    47144728        delete tria->matice; delete tria;
    47154729
    … …  
    47194733} /*}}}*/
    47204734/*FUNCTION Penta::GradjBbarPattyn {{{1*/
    4721 void  Penta::GradjBbarPattyn(Vec gradient){
     4735void  Penta::GradjBbarPattyn(Vec gradient,int control_index){
    47224736
    47234737        /*Gradient is computed on bed only (Bbar)*/
    … …  
    47294743        /*Collapse element to the base*/
    47304744        Tria* tria=(Tria*)SpawnTria(0,1,2);
    4731         tria->GradjBMacAyeal(gradient);    //We use MacAyeal as an estimate for now
     4745        tria->GradjBMacAyeal(gradient,control_index);    //We use MacAyeal as an estimate for now
    47324746        delete tria->matice; delete tria;
    47334747
    … …  
    47364750} /*}}}*/
    47374751/*FUNCTION Penta::GradjBbarStokes {{{1*/
    4738 void  Penta::GradjBbarStokes(Vec gradient){
     4752void  Penta::GradjBbarStokes(Vec gradient,int control_index){
    47394753
    47404754        /*Gradient is computed on bed only (Bbar)*/
    … …  
    47464760        /*Collapse element to the base*/
    47474761        Tria* tria=(Tria*)SpawnTria(0,1,2);
    4748         tria->GradjBMacAyeal(gradient);    //We use MacAyeal as an estimate for now
     4762        tria->GradjBMacAyeal(gradient,control_index);    //We use MacAyeal as an estimate for now
    47494763        delete tria->matice; delete tria;
    47504764

  • issm/trunk-jpl/src/c/objects/Elements/Penta.h

    r11260 r11313  
    141141                #ifdef _HAVE_CONTROL_
    142142                double DragCoefficientAbsGradient(bool process_units,int weight_index);
    143                 void   Gradj(Vec gradient,int control_type);
    144                 void   GradjDragMacAyeal(Vec gradient);
    145                 void   GradjDragPattyn(Vec gradient);
    146                 void   GradjDragStokes(Vec gradient);
    147                 void   GradjBbarMacAyeal(Vec gradient);
    148                 void   GradjBbarPattyn(Vec gradient);
    149                 void   GradjBbarStokes(Vec gradient);
    150                 void   ControlInputGetGradient(Vec gradient,int enum_type);
     143                void   GradientIndexing(int* indexing,int control_index);
     144                void   Gradj(Vec gradient,int control_type,int control_index);
     145                void   GradjDragMacAyeal(Vec gradient,int control_index);
     146                void   GradjDragPattyn(Vec gradient,int control_index);
     147                void   GradjDragStokes(Vec gradient,int control_index);
     148                void   GradjBbarMacAyeal(Vec gradient,int control_index);
     149                void   GradjBbarPattyn(Vec gradient,int control_index);
     150                void   GradjBbarStokes(Vec gradient,int control_index);
     151                void   ControlInputGetGradient(Vec gradient,int enum_type,int control_index);
    151152                void   ControlInputScaleGradient(int enum_type,double scale);
    152                 void   ControlInputSetGradient(double* gradient,int enum_type);
     153                void   ControlInputSetGradient(double* gradient,int enum_type,int control_index);
    153154                double RheologyBbarAbsGradient(bool process_units,int weight_index);
    154155                double ThicknessAbsMisfit(     bool process_units,int weight_index);

  • issm/trunk-jpl/src/c/objects/Elements/Tria.cpp

    r11292 r11313  
    32663266/*}}}*/
    32673267/*FUNCTION Tria::ControlInputGetGradient{{{1*/
    3268 void Tria::ControlInputGetGradient(Vec gradient,int enum_type){
     3268void Tria::ControlInputGetGradient(Vec gradient,int enum_type,int control_index){
    32693269
    32703270        int doflist1[NUMVERTICES];
    … …  
    32803280        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
    32813281
    3282         this->GetDofList1(&doflist1[0]);
     3282        GradientIndexing(&doflist1[0],control_index);
    32833283        ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]);
    32843284
    … …  
    33013301}/*}}}*/
    33023302/*FUNCTION Tria::ControlInputSetGradient{{{1*/
    3303 void Tria::ControlInputSetGradient(double* gradient,int enum_type){
     3303void Tria::ControlInputSetGradient(double* gradient,int enum_type,int control_index){
    33043304
    33053305        int    doflist1[NUMVERTICES];
    … …  
    33173317        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
    33183318
    3319         this->GetDofList1(&doflist1[0]);
     3319        GradientIndexing(&doflist1[0],control_index);
    33203320        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]];
    33213321        grad_input=new TriaP1Input(GradientEnum,grad_list);
    … …  
    33253325}/*}}}*/
    33263326/*FUNCTION Tria::Gradj {{{1*/
    3327 void  Tria::Gradj(Vec gradient,int control_type){
     3327void  Tria::Gradj(Vec gradient,int control_type,int control_index){
    33283328        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    33293329
    … …  
    33343334        switch(control_type){
    33353335                case FrictionCoefficientEnum:
    3336                         GradjDragMacAyeal(gradient);
     3336                        GradjDragMacAyeal(gradient,control_index);
    33373337                        break;
    33383338                case MaterialsRheologyBbarEnum:
    3339                         GradjBMacAyeal(gradient);
     3339                        GradjBMacAyeal(gradient,control_index);
    33403340                        break;
    33413341                case BalancethicknessThickeningRateEnum:
    3342                         GradjDhDtBalancedthickness(gradient);
     3342                        GradjDhDtBalancedthickness(gradient,control_index);
    33433343                        break;
    33443344                case VxEnum:
    3345                         GradjVxBalancedthickness(gradient);
     3345                        GradjVxBalancedthickness(gradient,control_index);
    33463346                        break;
    33473347                case VyEnum:
    3348                         GradjVyBalancedthickness(gradient);
     3348                        GradjVyBalancedthickness(gradient,control_index);
    33493349                        break;
    33503350                default:
    … …  
    33713371                        break;
    33723372                case DragCoefficientAbsGradientEnum:
    3373                         GradjDragGradient(gradient,resp);
     3373                        GradjDragGradient(gradient,resp,control_index);
    33743374                        break;
    33753375                case RheologyBbarAbsGradientEnum:
    3376                         GradjBGradient(gradient,resp);
     3376                        GradjBGradient(gradient,resp,control_index);
    33773377                        break;
    33783378                default:
    … …  
    33843384/*}}}*/
    33853385/*FUNCTION Tria::GradjBGradient{{{1*/
    3386 void  Tria::GradjBGradient(Vec gradient, int weight_index){
     3386void  Tria::GradjBGradient(Vec gradient,int weight_index,int control_index){
    33873387
    33883388        int        i,ig;
    … …  
    33973397        /*Retrieve all inputs we will be needing: */
    33983398        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3399         GetDofList1(&doflist1[0]);
     3399        GradientIndexing(&doflist1[0],control_index);
    34003400        Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    34013401        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
    … …  
    34243424/*}}}*/
    34253425/*FUNCTION Tria::GradjBMacAyeal{{{1*/
    3426 void  Tria::GradjBMacAyeal(Vec gradient){
     3426void  Tria::GradjBMacAyeal(Vec gradient,int control_index){
    34273427
    34283428        /*Intermediaries*/
    … …  
    34393439        /* Get node coordinates and dof list: */
    34403440        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3441         GetDofList1(&doflist[0]);
     3441        GradientIndexing(&doflist[0],control_index);
    34423442
    34433443        /*Retrieve all inputs*/
    3444         Input* thickness_input=inputs->GetInput(ThicknessEnum);            _assert_(thickness_input);
    3445         Input* vx_input=inputs->GetInput(VxEnum);                          _assert_(vx_input);
    3446         Input* vy_input=inputs->GetInput(VyEnum);                          _assert_(vy_input);
    3447         Input* adjointx_input=inputs->GetInput(AdjointxEnum);              _assert_(adjointx_input);
    3448         Input* adjointy_input=inputs->GetInput(AdjointyEnum);              _assert_(adjointy_input);
     3444        Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
     3445        Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
     3446        Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
     3447        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
     3448        Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
    34493449        Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    34503450
    … …  
    34813481/*}}}*/
    34823482/*FUNCTION Tria::GradjDragMacAyeal {{{1*/
    3483 void  Tria::GradjDragMacAyeal(Vec gradient){
     3483void  Tria::GradjDragMacAyeal(Vec gradient,int control_index){
    34843484
    34853485        int        i,ig;
    … …  
    35023502        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    35033503        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3504         GetDofList1(&doflist1[0]);
     3504        GradientIndexing(&doflist1[0],control_index);
    35053505
    35063506        /*Build frictoin element, needed later: */
    … …  
    35083508
    35093509        /*Retrieve all inputs we will be needing: */
    3510         Input* adjointx_input=inputs->GetInput(AdjointxEnum);               _assert_(adjointx_input);
    3511         Input* adjointy_input=inputs->GetInput(AdjointyEnum);               _assert_(adjointy_input);
    3512         Input* vx_input=inputs->GetInput(VxEnum);                           _assert_(vx_input);
    3513         Input* vy_input=inputs->GetInput(VyEnum);                           _assert_(vy_input);
     3510        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                   _assert_(adjointx_input);
     3511        Input* adjointy_input=inputs->GetInput(AdjointyEnum);                   _assert_(adjointy_input);
     3512        Input* vx_input=inputs->GetInput(VxEnum);                               _assert_(vx_input);
     3513        Input* vy_input=inputs->GetInput(VyEnum);                               _assert_(vy_input);
    35143514        Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    35153515
    … …  
    35683568/*}}}*/
    35693569/*FUNCTION Tria::GradjDragGradient{{{1*/
    3570 void  Tria::GradjDragGradient(Vec gradient, int weight_index){
     3570void  Tria::GradjDragGradient(Vec gradient, int weight_index,int control_index){
    35713571
    35723572        int        i,ig;
    … …  
    35823582        if(IsFloating())return;
    35833583        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3584         GetDofList1(&doflist1[0]);
     3584        GradientIndexing(&doflist1[0],control_index);
    35853585        Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    35863586        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                 _assert_(weights_input);
    … …  
    36123612/*}}}*/
    36133613/*FUNCTION Tria::GradjDhDtBalancedthickness{{{1*/
    3614 void  Tria::GradjDhDtBalancedthickness(Vec gradient){
     3614void  Tria::GradjDhDtBalancedthickness(Vec gradient,int control_index){
    36153615
    36163616        /*Intermediaries*/
    … …  
    36193619        double gradient_g[NUMVERTICES];
    36203620
    3621         GetDofList1(&doflist1[0]);
    3622 
    36233621        /*Compute Gradient*/
     3622        GradientIndexing(&doflist1[0],control_index);
    36243623        GetInputListOnVertices(&lambda[0],AdjointEnum);
    36253624        for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i];
    … …  
    36293628/*}}}*/
    36303629/*FUNCTION Tria::GradjVxBalancedthickness{{{1*/
    3631 void  Tria::GradjVxBalancedthickness(Vec gradient){
     3630void  Tria::GradjVxBalancedthickness(Vec gradient,int control_index){
    36323631
    36333632        /*Intermediaries*/
    … …  
    36433642        /* Get node coordinates and dof list: */
    36443643        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3645         GetDofList1(&doflist1[0]);
     3644        GradientIndexing(&doflist1[0],control_index);
    36463645
    36473646        /*Retrieve all inputs we will be needing: */
    … …  
    36723671/*}}}*/
    36733672/*FUNCTION Tria::GradjVyBalancedthickness{{{1*/
    3674 void  Tria::GradjVyBalancedthickness(Vec gradient){
     3673void  Tria::GradjVyBalancedthickness(Vec gradient,int control_index){
    36753674
    36763675        /*Intermediaries*/
    … …  
    36863685        /* Get node coordinates and dof list: */
    36873686        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3688         GetDofList1(&doflist1[0]);
     3687        GradientIndexing(&doflist1[0],control_index);
    36893688
    36903689        /*Retrieve all inputs we will be needing: */
    … …  
    37113710        /*Clean up and return*/
    37123711        delete gauss;
     3712}
     3713/*}}}*/
     3714/*FUNCTION Tria::GradientIndexing{{{1*/
     3715void  Tria::GradientIndexing(int* indexing,int control_index){
     3716
     3717        /*Get some parameters*/
     3718        int num_controls;
     3719        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     3720
     3721        /*get gradient indices*/
     3722        for(int i=0;i<NUMVERTICES;i++){
     3723                indexing[i]=num_controls*this->nodes[i]->GetVertexDof() + control_index;
     3724        }
     3725
    37133726}
    37143727/*}}}*/
    … …  
    42254238        /*Clean up and return*/
    42264239        delete gauss;
     4240        xfree((void**)&responses);
    42274241        return pe;
    42284242}

  • issm/trunk-jpl/src/c/objects/Elements/Tria.h

    r11260 r11313  
    140140                #ifdef _HAVE_CONTROL_
    141141                double DragCoefficientAbsGradient(bool process_units,int weight_index);
    142                 void   Gradj(Vec gradient,int control_type);
    143                 void   GradjBGradient(Vec gradient,int weight_index);
    144                 void   GradjBMacAyeal(Vec gradient);
    145                 void   GradjDragMacAyeal(Vec gradient);
    146                 void   GradjDragStokes(Vec gradient);
    147                 void   GradjDragGradient(Vec gradient,int weight_index);
    148                 void   GradjDhDtBalancedthickness(Vec gradient);
    149                 void   GradjVxBalancedthickness(Vec gradient);
    150                 void   GradjVyBalancedthickness(Vec gradient);
    151                 void   ControlInputGetGradient(Vec gradient,int enum_type);
     142                void   GradientIndexing(int* indexing,int control_index);
     143                void   Gradj(Vec gradient,int control_type,int control_index);
     144                void   GradjBGradient(Vec gradient,int weight_index,int control_index);
     145                void   GradjBMacAyeal(Vec gradient,int control_index);
     146                void   GradjDragMacAyeal(Vec gradient,int control_index);
     147                void   GradjDragStokes(Vec gradient,int control_index);
     148                void   GradjDragGradient(Vec gradient,int weight_index,int control_index);
     149                void   GradjDhDtBalancedthickness(Vec gradient,int control_index);
     150                void   GradjVxBalancedthickness(Vec gradient,int control_index);
     151                void   GradjVyBalancedthickness(Vec gradient,int control_index);
     152                void   ControlInputGetGradient(Vec gradient,int enum_type,int control_index);
    152153                void   ControlInputScaleGradient(int enum_type,double scale);
    153                 void   ControlInputSetGradient(double* gradient,int enum_type);
     154                void   ControlInputSetGradient(double* gradient,int enum_type,int control_index);
    154155                double RheologyBbarAbsGradient(bool process_units,int weight_index);
    155156                double ThicknessAbsMisfit(     bool process_units,int weight_index);

  • issm/trunk-jpl/src/c/solutions/controltao_core.cpp

    r11310 r11313  
    8585        adjointdiagnostic_core(user->femmodel);
    8686        //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum);
    87         Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,0);
     87        Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    8888        VecScale(gradient,-1.);
    8989        VecCopy(gradient,G);

  • issm/trunk-jpl/src/c/solutions/gradient_core.cpp

    r11310 r11313  
    1414
    1515void gradient_core(FemModel* femmodel,int step,bool orthogonalize){
    16        
    17         /*parameters: */
    18         bool    control_steady;
    19         int     num_controls;
    20         int    *control_type   = NULL;
    21         double *optscal_list   = NULL;
    22         double  optscal=1.e+100,norm_grad;
    2316
    2417        /*Intermediaries*/
    25         Vec new_gradient=NULL;
    26         Vec gradient=NULL;
    27         Vec old_gradient=NULL;
     18        double  norm_inf;
     19        double *norm_list    = NULL;
     20        Vec     new_gradient = NULL;
     21        Vec     gradient     = NULL;
     22        Vec     old_gradient = NULL;
    2823
    29         /*retrieve parameters:*/
    30         femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
    31         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    32         femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    33         femmodel->parameters->FindParam(&optscal_list,NULL,NULL,InversionGradientScalingEnum);
     24        /*Compute gradient*/
     25        _printf_(VerboseControl(),"   compute cost function gradient\n");
     26        Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
    3427
    35         /*Compute and norm gradient of all controls*/
    36         for (int i=0;i<num_controls;i++){
    37 
    38                 _printf_(VerboseControl(),"   compute gradient of J with respect to %s\n",EnumToStringx(control_type[i]));
    39                 Gradjx(&gradient,&norm_grad,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,i);
    40 
    41                 if (orthogonalize){
    42                         _printf_(VerboseControl(),"   orthogonalization\n");
    43                         ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]);
    44                         Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient);
    45                 }
    46                 else{
    47                         new_gradient=gradient;
    48                 }
    49 
    50                 /*Get scaling factor of current control:*/
    51                 VecNorm(new_gradient,NORM_INFINITY,&norm_grad);
    52                 if(norm_grad<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm of J with respect to %s is zero",EnumToStringx(control_type[i]));
    53                 if(isnan(norm_grad))_error_("||∂J/∂α||∞ = NaN  gradient norm of J with respect to %s is NaN" ,EnumToStringx(control_type[i]));
    54                 optscal=min(optscal,optscal_list[num_controls*step+i]/norm_grad);
    55 
    56                 /*plug back into inputs: */
    57                 ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],new_gradient);
    58                 VecFree(&new_gradient);
     28        if (orthogonalize){
     29                _printf_(VerboseControl(),"   orthogonalization\n");
     30                ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
     31                Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient);
     32        }
     33        else{
     34                new_gradient=gradient;
    5935        }
    6036
     37        /*Check that gradient is clean*/
     38        VecNorm(new_gradient,NORM_INFINITY,&norm_inf);
     39        if(norm_inf<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm is zero");
     40        if(isnan(norm_inf))_error_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
     41
     42        /*plug back into inputs: */
     43        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,new_gradient);
     44        VecFree(&new_gradient);
     45
    6146        /*Scale Gradients*/
    62         for (int i=0;i<num_controls;i++) ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],optscal);
     47        ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,norm_list,step);
    6348
    6449        /*Clean up and return*/
    65         xfree((void**)&control_type);
    66         xfree((void**)&optscal_list);
     50        xfree((void**)&norm_list);
    6751}

  • issm/trunk-jpl/src/m/solutions/control_core.m

    r11306 r11313  
    5050                eval(['femmodel=' adjointcore '(femmodel);']);
    5151
    52                 femmodel=gradient_core(femmodel,n,search_scalar==0);
     52                femmodel=gradient_core(femmodel,n-1,search_scalar==0);
    5353
    5454                %Return gradient if asked

  • issm/trunk-jpl/src/m/solutions/gradient_core.m

    r11310 r11313  
    2020end
    2121
    22 %recover parameters common to all solutions
    23 num_controls=femmodel.parameters.InversionNumControlParameters;
    24 control_type=femmodel.parameters.InversionControlParameters;
    25 control_steady=femmodel.parameters.ControlSteady;
    26 gradient_scaling_list=femmodel.parameters.InversionGradientScaling;
    27 gradient_scaling=10^100;
     22issmprintf(VerboseControl,['   compute cost function gradient']);
     23[grad norm_list]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    2824
    29 for i=1:num_controls,
    30 
    31         issmprintf(VerboseControl,['   compute gradient of J with respect to %s'],EnumToString(control_type(i)));
    32         [grad norm_grad]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,i-1);
    33 
    34         if orthogonalize,
    35                 issmprintf(VerboseControl,'%s',['   orthogonalization']);
    36                 old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i));
    37                 new_gradient=Orth(grad,old_gradient);
    38         else
    39                 new_gradient=grad;
    40         end
    41 
    42          %Get scaling factor of current control:
    43          norm_grad=norm(new_gradient,inf);
    44          if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm of J with respect to ' EnumToString(control_type(i))  ' is zero']); end
    45          if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i))  ' is NaN' ]); end
    46          gradient_scaling=min(gradient_scaling,gradient_scaling_list(step,i)/norm_grad);
    47 
    48         %plug back into inputs:
    49         [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);
    50 
     25if orthogonalize,
     26        issmprintf(VerboseControl,'%s',['   orthogonalization']);
     27        old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters);
     28        new_gradient=Orth(grad,old_gradient);
     29else
     30        new_gradient=grad;
    5131end
    5232
     33%Check that gradient is clean
     34norm_grad=norm(new_gradient,inf);
     35if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm is zero']); end
     36if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm is NaN' ]); end
     37
     38%plug back into inputs:
     39[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);
     40
    5341%Scale all gradients
    54 for i=1:num_controls,
    55         [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);
    56 end
     42[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);

  • issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.cpp

    r8910 r11313  
    1414        Materials  *materials  = NULL;
    1515        Parameters *parameters = NULL;
    16         int         control_type;
    1716        Vec         gradient   = NULL;
    1817
    … …  
    3231        FetchMatlabData((DataSet**)&materials,MATERIALS);
    3332        FetchMatlabData(&parameters,PARAMETERS);
    34         FetchMatlabData(&control_type,CONTROLTYPE);
    3533
    3634        /*configure: */
    … …  
    4139
    4240        /*!core code:*/
    43         ControlInputGetGradientx(&gradient,elements, nodes,vertices,loads, materials,parameters,control_type);
     41        ControlInputGetGradientx(&gradient,elements, nodes,vertices,loads, materials,parameters);
    4442
    4543        /*write output datasets: */
    … …  
    6159{
    6260        _printf_(true,"\n");
    63         _printf_(true,"   usage: [gradient] = %s(elements,nodes,vertices,loads, materials,parameters,control_type);\n",__FUNCT__);
     61        _printf_(true,"   usage: [gradient] = %s(elements,nodes,vertices,loads, materials,parameters);\n",__FUNCT__);
    6462        _printf_(true,"\n");
    6563}

  • issm/trunk-jpl/src/mex/ControlInputGetGradient/ControlInputGetGradient.h

    r6200 r11313  
    2323#define MATERIALS (mxArray*)prhs[4]
    2424#define PARAMETERS (mxArray*)prhs[5]
    25 #define CONTROLTYPE (mxArray*)prhs[6]
    2625
    2726/* serial output macros: */
    … …  
    3231#define NLHS  1
    3332#undef NRHS
    34 #define NRHS  7
     33#define NRHS  6
    3534
    3635#endif

  • issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.cpp

    r8910 r11313  
    88
    99/*input datasets: */
     10int       step;
    1011Elements   *elements   = NULL;
    1112Nodes      *nodes      = NULL;
    … …  
    1415Materials  *materials  = NULL;
    1516Parameters *parameters = NULL;
    16 int         control_type;
    17 double      scaling_factor;
     17double     *norm_list  = NULL;
    1818
    1919/*Boot module: */
    … …  
    3030FetchMatlabData((DataSet**)&materials,MATERIALSIN);
    3131FetchMatlabData(&parameters,PARAMETERSIN);
    32 FetchMatlabData(&control_type,CONTROLTYPE);
    33 FetchMatlabData(&scaling_factor,SCALE);
     32FetchMatlabData(&norm_list,NULL,NULL,NORMLIST);
     33FetchMatlabData(&step,STEP);
    3434
    3535/*configure: */
    … …  
    3939
    4040/*call "x" code layer*/
    41 ControlInputScaleGradientx(elements,nodes,vertices,loads, materials,parameters,control_type,scaling_factor);
     41ControlInputScaleGradientx(elements,nodes,vertices,loads, materials,parameters,norm_list,step);
    4242
    4343/*write output datasets: */
    … …  
    6464{
    6565        _printf_(true,"\n");
    66         _printf_(true,"   usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,control_type,scaling_factor);\n",__FUNCT__);
     66        _printf_(true,"   usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,norm_list,step);\n",__FUNCT__);
    6767        _printf_(true,"\n");
    6868}

  • issm/trunk-jpl/src/mex/ControlInputScaleGradient/ControlInputScaleGradient.h

    r6238 r11313  
    2424#define MATERIALSIN (mxArray*)prhs[4]
    2525#define PARAMETERSIN (mxArray*)prhs[5]
    26 #define CONTROLTYPE (mxArray*)prhs[6]
    27 #define SCALE (mxArray*)prhs[7]
     26#define NORMLIST (mxArray*)prhs[6]
     27#define STEP (mxArray*)prhs[7]
    2828
    2929/* serial output macros: */

  • issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.cpp

    r8910 r11313  
    1414Materials  *materials  = NULL;
    1515Parameters *parameters = NULL;
    16 int         control_type;
    1716double     *gradient   = NULL;
    1817
    … …  
    3029FetchMatlabData((DataSet**)&materials,MATERIALSIN);
    3130FetchMatlabData(&parameters,PARAMETERSIN);
    32 FetchMatlabData(&control_type,CONTROLTYPE);
    3331FetchMatlabData(&gradient,NULL,GRADIENT);
    3432
    … …  
    3937
    4038/*call "x" code layer*/
    41 ControlInputSetGradientx(elements,nodes,vertices,loads, materials,parameters,control_type,gradient);
     39ControlInputSetGradientx(elements,nodes,vertices,loads, materials,parameters,gradient);
    4240
    4341/*write output datasets: */
    … …  
    6462{
    6563        _printf_(true,"\n");
    66         _printf_(true,"   usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,control_type,gradient);\n",__FUNCT__);
     64        _printf_(true,"   usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,gradient);\n",__FUNCT__);
    6765        _printf_(true,"\n");
    6866}

  • issm/trunk-jpl/src/mex/ControlInputSetGradient/ControlInputSetGradient.h

    r6200 r11313  
    2424#define MATERIALSIN (mxArray*)prhs[4]
    2525#define PARAMETERSIN (mxArray*)prhs[5]
    26 #define CONTROLTYPE (mxArray*)prhs[6]
    27 #define GRADIENT (mxArray*)prhs[7]
     26#define GRADIENT (mxArray*)prhs[6]
    2827
    2928/* serial output macros: */
    … …  
    3938#define NLHS  6
    4039#undef NRHS
    41 #define NRHS  8
     40#define NRHS  7
    4241
    4342#endif

  • issm/trunk-jpl/src/mex/Gradj/Gradj.cpp

    r11310 r11313  
    88
    99        /*input datasets: */
    10         int         control_index;
    11         double      grad_norm;
     10        int         control_index,num_controls;
     11        double     *norm_list    = NULL;
    1212        Elements   *elements     = NULL;
    1313        Nodes      *nodes        = NULL;
    … …  
    3333        FetchMatlabData((DataSet**)&materials,MATERIALS);
    3434        FetchMatlabData(&parameters,PARAMETERS);
    35         FetchMatlabData(&control_index,CONTROLINDEX);
     35        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    3636
    3737        /*configure: */
    … …  
    4141
    4242        /*!Call core code: */
    43         Gradjx(&gradient,&grad_norm,elements,nodes, vertices,loads, materials,parameters, control_index);
     43        Gradjx(&gradient,&norm_list,elements,nodes, vertices,loads, materials,parameters);
    4444
    4545        /*write output : */
    46         WriteMatlabData(GRADNORM,grad_norm);
     46        WriteMatlabData(NORMLIST,norm_list,num_controls,1);
    4747        WriteMatlabData(GRADG,gradient);
    4848

  • issm/trunk-jpl/src/mex/Gradj/Gradj.h

    r11310 r11313  
    2424#define MATERIALS (mxArray*)prhs[4]
    2525#define PARAMETERS (mxArray*)prhs[5]
    26 #define CONTROLINDEX (mxArray*)prhs[6]
    2726
    2827/* serial output macros: */
    2928#define GRADG (mxArray**)&plhs[0]
    30 #define GRADNORM (mxArray**)&plhs[1]
     29#define NORMLIST (mxArray**)&plhs[1]
    3130
    3231/* serial arg counts: */
    … …  
    3433#define NLHS  2
    3534#undef NRHS
    36 #define NRHS  7
     35#define NRHS  6
    3736
    3837#endif  /* _GRADJ_H */
Note: See TracChangeset for help on using the changeset viewer.