Changeset 18058


Ignore:
Timestamp:
05/25/14 18:01:37 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on moving gradient to adjoint

Location:
issm/trunk-jpl/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

    r18057 r18058  
    974974}/*}}}*/
    975975void AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    976         _error_("not implemented yet");
     976
     977        /*return if floating*/
     978        if(element->IsFloating())return;
     979
     980        /*Intermediaries*/
     981        int      domaintype,dim;
     982        Element* basalelement;
     983
     984        /*Get basal element*/
     985        element->FindParam(&domaintype,DomainTypeEnum);
     986        switch(domaintype){
     987                case Domain2DhorizontalEnum:
     988                        basalelement = element;
     989                        dim          = 2;
     990                        break;
     991                case Domain2DverticalEnum:
     992                        if(!element->IsOnBase()) return;
     993                        basalelement = element->SpawnBasalElement();
     994                        dim          = 1;
     995                        break;
     996                case Domain3DEnum:
     997                        if(!element->IsOnBase()) return;
     998                        basalelement = element->SpawnBasalElement();
     999                        dim          = 2;
     1000                        break;
     1001                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     1002        }
     1003
     1004        /*Intermediaries*/
     1005        IssmDouble Jdet,weight;
     1006        IssmDouble dk[3];
     1007        IssmDouble *xyz_list= NULL;
     1008
     1009        /*Fetch number of vertices for this finite element*/
     1010        int numvertices = basalelement->GetNumberOfVertices();
     1011
     1012        /*Initialize some vectors*/
     1013        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
     1014        IssmDouble* ge            = xNew<IssmDouble>(numvertices);
     1015        int*        vertexpidlist = xNew<int>(numvertices);
     1016
     1017        /*Retrieve all inputs we will be needing: */
     1018        basalelement->GetVerticesCoordinates(&xyz_list);
     1019        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
     1020        Input* dragcoefficient_input = basalelement->GetInput(FrictionCoefficientEnum);                _assert_(dragcoefficient_input);
     1021        Input* weights_input         = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1022
     1023        /* Start  looping on the number of gaussian points: */
     1024        Gauss* gauss=basalelement->NewGauss(2);
     1025        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1026                gauss->GaussPoint(ig);
     1027
     1028                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1029                basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
     1030                weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
     1031
     1032                /*Build alpha_complement_list: */
     1033                dragcoefficient_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
     1034
     1035                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     1036                for(int i=0;i<numvertices;i++){
     1037                        if(dim==2){
     1038                                ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
     1039                        }
     1040                        else{
     1041                                ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
     1042                        }
     1043                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     1044                }
     1045        }
     1046        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1047
     1048        /*Clean up and return*/
     1049        xDelete<IssmDouble>(xyz_list);
     1050        xDelete<IssmDouble>(dbasis);
     1051        xDelete<IssmDouble>(ge);
     1052        xDelete<int>(vertexpidlist);
     1053        delete gauss;
     1054        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1055
    9771056}/*}}}*/
    9781057void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r18039 r18058  
    603603        return z;
    604604}/*}}}*/
     605void       Element::GradientIndexing(int* indexing,int control_index){/*{{{*/
     606
     607        /*Get number of controls*/
     608        int num_controls;
     609        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     610
     611        /*Get number of vertices*/
     612        int numvertices = this->GetNumberOfVertices();
     613
     614        /*get gradient indices*/
     615        for(int i=0;i<numvertices;i++){
     616                indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
     617        }
     618
     619}
     620/*}}}*/
    605621bool       Element::HasNodeOnBase(){/*{{{*/
    606622        return (this->inputs->Max(MeshVertexonbaseEnum)>0.);

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

    r18001 r18058  
    9494                IssmDouble GetYcoord(IssmDouble* xyz_list,Gauss* gauss);
    9595                IssmDouble GetZcoord(IssmDouble* xyz_list,Gauss* gauss);
     96                void       GradientIndexing(int* indexing,int control_index);
    9697                bool       HasNodeOnBase();
    9798                bool       HasNodeOnSurface();

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

    r18001 r18058  
    28952895
    28962896}/*}}}*/
    2897 void       Penta::GradientIndexing(int* indexing,int control_index){/*{{{*/
    2898 
    2899         /*Get some parameters*/
    2900         int num_controls;
    2901         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    2902 
    2903         /*get gradient indices*/
    2904         for(int i=0;i<NUMVERTICES;i++){
    2905                 indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
    2906         }
    2907 
    2908 }
    2909 /*}}}*/
    29102897void       Penta::Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){/*{{{*/
    29112898        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/

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

    r18001 r18058  
    131131                #endif
    132132
    133                 void   GradientIndexing(int* indexing,int control_index);
    134133                void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    135134                void   GradjDragSSA(Vector<IssmDouble>* gradient,int control_index);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.cpp

    r17972 r18058  
    116116
    117117}/*}}}*/
    118 void       Seg::GradientIndexing(int* indexing,int control_index){/*{{{*/
    119 
    120         /*Get some parameters*/
    121         int num_controls;
    122         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    123 
    124         /*get gradient indices*/
    125         for(int i=0;i<NUMVERTICES;i++){
    126                 indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
    127         }
    128 
    129 }
    130 /*}}}*/
    131118void       Seg::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    132119

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r18001 r18058  
    167167#endif
    168168
    169                 void       GradientIndexing(int* indexing,int control_index);
    170169                void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");};
    171170                void       GradjBGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};

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

    r18056 r18058  
    34033403}
    34043404/*}}}*/
    3405 void       Tria::GradientIndexing(int* indexing,int control_index){/*{{{*/
    3406 
    3407         /*Get some parameters*/
    3408         int num_controls;
    3409         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    3410 
    3411         /*get gradient indices*/
    3412         for(int i=0;i<NUMVERTICES;i++){
    3413                 indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
    3414         }
    3415 
    3416 }
    3417 /*}}}*/
    34183405void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    34193406

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

    r18003 r18058  
    136136                #endif
    137137
    138                 void       GradientIndexing(int* indexing,int control_index);
    139138                void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    140139                void       GradjBGradient(Vector<IssmDouble>* gradient,int control_index);
Note: See TracChangeset for help on using the changeset viewer.