Changeset 16860


Ignore:
Timestamp:
11/21/13 14:05:56 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: bal vel PVector

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

Legend:

Unmodified
Added
Removed

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

    r16782 r16860  
    6363}/*}}}*/
    6464ElementVector* BalancevelocityAnalysis::CreatePVector(Element* element){/*{{{*/
    65 _error_("not implemented yet");
     65
     66        /*Intermediaries*/
     67        int      meshtype;
     68        Element* basalelement;
     69
     70        /*Get basal element*/
     71        element->FindParam(&meshtype,MeshTypeEnum);
     72        switch(meshtype){
     73                case Mesh2DhorizontalEnum:
     74                        basalelement = element;
     75                        break;
     76                case Mesh3DEnum:
     77                        if(!element->IsOnBed()) return NULL;
     78                        basalelement = element->SpawnBasalElement();
     79                        break;
     80                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     81        }
     82
     83        /*Intermediaries */
     84        IssmDouble dhdt,mb,ms,Jdet;
     85        IssmDouble gamma,thickness;
     86        IssmDouble hnx,hny,dhnx[2],dhny[2];
     87        IssmDouble *xyz_list  = NULL;
     88
     89        /*Fetch number of nodes and dof for this finite element*/
     90        int numnodes = basalelement->GetNumberOfNodes();
     91
     92        /*Initialize Element vector*/
     93        ElementVector* pe    = basalelement->NewElementVector();
     94        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     95        IssmDouble*    dbasis = xNew<IssmDouble>(numnodes*2);
     96        IssmDouble*    H      = xNew<IssmDouble>(numnodes);
     97        IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
     98        IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
     99
     100        /*Retrieve all inputs and parameters*/
     101        basalelement->GetVerticesCoordinates(&xyz_list);
     102        Input* ms_input   = basalelement->GetInput(SurfaceforcingsMassBalanceEnum);     _assert_(ms_input);
     103        Input* mb_input   = basalelement->GetInput(BasalforcingsMeltingRateEnum);       _assert_(mb_input);
     104        Input* dhdt_input = basalelement->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
     105        Input* H_input    = basalelement->GetInput(ThicknessEnum);                      _assert_(H_input);
     106        IssmDouble h = basalelement->CharacteristicLength();
     107
     108        /*Get vector N for all nodes*/
     109        basalelement->GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
     110        basalelement->GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
     111        basalelement->GetInputListOnNodes(H,ThicknessEnum);
     112        for(int i=0;i<numnodes;i++){
     113                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
     114                Nx[i] = -H[i]*Nx[i]/norm;
     115                Ny[i] = -H[i]*Ny[i]/norm;
     116        }
     117
     118
     119        /* Start  looping on the number of gaussian points: */
     120        Gauss* gauss=basalelement->NewGauss(2);
     121        for(int ig=gauss->begin();ig<gauss->end();ig++){
     122                gauss->GaussPoint(ig);
     123
     124                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     125                basalelement->NodalFunctions(basis,gauss);
     126                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     127
     128                element->ValueP1DerivativesOnGauss(&dhnx[0],Nx,xyz_list,gauss);
     129                element->ValueP1DerivativesOnGauss(&dhny[0],Ny,xyz_list,gauss);
     130                element->ValueP1OnGauss(&hnx,Nx,gauss);
     131                element->ValueP1OnGauss(&hny,Ny,gauss);
     132
     133                ms_input->GetInputValue(&ms,gauss);
     134                mb_input->GetInputValue(&mb,gauss);
     135                dhdt_input->GetInputValue(&dhdt,gauss);
     136                H_input->GetInputValue(&thickness,gauss);
     137                if(thickness<50.) thickness=50.;
     138
     139                gamma=h/(2.*thickness+1.e-10);
     140
     141                for(int i=0;i<numnodes;i++){
     142                        pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*( basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[1])+hnx*dbasis[0*numnodes+i] + hny*dbasis[1*numnodes+i]));
     143                }
     144        }
     145
     146        /*Clean up and return*/
     147        xDelete<IssmDouble>(xyz_list);
     148        xDelete<IssmDouble>(basis);
     149        xDelete<IssmDouble>(dbasis);
     150        xDelete<IssmDouble>(H);
     151        xDelete<IssmDouble>(Nx);
     152        xDelete<IssmDouble>(Ny);
     153        delete gauss;
     154        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     155        return pe;
    66156}/*}}}*/
    67157void BalancevelocityAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16782 r16860  
    4949}/*}}}*/
    5050ElementVector* SmoothedSurfaceSlopeXAnalysis::CreatePVector(Element* element){/*{{{*/
    51 _error_("not implemented yet");
     51
     52        /*Intermediaries*/
     53        int      meshtype;
     54        Element* basalelement;
     55
     56        /*Get basal element*/
     57        element->FindParam(&meshtype,MeshTypeEnum);
     58        switch(meshtype){
     59                case Mesh2DhorizontalEnum:
     60                        basalelement = element;
     61                        break;
     62                case Mesh3DEnum:
     63                        if(!element->IsOnBed()) return NULL;
     64                        basalelement = element->SpawnBasalElement();
     65                        break;
     66                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     67        }
     68
     69        /*Intermediaries */
     70        int         input_enum;
     71        IssmDouble  Jdet,thickness,slope[2];
     72        IssmDouble  taud_x,norms,normv,vx,vy;
     73        IssmDouble *xyz_list  = NULL;
     74
     75        /*Fetch number of nodes and dof for this finite element*/
     76        int numnodes = basalelement->GetNumberOfNodes();
     77
     78        /*Initialize Element vector*/
     79        ElementVector* pe    = basalelement->NewElementVector();
     80        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     81
     82        /*Retrieve all inputs and parameters*/
     83        basalelement->GetVerticesCoordinates(&xyz_list);
     84        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     85        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
     86        Input* H_input       = basalelement->GetInput(ThicknessEnum); _assert_(H_input);
     87        Input* surface_input = basalelement->GetInput(SurfaceEnum);   _assert_(surface_input);
     88        Input* vx_input      = basalelement->GetInput(VxEnum);
     89        Input* vy_input      = basalelement->GetInput(VyEnum);
     90
     91        /* Start  looping on the number of gaussian points: */
     92        Gauss* gauss=basalelement->NewGauss(2);
     93        for(int ig=gauss->begin();ig<gauss->end();ig++){
     94                gauss->GaussPoint(ig);
     95
     96                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     97                basalelement->NodalFunctions(basis,gauss);
     98
     99                H_input->GetInputValue(&thickness,gauss);
     100                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     101                if(vx_input && vy_input){
     102                        vx_input->GetInputValue(&vx,gauss);
     103                        vy_input->GetInputValue(&vy,gauss);
     104                        norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
     105                        normv = sqrt(vx*vx + vy*vy);
     106                        if(normv>15./(365.*24.*3600.)) slope[0] = -vx/normv*norms;
     107                }
     108                taud_x = rho_ice*gravity*thickness*slope[0];
     109
     110                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_x*basis[i];
     111        }
     112
     113        /*Clean up and return*/
     114        xDelete<IssmDouble>(xyz_list);
     115        xDelete<IssmDouble>(basis);
     116        delete gauss;
     117        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     118        return pe;
    52119}/*}}}*/
    53120void SmoothedSurfaceSlopeXAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16782 r16860  
    4949}/*}}}*/
    5050ElementVector* SmoothedSurfaceSlopeYAnalysis::CreatePVector(Element* element){/*{{{*/
    51 _error_("not implemented yet");
     51
     52        /*Intermediaries*/
     53        int      meshtype;
     54        Element* basalelement;
     55
     56        /*Get basal element*/
     57        element->FindParam(&meshtype,MeshTypeEnum);
     58        switch(meshtype){
     59                case Mesh2DhorizontalEnum:
     60                        basalelement = element;
     61                        break;
     62                case Mesh3DEnum:
     63                        if(!element->IsOnBed()) return NULL;
     64                        basalelement = element->SpawnBasalElement();
     65                        break;
     66                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     67        }
     68
     69        /*Intermediaries */
     70        int         input_enum;
     71        IssmDouble  Jdet,thickness,slope[2];
     72        IssmDouble  taud_y,norms,normv,vx,vy;
     73        IssmDouble *xyz_list  = NULL;
     74
     75        /*Fetch number of nodes and dof for this finite element*/
     76        int numnodes = basalelement->GetNumberOfNodes();
     77
     78        /*Initialize Element vector*/
     79        ElementVector* pe    = basalelement->NewElementVector();
     80        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     81
     82        /*Retrieve all inputs and parameters*/
     83        basalelement->GetVerticesCoordinates(&xyz_list);
     84        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     85        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
     86        Input* H_input       = basalelement->GetInput(ThicknessEnum); _assert_(H_input);
     87        Input* surface_input = basalelement->GetInput(SurfaceEnum);   _assert_(surface_input);
     88        Input* vx_input      = basalelement->GetInput(VxEnum);
     89        Input* vy_input      = basalelement->GetInput(VyEnum);
     90
     91        /* Start  looping on the number of gaussian points: */
     92        Gauss* gauss=basalelement->NewGauss(2);
     93        for(int ig=gauss->begin();ig<gauss->end();ig++){
     94                gauss->GaussPoint(ig);
     95
     96                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     97                basalelement->NodalFunctions(basis,gauss);
     98
     99                H_input->GetInputValue(&thickness,gauss);
     100                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     101                if(vx_input && vy_input){
     102                        vx_input->GetInputValue(&vx,gauss);
     103                        vy_input->GetInputValue(&vy,gauss);
     104                        norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
     105                        normv = sqrt(vx*vx + vy*vy);
     106                        if(normv>15./(365.*24.*3600.)) slope[1] = -vy/normv*norms;
     107                }
     108                taud_y = rho_ice*gravity*thickness*slope[1];
     109
     110                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_y*basis[i];
     111        }
     112
     113        /*Clean up and return*/
     114        xDelete<IssmDouble>(xyz_list);
     115        xDelete<IssmDouble>(basis);
     116        delete gauss;
     117        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     118        return pe;
    52119}/*}}}*/
    53120void SmoothedSurfaceSlopeYAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16853 r16860  
    158158                virtual void   ResetCoordinateSystem()=0;
    159159                virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
     160                virtual void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss)=0;
     161                virtual void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss)=0;
    160162                virtual void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
    161163                virtual void   ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;

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

    r16854 r16860  
    34733473                        break;
    34743474        }
     3475}
     3476/*}}}*/
     3477/*FUNCTION Penta::ValueP1OnGauss{{{*/
     3478void Penta::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){
     3479        PentaRef::GetInputValue(pvalue,values,gauss);
     3480}
     3481/*}}}*/
     3482/*FUNCTION Penta::ValueP1DerivativesOnGauss{{{*/
     3483void Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){
     3484        PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
    34753485}
    34763486/*}}}*/

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

    r16853 r16860  
    132132                int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
    133133                IssmDouble TimeAdapt();
     134                void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
     135                void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
    134136                int    VertexConnectivity(int vertexindex);
    135137                void   VerticalSegmentIndices(int** pindices,int* pnumseg);

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

    r16853 r16860  
    138138                IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
    139139                int         PressureInterpolation(void){_error_("not implemented yet");};
     140                void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
     141                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    140142                int         VelocityInterpolation(void){_error_("not implemented yet");};
    141143                Input*      GetInput(int inputenum){_error_("not implemented yet");};

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

    r16857 r16860  
    29462946        delete gauss;
    29472947
     2948}
     2949/*}}}*/
     2950/*FUNCTION Tria::ValueP1OnGauss{{{*/
     2951void Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){
     2952        TriaRef::GetInputValue(pvalue,values,gauss);
     2953}
     2954/*}}}*/
     2955/*FUNCTION Tria::ValueP1DerivativesOnGauss{{{*/
     2956void Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){
     2957        TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
    29482958}
    29492959/*}}}*/

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

    r16853 r16860  
    139139                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    140140                IssmDouble  TimeAdapt();
     141                void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
     142                void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
    141143                int         VertexConnectivity(int vertexindex);
    142144                void   VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
Note: See TracChangeset for help on using the changeset viewer.