Changeset 8607


Ignore:
Timestamp:
06/10/11 16:54:02 (14 years ago)
Author:
Mathieu Morlighem
Message:

Now, a combination of responses can be requested, needs to be tested

Location:
issm/trunk/src
Files:
29 edited

Legend:

Unmodified
Added
Removed

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

    r8224 r8607  
    1111#include "../Responsex/Responsex.h"
    1212
    13 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,int response){
     13void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
    1414
    1515        /*Intermediary*/
    1616        int    i;
    17         double fit;
     17        int     num_responses;
    1818        double S;
    1919        Element* element=NULL;
     20        int*     responses=NULL;
    2021
    2122        /*output: */
    22         double J;
     23        double J,Jplus;
     24
     25       
     26        /*Recover parameters*/
     27        parameters->FindParam(&num_responses,NumResponsesEnum);
     28        parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     29
     30        /*Get response*/
     31        J=0;
     32        for(int i=0;i<num_responses;i++){
     33                Responsex(&Jplus,elements,nodes,vertices,loads,materials,parameters,EnumToStringx(responses[i]),false,i); //False means DO NOT process units
     34                J+=Jplus;
     35        }
     36
     37        /*REST TO BE DELETED*/
     38
     39        /*Add Regularization terms: */
    2340        double Jreg=0;
    2441        double Jreg_sum;
    25        
    26         /*Get response*/
    27         Responsex(&J,elements,nodes,vertices,loads,materials,parameters,EnumToStringx(response),false); //False means DO NOT process units
    28 
    29         /*Add Regularization terms: */
    3042        for (i=0;i<elements->Size();i++){
    3143                element=(Element*)elements->GetObjectByOffset(i);
    … …  
    3951
    4052        /*Assign output pointers: */
     53        xfree((void**)&responses);
    4154        *pJ=J;
    4255}

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

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

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

    r8303 r8607  
    7676                       
    7777                        //Responsex(responses_pointer,elements,nodes, vertices,loads,materials, parameters,root,process_units);
    78                         Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units);
     78                        Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units,0);//0 is the index for weights
    7979                       
    8080                        if(my_rank==0){
    … …  
    9595                       
    9696                        /*perfectly normal response function: */
    97                         Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units);
     97                        Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units,0);//0 is the weight index
    9898
    9999                        if(my_rank==0){

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

    r8224 r8607  
    1616#include "../modules.h"
    1717
    18 void Responsex(double* responses,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units){
     18void Responsex(double* responses,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units,int weight_index){
    1919
    2020        switch (StringToEnumx(response_descriptor)){
    2121               
    22                 case MinVelEnum: MinVelx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    23                 case MaxVelEnum: MaxVelx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    24                 case MinVxEnum:  MinVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    25                 case MaxVxEnum:  MaxVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    26                 case MaxAbsVxEnum: MaxAbsVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    27                 case MinVyEnum:  MinVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    28                 case MaxVyEnum:  MaxVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    29                 case MaxAbsVyEnum: MaxAbsVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    30                 case MinVzEnum: MinVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);  break;
    31                 case MaxVzEnum: MaxVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);  break;
    32                 case MaxAbsVzEnum: MaxAbsVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    33                 case MassFluxEnum: MassFluxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    34                 case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    35                 case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    36                 case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    37                 case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    38                 case SurfaceAverageVelMisfitEnum: SurfaceAverageVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
    39                 case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     22                case MinVelEnum:                 MinVelx(                  responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     23                case MaxVelEnum:                 MaxVelx(                  responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     24                case MinVxEnum:                  MinVxx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     25                case MaxVxEnum:                  MaxVxx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     26                case MaxAbsVxEnum:               MaxAbsVxx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     27                case MinVyEnum:                  MinVyx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     28                case MaxVyEnum:                  MaxVyx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     29                case MaxAbsVyEnum:               MaxAbsVyx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     30                case MinVzEnum:                  MinVzx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     31                case MaxVzEnum:                  MaxVzx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     32                case MaxAbsVzEnum:               MaxAbsVzx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     33                case MassFluxEnum:               MassFluxx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
     34                case SurfaceAbsVelMisfitEnum:    SurfaceAbsVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
     35                case SurfaceRelVelMisfitEnum:    SurfaceRelVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
     36                case SurfaceLogVelMisfitEnum:    SurfaceLogVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
     37                case SurfaceLogVxVyMisfitEnum:   SurfaceLogVxVyMisfitx(    responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
     38                case SurfaceAverageVelMisfitEnum:SurfaceAverageVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
     39                case ThicknessAbsMisfitEnum:     ThicknessAbsMisfitx(      responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
    4040                default: _error_(" response descriptor \"%s\" not supported yet!",response_descriptor); break;
    4141        }

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

    r5473 r8607  
    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,char* response_descriptor,bool process_units,int weight_index);
    1212
    1313#endif  /* _RESPONSESXX_H */

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

    r5284 r8607  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     12void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){
    1313
    1414        /*Intermediary*/
    … …  
    2323        for (i=0;i<elements->Size();i++){
    2424                element=(Element*)elements->GetObjectByOffset(i);
    25                 J+=element->SurfaceAbsVelMisfit(process_units);
     25                J+=element->SurfaceAbsVelMisfit(process_units,weight_index);
    2626        }
    2727

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

    r5284 r8607  
    1010
    1111/* local prototypes: */
    12 void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index);
    1313
    1414#endif

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

    r5414 r8607  
    1111#include "../SurfaceAreax/SurfaceAreax.h"
    1212
    13 void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     13void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){
    1414
    1515        /*Intermediary*/
    … …  
    2727        for (i=0;i<elements->Size();i++){
    2828                element=(Element*)elements->GetObjectByOffset(i);
    29                 J+=element->SurfaceAverageVelMisfit(process_units);
     29                J+=element->SurfaceAverageVelMisfit(process_units,weight_index);
    3030        }
    3131

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

    r5284 r8607  
    1010
    1111/* local prototypes: */
    12 void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index);
    1313
    1414#endif

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

    r5284 r8607  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     12void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){
    1313
    1414        /*Intermediary*/
    … …  
    2323        for (i=0;i<elements->Size();i++){
    2424                element=(Element*)elements->GetObjectByOffset(i);
    25                 J+=element->SurfaceLogVelMisfit(process_units);
     25                J+=element->SurfaceLogVelMisfit(process_units,weight_index);
    2626        }
    2727

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

    r5284 r8607  
    1010
    1111/* local prototypes: */
    12 void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index);
    1313
    1414#endif

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

    r5284 r8607  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     12void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){
    1313
    1414        /*Intermediary*/
    … …  
    2323        for (i=0;i<elements->Size();i++){
    2424                element=(Element*)elements->GetObjectByOffset(i);
    25                 J+=element->SurfaceLogVxVyMisfit(process_units);
     25                J+=element->SurfaceLogVxVyMisfit(process_units,weight_index);
    2626        }
    2727

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

    r5284 r8607  
    1010
    1111/* local prototypes: */
    12 void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index);
    1313
    1414#endif

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

    r5284 r8607  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     12void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){
    1313
    1414        /*Intermediary*/
    … …  
    2323        for (i=0;i<elements->Size();i++){
    2424                element=(Element*)elements->GetObjectByOffset(i);
    25                 J+=element->SurfaceRelVelMisfit(process_units);
     25                J+=element->SurfaceRelVelMisfit(process_units,weight_index);
    2626        }
    2727

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

    r5284 r8607  
    1010
    1111/* local prototypes: */
    12 void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index);
    1313
    1414#endif

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

    r5286 r8607  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
     12void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units, int weight_index){
    1313
    1414        /*Intermediary*/
    … …  
    2323        for (i=0;i<elements->Size();i++){
    2424                element=(Element*)elements->GetObjectByOffset(i);
    25                 J+=element->ThicknessAbsMisfit(process_units);
     25                J+=element->ThicknessAbsMisfit(process_units,weight_index);
    2626        }
    2727

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

    r5286 r8607  
    1010
    1111/* local prototypes: */
    12 void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units);
     12void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index);
    1313
    1414#endif

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

    r8365 r8607  
    4343                virtual void   GradjDrag(Vec gradient)=0;
    4444                virtual void   GradjB(Vec gradient)=0;
    45                 virtual double ThicknessAbsMisfit(bool process_units)=0;
    46                 virtual double SurfaceAbsVelMisfit(bool process_units)=0;
    47                 virtual double SurfaceRelVelMisfit(bool process_units)=0;
    48                 virtual double SurfaceLogVelMisfit(bool process_units)=0;
    49                 virtual double SurfaceLogVxVyMisfit(bool process_units)=0;
    50                 virtual double SurfaceAverageVelMisfit(bool process_units)=0;
     45                virtual double ThicknessAbsMisfit(bool process_units  ,int weight_index)=0;
     46                virtual double SurfaceAbsVelMisfit(bool process_units ,int weight_index)=0;
     47                virtual double SurfaceRelVelMisfit(bool process_units ,int weight_index)=0;
     48                virtual double SurfaceLogVelMisfit(bool process_units ,int weight_index)=0;
     49                virtual double SurfaceLogVxVyMisfit(bool process_units,int weight_index)=0;
     50                virtual double SurfaceAverageVelMisfit(bool process_units,int weight_index)=0;
    5151                virtual double RegularizeInversion(void)=0;
    5252                virtual double SurfaceArea(void)=0;

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

    r8601 r8607  
    68456845/*}}}*/
    68466846/*FUNCTION Penta::SurfaceAverageVelMisfit {{{1*/
    6847 double Penta::SurfaceAverageVelMisfit(bool process_units){
     6847double Penta::SurfaceAverageVelMisfit(bool process_units,int weight_index){
    68486848
    68496849        int    approximation;
    … …  
    68686868                 * and compute SurfaceAverageVelMisfit*/
    68696869                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    6870                 J=tria->SurfaceAverageVelMisfit(process_units);
     6870                J=tria->SurfaceAverageVelMisfit(process_units,weight_index);
    68716871                delete tria->matice; delete tria;
    68726872                return J;
    … …  
    68756875
    68766876                tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    6877                 J=tria->SurfaceAverageVelMisfit(process_units);
     6877                J=tria->SurfaceAverageVelMisfit(process_units,weight_index);
    68786878                delete tria->matice; delete tria;
    68796879                return J;
    … …  
    68826882/*}}}*/
    68836883/*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/
    6884 double Penta::SurfaceAbsVelMisfit(bool process_units){
     6884double Penta::SurfaceAbsVelMisfit(bool process_units,int weight_index){
    68856885
    68866886        int    approximation;
    … …  
    69056905                 * and compute SurfaceAbsVelMisfit*/
    69066906                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    6907                 J=tria->SurfaceAbsVelMisfit(process_units);
     6907                J=tria->SurfaceAbsVelMisfit(process_units,weight_index);
    69086908                delete tria->matice; delete tria;
    69096909                return J;
    … …  
    69126912
    69136913                tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    6914                 J=tria->SurfaceAbsVelMisfit(process_units);
     6914                J=tria->SurfaceAbsVelMisfit(process_units,weight_index);
    69156915                delete tria->matice; delete tria;
    69166916                return J;
    … …  
    69196919/*}}}*/
    69206920/*FUNCTION Penta::SurfaceLogVelMisfit {{{1*/
    6921 double Penta::SurfaceLogVelMisfit(bool process_units){
     6921double Penta::SurfaceLogVelMisfit(bool process_units,int weight_index){
    69226922
    69236923        int    approximation;
    … …  
    69426942                 * and compute SurfaceLogVelMisfit*/
    69436943                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    6944                 J=tria->SurfaceLogVelMisfit(process_units);
     6944                J=tria->SurfaceLogVelMisfit(process_units,weight_index);
    69456945                delete tria->matice; delete tria;
    69466946                return J;
    … …  
    69496949
    69506950                tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    6951                 J=tria->SurfaceLogVelMisfit(process_units);
     6951                J=tria->SurfaceLogVelMisfit(process_units,weight_index);
    69526952                delete tria->matice; delete tria;
    69536953                return J;
    … …  
    69566956/*}}}*/
    69576957/*FUNCTION Penta::SurfaceLogVxVyMisfit {{{1*/
    6958 double Penta::SurfaceLogVxVyMisfit(bool process_units){
     6958double Penta::SurfaceLogVxVyMisfit(bool process_units,int weight_index){
    69596959
    69606960        double J;
    … …  
    69816981                 * and compute SurfaceLogVxVyMisfit*/
    69826982                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    6983                 J=tria->SurfaceLogVxVyMisfit(process_units);
     6983                J=tria->SurfaceLogVxVyMisfit(process_units,weight_index);
    69846984                delete tria->matice; delete tria;
    69856985                return J;
    … …  
    69886988
    69896989                tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    6990                 J=tria->SurfaceLogVxVyMisfit(process_units);
     6990                J=tria->SurfaceLogVxVyMisfit(process_units,weight_index);
    69916991                delete tria->matice; delete tria;
    69926992                return J;
    … …  
    70197019/*}}}*/
    70207020/*FUNCTION Penta::SurfaceRelVelMisfit {{{1*/
    7021 double Penta::SurfaceRelVelMisfit(bool process_units){
     7021double Penta::SurfaceRelVelMisfit(bool process_units,int weight_index){
    70227022
    70237023        int    approximation;
    … …  
    70427042                 * and compute SurfaceRelVelMisfit*/
    70437043                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    7044                 J=tria->SurfaceRelVelMisfit(process_units);
     7044                J=tria->SurfaceRelVelMisfit(process_units,weight_index);
    70457045                delete tria->matice; delete tria;
    70467046                return J;
    … …  
    70497049
    70507050                tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    7051                 J=tria->SurfaceRelVelMisfit(process_units);
     7051                J=tria->SurfaceRelVelMisfit(process_units,weight_index);
    70527052                delete tria->matice; delete tria;
    70537053                return J;
    … …  
    71007100}
    71017101/*FUNCTION Penta::ThicknessAbsMisfit {{{1*/
    7102 double Penta::ThicknessAbsMisfit(bool process_units){
     7102double Penta::ThicknessAbsMisfit(bool process_units,int weight_index){
    71037103
    71047104        int    approximation;
    … …  
    71137113        _error_("Not implemented yet");
    71147114
    7115         tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    7116         J=tria->ThicknessAbsMisfit(process_units);
     7115        tria=(Tria*)SpawnTria(0,1,2);
     7116        J=tria->ThicknessAbsMisfit(process_units,weight_index);
    71177117        delete tria->matice; delete tria;
    71187118        return J;

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

    r8496 r8607  
    118118                void   MinVy(double* pminvy, bool process_units);
    119119                void   MinVz(double* pminvz, bool process_units);
    120                 double ThicknessAbsMisfit(bool process_units);
    121                 double SurfaceAbsVelMisfit(bool process_units);
    122                 double SurfaceRelVelMisfit(bool process_units);
    123                 double SurfaceLogVelMisfit(bool process_units);
    124                 double SurfaceLogVxVyMisfit(bool process_units);
    125                 double SurfaceAverageVelMisfit(bool process_units);
     120                double ThicknessAbsMisfit(     bool process_units,int weight_index);
     121                double SurfaceAbsVelMisfit(    bool process_units,int weight_index);
     122                double SurfaceRelVelMisfit(    bool process_units,int weight_index);
     123                double SurfaceLogVelMisfit(    bool process_units,int weight_index);
     124                double SurfaceLogVxVyMisfit(   bool process_units,int weight_index);
     125                double SurfaceAverageVelMisfit(bool process_units,int weight_index);
    126126                void   PatchFill(int* pcount, Patch* patch);
    127127                void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);

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

    r8603 r8607  
    17061706
    17071707                /*Get all parameters at gaussian point*/
    1708                 weights_input->GetParameterValue(&weight,gauss,0);
    17091708                vx_input->GetParameterValue(&vx,gauss);
    17101709                vy_input->GetParameterValue(&vy,gauss);
    … …  
    17141713
    17151714                /*Loop over all requested responses*/
    1716                 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    1717 
    1718                         case SurfaceAbsVelMisfitEnum:
    1719                                 /*
    1720                                  *      1  [           2              2 ]
    1721                                  * J = --- | (u - u   )  +  (v - v   )  |
    1722                                  *      2  [       obs            obs   ]
    1723                                  *
    1724                                  *        dJ
    1725                                  * DU = - -- = (u   - u )
    1726                                  *        du     obs
    1727                                  */
    1728                                 for (i=0;i<NUMVERTICES;i++){
    1729                                         dux=vxobs-vx;
    1730                                         duy=vyobs-vy;
    1731                                         pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1732                                         pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1733                                 }
    1734                                 break;
    1735                         case SurfaceRelVelMisfitEnum:
    1736                                 /*
    1737                                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    1738                                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    1739                                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    1740                                  *              obs                        obs                     
    1741                                  *
    1742                                  *        dJ     \bar{v}^2
    1743                                  * DU = - -- = ------------- (u   - u )
    1744                                  *        du   (u   + eps)^2    obs
    1745                                  *               obs
    1746                                  */
    1747                                 for (i=0;i<NUMVERTICES;i++){
    1748                                         scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
    1749                                         scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
    1750                                         dux=scalex*(vxobs-vx);
    1751                                         duy=scaley*(vyobs-vy);
    1752                                         pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1753                                         pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1754                                 }
    1755                                 break;
    1756                         case SurfaceLogVelMisfitEnum:
    1757                                 /*
    1758                                  *                 [        vel + eps     ] 2
    1759                                  * J = 4 \bar{v}^2 | log ( -----------  ) | 
    1760                                  *                 [       vel   + eps    ]
    1761                                  *                            obs
    1762                                  *
    1763                                  *        dJ                 2 * log(...)
    1764                                  * DU = - -- = - 4 \bar{v}^2 -------------  u
    1765                                  *        du                 vel^2 + eps
    1766                                  *           
    1767                                  */
    1768                                 for (i=0;i<NUMVERTICES;i++){
    1769                                         velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
    1770                                         obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
    1771                                         scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
    1772                                         dux=scale*vx;
    1773                                         duy=scale*vy;
    1774                                         pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1775                                         pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1776                                 }
    1777                                 break;
    1778                         case SurfaceAverageVelMisfitEnum:
    1779                                 /*
    1780                                  *      1                    2              2
    1781                                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    1782                                  *      S                obs            obs
    1783                                  *
    1784                                  *        dJ      1       1
    1785                                  * DU = - -- = - --- ----------- * 2 (u - u   )
    1786                                  *        du      S  2 sqrt(...)           obs
    1787                                  */
    1788                                 for (i=0;i<NUMVERTICES;i++){
    1789                                         scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
    1790                                         dux=scale*(vxobs-vx);
    1791                                         duy=scale*(vyobs-vy);
    1792                                         pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1793                                         pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1794                                 }
    1795                                 break;
    1796                         case SurfaceLogVxVyMisfitEnum:
    1797                                 /*
    1798                                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    1799                                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    1800                                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    1801                                  *                              obs                       obs
    1802                                  *        dJ                              1      u                             1
    1803                                  * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    1804                                  *        du                         |u| + eps  |u|                           u + eps
    1805                                  */
    1806                                 for (i=0;i<NUMVERTICES;i++){
    1807                                         dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    1808                                         duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    1809                                         pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1810                                         pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1811                                 }
    1812                                 break;
    1813                         default:
    1814                                 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
     1715                for(resp=0;resp<num_responses;resp++){
     1716
     1717                        weights_input->GetParameterValue(&weight,gauss,resp);
     1718
     1719                        switch(responses[resp]){
     1720                                case SurfaceAbsVelMisfitEnum:
     1721                                        /*
     1722                                         *      1  [           2              2 ]
     1723                                         * J = --- | (u - u   )  +  (v - v   )  |
     1724                                         *      2  [       obs            obs   ]
     1725                                         *
     1726                                         *        dJ
     1727                                         * DU = - -- = (u   - u )
     1728                                         *        du     obs
     1729                                         */
     1730                                        for (i=0;i<NUMVERTICES;i++){
     1731                                                dux=vxobs-vx;
     1732                                                duy=vyobs-vy;
     1733                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1734                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1735                                        }
     1736                                        break;
     1737                                case SurfaceRelVelMisfitEnum:
     1738                                        /*
     1739                                         *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     1740                                         * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     1741                                         *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     1742                                         *              obs                        obs                     
     1743                                         *
     1744                                         *        dJ     \bar{v}^2
     1745                                         * DU = - -- = ------------- (u   - u )
     1746                                         *        du   (u   + eps)^2    obs
     1747                                         *               obs
     1748                                         */
     1749                                        for (i=0;i<NUMVERTICES;i++){
     1750                                                scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
     1751                                                scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     1752                                                dux=scalex*(vxobs-vx);
     1753                                                duy=scaley*(vyobs-vy);
     1754                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1755                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1756                                        }
     1757                                        break;
     1758                                case SurfaceLogVelMisfitEnum:
     1759                                        /*
     1760                                         *                 [        vel + eps     ] 2
     1761                                         * J = 4 \bar{v}^2 | log ( -----------  ) | 
     1762                                         *                 [       vel   + eps    ]
     1763                                         *                            obs
     1764                                         *
     1765                                         *        dJ                 2 * log(...)
     1766                                         * DU = - -- = - 4 \bar{v}^2 -------------  u
     1767                                         *        du                 vel^2 + eps
     1768                                         *           
     1769                                         */
     1770                                        for (i=0;i<NUMVERTICES;i++){
     1771                                                velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
     1772                                                obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
     1773                                                scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
     1774                                                dux=scale*vx;
     1775                                                duy=scale*vy;
     1776                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1777                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1778                                        }
     1779                                        break;
     1780                                case SurfaceAverageVelMisfitEnum:
     1781                                        /*
     1782                                         *      1                    2              2
     1783                                         * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     1784                                         *      S                obs            obs
     1785                                         *
     1786                                         *        dJ      1       1
     1787                                         * DU = - -- = - --- ----------- * 2 (u - u   )
     1788                                         *        du      S  2 sqrt(...)           obs
     1789                                         */
     1790                                        for (i=0;i<NUMVERTICES;i++){
     1791                                                scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
     1792                                                dux=scale*(vxobs-vx);
     1793                                                duy=scale*(vyobs-vy);
     1794                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1795                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1796                                        }
     1797                                        break;
     1798                                case SurfaceLogVxVyMisfitEnum:
     1799                                        /*
     1800                                         *      1            [        |u| + eps     2          |v| + eps     2  ]
     1801                                         * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     1802                                         *      2            [       |u    |+ eps              |v    |+ eps     ]
     1803                                         *                              obs                       obs
     1804                                         *        dJ                              1      u                             1
     1805                                         * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     1806                                         *        du                         |u| + eps  |u|                           u + eps
     1807                                         */
     1808                                        for (i=0;i<NUMVERTICES;i++){
     1809                                                dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     1810                                                duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     1811                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1812                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1813                                        }
     1814                                        break;
     1815                                default:
     1816                                        _error_("response %s not supported yet",EnumToStringx(responses[resp]));
     1817                        }
    18151818                }
    18161819        }
    … …  
    18701873
    18711874                /*Get all parameters at gaussian point*/
    1872                 weights_input->GetParameterValue(&weight,gauss,0);
    18731875                vx_input->GetParameterValue(&vx,gauss);
    18741876                vy_input->GetParameterValue(&vy,gauss);
    … …  
    18781880
    18791881                /*Loop over all requested responses*/
    1880                 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    1881 
    1882                         case SurfaceAbsVelMisfitEnum:
    1883                                 /*
    1884                                  *      1  [           2              2 ]
    1885                                  * J = --- | (u - u   )  +  (v - v   )  |
    1886                                  *      2  [       obs            obs   ]
    1887                                  *
    1888                                  *        dJ
    1889                                  * DU = - -- = (u   - u )
    1890                                  *        du     obs
    1891                                  */
    1892                                 for (i=0;i<NUMVERTICES;i++){
    1893                                         dux=vxobs-vx;
    1894                                         duy=vyobs-vy;
    1895                                         pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1896                                         pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1897                                 }
    1898                                 break;
    1899                         case SurfaceRelVelMisfitEnum:
    1900                                 /*
    1901                                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    1902                                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    1903                                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    1904                                  *              obs                        obs                     
    1905                                  *
    1906                                  *        dJ     \bar{v}^2
    1907                                  * DU = - -- = ------------- (u   - u )
    1908                                  *        du   (u   + eps)^2    obs
    1909                                  *               obs
    1910                                  */
    1911                                 for (i=0;i<NUMVERTICES;i++){
    1912                                         scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
    1913                                         scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
    1914                                         dux=scalex*(vxobs-vx);
    1915                                         duy=scaley*(vyobs-vy);
    1916                                         pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1917                                         pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1918                                 }
    1919                                 break;
    1920                         case SurfaceLogVelMisfitEnum:
    1921                                 /*
    1922                                  *                 [        vel + eps     ] 2
    1923                                  * J = 4 \bar{v}^2 | log ( -----------  ) | 
    1924                                  *                 [       vel   + eps    ]
    1925                                  *                            obs
    1926                                  *
    1927                                  *        dJ                 2 * log(...)
    1928                                  * DU = - -- = - 4 \bar{v}^2 -------------  u
    1929                                  *        du                 vel^2 + eps
    1930                                  *           
    1931                                  */
    1932                                 for (i=0;i<NUMVERTICES;i++){
    1933                                         velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
    1934                                         obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
    1935                                         scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
    1936                                         dux=scale*vx;
    1937                                         duy=scale*vy;
    1938                                         pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1939                                         pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1940                                 }
    1941                                 break;
    1942                         case SurfaceAverageVelMisfitEnum:
    1943                                 /*
    1944                                  *      1                    2              2
    1945                                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    1946                                  *      S                obs            obs
    1947                                  *
    1948                                  *        dJ      1       1
    1949                                  * DU = - -- = - --- ----------- * 2 (u - u   )
    1950                                  *        du      S  2 sqrt(...)           obs
    1951                                  */
    1952                                 for (i=0;i<NUMVERTICES;i++){
    1953                                         scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
    1954                                         dux=scale*(vxobs-vx);
    1955                                         duy=scale*(vyobs-vy);
    1956                                         pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1957                                         pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1958                                 }
    1959                                 break;
    1960                         case SurfaceLogVxVyMisfitEnum:
    1961                                 /*
    1962                                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    1963                                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    1964                                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    1965                                  *                              obs                       obs
    1966                                  *        dJ                              1      u                             1
    1967                                  * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    1968                                  *        du                         |u| + eps  |u|                           u + eps
    1969                                  */
    1970                                 for (i=0;i<NUMVERTICES;i++){
    1971                                         dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    1972                                         duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    1973                                         pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1974                                         pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
    1975                                 }
    1976                                 break;
    1977                         default:
    1978                                 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
     1882                for(resp=0;resp<num_responses;resp++){
     1883
     1884                        weights_input->GetParameterValue(&weight,gauss,resp);
     1885
     1886                        switch(responses[resp]){
     1887
     1888                                case SurfaceAbsVelMisfitEnum:
     1889                                        /*
     1890                                         *      1  [           2              2 ]
     1891                                         * J = --- | (u - u   )  +  (v - v   )  |
     1892                                         *      2  [       obs            obs   ]
     1893                                         *
     1894                                         *        dJ
     1895                                         * DU = - -- = (u   - u )
     1896                                         *        du     obs
     1897                                         */
     1898                                        for (i=0;i<NUMVERTICES;i++){
     1899                                                dux=vxobs-vx;
     1900                                                duy=vyobs-vy;
     1901                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1902                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1903                                        }
     1904                                        break;
     1905                                case SurfaceRelVelMisfitEnum:
     1906                                        /*
     1907                                         *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     1908                                         * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     1909                                         *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     1910                                         *              obs                        obs                     
     1911                                         *
     1912                                         *        dJ     \bar{v}^2
     1913                                         * DU = - -- = ------------- (u   - u )
     1914                                         *        du   (u   + eps)^2    obs
     1915                                         *               obs
     1916                                         */
     1917                                        for (i=0;i<NUMVERTICES;i++){
     1918                                                scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
     1919                                                scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     1920                                                dux=scalex*(vxobs-vx);
     1921                                                duy=scaley*(vyobs-vy);
     1922                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1923                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1924                                        }
     1925                                        break;
     1926                                case SurfaceLogVelMisfitEnum:
     1927                                        /*
     1928                                         *                 [        vel + eps     ] 2
     1929                                         * J = 4 \bar{v}^2 | log ( -----------  ) | 
     1930                                         *                 [       vel   + eps    ]
     1931                                         *                            obs
     1932                                         *
     1933                                         *        dJ                 2 * log(...)
     1934                                         * DU = - -- = - 4 \bar{v}^2 -------------  u
     1935                                         *        du                 vel^2 + eps
     1936                                         *           
     1937                                         */
     1938                                        for (i=0;i<NUMVERTICES;i++){
     1939                                                velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
     1940                                                obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
     1941                                                scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
     1942                                                dux=scale*vx;
     1943                                                duy=scale*vy;
     1944                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1945                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1946                                        }
     1947                                        break;
     1948                                case SurfaceAverageVelMisfitEnum:
     1949                                        /*
     1950                                         *      1                    2              2
     1951                                         * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     1952                                         *      S                obs            obs
     1953                                         *
     1954                                         *        dJ      1       1
     1955                                         * DU = - -- = - --- ----------- * 2 (u - u   )
     1956                                         *        du      S  2 sqrt(...)           obs
     1957                                         */
     1958                                        for (i=0;i<NUMVERTICES;i++){
     1959                                                scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
     1960                                                dux=scale*(vxobs-vx);
     1961                                                duy=scale*(vyobs-vy);
     1962                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1963                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1964                                        }
     1965                                        break;
     1966                                case SurfaceLogVxVyMisfitEnum:
     1967                                        /*
     1968                                         *      1            [        |u| + eps     2          |v| + eps     2  ]
     1969                                         * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     1970                                         *      2            [       |u    |+ eps              |v    |+ eps     ]
     1971                                         *                              obs                       obs
     1972                                         *        dJ                              1      u                             1
     1973                                         * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     1974                                         *        du                         |u| + eps  |u|                           u + eps
     1975                                         */
     1976                                        for (i=0;i<NUMVERTICES;i++){
     1977                                                dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     1978                                                duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     1979                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
     1980                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1981                                        }
     1982                                        break;
     1983                                default:
     1984                                        _error_("response %s not supported yet",EnumToStringx(responses[resp]));
     1985                        }
    19791986                }
    19801987        }
    … …  
    48204827/*}}}*/
    48214828/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/
    4822 double Tria::SurfaceAbsVelMisfit(bool process_units){
     4829double Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){
    48234830
    48244831        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    48544861
    48554862                /*Get all parameters at gaussian point*/
    4856                 weights_input->GetParameterValue(&weight,gauss,0);
     4863                weights_input->GetParameterValue(&weight,gauss,weight_index);
    48574864                vx_input->GetParameterValue(&vx,gauss);
    48584865                vy_input->GetParameterValue(&vy,gauss);
    … …  
    49104917/*}}}*/
    49114918/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/
    4912 double Tria::SurfaceAverageVelMisfit(bool process_units){
     4919double Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){
    49134920
    49144921        const int    numdof=2*NUMVERTICES;
    … …  
    49454952
    49464953                /*Get all parameters at gaussian point*/
    4947                 weights_input->GetParameterValue(&weight,gauss,0);
     4954                weights_input->GetParameterValue(&weight,gauss,weight_index);
    49484955                vx_input->GetParameterValue(&vx,gauss);
    49494956                vy_input->GetParameterValue(&vy,gauss);
    … …  
    49714978/*}}}*/
    49724979/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/
    4973 double Tria::SurfaceLogVelMisfit(bool process_units){
     4980double Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){
    49744981
    49754982        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    50085015
    50095016                /*Get all parameters at gaussian point*/
    5010                 weights_input->GetParameterValue(&weight,gauss,0);
     5017                weights_input->GetParameterValue(&weight,gauss,weight_index);
    50115018                vx_input->GetParameterValue(&vx,gauss);
    50125019                vy_input->GetParameterValue(&vy,gauss);
    … …  
    50365043/*}}}*/
    50375044/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
    5038 double Tria::SurfaceLogVxVyMisfit(bool process_units){
     5045double Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){
    50395046
    50405047        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    50735080
    50745081                /*Get all parameters at gaussian point*/
    5075                 weights_input->GetParameterValue(&weight,gauss,0);
     5082                weights_input->GetParameterValue(&weight,gauss,weight_index);
    50765083                vx_input->GetParameterValue(&vx,gauss);
    50775084                vy_input->GetParameterValue(&vy,gauss);
    … …  
    51265133/*}}}*/
    51275134/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
    5128 double Tria::SurfaceRelVelMisfit(bool process_units){
     5135double Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){
    51295136        const int  numdof=2*NUMVERTICES;
    51305137
    … …  
    51625169
    51635170                /*Get all parameters at gaussian point*/
    5164                 weights_input->GetParameterValue(&weight,gauss,0);
     5171                weights_input->GetParameterValue(&weight,gauss,weight_index);
    51655172                vx_input->GetParameterValue(&vx,gauss);
    51665173                vy_input->GetParameterValue(&vy,gauss);
    … …  
    51905197/*}}}*/
    51915198/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
    5192 double Tria::ThicknessAbsMisfit(bool process_units){
     5199double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){
    51935200
    51945201        /* Constants */
    … …  
    52305237                thickness_input->GetParameterDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
    52315238                thicknessobs_input->GetParameterValue(&thicknessobs,gauss);
    5232                 weights_input->GetParameterValue(&weight,gauss,0);
     5239                weights_input->GetParameterValue(&weight,gauss,weight_index);
    52335240
    52345241                /*compute ThicknessAbsMisfit*/

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

    r8592 r8607  
    123123                void   MinVy(double* pminvy, bool process_units);
    124124                void   MinVz(double* pminvz, bool process_units);
    125                 double ThicknessAbsMisfit(bool process_units);
    126                 double SurfaceAbsVelMisfit(bool process_units);
    127                 double SurfaceRelVelMisfit(bool process_units);
    128                 double SurfaceLogVelMisfit(bool process_units);
    129                 double SurfaceLogVxVyMisfit(bool process_units);
    130                 double SurfaceAverageVelMisfit(bool process_units);
     125                double ThicknessAbsMisfit(     bool process_units,int weight_index);
     126                double SurfaceAbsVelMisfit(    bool process_units,int weight_index);
     127                double SurfaceRelVelMisfit(    bool process_units,int weight_index);
     128                double SurfaceLogVelMisfit(    bool process_units,int weight_index);
     129                double SurfaceLogVxVyMisfit(   bool process_units,int weight_index);
     130                double SurfaceAverageVelMisfit(bool process_units,int weight_index);
    131131                void   PatchFill(int* pcount, Patch* patch);
    132132                void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);

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

    r8429 r8607  
    113113        //VecScale(gradient,-1.);
    114114        VecCopy(gradient,G);
    115         CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum);
     115        CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    116116
    117117        //printf("X\n");

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

    r8601 r8607  
    2525       
    2626        /*output: */
    27         double J=0,Jplus;
     27        double J;
    2828       
    2929        /*parameters: */
    30         int        num_responses;
    3130        int        solution_type,analysis_type;
    3231        bool       isstokes       = false;
    3332        bool       conserve_loads = true;
    34         int       *responses      = NULL;
    3533        FemModel  *femmodel       = NULL;
    3634
    … …  
    3937
    4038        /*Recover parameters: */
    41         femmodel->parameters->FindParam(&num_responses,NumResponsesEnum);
    42         femmodel->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
    4339        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    4440        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    7470
    7571        /*Compute misfit for this velocity field.*/
    76         for(int i=0;i<num_responses;i++){
    77                 CostFunctionx(&Jplus, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,responses[i]);
    78                 J+=Jplus;
    79         }
     72        CostFunctionx(&J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    8073
    8174        /*Free ressources:*/
    82         xfree((void**)&responses);
    8375        return J;
    8476}

  • issm/trunk/src/m/solutions/objectivefunctionC.m

    r8602 r8607  
    66
    77%recover some parameters
    8 num_responses = femmodel.parameters.NumResponses;
    9 responses     = femmodel.parameters.StepResponses;
    108analysis_type = femmodel.parameters.AnalysisType;
    119solution_type = femmodel.parameters.SolutionType;
    … …  
    3634
    3735%Compute misfit for this velocity field
    38 for i=1:num_responses
    39         J=J+CostFunction(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials, femmodel.parameters,responses(i));
    40 end
     36J=CostFunction(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials, femmodel.parameters);

  • issm/trunk/src/mex/CostFunction/CostFunction.cpp

    r6412 r8607  
    1414        Materials  *materials  = NULL;
    1515        Parameters *parameters = NULL;
    16         int         response;
    1716
    1817        /* output datasets: */
    … …  
    3231        FetchData((DataSet**)&materials,MATERIALS);
    3332        FetchParams(&parameters,PARAMETERS);
    34         FetchData(&response,RESPONSE);
    3533
    3634        /*configure: */
    … …  
    4038
    4139        /*!Call core code: */
    42         CostFunctionx(&J, elements,nodes,vertices, loads,materials,parameters,response);
     40        CostFunctionx(&J, elements,nodes,vertices, loads,materials,parameters);
    4341
    4442        /*write output : */
    … …  
    6058{
    6159        _printf_(true,"\n");
    62         _printf_(true,"   usage: [J] = %s(elements,nodes,vertices,loads, materials, parameters,response);\n",__FUNCT__);
     60        _printf_(true,"   usage: [J] = %s(elements,nodes,vertices,loads, materials, parameters);\n",__FUNCT__);
    6361        _printf_(true,"\n");
    6462}

  • issm/trunk/src/mex/CostFunction/CostFunction.h

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

  • issm/trunk/src/mex/Response/Response.cpp

    r6716 r8607  
    1616        char       *response   = NULL;
    1717        bool        process_units;
     18        int         weight_index;
    1819
    1920        /* output datasets: */
    … …  
    3536        FetchData(&response,RESPONSE);
    3637        FetchData(&process_units,PROCESSUNITS);
     38        FetchData(&weight_index,WEIGHTINDEX);
    3739
    3840        /*configure: */
    … …  
    4244
    4345        /*!Call core code: */
    44         Responsex(&resp, elements,nodes,vertices, loads,materials,parameters,response,process_units);
     46        Responsex(&resp, elements,nodes,vertices, loads,materials,parameters,response,process_units,weight_index);
    4547
    4648        /*write output : */

  • issm/trunk/src/mex/Response/Response.h

    r6689 r8607  
    2626#define RESPONSE (mxArray*)prhs[6]
    2727#define PROCESSUNITS (mxArray*)prhs[7]
     28#define WEIGHTINDEX (mxArray*)prhs[8]
    2829
    2930/* serial output macros: */
    … …  
    3435#define NLHS  1
    3536#undef NRHS
    36 #define NRHS  8
     37#define NRHS  9
    3738
    3839#endif  /* _RESPONSE_H */
Note: See TracChangeset for help on using the changeset viewer.