Changeset 16862


Ignore:
Timestamp:
11/21/13 15:03:46 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: Adjoint BalH

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

Legend:

Unmodified
Added
Removed

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

    r16782 r16862  
    3030}/*}}}*/
    3131ElementVector* AdjointBalancethicknessAnalysis::CreatePVector(Element* element){/*{{{*/
    32 _error_("not implemented yet");
     32
     33        /*Intermediaries*/
     34        int      meshtype;
     35        Element* basalelement;
     36
     37        /*Get basal element*/
     38        element->FindParam(&meshtype,MeshTypeEnum);
     39        switch(meshtype){
     40                case Mesh2DhorizontalEnum:
     41                        basalelement = element;
     42                        break;
     43                case Mesh3DEnum:
     44                        if(!element->IsOnBed()) return NULL;
     45                        basalelement = element->SpawnBasalElement();
     46                        break;
     47                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     48        }
     49
     50        /*Intermediaries */
     51        int         num_responses,i;
     52        IssmDouble  dH[2];
     53        IssmDouble  vx,vy,vel,Jdet;
     54        IssmDouble  thickness,thicknessobs,weight;
     55        int        *responses = NULL;
     56        IssmDouble *xyz_list  = NULL;
     57
     58        /*Fetch number of nodes and dof for this finite element*/
     59        int numnodes = basalelement->GetNumberOfNodes();
     60
     61        /*Initialize Element vector and vectors*/
     62        ElementVector* pe     = basalelement->NewElementVector(SSAApproximationEnum);
     63        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     64        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     65
     66        /*Retrieve all inputs and parameters*/
     67        basalelement->GetVerticesCoordinates(&xyz_list);
     68        basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     69        basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum);
     70        Input* thickness_input    = basalelement->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     71        Input* thicknessobs_input = basalelement->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
     72        Input* weights_input      = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     73        Input* vx_input           = basalelement->GetInput(VxEnum);                                 _assert_(vx_input);
     74        Input* vy_input           = basalelement->GetInput(VyEnum);                                 _assert_(vy_input);
     75
     76        /* Start  looping on the number of gaussian points: */
     77        Gauss* gauss=basalelement->NewGauss(2);
     78        for(int ig=gauss->begin();ig<gauss->end();ig++){
     79                gauss->GaussPoint(ig);
     80
     81                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     82                basalelement->NodalFunctions(basis,gauss);
     83                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     84
     85                thickness_input->GetInputValue(&thickness, gauss);
     86                thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
     87                thicknessobs_input->GetInputValue(&thicknessobs, gauss);
     88
     89                /*Loop over all requested responses*/
     90                for(int resp=0;resp<num_responses;resp++){
     91                        weights_input->GetInputValue(&weight,gauss,responses[resp]);
     92
     93                        switch(responses[resp]){
     94                                case ThicknessAbsMisfitEnum:
     95                                        for(i=0;i<numnodes;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
     96                                        break;
     97                                case ThicknessAbsGradientEnum:
     98                                        for(i=0;i<numnodes;i++) pe->values[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight;
     99                                        for(i=0;i<numnodes;i++) pe->values[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight;
     100                                        break;
     101                                case ThicknessAlongGradientEnum:
     102                                        vx_input->GetInputValue(&vx,gauss);
     103                                        vy_input->GetInputValue(&vy,gauss);
     104                                        vel = sqrt(vx*vx+vy*vy);
     105                                        vx  = vx/(vel+1.e-9);
     106                                        vy  = vy/(vel+1.e-9);
     107                                        for(i=0;i<numnodes;i++) pe->values[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight;
     108                                        break;
     109                                case ThicknessAcrossGradientEnum:
     110                                        vx_input->GetInputValue(&vx,gauss);
     111                                        vy_input->GetInputValue(&vy,gauss);
     112                                        vel = sqrt(vx*vx+vy*vy);
     113                                        vx  = vx/(vel+1.e-9);
     114                                        vy  = vy/(vel+1.e-9);
     115                                        for(i=0;i<numnodes;i++) pe->values[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight;
     116                                        break;
     117                                default:
     118                                        _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
     119                        }
     120                }
     121        }
     122
     123        /*Clean up and return*/
     124        xDelete<IssmDouble>(xyz_list);
     125        xDelete<IssmDouble>(basis);
     126        xDelete<IssmDouble>(dbasis);
     127        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     128        delete gauss;
     129        return pe;
    33130}/*}}}*/
    34131void AdjointBalancethicknessAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16782 r16862  
    4141}/*}}}*/
    4242ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/
    43 _error_("not implemented yet");
     43        return NULL;
    4444}/*}}}*/
    4545void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16831 r16862  
    121121                case Mesh2DhorizontalEnum:
    122122                        basalelement = element;
     123                        break;
     124                case Mesh2DverticalEnum:
     125                        if(!element->IsOnBed()) return NULL;
     126                        basalelement = element->SpawnBasalElement();
    123127                        break;
    124128                case Mesh3DEnum:

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

    r16859 r16862  
    20542054        /*Retrieve all inputs and parameters*/
    20552055        element->GetVerticesCoordinates(&xyz_list);
     2056        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     2057        IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
    20562058        Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
    20572059        Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    2058         Input*      loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
    2059         IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
    2060         IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
     2060        Input*      loadingforcez_input=NULL;
     2061        if(dim==3){
     2062                loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     2063        }
    20612064
    20622065        /* Start  looping on the number of gaussian points: */

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

    r16847 r16862  
    514514}
    515515/*}}}*/
     516/*FUNCTION Seg::DeleteMaterials{{{*/
     517void Seg::DeleteMaterials(void){
     518        delete this->material;
     519}
     520/*}}}*/
     521/*FUNCTION Seg::GetInput(int inputenum) {{{*/
     522Input* Seg::GetInput(int inputenum){
     523        return inputs->GetInput(inputenum);
     524}
     525/*}}}*/
     526/*FUNCTION Seg::GetNumberOfNodes;{{{*/
     527int Seg::GetNumberOfNodes(void){
     528        return this->NumberofNodes();
     529}
     530/*}}}*/
     531/*FUNCTION Seg::GetVerticesCoordinates(IssmDouble** pxyz_list){{{*/
     532void Seg::GetVerticesCoordinates(IssmDouble** pxyz_list){
     533
     534        IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES*3);
     535        ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES);
     536
     537        /*Assign output pointer*/
     538        *pxyz_list = xyz_list;
     539
     540}/*}}}*/
     541/*FUNCTION Seg::JacobianDeterminant{{{*/
     542void Seg::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
     543
     544        _assert_(gauss->Enum()==GaussSegEnum);
     545        this->GetJacobianDeterminant(pJdet,xyz_list,(GaussSeg*)gauss);
     546
     547}
     548/*}}}*/
     549/*FUNCTION Seg::NewGauss(int order){{{*/
     550Gauss* Seg::NewGauss(int order){
     551        return new GaussSeg(order);
     552}
     553/*}}}*/
     554/*FUNCTION Seg::NewElementVector{{{*/
     555ElementVector* Seg::NewElementVector(int approximation_enum){
     556        return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
     557}
     558/*}}}*/
     559/*FUNCTION Seg::NewElementMatrix{{{*/
     560ElementMatrix* Seg::NewElementMatrix(int approximation_enum){
     561        return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
     562}
     563/*}}}*/
     564/*FUNCTION Seg::NodalFunctions{{{*/
     565void Seg::NodalFunctions(IssmDouble* basis, Gauss* gauss){
     566
     567        _assert_(gauss->Enum()==GaussSegEnum);
     568        this->GetNodalFunctions(basis,(GaussSeg*)gauss);
     569
     570}
     571/*}}}*/
     572/*FUNCTION Seg::NodalFunctionsDerivatives{{{*/
     573void Seg::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
     574
     575        _assert_(gauss->Enum()==GaussSegEnum);
     576        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussSeg*)gauss);
     577
     578}
     579/*}}}*/

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

    r16860 r16862  
    7777                void        ComputeStressTensor(){_error_("not implemented yet");};
    7878                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    79                 void        DeleteMaterials(void){_error_("not implemented yet");};
     79                void        DeleteMaterials(void);
    8080                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    8181                void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
    … …  
    104104                void        GetNodesSidList(int* sidlist){_error_("not implemented yet");};
    105105                void        GetNodesLidList(int* lidlist){_error_("not implemented yet");};
    106                 int         GetNumberOfNodes(void){_error_("not implemented yet");};
     106                int         GetNumberOfNodes(void);
    107107                int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
    108108                int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
    109109                int         GetNumberOfVertices(void){_error_("not implemented yet");};
    110                 void        GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");};
     110                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
    111111                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
    112112                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
    … …  
    117117                bool        IsFloating(){_error_("not implemented yet");};
    118118                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    119                 void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     119                void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    120120                void        JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    121121                void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    … …  
    123123                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    124124                IssmDouble  MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
    125                 void        NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
     125                void        NodalFunctions(IssmDouble* basis,Gauss* gauss);
    126126                void        NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    127127                void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    128                 void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     128                void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    129129                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    130130                bool        NoIceInElement(){_error_("not implemented yet");};
    … …  
    141141                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    142142                int         VelocityInterpolation(void){_error_("not implemented yet");};
    143                 Input*      GetInput(int inputenum){_error_("not implemented yet");};
     143                Input*      GetInput(int inputenum);
    144144                Input*      GetMaterialInput(int inputenum){_error_("not implemented yet");};
    145145                void        GetInputListOnVertices(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");};
    … …  
    152152                IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
    153153                Gauss*      NewGauss(void){_error_("not implemented yet");};
    154                 Gauss*      NewGauss(int order){_error_("not implemented yet");};
     154                Gauss*      NewGauss(int order);
    155155      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    156156      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
    … …  
    158158                Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
    159159                Gauss*      NewGaussTop(int order){_error_("not implemented yet");};
    160                 ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");};
    161                 ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");};
     160                ElementVector* NewElementVector(int approximation_enum);
     161                ElementMatrix* NewElementMatrix(int approximation_enum);
    162162                int         VertexConnectivity(int vertexindex){_error_("not implemented yet");};
    163163                void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};

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

    r16860 r16862  
    22492249}
    22502250/*}}}*/
     2251/*FUNCTION Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){{{*/
     2252Gauss* Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){
     2253
     2254        IssmDouble  area_coordinates[4][3];
     2255        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     2256        return new GaussTria(area_coordinates,order_vert);
     2257}
     2258/*}}}*/
    22512259/*FUNCTION Tria::NewElementVector{{{*/
    22522260ElementVector* Tria::NewElementVector(int approximation_enum){
    … …  
    26772685Element*  Tria::SpawnBasalElement(void){
    26782686
     2687        int index1,index2;
    26792688        int meshtype;
    26802689
    2681         /*go into parameters and get the analysis_counter: */
    26822690        this->parameters->FindParam(&meshtype,MeshTypeEnum);
    2683 
    2684         if(meshtype==Mesh2DhorizontalEnum){
    2685                 return this;
    2686         }
    2687         else _error_("not implemented yet");
     2691        switch(meshtype){
     2692                case Mesh2DhorizontalEnum:
     2693                        return this;
     2694                case Mesh2DverticalEnum:
     2695                        _assert_(HasEdgeOnBed());
     2696                        this->EdgeOnBedIndices(&index1,&index2);
     2697                        return SpawnSeg(index1,index2);
     2698                default:
     2699                        _error_("not implemented yet");
     2700        }
    26882701}
    26892702/*}}}*/

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

    r16860 r16862  
    296296                Gauss*         NewGauss(int order);
    297297      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
    298       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
     298      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
    299299                Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
    300300                Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
Note: See TracChangeset for help on using the changeset viewer.