Changeset 5281


Ignore:
Timestamp:
08/16/10 11:37:04 (15 years ago)
Author:
Mathieu Morlighem
Message:

moved md.fit to md.cm_responses and 0-4 to enums

Location:
issm/trunk/src/c
Files:
1 deleted
24 edited

Legend:

Unmodified
Added
Removed

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

    r5251 r5281  
    170170        BedSlopeYEnum,
    171171        BoundaryEnum,
     172        CmResponseEnum,
     173        CmResponsesEnum,
    172174        CmMaxDmpSlopeEnum,
    173175        CmMaxDmpValueEnum,
    … …  
    190192        EpsVelEnum,
    191193        FillEnum,
    192         FitEnum,
    193194        FractionIncrementEnum,
    194195        FrictionEnum,
    … …  
    197198        KflagEnum,
    198199        MassFluxEnum,
     200        SurfaceAbsVelMisfitEnum,
     201        SurfaceRelVelMisfitEnum,
     202        SurfaceLogVelMisfitEnum,
     203        SurfaceLogVxVyMisfitEnum,
     204        SurfaceAverageVelMisfitEnum,
    199205        MaxPenetrationEnum,
    200206        MeanVelEnum,

  • issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

    r5251 r5281  
    144144                case BedSlopeYEnum : return "BedSlopeY";
    145145                case BoundaryEnum : return "Boundary";
     146                case CmResponseEnum : return "CmResponse";
     147                case CmResponsesEnum : return "CmResponses";
    146148                case CmMaxDmpSlopeEnum : return "CmMaxDmpSlope";
    147149                case CmMaxDmpValueEnum : return "CmMaxDmpValue";
    … …  
    164166                case EpsVelEnum : return "EpsVel";
    165167                case FillEnum : return "Fill";
    166                 case FitEnum : return "Fit";
    167168                case FractionIncrementEnum : return "FractionIncrement";
    168169                case FrictionEnum : return "Friction";
    … …  
    171172                case KflagEnum : return "Kflag";
    172173                case MassFluxEnum : return "MassFlux";
     174                case SurfaceAbsVelMisfitEnum : return "SurfaceAbsVelMisfit";
     175                case SurfaceRelVelMisfitEnum : return "SurfaceRelVelMisfit";
     176                case SurfaceLogVelMisfitEnum : return "SurfaceLogVelMisfit";
     177                case SurfaceLogVxVyMisfitEnum : return "SurfaceLogVxVyMisfit";
     178                case SurfaceAverageVelMisfitEnum : return "SurfaceAverageVelMisfit";
    173179                case MaxPenetrationEnum : return "MaxPenetration";
    174180                case MeanVelEnum : return "MeanVel";

  • issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

    r5251 r5281  
    142142        else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    143143        else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
     144        else if (strcmp(name,"CmResponse")==0) return CmResponseEnum;
     145        else if (strcmp(name,"CmResponses")==0) return CmResponsesEnum;
    144146        else if (strcmp(name,"CmMaxDmpSlope")==0) return CmMaxDmpSlopeEnum;
    145147        else if (strcmp(name,"CmMaxDmpValue")==0) return CmMaxDmpValueEnum;
    … …  
    162164        else if (strcmp(name,"EpsVel")==0) return EpsVelEnum;
    163165        else if (strcmp(name,"Fill")==0) return FillEnum;
    164         else if (strcmp(name,"Fit")==0) return FitEnum;
    165166        else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    166167        else if (strcmp(name,"Friction")==0) return FrictionEnum;
    … …  
    169170        else if (strcmp(name,"Kflag")==0) return KflagEnum;
    170171        else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
     172        else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
     173        else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
     174        else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
     175        else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     176        else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    171177        else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum;
    172178        else if (strcmp(name,"MeanVel")==0) return MeanVelEnum;
    … …  
    305311        else if (strcmp(name,"WaitOnLock")==0) return WaitOnLockEnum;
    306312        else if (strcmp(name,"Yts")==0) return YtsEnum;
    307         else ISSMERROR("Enum %s not found",name);
    308 
    309 }

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

    r5250 r5281  
    433433                                        ./modules/SurfaceAreax/SurfaceAreax.h\
    434434                                        ./modules/SurfaceAreax/SurfaceAreax.cpp\
    435                                         ./modules/Misfitx/Misfitx.h\
    436                                         ./modules/Misfitx/Misfitx.cpp\
     435                                        ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h\
     436                                        ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp\
     437                                        ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h\
     438                                        ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp\
     439                                        ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h\
     440                                        ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp\
     441                                        ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h\
     442                                        ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp\
     443                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\
     444                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\
    437445                                        ./modules/CostFunctionx/CostFunctionx.h\
    438446                                        ./modules/CostFunctionx/CostFunctionx.cpp\
    … …  
    968976                                        ./modules/SurfaceAreax/SurfaceAreax.h\
    969977                                        ./modules/SurfaceAreax/SurfaceAreax.cpp\
    970                                         ./modules/Misfitx/Misfitx.h\
    971                                         ./modules/Misfitx/Misfitx.cpp\
     978                                        ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h\
     979                                        ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp\
     980                                        ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h\
     981                                        ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp\
     982                                        ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h\
     983                                        ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp\
     984                                        ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h\
     985                                        ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp\
     986                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\
     987                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\
    972988                                        ./modules/CostFunctionx/CostFunctionx.h\
    973989                                        ./modules/CostFunctionx/CostFunctionx.cpp\

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

    r4974 r5281  
    99#include "../../toolkits/toolkits.h"
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    11 #include "../SurfaceAreax/SurfaceAreax.h"
     11#include "../Responsex/Responsex.h"
    1212
    13 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     13void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,int response){
    1414
    1515        /*Intermediary*/
    … …  
    2020
    2121        /*output: */
    22         double J=0;
    23         double J_sum;
     22        double J;
     23        double Jreg=0;
     24        double Jreg_sum;
     25       
     26        /*Get response*/
     27        Responsex(&J,elements,nodes,vertices,loads,materials,parameters,(const char*)EnumToString(response),false); //False means DO NOT process units
    2428
    25         /*Compute surface area: */
    26         SurfaceAreax(&S,elements,nodes,vertices, loads,materials,parameters);
    27 
    28         /*add surface area to elements :*/
    29         for(i=0;i<elements->Size();i++){
    30                 Element* element=(Element*)elements->GetObjectByOffset(i);
    31                 element->InputUpdateFromVector(&S,SurfaceAreaEnum,ConstantEnum);
    32         }
    33        
    34         /*Compute cost function: */
     29        /*Add Regularization terms: */
    3530        for (i=0;i<elements->Size();i++){
    3631                element=(Element*)elements->GetObjectByOffset(i);
    37                 J+=element->CostFunction(process_units);
     32                Jreg+=element->RegularizeInversion();
    3833        }
    3934
    4035        /*Sum all J from all cpus of the cluster:*/
    41         MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD );
    42         MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
    43         J=J_sum;
     36        MPI_Reduce (&Jreg,&Jreg_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD );
     37        MPI_Bcast(&Jreg_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
     38        J=J+Jreg_sum;
    4439
    4540        /*Assign output pointers: */

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

    r4974 r5281  
    1010
    1111/* local prototypes: */
    12 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int response);
    1313
    14 #endif  /* _MISFITX_H */
    15 
     14#endif

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

    r5250 r5281  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* descriptor,bool process_units){
     12void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* descriptor,bool process_units){
    1313
    1414        int i,j;

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

    r5250 r5281  
    1010
    1111/* local prototypes: */
    12 void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* descriptor,bool process_units);
     12void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* descriptor,bool process_units);
    1313
    1414

  • issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

    r5247 r5281  
    6060
    6161                /*Now, recover fit, optscal and maxiter as vectors: */
    62                 IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit");
     62                IoModelFetchData(&iomodel->cm_responses,NULL,NULL,iomodel_handle,"cm_responses");
    6363                IoModelFetchData(&iomodel->cm_jump,NULL,NULL,iomodel_handle,"cm_jump");
    6464                IoModelFetchData(&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal");
    6565                IoModelFetchData(&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter");
    6666
    67                 parameters->AddObject(new DoubleVecParam(FitEnum,iomodel->fit,iomodel->nsteps));
     67                parameters->AddObject(new DoubleVecParam(CmResponsesEnum,iomodel->cm_responses,iomodel->nsteps));
    6868                parameters->AddObject(new DoubleVecParam(CmJumpEnum,iomodel->cm_jump,iomodel->nsteps));
    6969                parameters->AddObject(new DoubleVecParam(OptScalEnum,iomodel->optscal,iomodel->nsteps));
    7070                parameters->AddObject(new DoubleVecParam(MaxIterEnum,iomodel->maxiter,iomodel->nsteps));
    7171
    72                 xfree((void**)&iomodel->fit);
     72                xfree((void**)&iomodel->cm_responses);
    7373                xfree((void**)&iomodel->cm_jump);
    7474                xfree((void**)&iomodel->optscal);

  • issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp

    r5240 r5281  
    207207                                qmu_mass_flux_present=true;
    208208                        }
    209 
    210                         if (strcmp(descriptor,"misfit")==0){
    211 
    212                                 /*We need the fit: */
    213                                 IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit");
    214                                 parameters->SetParam(iomodel->fit,iomodel->nsteps,FitEnum);
    215                                 xfree((void**)&iomodel->fit);
    216 
    217                         }
    218209                }
    219210               

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

    r5250 r5281  
    1616#include "../modules.h"
    1717
    18 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units){
     18void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* response_descriptor,bool process_units){
    1919
    2020        if(strcmp(response_descriptor,"MinVel")==0){
    … …  
    5151                MaxAbsVzx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
    5252        }
    53         else if(strcmp(response_descriptor,"Misfit")==0){
    54                 Misfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
     53        else if(strcmp(response_descriptor,"SurfaceAbsVelMisfit")==0){
     54                SurfaceAbsVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
    5555        }
    56         else if(strlen(response_descriptor)>=8){
    57                 if(strncmp(response_descriptor,"MassFlux",8)==0) MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units);
     56        else if(strcmp(response_descriptor,"SurfaceRelVelMisfit")==0){
     57                SurfaceRelVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
     58        }
     59        else if(strcmp(response_descriptor,"SurfaceLogVelMisfit")==0){
     60                SurfaceLogVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
     61        }
     62        else if(strcmp(response_descriptor,"SurfaceLogVxVyMisfit")==0){
     63                SurfaceLogVxVyMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
     64        }
     65        else if(strcmp(response_descriptor,"SurfaceAverageVelMisfit")==0){
     66                SurfaceAverageVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
     67        }
     68        else if(strncmp(response_descriptor,"MassFlux",8)==0){
     69                MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units);
    5870        }
    5971        else{
    60                 ISSMERROR("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!");
     72                ISSMERROR(" response descriptor \"%s\" not supported yet!",response_descriptor);
    6173        }
    6274}

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

    r5250 r5281  
    99#include "../../Container/Container.h"
    1010
    11 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units);
     11void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* response_descriptor,bool process_units);
    1212
    1313#endif  /* _RESPONSESXX_H */

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

    r5250 r5281  
    5555#include "./MinVyx/MinVyx.h"
    5656#include "./MinVzx/MinVzx.h"
    57 #include "./Misfitx/Misfitx.h"
     57#include "./SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h"
     58#include "./SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h"
     59#include "./SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h"
     60#include "./SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h"
     61#include "./SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h"
    5862#include "./ModelProcessorx/ModelProcessorx.h"
    5963#include "./NodeConnectivityx/NodeConnectivityx.h"

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

    r4974 r5281  
    4242                virtual void   GradjDrag(Vec gradient)=0;
    4343                virtual void   GradjB(Vec gradient)=0;
    44                 virtual double Misfit(bool process_units)=0;
    45                 virtual double CostFunction(bool process_units)=0;
     44                virtual double SurfaceAbsVelMisfit(bool process_units)=0;
     45                virtual double SurfaceRelVelMisfit(bool process_units)=0;
     46                virtual double SurfaceLogVelMisfit(bool process_units)=0;
     47                virtual double SurfaceLogVxVyMisfit(bool process_units)=0;
     48                virtual double SurfaceAverageVelMisfit(bool process_units)=0;
     49                virtual double RegularizeInversion(void)=0;
    4650                virtual double SurfaceArea(void)=0;
    4751                virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0;

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

    r5263 r5281  
    666666}
    667667/*}}}*/
    668 /*FUNCTION Penta::CostFunction {{{1*/
    669 double Penta::CostFunction(bool process_units){
     668/*FUNCTION Penta::RegularizeInversion {{{1*/
     669double Penta::RegularizeInversion(void){
    670670
    671671        double J;
    … …  
    696696
    697697                /*This element should be collapsed into a tria element at its base. Create this tria element,
    698                  * and compute CostFunction*/
     698                 * and compute RegularizeInversion*/
    699699
    700700                /*Depth Average B*/
    … …  
    702702
    703703                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    704                 J=tria->CostFunction(process_units);
     704                J=tria->RegularizeInversion();
    705705                delete tria->matice; delete tria;
    706706
    … …  
    720720
    721721                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    722                 J=tria->CostFunction(process_units);
     722                J=tria->RegularizeInversion();
    723723                delete tria->matice; delete tria;
    724724
    … …  
    15971597}
    15981598/*}}}*/
    1599 /*FUNCTION Penta::Misfit {{{1*/
    1600 double Penta::Misfit(bool process_units){
     1599/*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/
     1600double Penta::SurfaceAbsVelMisfit(bool process_units){
    16011601
    16021602        double J;
    … …  
    16271627
    16281628                /*This element should be collapsed into a tria element at its base. Create this tria element,
    1629                  * and compute Misfit*/
     1629                 * and compute SurfaceAbsVelMisfit*/
    16301630                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    1631                 J=tria->Misfit(process_units);
     1631                J=tria->SurfaceAbsVelMisfit(process_units);
    16321632                delete tria->matice; delete tria;
    16331633                return J;
    … …  
    16361636
    16371637                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    1638                 J=tria->Misfit(process_units);
     1638                J=tria->SurfaceAbsVelMisfit(process_units);
     1639                delete tria->matice; delete tria;
     1640                return J;
     1641        }
     1642}
     1643/*}}}*/
     1644/*FUNCTION Penta::SurfaceRelVelMisfit {{{1*/
     1645double Penta::SurfaceRelVelMisfit(bool process_units){
     1646
     1647        double J;
     1648        Tria* tria=NULL;
     1649
     1650        /*inputs: */
     1651        bool onwater;
     1652        bool onsurface;
     1653        bool onbed;
     1654        int  approximation;
     1655
     1656        /*retrieve inputs :*/
     1657        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1658        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1659        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     1660        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1661
     1662        /*If on water, return 0: */
     1663        if(onwater)return 0;
     1664
     1665        /*Bail out if this element if:
     1666         * -> Non MacAyeal and not on the surface
     1667         * -> MacAyeal (2d model) and not on bed) */
     1668        if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){
     1669                return 0;
     1670        }
     1671        else if (MacAyealApproximationEnum){
     1672
     1673                /*This element should be collapsed into a tria element at its base. Create this tria element,
     1674                 * and compute SurfaceRelVelMisfit*/
     1675                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     1676                J=tria->SurfaceRelVelMisfit(process_units);
     1677                delete tria->matice; delete tria;
     1678                return J;
     1679        }
     1680        else{
     1681
     1682                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     1683                J=tria->SurfaceRelVelMisfit(process_units);
     1684                delete tria->matice; delete tria;
     1685                return J;
     1686        }
     1687}
     1688/*}}}*/
     1689/*FUNCTION Penta::SurfaceLogVelMisfit {{{1*/
     1690double Penta::SurfaceLogVelMisfit(bool process_units){
     1691
     1692        double J;
     1693        Tria* tria=NULL;
     1694
     1695        /*inputs: */
     1696        bool onwater;
     1697        bool onsurface;
     1698        bool onbed;
     1699        int  approximation;
     1700
     1701        /*retrieve inputs :*/
     1702        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1703        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1704        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     1705        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1706
     1707        /*If on water, return 0: */
     1708        if(onwater)return 0;
     1709
     1710        /*Bail out if this element if:
     1711         * -> Non MacAyeal and not on the surface
     1712         * -> MacAyeal (2d model) and not on bed) */
     1713        if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){
     1714                return 0;
     1715        }
     1716        else if (MacAyealApproximationEnum){
     1717
     1718                /*This element should be collapsed into a tria element at its base. Create this tria element,
     1719                 * and compute SurfaceLogVelMisfit*/
     1720                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     1721                J=tria->SurfaceLogVelMisfit(process_units);
     1722                delete tria->matice; delete tria;
     1723                return J;
     1724        }
     1725        else{
     1726
     1727                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     1728                J=tria->SurfaceLogVelMisfit(process_units);
     1729                delete tria->matice; delete tria;
     1730                return J;
     1731        }
     1732}
     1733/*}}}*/
     1734/*FUNCTION Penta::SurfaceLogVxVyMisfit {{{1*/
     1735double Penta::SurfaceLogVxVyMisfit(bool process_units){
     1736
     1737        double J;
     1738        Tria* tria=NULL;
     1739
     1740        /*inputs: */
     1741        bool onwater;
     1742        bool onsurface;
     1743        bool onbed;
     1744        int  approximation;
     1745
     1746        /*retrieve inputs :*/
     1747        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1748        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1749        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     1750        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1751
     1752        /*If on water, return 0: */
     1753        if(onwater)return 0;
     1754
     1755        /*Bail out if this element if:
     1756         * -> Non MacAyeal and not on the surface
     1757         * -> MacAyeal (2d model) and not on bed) */
     1758        if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){
     1759                return 0;
     1760        }
     1761        else if (MacAyealApproximationEnum){
     1762
     1763                /*This element should be collapsed into a tria element at its base. Create this tria element,
     1764                 * and compute SurfaceLogVxVyMisfit*/
     1765                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     1766                J=tria->SurfaceLogVxVyMisfit(process_units);
     1767                delete tria->matice; delete tria;
     1768                return J;
     1769        }
     1770        else{
     1771
     1772                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     1773                J=tria->SurfaceLogVxVyMisfit(process_units);
     1774                delete tria->matice; delete tria;
     1775                return J;
     1776        }
     1777}
     1778/*}}}*/
     1779/*FUNCTION Penta::SurfaceAverageVelMisfit {{{1*/
     1780double Penta::SurfaceAverageVelMisfit(bool process_units){
     1781
     1782        double J;
     1783        Tria* tria=NULL;
     1784
     1785        /*inputs: */
     1786        bool onwater;
     1787        bool onsurface;
     1788        bool onbed;
     1789        int  approximation;
     1790
     1791        /*retrieve inputs :*/
     1792        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1793        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1794        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     1795        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1796
     1797        /*If on water, return 0: */
     1798        if(onwater)return 0;
     1799
     1800        /*Bail out if this element if:
     1801         * -> Non MacAyeal and not on the surface
     1802         * -> MacAyeal (2d model) and not on bed) */
     1803        if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){
     1804                return 0;
     1805        }
     1806        else if (MacAyealApproximationEnum){
     1807
     1808                /*This element should be collapsed into a tria element at its base. Create this tria element,
     1809                 * and compute SurfaceAverageVelMisfit*/
     1810                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     1811                J=tria->SurfaceAverageVelMisfit(process_units);
     1812                delete tria->matice; delete tria;
     1813                return J;
     1814        }
     1815        else{
     1816
     1817                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     1818                J=tria->SurfaceAverageVelMisfit(process_units);
    16391819                delete tria->matice; delete tria;
    16401820                return J;
    … …  
    59006080                                name==TemperatureEnum ||
    59016081                                name==ControlParameterEnum ||
    5902                                 name==FitEnum ||
     6082                                name==CmResponseEnum ||
    59036083                                name==DragCoefficientEnum ||
    59046084                                name==GradientEnum ||

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

    r5258 r5281  
    7171                void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    7272                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    73                 double CostFunction(bool process_units);
     73                double RegularizeInversion(void);
    7474                void   CreateKMatrix(Mat Kgg);
    7575                void   CreatePVector(Vec pg);
    … …  
    104104                void   MinVy(double* pminvy, bool process_units);
    105105                void   MinVz(double* pminvz, bool process_units);
    106                 double Misfit(bool process_units);
     106                double SurfaceAbsVelMisfit(bool process_units);
     107                double SurfaceRelVelMisfit(bool process_units);
     108                double SurfaceLogVelMisfit(bool process_units);
     109                double SurfaceLogVxVyMisfit(bool process_units);
     110                double SurfaceAverageVelMisfit(bool process_units);
    107111                void   PatchFill(int* pcount, Patch* patch);
    108112                void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);

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

    r5263 r5281  
    548548}
    549549/*}}}*/
    550 /*FUNCTION Tria::CostFunction {{{1*/
    551 double Tria::CostFunction(bool process_units){
     550/*FUNCTION Tria::RegularizeInversion {{{1*/
     551double Tria::RegularizeInversion(){
    552552
    553553        int i;
    554554
    555555        /* output: */
    556         double Jelem;
     556        double Jelem=0;
    557557
    558558        /* node data: */
    … …  
    602602        /* Get node coordinates and dof list: */
    603603        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    604 
    605         /*First, get Misfit*/
    606         Jelem=Misfit(process_units);
    607604
    608605        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    18601857}
    18611858/*}}}*/
    1862 /*FUNCTION Tria::Misfit {{{1*/
    1863 double Tria::Misfit(bool process_units){
     1859/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/
     1860double Tria::SurfaceAbsVelMisfit(bool process_units){
    18641861
    18651862        int i;
    … …  
    19001897        double Jdet;
    19011898
    1902         /*relative and logarithmic control method :*/
    1903         double scalex=1;
    1904         double scaley=1;
    1905         double S=0;
    1906         int    fit=-1;
    1907 
    19081899        /*inputs: */
    19091900        bool onwater;
    … …  
    19211912
    19221913        /* Recover input data: */
    1923         inputs->GetParameterValue(&fit,FitEnum);
    1924         if(fit==3){
    1925                 inputs->GetParameterValue(&S,SurfaceAreaEnum);
    1926         }
    19271914        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
    19281915        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
    … …  
    19351922        this->parameters->FindParam(&epsvel,EpsVelEnum);
    19361923       
    1937         /* Compute Misfit at the 3 nodes
     1924        /* Compute SurfaceAbsVelMisfit at the 3 nodes
    19381925         * Here we integrate linearized functions:
    19391926         *               
    … …  
    19441931         *       the vertex i
    19451932         */
    1946         if(fit==0){
    1947                 /*We are using an absolute misfit:
    1948                  *
    1949                  *      1  [           2              2 ]
    1950                  * J = --- | (u - u   )  +  (v - v   )  |
    1951                  *      2  [       obs            obs   ]
    1952                  *
    1953                  */
    1954                 for (i=0;i<numgrids;i++){
    1955                         misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
    1956                 }
    1957                 /*Process units: */
    1958                 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters);
    1959 
    1960         }
    1961         else if(fit==1){
    1962                 /*We are using a relative misfit:
    1963                  *                       
    1964                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    1965                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    1966                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    1967                  *              obs                        obs                     
    1968                  */
    1969                 for (i=0;i<numgrids;i++){
    1970                         scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
    1971                         scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
    1972                         if(obs_vx_list[i]==0)scalex=0;
    1973                         if(obs_vy_list[i]==0)scaley=0;
    1974                         misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2));
    1975                 }
    1976 
    1977                 /*Process units: */
    1978                 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters);
    1979 
    1980         }
    1981         else if(fit==2){
    1982                 /*We are using a logarithmic misfit:
    1983                 *                       
    1984                 *                 [        vel + eps     ] 2
    1985                 * J = 4 \bar{v}^2 | log ( -----------  ) | 
    1986                 *                 [       vel   + eps    ]
    1987                 *                            obs
    1988                 */
    1989                 for (i=0;i<numgrids;i++){
    1990                         velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
    1991                         obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
    1992                         misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
    1993                 }
    1994                
    1995                 /*Process units: */
    1996                 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters);
    1997 
    1998         }
    1999         else if(fit==3){
    2000                 /*We are using a spacially average absolute misfit:
    2001                  *
    2002                  *      1                    2              2
    2003                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    2004                  *      S                obs            obs
    2005                  */
    2006                 for (i=0;i<numgrids;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
    2007                
    2008                 /*Process units: */
    2009                 if(process_units)NodalValuesUnitConversion(&misfit_square_list[0],numgrids,MisfitEnum,this->parameters);
    2010        
    2011                 /*Take the square root, and scale by surface: */
    2012                 for (i=0;i<numgrids;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
    2013         }
    2014         else if(fit==4){
    2015                 /*We are using an logarithmic 2 misfit:
    2016                  *
    2017                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    2018                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    2019                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    2020                  *                              obs                       obs
    2021                  */
    2022                 for (i=0;i<numgrids;i++){
    2023                         misfit_list[i]=0.5*pow(meanvel,(double)2)*(
    2024                           pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
    2025                           pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) );
    2026                 }
    2027                
    2028                 /*Process units: */
    2029                 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters);
    2030 
    2031         }
    2032         else{
    2033                 /*Not supported yet! : */
    2034                 ISSMERROR("%s%g","unsupported type of fit: ",fit);
    2035         }
     1933
     1934        /*We are using an absolute misfit:
     1935         *
     1936         *      1  [           2              2 ]
     1937         * J = --- | (u - u   )  +  (v - v   )  |
     1938         *      2  [       obs            obs   ]
     1939         *
     1940         */
     1941        for (i=0;i<numgrids;i++){
     1942                misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
     1943        }
     1944        /*Process units: */
     1945        if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceAbsVelMisfitEnum,this->parameters);
    20361946
    20371947        /*Apply weights to misfits*/
    … …  
    20571967                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
    20581968
    2059                 /*compute Misfit*/
     1969                /*compute SurfaceAbsVelMisfit*/
    20601970                Jelem+=misfit*Jdet*gauss_weight;
    20611971        }
    … …  
    20661976        xfree((void**)&gauss_weights);
    20671977               
     1978        /*Return: */
     1979        return Jelem;
     1980}
     1981/*}}}*/
     1982/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
     1983double Tria::SurfaceRelVelMisfit(bool process_units){
     1984
     1985        int i;
     1986
     1987        /* output: */
     1988        double Jelem=0;
     1989
     1990        /* node data: */
     1991        const int    numgrids=3;
     1992        const int    numdof=2*numgrids;
     1993        const int    NDOF2=2;
     1994        double       xyz_list[numgrids][3];
     1995
     1996        /* grid data: */
     1997        double vx_list[numgrids];
     1998        double vy_list[numgrids];
     1999        double obs_vx_list[numgrids];
     2000        double obs_vy_list[numgrids];
     2001        double misfit_square_list[numgrids];
     2002        double misfit_list[numgrids];
     2003        double weights_list[numgrids];
     2004
     2005        /* gaussian points: */
     2006        int     num_gauss,ig;
     2007        double* first_gauss_area_coord  =  NULL;
     2008        double* second_gauss_area_coord =  NULL;
     2009        double* third_gauss_area_coord  =  NULL;
     2010        double* gauss_weights           =  NULL;
     2011        double  gauss_weight;
     2012        double  gauss_l1l2l3[3];
     2013        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2014
     2015        /* parameters: */
     2016        double  velocity_mag,obs_velocity_mag;
     2017        double  misfit;
     2018
     2019        /* Jacobian: */
     2020        double Jdet;
     2021
     2022        /*relative and logarithmic control method :*/
     2023        double scalex=1;
     2024        double scaley=1;
     2025
     2026        /*inputs: */
     2027        bool onwater;
     2028
     2029        double  meanvel, epsvel;
     2030
     2031        /*retrieve inputs :*/
     2032        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2033
     2034        /*If on water, return 0: */
     2035        if(onwater)return 0;
     2036
     2037        /* Get node coordinates and dof list: */
     2038        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2039
     2040        /* Recover input data: */
     2041        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     2042        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     2043        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     2044        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     2045        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     2046
     2047        /*retrieve some parameters: */
     2048        this->parameters->FindParam(&meanvel,MeanVelEnum);
     2049        this->parameters->FindParam(&epsvel,EpsVelEnum);
     2050
     2051        /* Compute SurfaceRelVelMisfit at the 3 nodes
     2052         * Here we integrate linearized functions:
     2053         *               
     2054         * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
     2055         *
     2056         * where J_i are the misfits at the 3 nodes of the triangle
     2057         *       Phi_i is the nodal function (P1) with respect to
     2058         *       the vertex i
     2059         */
     2060
     2061        /*We are using a relative misfit:
     2062         *                       
     2063         *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     2064         * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     2065         *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     2066         *              obs                        obs                     
     2067         */
     2068        for (i=0;i<numgrids;i++){
     2069                scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
     2070                scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
     2071                if(obs_vx_list[i]==0)scalex=0;
     2072                if(obs_vy_list[i]==0)scaley=0;
     2073                misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2));
     2074        }
     2075
     2076        /*Process units: */
     2077        if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceRelVelMisfitEnum,this->parameters);
     2078
     2079        /*Apply weights to misfits*/
     2080        for (i=0;i<numgrids;i++){
     2081                misfit_list[i]=weights_list[i]*misfit_list[i];
     2082        }
     2083
     2084        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2085        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2086
     2087        /* Start  looping on the number of gaussian points: */
     2088        for (ig=0; ig<num_gauss; ig++){
     2089                /*Pick up the gaussian point: */
     2090                gauss_weight=*(gauss_weights+ig);
     2091                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2092                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2093                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2094
     2095                /* Get Jacobian determinant: */
     2096                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2097
     2098                /*Compute misfit at gaussian point: */
     2099                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
     2100
     2101                /*compute SurfaceRelVelMisfit*/
     2102                Jelem+=misfit*Jdet*gauss_weight;
     2103        }
     2104
     2105        xfree((void**)&first_gauss_area_coord);
     2106        xfree((void**)&second_gauss_area_coord);
     2107        xfree((void**)&third_gauss_area_coord);
     2108        xfree((void**)&gauss_weights);
     2109
     2110        /*Return: */
     2111        return Jelem;
     2112}
     2113/*}}}*/
     2114/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/
     2115double Tria::SurfaceLogVelMisfit(bool process_units){
     2116
     2117        int i;
     2118
     2119        /* output: */
     2120        double Jelem=0;
     2121
     2122        /* node data: */
     2123        const int    numgrids=3;
     2124        const int    numdof=2*numgrids;
     2125        const int    NDOF2=2;
     2126        double       xyz_list[numgrids][3];
     2127
     2128        /* grid data: */
     2129        double vx_list[numgrids];
     2130        double vy_list[numgrids];
     2131        double obs_vx_list[numgrids];
     2132        double obs_vy_list[numgrids];
     2133        double misfit_square_list[numgrids];
     2134        double misfit_list[numgrids];
     2135        double weights_list[numgrids];
     2136
     2137        /* gaussian points: */
     2138        int     num_gauss,ig;
     2139        double* first_gauss_area_coord  =  NULL;
     2140        double* second_gauss_area_coord =  NULL;
     2141        double* third_gauss_area_coord  =  NULL;
     2142        double* gauss_weights           =  NULL;
     2143        double  gauss_weight;
     2144        double  gauss_l1l2l3[3];
     2145        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2146
     2147        /* parameters: */
     2148        double  velocity_mag,obs_velocity_mag;
     2149        double  misfit;
     2150
     2151        /* Jacobian: */
     2152        double Jdet;
     2153
     2154        /*relative and logarithmic control method :*/
     2155        double scalex=1;
     2156        double scaley=1;
     2157
     2158        /*inputs: */
     2159        bool onwater;
     2160       
     2161        double  meanvel, epsvel;
     2162
     2163        /*retrieve inputs :*/
     2164        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2165
     2166        /*If on water, return 0: */
     2167        if(onwater)return 0;
     2168
     2169        /* Get node coordinates and dof list: */
     2170        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2171
     2172        /* Recover input data: */
     2173        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     2174        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     2175        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     2176        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     2177        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     2178
     2179        /*retrieve some parameters: */
     2180        this->parameters->FindParam(&meanvel,MeanVelEnum);
     2181        this->parameters->FindParam(&epsvel,EpsVelEnum);
     2182       
     2183        /* Compute SurfaceLogVelMisfit at the 3 nodes
     2184         * Here we integrate linearized functions:
     2185         *               
     2186         * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
     2187         *
     2188         * where J_i are the misfits at the 3 nodes of the triangle
     2189         *       Phi_i is the nodal function (P1) with respect to
     2190         *       the vertex i
     2191         */
     2192
     2193        /*We are using a logarithmic misfit:
     2194         *                       
     2195         *                 [        vel + eps     ] 2
     2196         * J = 4 \bar{v}^2 | log ( -----------  ) | 
     2197         *                 [       vel   + eps    ]
     2198         *                            obs
     2199         */
     2200        for (i=0;i<numgrids;i++){
     2201                velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
     2202                obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
     2203                misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
     2204        }
     2205
     2206        /*Process units: */
     2207        if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceLogVelMisfitEnum,this->parameters);
     2208
     2209        /*Apply weights to misfits*/
     2210        for (i=0;i<numgrids;i++){
     2211                misfit_list[i]=weights_list[i]*misfit_list[i];
     2212        }
     2213
     2214        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2215        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2216
     2217        /* Start  looping on the number of gaussian points: */
     2218        for (ig=0; ig<num_gauss; ig++){
     2219                /*Pick up the gaussian point: */
     2220                gauss_weight=*(gauss_weights+ig);
     2221                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2222                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2223                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2224
     2225                /* Get Jacobian determinant: */
     2226                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2227
     2228                /*Compute misfit at gaussian point: */
     2229                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
     2230
     2231                /*compute SurfaceLogVelMisfit*/
     2232                Jelem+=misfit*Jdet*gauss_weight;
     2233        }
     2234
     2235        xfree((void**)&first_gauss_area_coord);
     2236        xfree((void**)&second_gauss_area_coord);
     2237        xfree((void**)&third_gauss_area_coord);
     2238        xfree((void**)&gauss_weights);
     2239               
     2240        /*Return: */
     2241        return Jelem;
     2242}
     2243/*}}}*/
     2244/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
     2245double Tria::SurfaceLogVxVyMisfit(bool process_units){
     2246
     2247        int i;
     2248
     2249        /* output: */
     2250        double Jelem=0;
     2251
     2252        /* node data: */
     2253        const int    numgrids=3;
     2254        const int    numdof=2*numgrids;
     2255        const int    NDOF2=2;
     2256        double       xyz_list[numgrids][3];
     2257
     2258        /* grid data: */
     2259        double vx_list[numgrids];
     2260        double vy_list[numgrids];
     2261        double obs_vx_list[numgrids];
     2262        double obs_vy_list[numgrids];
     2263        double misfit_square_list[numgrids];
     2264        double misfit_list[numgrids];
     2265        double weights_list[numgrids];
     2266
     2267        /* gaussian points: */
     2268        int     num_gauss,ig;
     2269        double* first_gauss_area_coord  =  NULL;
     2270        double* second_gauss_area_coord =  NULL;
     2271        double* third_gauss_area_coord  =  NULL;
     2272        double* gauss_weights           =  NULL;
     2273        double  gauss_weight;
     2274        double  gauss_l1l2l3[3];
     2275        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2276
     2277        /* parameters: */
     2278        double  velocity_mag,obs_velocity_mag;
     2279        double  misfit;
     2280
     2281        /* Jacobian: */
     2282        double Jdet;
     2283
     2284        /*relative and logarithmic control method :*/
     2285        double scalex=1;
     2286        double scaley=1;
     2287        double S=0;
     2288        int    fit=-1;
     2289
     2290        /*inputs: */
     2291        bool onwater;
     2292       
     2293        double  meanvel, epsvel;
     2294
     2295        /*retrieve inputs :*/
     2296        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2297
     2298        /*If on water, return 0: */
     2299        if(onwater)return 0;
     2300
     2301        /* Get node coordinates and dof list: */
     2302        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2303
     2304        /* Recover input data: */
     2305        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     2306        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     2307        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     2308        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     2309        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     2310
     2311        /*retrieve some parameters: */
     2312        this->parameters->FindParam(&meanvel,MeanVelEnum);
     2313        this->parameters->FindParam(&epsvel,EpsVelEnum);
     2314       
     2315        /* Compute SurfaceLogVxVyMisfit at the 3 nodes
     2316         * Here we integrate linearized functions:
     2317         *               
     2318         * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
     2319         *
     2320         * where J_i are the misfits at the 3 nodes of the triangle
     2321         *       Phi_i is the nodal function (P1) with respect to
     2322         *       the vertex i
     2323         */
     2324
     2325        /*We are using an logarithmic 2 misfit:
     2326         *
     2327         *      1            [        |u| + eps     2          |v| + eps     2  ]
     2328         * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     2329         *      2            [       |u    |+ eps              |v    |+ eps     ]
     2330         *                              obs                       obs
     2331         */
     2332        for (i=0;i<numgrids;i++){
     2333                misfit_list[i]=0.5*pow(meanvel,(double)2)*(
     2334                                        pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
     2335                                        pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) );
     2336        }
     2337
     2338        /*Process units: */
     2339        if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceLogVxVyMisfitEnum,this->parameters);
     2340
     2341        /*Apply weights to misfits*/
     2342        for (i=0;i<numgrids;i++){
     2343                misfit_list[i]=weights_list[i]*misfit_list[i];
     2344        }
     2345
     2346        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2347        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2348
     2349        /* Start  looping on the number of gaussian points: */
     2350        for (ig=0; ig<num_gauss; ig++){
     2351                /*Pick up the gaussian point: */
     2352                gauss_weight=*(gauss_weights+ig);
     2353                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2354                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2355                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2356
     2357                /* Get Jacobian determinant: */
     2358                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2359
     2360                /*Compute misfit at gaussian point: */
     2361                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
     2362
     2363                /*compute SurfaceLogVxVyMisfit*/
     2364                Jelem+=misfit*Jdet*gauss_weight;
     2365        }
     2366
     2367        xfree((void**)&first_gauss_area_coord);
     2368        xfree((void**)&second_gauss_area_coord);
     2369        xfree((void**)&third_gauss_area_coord);
     2370        xfree((void**)&gauss_weights);
     2371               
     2372        /*Return: */
     2373        return Jelem;
     2374}
     2375/*}}}*/
     2376/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/
     2377double Tria::SurfaceAverageVelMisfit(bool process_units){
     2378
     2379        int i;
     2380
     2381        /* output: */
     2382        double Jelem=0;
     2383
     2384        /* node data: */
     2385        const int    numgrids=3;
     2386        const int    numdof=2*numgrids;
     2387        const int    NDOF2=2;
     2388        double       xyz_list[numgrids][3];
     2389
     2390        /* grid data: */
     2391        double vx_list[numgrids];
     2392        double vy_list[numgrids];
     2393        double obs_vx_list[numgrids];
     2394        double obs_vy_list[numgrids];
     2395        double misfit_square_list[numgrids];
     2396        double misfit_list[numgrids];
     2397        double weights_list[numgrids];
     2398
     2399        /* gaussian points: */
     2400        int     num_gauss,ig;
     2401        double* first_gauss_area_coord  =  NULL;
     2402        double* second_gauss_area_coord =  NULL;
     2403        double* third_gauss_area_coord  =  NULL;
     2404        double* gauss_weights           =  NULL;
     2405        double  gauss_weight;
     2406        double  gauss_l1l2l3[3];
     2407        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2408
     2409        /* parameters: */
     2410        double  velocity_mag,obs_velocity_mag;
     2411        double  misfit;
     2412
     2413        /* Jacobian: */
     2414        double Jdet;
     2415
     2416        /*relative and logarithmic control method :*/
     2417        double scalex=1;
     2418        double scaley=1;
     2419        double S=0;
     2420        int    fit=-1;
     2421
     2422        /*inputs: */
     2423        bool onwater;
     2424
     2425        double  meanvel, epsvel;
     2426
     2427        /*retrieve inputs :*/
     2428        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2429
     2430        /*If on water, return 0: */
     2431        if(onwater)return 0;
     2432
     2433        /* Get node coordinates and dof list: */
     2434        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2435
     2436        /* Recover input data: */
     2437        inputs->GetParameterValue(&S,SurfaceAreaEnum);
     2438        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     2439        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     2440        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     2441        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     2442        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     2443
     2444        /*retrieve some parameters: */
     2445        this->parameters->FindParam(&meanvel,MeanVelEnum);
     2446        this->parameters->FindParam(&epsvel,EpsVelEnum);
     2447
     2448        /* Compute SurfaceAverageVelMisfit at the 3 nodes
     2449         * Here we integrate linearized functions:
     2450         *               
     2451         * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
     2452         *
     2453         * where J_i are the misfits at the 3 nodes of the triangle
     2454         *       Phi_i is the nodal function (P1) with respect to
     2455         *       the vertex i
     2456         */
     2457
     2458        /*We are using a spacially average absolute misfit:
     2459         *
     2460         *      1                    2              2
     2461         * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     2462         *      S                obs            obs
     2463         */
     2464        for (i=0;i<numgrids;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
     2465
     2466        /*Process units: */
     2467        if(process_units)NodalValuesUnitConversion(&misfit_square_list[0],numgrids,SurfaceAverageVelMisfitEnum,this->parameters);
     2468
     2469        /*Take the square root, and scale by surface: */
     2470        for (i=0;i<numgrids;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
     2471
     2472        /*Apply weights to misfits*/
     2473        for (i=0;i<numgrids;i++){
     2474                misfit_list[i]=weights_list[i]*misfit_list[i];
     2475        }
     2476
     2477        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2478        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2479
     2480        /* Start  looping on the number of gaussian points: */
     2481        for (ig=0; ig<num_gauss; ig++){
     2482                /*Pick up the gaussian point: */
     2483                gauss_weight=*(gauss_weights+ig);
     2484                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2485                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2486                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2487
     2488                /* Get Jacobian determinant: */
     2489                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2490
     2491                /*Compute misfit at gaussian point: */
     2492                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
     2493
     2494                /*compute SurfaceAverageVelMisfit*/
     2495                Jelem+=misfit*Jdet*gauss_weight;
     2496        }
     2497
     2498        xfree((void**)&first_gauss_area_coord);
     2499        xfree((void**)&second_gauss_area_coord);
     2500        xfree((void**)&third_gauss_area_coord);
     2501        xfree((void**)&gauss_weights);
     2502
    20682503        /*Return: */
    20692504        return Jelem;
    … …  
    21562591        /*inputs: */
    21572592        bool onwater;
    2158         int fit;
    21592593
    21602594        /*retrieve inputs :*/
    2161         inputs->GetParameterValue(&fit,FitEnum);
    21622595        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    2163 
    2164         /*If fit!=3, do not compute surface: */
    2165         if(fit!=3)return 0;
    21662596
    21672597        /*If on water, return 0: */
    … …  
    44714901        double scale=0;
    44724902        double S=0;
    4473         int    fit=-1;
     4903        int    response;
    44744904
    44754905        /*retrieve some parameters: */
    … …  
    44904920        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
    44914921
    4492         inputs->GetParameterValue(&fit,FitEnum);
    4493         if(fit==3){
     4922        inputs->GetParameterValue(&response,CmResponseEnum);
     4923        if(response==SurfaceAverageVelMisfitEnum){
    44944924                inputs->GetParameterValue(&S,SurfaceAreaEnum);
    44954925        }
    … …  
    45084938         *       the vertex i
    45094939         */
    4510         if(fit==0){
     4940        if(response==SurfaceAbsVelMisfitEnum){
    45114941                /*We are using an absolute misfit:
    45124942                 *
    … …  
    45244954                }
    45254955        }
    4526         else if(fit==1){
     4956        else if(response==SurfaceRelVelMisfitEnum){
    45274957                /*We are using a relative misfit:
    45284958                 *                       
    … …  
    45464976                }
    45474977        }
    4548         else if(fit==2){
     4978        else if(response==SurfaceLogVelMisfitEnum){
    45494979                /*We are using a logarithmic misfit:
    45504980                 *                       
    … …  
    45674997                }
    45684998        }
    4569         else if(fit==3){
     4999        else if(response==SurfaceAverageVelMisfitEnum){
    45705000                /*We are using an spacially average absolute misfit:
    45715001                 *
    … …  
    45845014                }
    45855015        }
    4586         else if(fit==4){
     5016        else if(response==SurfaceLogVxVyMisfitEnum){
    45875017                /*We are using an logarithmic 2 misfit:
    45885018                 *
    … …  
    46045034        else{
    46055035                /*Not supported yet! : */
    4606                 ISSMERROR("%s%g","unsupported type of fit: ",fit);
     5036                ISSMERROR("response %s not supported yet",EnumToString(response));
    46075037        }
    46085038
    … …  
    47105140        double S=0;
    47115141        int    fit=-1;
     5142        int    response;
    47125143
    47135144        /*retrieve some parameters: */
    … …  
    47285159        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
    47295160
    4730         inputs->GetParameterValue(&fit,FitEnum);
    4731         if(fit==3){
     5161        inputs->GetParameterValue(&response,CmResponseEnum);
     5162        if(response==SurfaceAverageVelMisfitEnum){
    47325163                inputs->GetParameterValue(&S,SurfaceAreaEnum);
    47335164        }
    … …  
    47465177         *       the vertex i
    47475178         */
    4748         if(fit==0){
     5179        if(response==SurfaceAbsVelMisfitEnum){
    47495180                /*We are using an absolute misfit:
    47505181                 *
    … …  
    47625193                }
    47635194        }
    4764         else if(fit==1){
     5195        else if(response==SurfaceRelVelMisfitEnum){
    47655196                /*We are using a relative misfit:
    47665197                 *                       
    … …  
    47845215                }
    47855216        }
    4786         else if(fit==2){
     5217        else if(response==SurfaceLogVelMisfitEnum){
    47875218                /*We are using a logarithmic misfit:
    47885219                 *                       
    … …  
    48055236                }
    48065237        }
    4807         else if(fit==3){
     5238        else if(response==SurfaceAverageVelMisfitEnum){
    48085239                /*We are using an spacially average absolute misfit:
    48095240                 *
    … …  
    48225253                }
    48235254        }
    4824         else if(fit==4){
     5255        else if(response==SurfaceLogVxVyMisfitEnum){
    48255256                /*We are using an logarithmic 2 misfit:
    48265257                 *
    … …  
    48425273        else{
    48435274                /*Not supported yet! : */
    4844                 ISSMERROR("%s%g","unsupported type of fit: ",fit);
     5275                ISSMERROR("response %s not supported yet",EnumToString(response));
    48455276        }
    48465277
    … …  
    62986729                                name==VxObsEnum ||
    62996730                                name==VyObsEnum ||
    6300                                 name==FitEnum ||
     6731                                name==CmResponseEnum ||
    63016732                                name==DragCoefficientEnum ||
    63026733                                name==GradientEnum ||

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

    r5258 r5281  
    6868                void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    6969                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    70                 double CostFunction(bool process_units);
     70                double RegularizeInversion(void);
    7171                void   CreateKMatrix(Mat Kgg);
    7272                void   CreatePVector(Vec pg);
    … …  
    103103                void   MinVy(double* pminvy, bool process_units);
    104104                void   MinVz(double* pminvz, bool process_units);
    105                 double Misfit(bool process_units);
     105                double SurfaceAbsVelMisfit(bool process_units);
     106                double SurfaceRelVelMisfit(bool process_units);
     107                double SurfaceLogVelMisfit(bool process_units);
     108                double SurfaceLogVxVyMisfit(bool process_units);
     109                double SurfaceAverageVelMisfit(bool process_units);
    106110                void   PatchFill(int* pcount, Patch* patch);
    107111                void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);

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

    r5251 r5281  
    8686        xfree((void**)&this->rheology_B);
    8787        xfree((void**)&this->rheology_n);
    88         xfree((void**)&this->fit);
     88        xfree((void**)&this->cm_responses);
    8989        xfree((void**)&this->weights);
    9090        xfree((void**)&this->cm_jump);
    … …  
    302302
    303303        /*!solution parameters: */
    304         this->fit=NULL;
     304        this->cm_responses=NULL;
    305305        this->weights=NULL;
    306306        this->cm_jump=NULL;

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

    r5251 r5281  
    129129
    130130                /*solution parameters: */
    131                 double* fit;
     131                double* cm_responses;
    132132                double* weights;
    133133                double* cm_jump;

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

    r5247 r5281  
    2626        int     solution_type;
    2727
    28         double* fit=NULL;
     28        double* responses=NULL;
    2929        double* optscal=NULL;
    3030        double* maxiter=NULL;
    … …  
    4646        femmodel->parameters->FindParam(&nsteps,NStepsEnum);
    4747        femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
    48         femmodel->parameters->FindParam(&fit,NULL,FitEnum);
     48        femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    4949        femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
    5050        femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum);
    … …  
    6565        }
    6666
    67         /*Initialize misfit: */
     67        /*Initialize responses: */
    6868        J=(double*)xmalloc(nsteps*sizeof(double));
    6969               
    … …  
    7575        for(n=0;n<nsteps;n++){
    7676
     77                /*Display info*/
    7778                _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    78                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)fit[n],FitEnum);
     79                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n],CmResponseEnum);
    7980               
    8081                /*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */
    … …  
    9798                InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true);
    9899               
    99                 if(controlconvergence(J,fit,eps_cm,n)) break;
     100                if(controlconvergence(J,responses,eps_cm,n)) break;
    100101
    101102                /*Temporary saving every 5 control steps: */
    … …  
    117118        cleanup_and_return:
    118119        /*Free ressources: */
    119         xfree((void**)&fit);
     120        xfree((void**)&responses);
    120121        xfree((void**)&optscal);
    121122        xfree((void**)&maxiter);

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

    r4778 r5281  
    1818
    1919
    20 bool controlconvergence(double* J, double* fit, double eps_cm, int n){
     20bool controlconvergence(double* J,double* responses, double eps_cm, int n){
    2121
    2222        int i;
    … …  
    2929                //go through the previous misfits(starting from n-2)
    3030                while(i>=0){
    31                         if (fit[i]==fit[n]){
     31                        if (responses[i]==responses[n]){
    3232                                //convergence test only if we have the same misfits
    3333                                if ((J[i]-J[n])/J[n] <= eps_cm){

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

    r5251 r5281  
    2828       
    2929        /*parameters: */
    30         FemModel *femmodel       = NULL;
    31         int       n;
    32         double   *optscal        = NULL;
    33         double   *fit            = NULL;
    34         int       control_type;
    35         int       solution_type;
    36         int       analysis_type;
    37         bool      isstokes       = false;
    38         bool      conserve_loads = true;
    39         bool      process_units  = false;
    40        
     30        FemModel  *femmodel       = NULL;
     31        int        n;
     32        double    *optscal        = NULL;
     33        double    *responses      =NULL;
     34        int        control_type;
     35        int        solution_type;
     36        int        analysis_type;
     37        bool       isstokes       = false;
     38        bool       conserve_loads = true;
     39
    4140        /*Recover finite element model: */
    4241        femmodel=optargs->femmodel;
    … …  
    4645
    4746        femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
    48         femmodel->parameters->FindParam(&fit,NULL,FitEnum);
     47        femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    4948        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    5049        femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
    … …  
    8281
    8382        /*Compute misfit for this velocity field.*/
    84         InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)fit[n],FitEnum);
    85         CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,process_units);
     83        CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n]);
    8684
    8785        /*Free ressources:*/
    88         xfree((void**)&fit);
     86        xfree((void**)&responses);
    8987        xfree((void**)&optscal);
    9088        return J;

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

    r5247 r5281  
    3232//convergence:
    3333void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters);
    34 bool controlconvergence(double* J, double* fit, double eps_cm, int n);
     34bool controlconvergence(double* J, double* response, double eps_cm, int n);
    3535bool steadystateconvergence(FemModel* femmodel);
    3636
Note: See TracChangeset for help on using the changeset viewer.