Changeset 3169


Ignore:
Timestamp:
03/03/10 16:21:37 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added Cost function

Location:
issm/trunk/src/c
Files:
3 added
17 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/DataSet/DataSet.cpp

    r3045 r3169  
    14861486}               
    14871487               
    1488 void  DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type,int real){
     1488void  DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){
    14891489
    14901490        double J=0;;
    … …  
    14981498
    14991499                        element=(Element*)(*object);
    1500                         J+=element->Misfit(inputs,analysis_type,sub_analysis_type,real);
     1500                        J+=element->Misfit(inputs,analysis_type,sub_analysis_type);
    15011501
    15021502                }
    … …  
    15081508}
    15091509
     1510void  DataSet::CostFunction(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){
     1511
     1512        double J=0;;
     1513
     1514        vector<Object*>::iterator object;
     1515        Element* element=NULL;
     1516
     1517        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1518
     1519                if(EnumIsElement((*object)->Enum())){
     1520
     1521                        element=(Element*)(*object);
     1522                        J+=element->CostFunction(inputs,analysis_type,sub_analysis_type);
     1523
     1524                }
     1525        }
     1526
     1527        /*Assign output pointers:*/
     1528        *pJ=J;
     1529
     1530}
     1531
    15101532void  DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){
    15111533

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

    r3045 r3169  
    8181                void  Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type);
    8282                void  Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
    83                 void  Misfit(double* pJ, void* inputs,int analysis_type,int sub_analysis_type,int real);
     83                void  Misfit(double* pJ, void* inputs,int analysis_type,int sub_analysis_type);
     84                void  CostFunction(double* pJ, void* inputs,int analysis_type,int sub_analysis_type);
    8485                void  FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname);
    8586                int   DeleteObject(Object* object);

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

    r3159 r3169  
    235235                                        ./Misfitx/Misfitx.h\
    236236                                        ./Misfitx/Misfitx.cpp\
     237                                        ./CostFunctionx/CostFunctionx.h\
     238                                        ./CostFunctionx/CostFunctionx.cpp\
    237239                                        ./Orthx/Orthx.h\
    238240                                        ./Orthx/Orthx.cpp\
    … …  
    561563                                        ./Misfitx/Misfitx.h\
    562564                                        ./Misfitx/Misfitx.cpp\
     565                                        ./CostFunctionx/CostFunctionx.h\
     566                                        ./CostFunctionx/CostFunctionx.cpp\
    563567                                        ./Orthx/Orthx.h\
    564568                                        ./Orthx/Orthx.cpp\

  • issm/trunk/src/c/Misfitx/Misfitx.cpp

    r3045 r3169  
    1414
    1515void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,DataSet* parameters,
    16                         ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int real){
     16                        ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1717       
    1818        /*output: */
    … …  
    2525
    2626        /*Compute gradients: */
    27         elements->Misfit(&J,inputs,analysis_type,sub_analysis_type,real);
     27        elements->Misfit(&J,inputs,analysis_type,sub_analysis_type);
    2828
    2929        /*Sum all J from all cpus of the cluster:*/

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

    r3045 r3169  
    1010/* local prototypes: */
    1111void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, DataSet* parameters,
    12                         ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int real);
     12                        ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1313
    1414#endif  /* _MISFITX_H */

  • issm/trunk/src/c/Qmux/DakotaResponses.cpp

    r3057 r3169  
    264264
    265265                        /*Compute misfit: */
    266                         Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type,1);
     266                        Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type);
    267267                       
    268268

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

    r3127 r3169  
    5353#include "./Orthx/Orthx.h"
    5454#include "./Misfitx/Misfitx.h"
     55#include "./CostFunctionx/CostFunctionx.h"
    5556#include "./ControlConstrainx/ControlConstrainx.h"
    5657#include "./FieldDepthAveragex/FieldDepthAveragex.h"

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

    r3159 r3169  
    233233}
    234234/*}}}*/
     235/*FUNCTION Beam CostFunction{{{1*/
     236#undef __FUNCT__
     237#define __FUNCT__ "Beam::CostFunction"
     238double Beam::CostFunction(void*,int,int){
     239        throw ErrorException(__FUNCT__," not supported yet!");
     240}
     241/*}}}*/
    235242/*FUNCTION Beam CreateKMatrix{{{1*/
    236243#undef __FUNCT__
    … …  
    666673#undef __FUNCT__
    667674#define __FUNCT__ "Beam::Misfit"
    668 double Beam::Misfit(void*,int,int,int){
     675double Beam::Misfit(void*,int,int){
    669676        throw ErrorException(__FUNCT__," not supported yet!");
    670677}

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

    r3041 r3169  
    8787                void  GradjDrag(_p_Vec*, void*, int,int );
    8888                void  GradjB(_p_Vec*, void*, int,int );
    89                 double Misfit(void*,int,int,int);
     89                double Misfit(void*,int,int);
     90                double CostFunction(void*,int,int);
    9091                void  GetNodalFunctions(double* l1l2, double gauss_coord);
    9192                void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);

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

    r3041 r3169  
    3838                virtual void  GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type)=0;
    3939                virtual void  GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type)=0;
    40                 virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real)=0;
     40                virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type)=0;
     41                virtual double CostFunction(void* inputs,int analysis_type,int sub_analysis_type)=0;
    4142                virtual void  ComputePressure(Vec p_g)=0;
    4243                virtual double MassFlux(double* segment,double* ug)=0;

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

    r3161 r3169  
    282282Object* Penta::copy() {
    283283        return new Penta(*this);
     284}
     285/*}}}*/
     286/*FUNCTION CostFunction {{{1*/
     287#undef __FUNCT__
     288#define __FUNCT__ "Penta::CostFunction"
     289double Penta::CostFunction(void* inputs,int analysis_type,int sub_analysis_type){
     290
     291        double J;
     292        Tria* tria=NULL;
     293
     294        /*If on water, return 0: */
     295        if(onwater)return 0;
     296
     297        /*Bail out if this element if:
     298         * -> Non collapsed and not on the surface
     299         * -> collapsed (2d model) and not on bed) */
     300        if ((!collapse && !onsurface) || (collapse && !onbed)){
     301                return 0;
     302        }
     303        else if (collapse){
     304
     305                /*This element should be collapsed into a tria element at its base. Create this tria element,
     306                 * and compute CostFunction*/
     307                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     308                J=tria->CostFunction(inputs,analysis_type,sub_analysis_type);
     309                delete tria;
     310                return J;
     311        }
     312        else{
     313
     314                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     315                J=tria->CostFunction(inputs,analysis_type,sub_analysis_type);
     316                delete tria;
     317                return J;
     318        }
    284319}
    285320/*}}}*/
    … …  
    39153950#undef __FUNCT__
    39163951#define __FUNCT__ "Penta::Misfit"
    3917 double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real){
     3952double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){
    39183953
    39193954        double J;
    … …  
    39343969                 * and compute Misfit*/
    39353970                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    3936                 J=tria->Misfit(inputs,analysis_type,sub_analysis_type,real);
     3971                J=tria->Misfit(inputs,analysis_type,sub_analysis_type);
    39373972                delete tria;
    39383973                return J;
    … …  
    39413976
    39423977                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    3943                 J=tria->Misfit(inputs,analysis_type,sub_analysis_type,real);
     3978                J=tria->Misfit(inputs,analysis_type,sub_analysis_type);
    39443979                delete tria;
    39453980                return J;

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

    r3041 r3169  
    8888                void  GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
    8989                void  GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
    90                 double Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real);
     90                double Misfit(void* inputs,int analysis_type,int sub_analysis_type);
     91                double CostFunction(void* inputs,int analysis_type,int sub_analysis_type);
    9192               
    9293                void          GetThicknessList(double* thickness_list);

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

    r3041 r3169  
    216216}
    217217/*}}}*/
     218/*FUNCTION Sing::CostFunction {{{1*/
     219#undef __FUNCT__
     220#define __FUNCT__ "Sing::CostFunction"
     221double Sing::CostFunction(void*, int,int){
     222        throw ErrorException(__FUNCT__," not supported yet!");
     223}
     224/*}}}*/
    218225/*FUNCTION Sing::CreateKMatrix {{{1*/
    219226#undef __FUNCT__
    … …  
    515522#undef __FUNCT__
    516523#define __FUNCT__ "Sing::Misfit"
    517 double Sing::Misfit(void*, int,int,int){
     524double Sing::Misfit(void*, int,int){
    518525        throw ErrorException(__FUNCT__," not supported yet!");
    519526}

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

    r3041 r3169  
    8282                void  GradjDrag(_p_Vec*, void*, int,int);
    8383                void  GradjB(_p_Vec*, void*, int,int);
    84                 double Misfit(void*,int,int,int);
     84                double Misfit(void*,int,int);
     85                double CostFunction(void*,int,int);
    8586                double MassFlux(double* segment,double* ug);
    8687

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

    r3160 r3169  
    268268}
    269269/*}}}*/
     270/*FUNCTION CostFunction {{{1*/
     271#undef __FUNCT__
     272#define __FUNCT__ "Tria::CostFunction"
     273double Tria::CostFunction(void* vinputs,int analysis_type,int sub_analysis_type){
     274
     275        int i;
     276
     277        /* output: */
     278        double Jelem;
     279
     280        /* node data: */
     281        const int    numgrids=3;
     282        const int    numdof=2*numgrids;
     283        const int    NDOF2=2;
     284        int          dofs1[1]={0};
     285        int          dofs2[2]={0,1};
     286        double       xyz_list[numgrids][3];
     287
     288        /* grid data: */
     289        double B[numgrids];
     290
     291        /* gaussian points: */
     292        int     num_gauss,ig;
     293        double* first_gauss_area_coord  =  NULL;
     294        double* second_gauss_area_coord =  NULL;
     295        double* third_gauss_area_coord  =  NULL;
     296        double* gauss_weights           =  NULL;
     297        double  gauss_weight;
     298        double  gauss_l1l2l3[3];
     299        double  k_gauss;
     300        double  B_gauss;
     301
     302        /* parameters: */
     303        double  dk[NDOF2];
     304        double  dB[NDOF2];
     305
     306        /* Jacobian: */
     307        double Jdet;
     308
     309        ParameterInputs* inputs=NULL;
     310
     311        /*If on water, return 0: */
     312        if(onwater)return 0;
     313
     314        /*recover pointers: */
     315        inputs=(ParameterInputs*)vinputs;
     316
     317        /* Get node coordinates and dof list: */
     318        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     319
     320        /*First, get Misfit*/
     321        Jelem=Misfit(inputs,analysis_type,sub_analysis_type);
     322
     323          /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     324          GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     325
     326#ifdef _ISSM_DEBUG_
     327        for (i=0;i<num_gauss;i++){
     328                printf("Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(gauss_weights+i));
     329        }
     330#endif
     331
     332        /* Start  looping on the number of gaussian points: */
     333        for (ig=0; ig<num_gauss; ig++){
     334                /*Pick up the gaussian point: */
     335                gauss_weight=*(gauss_weights+ig);
     336                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     337                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     338                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     339
     340                /* Get Jacobian determinant: */
     341                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     342
     343                /*Add Tikhonov regularization term to misfit*/
     344                if (strcmp(numpar->control_type,"drag")==0){
     345                        if (!shelf){
     346
     347                                GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
     348                                Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
     349
     350                        }
     351                }
     352                else if (strcmp(numpar->control_type,"B")==0){
     353
     354                        if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){
     355                                throw ErrorException(__FUNCT__,"parameter B not found in input");
     356                        }
     357                        GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
     358                        Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
     359
     360                }
     361                else{
     362                        throw ErrorException(__FUNCT__,exprintf("%s%s","unsupported control type: ",numpar->control_type));
     363                }
     364
     365        }
     366
     367        xfree((void**)&first_gauss_area_coord);
     368        xfree((void**)&second_gauss_area_coord);
     369        xfree((void**)&third_gauss_area_coord);
     370        xfree((void**)&gauss_weights);
     371
     372        /*Return: */
     373        return Jelem;
     374}
     375/*}}}*/
    270376/*FUNCTION CreateKMatrix {{{1*/
    271377#undef __FUNCT__
    … …  
    43184424#undef __FUNCT__
    43194425#define __FUNCT__ "Tria::Misfit"
    4320 double Tria::Misfit(void* vinputs,int analysis_type,int sub_analysis_type,int real){
     4426double Tria::Misfit(void* vinputs,int analysis_type,int sub_analysis_type){
    43214427
    43224428        int i;
    … …  
    43434449        double relative_list[numgrids];
    43444450        double logarithmic_list[numgrids];
    4345         double B[numgrids];
    43464451
    43474452        /* gaussian points: */
    … …  
    43534458        double  gauss_weight;
    43544459        double  gauss_l1l2l3[3];
    4355         double  k_gauss;
    4356         double  B_gauss;
    43574460
    43584461        /* parameters: */
    43594462        double  velocity_mag,obs_velocity_mag;
    43604463        double  absolute,relative,logarithmic;
    4361         double  dk[NDOF2];
    4362         double  dB[NDOF2];
    43634464
    43644465        /* Jacobian: */
    … …  
    44314532        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    44324533
    4433 #ifdef _ISSM_DEBUG_
    4434         for (i=0;i<num_gauss;i++){
    4435                 printf("Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(gauss_weights+i));
    4436         }
    4437 #endif
    4438 
    44394534        /* Start  looping on the number of gaussian points: */
    44404535        for (ig=0; ig<num_gauss; ig++){
    … …  
    44474542                /* Get Jacobian determinant: */
    44484543                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    4449 #ifdef _ISSM_DEBUG_
    4450                 printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
    4451 #endif
    4452 
    4453                 /*Add dampening terms to misfit*/
    4454                 if(!real){
    4455                         if (strcmp(numpar->control_type,"drag")==0){
    4456                                 if (!shelf){
    4457 
    4458                                         //noise dampening
    4459                                         GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
    4460                                         Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
    4461 
    4462                                 }
    4463                         }
    4464                         else if (strcmp(numpar->control_type,"B")==0){
    4465                                 if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){
    4466                                         throw ErrorException(__FUNCT__,"parameter B not found in input");
    4467                                 }
    4468                                 //noise dampening
    4469                                 GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
    4470                                 Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
    4471 
    4472                                 //min dampening
    4473                                 GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3);
    4474                                 if(B_gauss<numpar->cm_mindmp_value){
    4475                                         Jelem+=numpar->cm_mindmp_slope*B_gauss*Jdet*gauss_weight;
    4476                                 }
    4477 
    4478                                 //max dampening
    4479                                 if(B_gauss>numpar->cm_maxdmp_value){
    4480                                         Jelem+=numpar->cm_maxdmp_slope*B_gauss*Jdet*gauss_weight;
    4481                                 }
    4482                         }
    4483                         else{
    4484                                 throw ErrorException(__FUNCT__,exprintf("%s%s","unsupported control type: ",numpar->control_type));
    4485                         }
    4486                 }
    44874544
    44884545                /*Differents misfits are allowed: */
    … …  
    45094566                }
    45104567                else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!"));
    4511                        
    4512 
    4513         }
    4514 cleanup_and_return:
     4568
     4569        }
     4570
    45154571        xfree((void**)&first_gauss_area_coord);
    45164572        xfree((void**)&second_gauss_area_coord);

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

    r3041 r3169  
    9999                void  SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
    100100                void  GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
    101                 double Misfit(void* inputs,int analysis_type,int sub_analysis_type,int real);
     101                double Misfit(void* inputs,int analysis_type,int sub_analysis_type);
     102                double CostFunction(void* inputs,int analysis_type,int sub_analysis_type);
    102103
    103104                void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);

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

    r3045 r3169  
    9797        /*Compute misfit for this velocity field.*/
    9898        inputs->Add("fit",fit[n]);
    99         Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type,0);
     99        CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, femmodel->parameters,inputs,analysis_type,sub_analysis_type);
    100100
    101101        /*Free ressources:*/
Note: See TracChangeset for help on using the changeset viewer.