Changeset 4048


Ignore:
Timestamp:
06/16/10 18:20:52 (15 years ago)
Author:
Eric.Larour
Message:

Serious cleanup of the control solutoin

Location:
issm/trunk/src/c
Files:
7 added
36 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r4047 r4048  
    393393                                        ./modules/ScaleInputx/ScaleInputx.h\
    394394                                        ./modules/ScaleInputx/ScaleInputx.cpp\
     395                                        ./modules/AXPYInputx/AXPYInputx.h\
     396                                        ./modules/AXPYInputx/AXPYInputx.cpp\
    395397                                        ./modules/ControlConstrainx/ControlConstrainx.h\
    396398                                        ./modules/ControlConstrainx/ControlConstrainx.cpp\
    … …  
    891893                                        ./modules/ScaleInputx/ScaleInputx.h\
    892894                                        ./modules/ScaleInputx/ScaleInputx.cpp\
     895                                        ./modules/AXPYInputx/AXPYInputx.h\
     896                                        ./modules/AXPYInputx/AXPYInputx.cpp\
    893897                                        ./modules/ControlConstrainx/ControlConstrainx.h\
    894898                                        ./modules/ControlConstrainx/ControlConstrainx.cpp\

  • issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.cpp

    r3913 r4048  
    1 /*!\file ControlConstrainx
    2  * \brief constrain optimization parameters during inverse method
     1/*!\file ControlContrainInputx
     2 * \brief: Y=Y+aX operation on inputs.
    33 */
    44
    5 #include "./ControlConstrainx.h"
    6 
     5#include "./ControlContrainInputx.h"
    76#include "../../shared/shared.h"
    87#include "../../include/include.h"
    … …  
    109#include "../../EnumDefinitions/EnumDefinitions.h"
    1110
    12 void ControlConstrainx( double* p_g, int numdofnodes, double cm_min, double cm_max,int control_type){
     11void ControlContrainInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int control_type,double cm_min, double cm_max){
    1312
    14         int i;
    15        
    16         if(isnan(cm_min) & isnan(cm_max)){
    17                 /*do nothing*/
     13        /*intermediary:*/
     14        int      i;
     15
     16    /*some early returns: */
     17        if(isnan(cm_min) & isnan(cm_max))return;
     18
     19        /*First, get elements and nodes configured: */
     20        elements->Configure(elements,loads, nodes,vertices, materials,parameters);
     21
     22        /*Go through elemnets, and ask to carry out the ControlContrain operation on inputs: */
     23        for(i=0;i<elements->Size();i++){
     24                Element* element=(Element*)elements->GetObjectByOffset(i);
     25                element->ControlContrainInput(control_type,cm_min,cm_max);
    1826        }
    19         else{
    20                 for(i=0;i<numdofnodes;i++){
    21                         if(isnan(p_g[i])){
    22                                 ISSMERROR("NaN found in parameter p_g[%i]",i);
    23                         }
    24                         if(!isnan(cm_min)){
    25                                 if (p_g[i]<cm_min)p_g[i]=cm_min;
    26                         }
    27                         if(!isnan(cm_max)){
    28                                 if (p_g[i]>cm_max)p_g[i]=cm_max;
    29                         }
    30                 }
    31         }
    32                
    3327}

  • issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.h

    r3913 r4048  
    1 /*!\file:  ControlConstrainx.h
    2  * \brief constrain optimization parameters during inverse method
     1/*!\file:  ControlContrainInputx.h
     2 * \brief header file for field extrusion
    33 */
    44
    5 #ifndef _CONTROLCONSTRAINX_H
    6 #define _CONTROLCONSTRAINX_H
     5#ifndef _CONTROLCONTRAININPUTX_H
     6#define _CONTROLCONTRAININPUTX_H
    77
    88#include "../../DataSet/DataSet.h"
    99
    1010/* local prototypes: */
    11 void ControlConstrainx( double* p_g, int numberofnodes, double cm_min, double cm_max,int control_type);
     11void ControlContrainInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int control_type,double cm_min, double cm_max);
    1212
    13 #endif  /* _CONTROLCONSTRAINX_H */
     13#endif  /* _ControlContrainINPUTX_H */
    1414

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

    r4047 r4048  
    2424        parameters->FindParam(&dim,DimEnum);
    2525
    26 
    2726        /*Compute gradients: */
    2827        for (i=0;i<elements->Size();i++){

  • issm/trunk/src/c/modules/Orthx/Orthx.cpp

    r3519 r4048  
    66
    77void    Orthx( Vec* pnewgradj, Vec gradj, Vec oldgradj){
     8       
     9        GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,  femmodel->parameters, analysis_type,sub_analysis_type);
    810
    911        /*output: */

  • issm/trunk/src/c/modules/modules.h

    r4047 r4048  
    8080#include "./DuplicateInputx/DuplicateInputx.h"
    8181#include "./ScaleInputx/ScaleInputx.h"
     82#include "./AXPYInputx/AXPYInputx.h"
     83#include "./ControlContrainx/ControlContrainx.h"
    8284
    8385#endif

  • issm/trunk/src/c/objects/Elements/Beam.cpp

    r4047 r4048  
    10021002}
    10031003/*}}}*/
     1004/*FUNCTION Beam::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
     1005void  Beam::AXPY(int YEnum, double scalar, int XEnum){
     1006
     1007        Input* xinput=NULL;
     1008        Input* yinput=NULL;
     1009
     1010        /*Find x and y inputs: */
     1011        xinput=(Input*)this->inputs->GetInput(XEnum);
     1012        yinput=(Input*)this->inputs->GetInput(YEnum);
     1013
     1014        /*some checks: */
     1015        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
     1016        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
     1017
     1018        /*Scale: */
     1019        yinput->AXPY(xinput,scalar);
     1020}
     1021/*}}}*/
     1022/*FUNCTION Beam::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
     1023void  Beam::ControlConstrain(int control_type, double cm_min, double cm_max){
     1024
     1025        Input* input=NULL;
     1026
     1027        /*Find input: */
     1028        input=(Input*)this->inputs->GetInput(control_type);
     1029       
     1030        /*Do nothing if we  don't find it: */
     1031        if(!input)return;
     1032
     1033        /*Constrain input using cm_min and cm_max: */
     1034        input->Constrain(cm_min,cm_max);
     1035
     1036}
     1037/*}}}*/
     1038/*FUNCTION Beam::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
     1039void  Beam::GetVectorFromInputs(Vec vector,int NameEnum){
     1040
     1041        const int numvertices=2;
     1042        int doflist1[numvertices];
     1043
     1044        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
     1045        for(i=0;i<this->inputs->Size();i++){
     1046                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
     1047                if(input->EnumType==NameEnum){
     1048                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     1049                        this->GetDofList1(&doflist1[0]);
     1050                        input->GetVectorFromInputs(vector,&doflist1[0]);
     1051                        break;
     1052                }
     1053        }
     1054}
     1055/*}}}*/

  • issm/trunk/src/c/objects/Elements/Beam.h

    r4047 r4048  
    9696                void  DuplicateInput(int original_enum,int new_enum);
    9797                void  ScaleInput(int enum_type,double scale_factor);
     98                void  AXPY(int YEnum, double scalar, int XEnum);
     99                void  ControlConstrain(int control_type,double cm_min, double cm_max);
    98100
    99101                /*}}}*/
    … …  
    114116                void  GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord);
    115117                double MassFlux(double* segment);
     118                void GetVectorFromInputs(Vec vector,int NameEnum);
     119
    116120                /*}}}*/
    117121

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

    r4047 r4048  
    6666                virtual void   DuplicateInput(int original_enum,int new_enum)=0;
    6767                virtual void   ScaleInput(int enum_type,double scale_factor)=0;
    68 
    69 
     68                virtual void   GetVectorFromInputs(Vec vector,int NameEnum)=0;
     69                virtual void   AXPY(int YEnum, double scalar, int XEnum)=0;
     70                virtual void   ControlConstrain(int control_type,double cm_min, double cm_max)=0;
    7071                /*Implementation: */
    7172

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

    r4047 r4048  
    53525352}
    53535353/*}}}*/
     5354/*FUNCTION Penta::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
     5355void  Penta::AXPY(int YEnum, double scalar, int XEnum){
     5356
     5357        Input* xinput=NULL;
     5358        Input* yinput=NULL;
     5359
     5360        /*Find x and y inputs: */
     5361        xinput=(Input*)this->inputs->GetInput(XEnum);
     5362        yinput=(Input*)this->inputs->GetInput(YEnum);
     5363
     5364        /*some checks: */
     5365        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
     5366        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
     5367
     5368        /*Scale: */
     5369        yinput->AXPY(xinput,scalar);
     5370}
     5371/*}}}*/
     5372/*FUNCTION Penta::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
     5373void  Penta::ControlConstrain(int control_type, double cm_min, double cm_max){
     5374
     5375        Input* input=NULL;
     5376
     5377        /*Find input: */
     5378        input=(Input*)this->inputs->GetInput(control_type);
     5379       
     5380        /*Do nothing if we  don't find it: */
     5381        if(!input)return;
     5382
     5383        /*Constrain input using cm_min and cm_max: */
     5384        input->Constrain(cm_min,cm_max);
     5385
     5386}
     5387/*}}}*/
     5388/*FUNCTION Penta::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
     5389void  Penta::GetVectorFromInputs(Vec vector,int NameEnum){
     5390
     5391        const int numvertices=6;
     5392        int doflist1[numvertices];
     5393
     5394        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
     5395        for(i=0;i<this->inputs->Size();i++){
     5396                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
     5397                if(input->EnumType==NameEnum){
     5398                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     5399                        this->GetDofList1(&doflist1[0]);
     5400                        input->GetVectorFromInputs(vector,&doflist1[0]);
     5401                        break;
     5402                }
     5403        }
     5404}
     5405/*}}}*/

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

    r4047 r4048  
    163163                void  DuplicateInput(int original_enum,int new_enum);
    164164                void  ScaleInput(int enum_type,double scale_factor);
     165                void  AXPY(int YEnum, double scalar, int XEnum);
     166                void  ControlConstrain(int control_type,double cm_min, double cm_max);
     167                void  GetVectorFromInputs(Vec vector,int NameEnum);
    165168
    166169

  • issm/trunk/src/c/objects/Elements/Sing.cpp

    r4047 r4048  
    708708}
    709709/*}}}*/
     710/*FUNCTION Sing::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
     711void  Sing::AXPY(int YEnum, double scalar, int XEnum){
     712
     713        Input* xinput=NULL;
     714        Input* yinput=NULL;
     715
     716        /*Find x and y inputs: */
     717        xinput=(Input*)this->inputs->GetInput(XEnum);
     718        yinput=(Input*)this->inputs->GetInput(YEnum);
     719
     720        /*some checks: */
     721        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
     722        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
     723
     724        /*Scale: */
     725        yinput->AXPY(xinput,scalar);
     726}
     727/*}}}*/
     728/*FUNCTION Sing::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
     729void  Sing::ControlConstrain(int control_type, double cm_min, double cm_max){
     730
     731        Input* input=NULL;
     732
     733        /*Find input: */
     734        input=(Input*)this->inputs->GetInput(control_type);
     735       
     736        /*Do nothing if we  don't find it: */
     737        if(!input)return;
     738
     739        /*Constrain input using cm_min and cm_max: */
     740        input->Constrain(cm_min,cm_max);
     741
     742}
     743/*}}}*/
     744/*FUNCTION Sing::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
     745void  Sing::GetVectorFromInputs(Vec vector,int NameEnum){
     746
     747        const int numvertices=1;
     748        int doflist1[numvertices];
     749
     750        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
     751        for(i=0;i<this->inputs->Size();i++){
     752                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
     753                if(input->EnumType==NameEnum){
     754                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     755                        this->GetDofList1(&doflist1[0]);
     756                        input->GetVectorFromInputs(vector,&doflist1[0]);
     757                        break;
     758                }
     759        }
     760}
     761/*}}}*/

  • issm/trunk/src/c/objects/Elements/Sing.h

    r4047 r4048  
    9595                void  DuplicateInput(int original_enum,int new_enum);
    9696                void  ScaleInput(int enum_type,double scale_factor);
     97                void  AXPY(int YEnum, double scalar, int XEnum);
     98                void  ControlConstrain(int control_type,double cm_min, double cm_max);
    9799
    98100
    … …  
    111113                double CostFunction(void);
    112114                double MassFlux(double* segment);
     115                void  GetVectorFromInputs(Vec vector,int NameEnum);
    113116                /*}}}*/
    114117

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

    r4047 r4048  
    54505450}
    54515451/*}}}*/
     5452/*FUNCTION Tria::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
     5453void  Tria::AXPY(int YEnum, double scalar, int XEnum){
     5454
     5455        Input* xinput=NULL;
     5456        Input* yinput=NULL;
     5457
     5458        /*Find x and y inputs: */
     5459        xinput=(Input*)this->inputs->GetInput(XEnum);
     5460        yinput=(Input*)this->inputs->GetInput(YEnum);
     5461
     5462        /*some checks: */
     5463        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
     5464        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
     5465
     5466        /*Scale: */
     5467        yinput->AXPY(xinput,scalar);
     5468}
     5469/*}}}*/
     5470/*FUNCTION Tria::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
     5471void  Tria::ControlConstrain(int control_type, double cm_min, double cm_max){
     5472
     5473        Input* input=NULL;
     5474
     5475        /*Find input: */
     5476        input=(Input*)this->inputs->GetInput(control_type);
     5477       
     5478        /*Do nothing if we  don't find it: */
     5479        if(!input)return;
     5480
     5481        /*Constrain input using cm_min and cm_max: */
     5482        input->Constrain(cm_min,cm_max);
     5483
     5484}
     5485/*}}}*/
     5486/*FUNCTION Tria::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
     5487void  Tria::GetVectorFromInputs(Vec vector,int NameEnum){
     5488
     5489        const int numvertices=3;
     5490        int doflist1[numvertices];
     5491
     5492        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
     5493        for(i=0;i<this->inputs->Size();i++){
     5494                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
     5495                if(input->EnumType==NameEnum){
     5496                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     5497                        this->GetDofList1(&doflist1[0]);
     5498                        input->GetVectorFromInputs(vector,&doflist1[0]);
     5499                        break;
     5500                }
     5501        }
     5502}
     5503/*}}}*/

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

    r4047 r4048  
    141141                void  DuplicateInput(int original_enum,int new_enum);
    142142                void  ScaleInput(int enum_type,double scale_factor);
     143                void  AXPY(int YEnum, double scalar, int XEnum);
     144                void  ControlConstrain(int control_type,double cm_min, double cm_max);
     145                void  GetVectorFromInputs(Vec vector,int NameEnum);
    143146
    144147

  • issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp

    r4047 r4048  
    259259}
    260260/*}}}*/
     261/*FUNCTION BeamVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
     262void BeamVertexInput::AXPY(Input* xinput,double scalar){
     263
     264        int i;
     265        const int numgrids=2;
     266        BeamVertexInput*  xbeamvertexinput=NULL;
     267
     268        /*xinput is of the same type, so cast it: */
     269        xbeamvertexinput=(BeamVertexInput)xinput;
     270
     271        /*Carry out the AXPY operation:*/
     272        for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xbeamvertexinput->values[i];
     273
     274}
     275/*}}}*/
     276/*FUNCTION BeamVertexInput::Constrain(double cm_min, double cm_max){{{1*/
     277void BeamVertexInput::Constrain(double cm_min, double cm_max){
     278
     279        int i;
     280        const int numgrids=2;
     281               
     282        if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
     283        if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
     284
     285}
     286/*}}}*/
     287/*FUNCTION BeamVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     288void BeamVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
     289
     290        const int numvertices=2;
     291        VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
     292
     293
     294}
     295/*}}}*/

  • issm/trunk/src/c/objects/Inputs/BeamVertexInput.h

    r4047 r4048  
    7575                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    7676                void Scale(double scale_factor);
     77                void AXPY(Input* xinput,double scalar);
     78                void Constrain(double cm_min, double cm_max);
     79                void GetVectorFromInputs(Vec vector,int* doflist);
    7780                /*}}}*/
    7881

  • issm/trunk/src/c/objects/Inputs/BoolInput.cpp

    r4047 r4048  
    228228}
    229229/*}}}*/
     230/*FUNCTION BoolInput::AXPY(Input* xinput,double scalar);{{{1*/
     231void BoolInput::AXPY(Input* xinput,double scalar){
     232
     233        BoolInput*  xboolinput=NULL;
     234
     235        /*xinput is of the same type, so cast it: */
     236        xboolinput=(BoolInput)xinput;
     237
     238        /*Carry out the AXPY operation:*/
     239        this->value=this->value+scalar*xboolinput->value;
     240
     241}
     242/*}}}*/
     243/*FUNCTION BoolInput::Constrain(double cm_min, double cm_max){{{1*/
     244void BoolInput::Constrain(double cm_min, double cm_max){
     245
     246        if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
     247        if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
     248
     249}
     250/*}}}*/
     251/*FUNCTION BoolInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     252void BoolInput::GetVectorFromInputs(Vec vector,int* doflist){
     253
     254        ISSMERROR(" not supporte yet!");
     255
     256}
     257/*}}}*/

  • issm/trunk/src/c/objects/Inputs/BoolInput.h

    r4047 r4048  
    7575                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    7676                void Scale(double scale_factor);
     77                void AXPY(Input* xinput,double scalar);
     78                void Constrain(double cm_min, double cm_max);
     79                void GetVectorFromInputs(Vec vector,int* doflist);
    7780                /*}}}*/
    7881

  • issm/trunk/src/c/objects/Inputs/DoubleInput.cpp

    r4047 r4048  
    239239}
    240240/*}}}*/
     241/*FUNCTION DoubleInput::AXPY(Input* xinput,double scalar);{{{1*/
     242void DoubleInput::AXPY(Input* xinput,double scalar){
     243
     244        DoubleInput*  xdoubleinput=NULL;
     245
     246        /*xinput is of the same type, so cast it: */
     247        xdoubleinput=(DoubleInput)xinput;
     248
     249        /*Carry out the AXPY operation:*/
     250        this->value=this->value+scalar*xdoubleinput->value;
     251
     252}
     253/*}}}*/
     254/*FUNCTION DoubleInput::Constrain(double cm_min, double cm_max){{{1*/
     255void DoubleInput::Constrain(double cm_min, double cm_max){
     256
     257        if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
     258        if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
     259
     260}
     261/*}}}*/
     262/*FUNCTION DoubleInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     263void DoubleInput::GetVectorFromInputs(Vec vector,int* doflist){
     264
     265        ISSMERROR(" not supporte yet!");
     266
     267}
     268/*}}}*/

  • issm/trunk/src/c/objects/Inputs/DoubleInput.h

    r4047 r4048  
    7676                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    7777                void Scale(double scale_factor);
     78                void AXPY(Input* xinput,double scalar);
     79                void Constrain(double cm_min, double cm_max);
     80                void GetVectorFromInputs(Vec vector,int* doflist);
    7881                /*}}}*/
    7982

  • issm/trunk/src/c/objects/Inputs/Input.h

    r4047 r4048  
    5050                virtual Input* SpawnTriaInput(int* indices)=0;
    5151                virtual Result* SpawnResult(int step, double time)=0;
    52                 virtual void   SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0;
    53                 virtual void   Scale(double scale_factor)=0;
     52                virtual void SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0;
     53                virtual void Scale(double scale_factor)=0;
     54                virtual void AXPY(Input* xinput,double scalar)=0;
     55                virtual void Constrain(double cm_min, double cm_max)=0;
     56                virtual void GetVectorFromInputs(Vec vector,int* doflist)=0;
    5457                /*}}}*/
    5558

  • issm/trunk/src/c/objects/Inputs/IntInput.cpp

    r4047 r4048  
    226226}
    227227/*}}}*/
     228/*FUNCTION IntInput::AXPY(Input* xinput,int scalar);{{{1*/
     229void IntInput::AXPY(Input* xinput,int scalar){
     230
     231        IntInput*  xintinput=NULL;
     232
     233        /*xinput is of the same type, so cast it: */
     234        xintinput=(IntInput)xinput;
     235
     236        /*Carry out the AXPY operation:*/
     237        this->value=this->value+scalar*xintinput->value;
     238
     239}
     240/*}}}*/
     241/*FUNCTION IntInput::Constrain(double cm_min, double cm_max){{{1*/
     242void IntInput::Constrain(double cm_min, double cm_max){
     243
     244        if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
     245        if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
     246
     247}
     248/*}}}*/
     249/*FUNCTION IntInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     250void IntInput::GetVectorFromInputs(Vec vector,int* doflist){
     251
     252        ISSMERROR(" not supporte yet!");
     253
     254}
     255/*}}}*/

  • issm/trunk/src/c/objects/Inputs/IntInput.h

    r4047 r4048  
    7575                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    7676                void Scale(double scale_factor);
     77                void AXPY(Input* xinput,double scalar);
     78                void Constrain(double cm_min, double cm_max);
     79                void GetVectorFromInputs(Vec vector,int* doflist);
    7780                /*}}}*/
    7881

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp

    r4047 r4048  
    908908}
    909909/*}}}*/
     910/*FUNCTION PentaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
     911void PentaVertexInput::AXPY(Input* xinput,double scalar){
     912
     913        int i;
     914        const int numgrids=6;
     915        PentaVertexInput*  xpentavertexinput=NULL;
     916
     917        /*xinput is of the same type, so cast it: */
     918        xpentavertexinput=(PentaVertexInput)xinput;
     919
     920        /*Carry out the AXPY operation:*/
     921        for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
     922
     923}
     924/*}}}*/
     925/*FUNCTION PentaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
     926void PentaVertexInput::Constrain(double cm_min, double cm_max){
     927
     928        int i;
     929        const int numgrids=6;
     930               
     931        if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
     932        if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
     933
     934}
     935/*}}}*/
     936/*FUNCTION PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     937void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
     938
     939        const int numvertices=6;
     940        VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
     941
     942
     943}
     944/*}}}*/

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.h

    r4047 r4048  
    8484                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    8585                void Scale(double scale_factor);
     86                void AXPY(Input* xinput,double scalar);
     87                void Constrain(double cm_min, double cm_max);
     88                void GetVectorFromInputs(Vec vector,int* doflist);
    8689                /*}}}*/
    8790

  • issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp

    r4047 r4048  
    229229}
    230230/*}}}*/
     231/*FUNCTION SingVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
     232void SingVertexInput::AXPY(Input* xinput,double scalar){
     233
     234        SingVertexInput*  xsingvertexinput=NULL;
     235
     236        /*xinput is of the same type, so cast it: */
     237        xsingvertexinput=(SingVertexInput)xinput;
     238
     239        /*Carry out the AXPY operation:*/
     240        this->value=this->value+scalar*xngvertexinput->value;
     241
     242}
     243/*}}}*/
     244/*FUNCTION SingVertexInput::Constrain(double cm_min, double cm_max){{{1*/
     245void SingVertexInput::Constrain(double cm_min, double cm_max){
     246
     247        if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
     248        if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
     249
     250}
     251/*}}}*/
     252/*FUNCTION SingVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     253void SingVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
     254
     255        const int numvertices=1;
     256        VecSetValues(vector,numvertices,doflist,(const double*)&this->value,ADD_VALUES);
     257
     258
     259}
     260/*}}}*/

  • issm/trunk/src/c/objects/Inputs/SingVertexInput.h

    r4047 r4048  
    7474                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    7575                void Scale(double scale_factor);
     76                void AXPY(Input* xinput,double scalar);
     77                void Constrain(double cm_min, double cm_max);
     78                void GetVectorFromInputs(Vec vector,int* doflist);
    7679                /*}}}*/
    7780

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r4047 r4048  
    482482}
    483483/*}}}*/
     484/*FUNCTION TriaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
     485void TriaVertexInput::AXPY(Input* xinput,double scalar){
     486
     487        int i;
     488        const int numgrids=3;
     489        TriaVertexInput*  xtriavertexinput=NULL;
     490
     491        /*xinput is of the same type, so cast it: */
     492        xtriavertexinput=(TriaVertexInput)xinput;
     493
     494        /*Carry out the AXPY operation:*/
     495        for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xtriavertexinput->values[i];
     496
     497}
     498/*}}}*/
     499/*FUNCTION TriaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
     500void TriaVertexInput::Constrain(double cm_min, double cm_max){
     501
     502        int i;
     503        const int numgrids=3;
     504               
     505        if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
     506        if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
     507
     508}
     509/*}}}*/
     510/*FUNCTION TriaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
     511void TriaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
     512
     513        const int numvertices=3;
     514        VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
     515
     516
     517}
     518/*}}}*/

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

    r4047 r4048  
    8181                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
    8282                void Scale(double scale_factor);
     83                void AXPY(Input* xinput,double scalar);
     84                void Constrain(double cm_min, double cm_max);
     85                void GetVectorFromInputs(Vec vector,int* doflist);
    8386                /*}}}*/
    8487

  • issm/trunk/src/c/objects/OptArgs.h

    r3681 r4048  
    2525class Model;
    2626struct OptArgs{
    27         Model* model;
    28         double* param_g;
    29         double* grad_g;
     27        FemModel* femmodel;
    3028        int n;
    3129};

  • issm/trunk/src/c/solutions/ControlRestart.cpp

    r3938 r4048  
    1 
    21/*!\file: ControlRestart.cpp
    3  * \brief: core of the control solution
     2 * \brief: save as much as possible, to be able to restart the control_core solution
    43 */
    54
    … …  
    87#include "../EnumDefinitions/EnumDefinitions.h"
    98
    10 void ControlRestart(Model* model,double* param_g){
     9void ControlRestart(FemModel* femmodel){
    1110
    12         extern int my_rank;
    13 
    14         /*output: */
    15         Results* results=NULL;
    1611        char*    outputfilename=NULL;
    17 
    18         /*Intermediary: */
    19         int      i;
    20         int      numberofnodes;
    21         double* param_g_copy;
    22 
    23         /*Recover parameters used throughout the solution:*/
    24         model->FindParam(&numberofnodes,NumberOfNodesEnum);
     12       
     13        /*retrieve output file name: */
    2514        model->FindParam(&outputfilename,OutputFileNameEnum);
    2615
    27         /*Plug COPYS of the results into output dataset:
    28          * only the pointer is given to temporary_results and at the
    29          * end of ProcessResultsx the pointer is deleted. That would
    30          * destroy param_g*/
     16        /*we essentially want J and the parameter: */
     17        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
     18        femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
     19        femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
    3120
    32         param_g_copy=(double*)xcalloc(numberofnodes,sizeof(double));
    33         for(i=0;i<numberofnodes;i++) param_g_copy[i]=param_g[i];
    34 
    35         results=new Results();
    36         results->AddObject(new Result(results->Size()+1,0,1,"param_g",param_g_copy,numberofnodes));
    37         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(DiagnosticAnalysisEnum)));
    38 
    39         //Write results on disk
    40         OutputResults(results,outputfilename);
     21        /*write to disk: */
     22        OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
    4123       
    4224        /*Free ressources:*/
    43         delete results;
    4425        xfree((void**)&outputfilename);
    45         xfree((void**)&param_g_copy);
    4626}

  • issm/trunk/src/c/solutions/control_core.cpp

    r4047 r4048  
    1414void control_core(FemModel* femmodel){
    1515
    16         Vec     u_g=NULL;
    17         Vec     t_g=NULL;
    18         Vec     m_g=NULL;
    19         double  search_scalar;
     16        int     i,n;
     17       
     18        /*parameters: */
     19        int     verbose=0;
    2020        int     control_type;
     21        int     control_steady;
     22        int     nsteps;
     23        double  eps_cm;
     24        double  tolx;
     25        double  cm_min;
     26        double  cm_max;
     27        int     cm_gradient;
     28
    2129        double* fit=NULL;
    2230        double* optscal=NULL;
    2331        double* maxiter=NULL;
    2432        double* cm_jump=NULL;
    25         double  eps_cm;
    26         double  tolx;
    27         double* param_g=NULL;
    28         Vec     grad_g=NULL;
    29         Vec     new_grad_g=NULL;
    30         Vec     grad_g_old=NULL;
    31         double* grad_g_double=NULL;
    32         double  cm_min;
    33         double  cm_max;
    34         int     cm_gradient;
    35         int     nsteps,n,i;
    36         double* J=NULL;
     33
     34               
     35        /*intermediary: */
     36        double  search_scalar;
    3737        OptArgs optargs;
    3838        OptPars optpars;
    39         Param*  param=NULL;
    4039
    41         /*flags: */
    42         int analysis_type;
    43         int sub_analysis_type;
    44         int     control_steady;
    45         int verbose=0;
    46         int converged=0;
    47         int numberofnodes;
     40        /*output: */
     41        double* J=NULL;
    4842
    4943        /*Process models*/
    … …  
    6256        femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
    6357        femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
    64         femmodel->parameters->FindParam(&param_g,NULL,NULL,ControlParameterEnum);
    65         femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    66         femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
    67         femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
    6858        femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
    6959
    7060        /*Initialize misfit: */
    7161        J=(double*)xmalloc(nsteps*sizeof(double));
    72 
     62               
     63        /*Initialize some of the BrentSearch arguments: */
     64        optargs.femmodel=femmodel;
     65        optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ;
     66       
    7367        /*Start looping: */
    7468        for(n=0;n<nsteps;n++){
    7569
    7670                _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    77                 inputs->Add(control_type,param_g,1,numberofnodes);
    7871                femmodel->UpdateInputsFromConstant(fit[n],FitEnum);
    7972               
    80                 /*In case we are running a steady state control method, compute new temperature field using new parameter
    81                  * distribution: */
     73                /*In case we are running a steady state control method, compute new temperature field using new parameter * distribution: */
    8274                if (control_steady) steadystate_core(model);
    8375       
    … …  
    8779                /*Return gradient if asked: */
    8880                if (cm_gradient){
    89                         /*Transfer gradient from input to results: */
    9081                        InputToResultx(femodel->elements,femodel->nodes,femodel->vertices,femodel->loads,femodel->materials,femodel->parameters,GradientEnum);
    91                         return;
     82                        goto cleanup_and_return;
    9283                }
    9384
    94                 /*Normalize if last gradient not satisfying (search_scalar==0)*/
    95                 if (n>0 && search_scalar==0){
    96                         _printf_("%s","      orthogonalization...");
    97                         Orthx(&new_grad_g,grad_g,grad_g_old);
    98                         _printf_("%s\n"," done.");
    99                 }
    100                 else{
    101                         _printf_("%s","      normalizing directions...");
    102                         Orthx(&new_grad_g,grad_g,NULL);
    103                         _printf_("%s\n"," done.");
    104                 }
    105                 VecFree(&grad_g); VecFree(&grad_g_old);
    106                 grad_g_old=new_grad_g;
    107                 VecToMPISerial(&grad_g_double,new_grad_g);
    108 
    10985                _printf_("%s\n","      optimizing along gradient direction");
    110                 optargs.model=model;
    111                 optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.n=n;
    112                 optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];optpars.cm_jump=cm_jump[n];
     86                optargs.n=n; optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
    11387                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
    11488
    11589                _printf_("%s","      updating parameter using optimized search scalar...");
    116                 for(i=0;i<numberofnodes;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
    117                 _printf_("%s\n"," done.");
     90                AXPYInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],ControlParameterEnum);
    11891
    11992                _printf_("%s","      constraining the new distribution...");   
    120                 ControlConstrainx(param_g,numberofnodes,cm_min,cm_max,control_type);
    121                 _printf_("%s\n"," done.");
     93                ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
     94               
     95                _printf_("%s","      save new parameter...");
     96                DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum);
    12297               
    12398                _printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
     99               
     100                if(controlconvergence)goto cleanup_and_return;
    124101
    125                 /*some freeing:*/
    126                 xfree((void**)&grad_g_double);
    127 
    128                 /*Has convergence been reached?*/
    129                 if (!isnan(eps_cm)){
    130                         i=n-2;
    131                         //go through the previous misfits(starting from n-2)
    132                         while(i>=0){
    133                                 if (fit[i]==fit[n]){
    134                                         //convergence test only if we have the same misfits
    135                                         if ((J[i]-J[n])/J[n] <= eps_cm){
    136                                                 //convergence if convergence criteria fullfilled
    137                                                 converged=1;
    138                                                 _printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],"<",eps_cm);
    139                                         }
    140                                         else{
    141                                                 _printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],">",eps_cm);
    142                                         }
    143                                         break;
    144                                 }
    145                                 i=i-1;
    146                         }
    147                 }
    148                 //stop if convergence has been reached
    149                 if(converged) break;
    150 
    151                 //some temporary saving
     102                /*Temporary saving every 5 control steps: */
    152103                if (((n+1)%5)==0){
    153104                        _printf_("%s","      saving temporary results...");
    154                         ControlRestart(model,param_g);
    155                         _printf_("%s\n"," done.");
     105                        ControlRestart(femmodel);
    156106                }
    157107        }
    158108
    159         /*Write results to disk: */
    160109        _printf_("%s","      preparing final velocity solution");
    161         /*Launch diagnostic with the last parameter distribution*/
    162         if (control_steady){
    163                 model->UpdateInputsFromVector(param_g,control_type,VertexEnum);
    164                 steadystate_results=steadystate_core(model);
     110        if (control_steady) steadystate_core(femmodel);
     111        else diagnostic_core(femmodel);
    165112
    166                 //extract u_g ,t_g and m_g from steadystate results, and erase diagnostic_results;
    167                 steadystate_results->FindResult(&u_g,"u_g");
    168                 steadystate_results->FindResult(&m_g,"m_g");
    169                 steadystate_results->FindResult(&t_g,"t_g");
    170                 delete steadystate_results;
    171         }
    172         else{
    173                 model->UpdateInputsFromVector(param_g,control_type,VertexEnum);
    174                 diagnostic_results=diagnostic_core(model);
     113        /*some results not computed by steadystate_core or diagnostic_core: */
     114        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
     115        femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
     116        femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
    175117
    176                 //extract u_g from diagnostic_results, and erase diagnostic_results;
    177                 diagnostic_results->FindResult(&u_g,"u_g");
    178                 delete diagnostic_results;
    179         }
    180 
    181         /*Plug results into output dataset: */
    182         results->AddObject(new Result(results->Size()+1,0,1,"u_g",u_g));
    183         results->AddObject(new Result(results->Size()+1,0,1,"param_g",param_g,numberofnodes));
    184         results->AddObject(new Result(results->Size()+1,0,1,"J",J,nsteps));
    185         if (control_steady){
    186                 results->AddObject(new Result(results->Size()+1,0,1,"t_g",t_g));
    187                 results->AddObject(new Result(results->Size()+1,0,1,"m_g",m_g));
    188         }
     118        cleanup_and_return:
    189119       
    190         /*Add analysis_type and control_type to results: */
    191         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(analysis_type)));
    192         results->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
    193 
    194120        /*Free ressources: */
    195121        xfree((void**)&control_type);
    … …  
    198124        xfree((void**)&maxiter);
    199125        xfree((void**)&cm_jump);
    200         VecFree(&new_grad_g); //do not VecFree grad_g and grad_g_old, they point to new_grad_g
    201         xfree((void**)&grad_g_double);
    202         xfree((void**)&param_g);
    203         VecFree(&u_g);
    204         VecFree(&t_g);
    205         VecFree(&m_g);
    206126        xfree((void**)&J);
    207 
    208         //return:
    209         return results;
    210127}

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

    r4047 r4048  
    1313
    1414
    15 void gradient_core(FemModel* femmodel){
    16        
     15void gradient_core(FemModel* femmodel,int step, double search_scalar){ //step defaults to 0, search_scalar defaults to 0
    1716       
    1817        /*parameters: */
    … …  
    2019        int control_type;
    2120        int verbose;
     21        Vec new_gradient=NULL;
     22        Vec gradient=NULL;
     23        Vec old_gradient=NULL;
    2224
    2325        /*retrieve parameters:*/
    … …  
    3032       
    3133        if(verbose)_printf_("%s\n","      compute gradient:");
    32         Gradjx( femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type);
     34        Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type,AdjointEnum);
     35       
     36        if(verbose)_printf_("%s\n","      retrieve old gradient:");
     37        GetVectorFromInputsx(&old_gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, OldGradientEnum,VertexEnum);
    3338
    3439        if(control_steady)diagnostic_core(model);
     40       
     41        if (step>0 && search_scalar==0){
     42                _printf_("%s","      orthogonalization...");
     43                Orthx(&new_gradient,gradient,old_gradient);
     44        }
     45        else{
     46                _printf_("%s","      normalizing directions...");
     47                Orthx(&new_gradient,gradient,NULL);
     48        }
     49        _printf_("%s\n"," done.");
    3550
     51       
     52        /*point gradient and old_gradient to new_gradient: */
     53        VecFree(&gradient); gradient=new_gradient;
     54        VecFree(&old_gradient); old_gradient=new_gradient;
     55
     56        /*plug back into inputs: */
     57        UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmoel->parameters,gradient,GradientEnum,VertexEnum);
     58        UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmoel->parameters,old_gradient,OldGradientEnum,VertexEnum);
     59
     60        /*Free ressources and return:*/
     61        VecFree(&new_gradient);
    3662}

  • issm/trunk/src/c/solutions/objectivefunctionC.cpp

    r3938 r4048  
    2121       
    2222        /*parameters: */
    23         Model*    model=NULL;
    2423        FemModel* femmodel=NULL;
    25         Results* diagnostic_results=NULL;
    26         double* param_g=NULL;
    27         double* grad_g=NULL;
    28         int numberofdofspernode;
    2924        int n;
    30 
    31         /*intermediary:*/
    32         int gsize;
    3325        double* optscal=NULL;
    3426        double* fit=NULL;
    3527        double  cm_min;
    3628        double  cm_max;
    37         int   control_type;
    38         double* param_g_copy=NULL;
     29        int     control_type;
    3930        int     analysis_type;
    4031        int     sub_analysis_type;
    4132        int     control_steady;
    42         Vec     u_g=NULL;
    43         Vec     u_g_full=NULL;
    44         double* u_g_double=NULL;
    45         double* vx=NULL;
    46         double* vy=NULL;
    47         double* vz=NULL;
    48         int     numberofnodes;
    4933               
    50         /*steadystate: */
    51         int     dt=0;
    52         int     isstokes=0;
    53         Results* results_steadystate=NULL;
    54         int dofs01[2]={0,1};
    55         double* dofset=NULL;
    56 
    57         /*Recover active model: */
    58         model=optargs->model;
    59         femmodel=model->GetActiveFormulation();
     34        /*Recover finite element model: */
     35        femmodel=optargs->femmodel;
    6036
    6137        /*Recover parameters: */
    62         param_g=optargs->param_g;
    63         grad_g=optargs->grad_g;
    6438        n=optargs->n;
    6539
    66         gsize=femmodel->nodesets->GetGSize();
    6740        femmodel->parameters->FindParam(&optscal,NULL,NULL,OptScalEnum);
    6841        femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum);
    … …  
    7346        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    7447        femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
    75         femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    76         femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
    77         femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
    7848
    79         /*First copy param_g so we don't modify it: */
    80         param_g_copy=(double*)xmalloc(numberofnodes*sizeof(double));
    81         memcpy(param_g_copy,param_g,numberofnodes*sizeof(double));
    82 
    83         /*First, update param_g using search_scalar: */
    84         for(i=0;i<numberofnodes;i++)param_g_copy[i]=param_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
     49        /*Use ControlParameterEnum input to  reinitialize our input parameter: */
     50        DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type);
     51       
     52        /*Use search scalar to shoot parameter in the gradient direction: */
     53        AXPYInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],ControlParameterEnum);
    8554
    8655        /*Constrain:*/
    87         ControlConstrainx(param_g_copy,numberofnodes,cm_min,cm_max,control_type);
     56        ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
    8857
    89         /*Add new parameter to inputs: */
    90         UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,param_g_copy,control_type,VertexEnum);
    91 
    92         /*Run diagnostic with updated parameters.*/
    93         if(!control_steady){
    94                 diagnostic_core_nonlinear(&u_g,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,sub_analysis_type);
    95                 UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,u_g,DiagnosticAnalysisEnum,sub_analysis_type);
    96                 VecFree(&u_g);
    97         }
    98         else{
    99                 //We need a 3D velocity!! (vz is required for the next thermal run)
    100                 diagnostic_results=     diagnostic_core(model);
    101 
    102                 //extract u_g and add it to input (3d velocity needed by thermal_core)
    103                 diagnostic_results->FindResult(&u_g,"u_g");
    104                
    105                 SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz);
    106                 UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vx,VxEnum,VertexEnum);
    107                 UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vy,VyEnum,VertexEnum);
    108                 UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vz,VzEnum,VertexEnum);
    109                
    110                 delete diagnostic_results;
    111         }
     58        /*Run diagnostic with updated inputs: */
     59        if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,true,femmodel,DiagnosticAnalysisEnum,sub_analysis_type);  //true means we conserve loads at each diagnostic run
     60        else                diagnostic_core(femmodel);  //We need a 3D velocity!! (vz is required for the next thermal run)
    11261
    11362        /*Compute misfit for this velocity field.*/
    … …  
    11867        xfree((void**)&fit);
    11968        xfree((void**)&optscal);
    120         xfree((void**)&param_g_copy);
    121         xfree((void**)&dofset);
    122         xfree((void**)&vx);
    123         xfree((void**)&vy);
    124         xfree((void**)&vz);
    125 
    12669        return J;
    12770}

  • issm/trunk/src/c/solutions/solutions.h

    r4047 r4048  
    2525Results* transient_core_3d(FemModel* model);
    2626void adjoint_core(FemModel* model);
    27 void gradient_core(FemModel* model);
     27void gradient_core(FemModel* model,int step=0, double search_scalar=0);
    2828void diagnostic_core(FemModel* model);
    2929void thermal_core(FemModel* model);
    … …  
    3333//int GradJOrth(WorkspaceParams* workspaceparams);
    3434void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters);
     35int controlconvergence(double* J, double* fit, double eps_cm, int n);
    3536
    3637int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);
    … …  
    4748
    4849void ControlInitialization(FemModel* model);
    49 void ControlRestart(FemModel* model,double* param_g);
     50void ControlRestart(FemModel* femmodel);
    5051
    5152void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
Note: See TracChangeset for help on using the changeset viewer.