Changeset 16849


Ignore:
Timestamp:
11/21/13 09:48:28 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: PVector SIA 2D

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

Legend:

Unmodified
Added
Removed

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

    r16819 r16849  
    139139}/*}}}*/
    140140ElementVector* StressbalanceSIAAnalysis::CreatePVector(Element* element){/*{{{*/
    141 _error_("not implemented yet");
     141
     142        int meshtype;
     143        element->FindParam(&meshtype,MeshTypeEnum);
     144        switch(meshtype){
     145                case Mesh2DhorizontalEnum:
     146                        return CreatePVector2D(element);
     147                case Mesh3DEnum:
     148                        return CreatePVector3D(element);
     149                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     150        }
     151}/*}}}*/
     152ElementVector* StressbalanceSIAAnalysis::CreatePVector2D(Element* element){/*{{{*/
     153
     154        /*Intermediaries */
     155        IssmDouble ub,vb,slope2,drag,thickness,connectivity;
     156        IssmDouble slope[2];
     157
     158        /*Fetch number vertices for this element*/
     159        int numvertices = element->GetNumberOfVertices();
     160
     161        /*Initialize Element vector*/
     162        ElementVector* pe=element->NewElementVector();
     163
     164        /*Retrieve all inputs and parameters*/
     165        IssmDouble  rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
     166        IssmDouble  gravity    = element->GetMaterialParameter(ConstantsGEnum);
     167        IssmDouble  n          = element->GetMaterialParameter(MaterialsRheologyNEnum);
     168        IssmDouble  B          = element->GetMaterialParameter(MaterialsRheologyBbarEnum);
     169        Input* slopex_input    = element->GetInput(SurfaceSlopeXEnum);        _assert_(slopex_input);
     170        Input* slopey_input    = element->GetInput(SurfaceSlopeYEnum);        _assert_(slopey_input);
     171        Input* thickness_input = element->GetInput(ThicknessEnum);            _assert_(thickness_input);
     172        Input* drag_input      = element->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
     173
     174        Gauss* gauss=element->NewGauss();
     175        for(int iv=0;iv<numvertices;iv++){
     176                gauss->GaussVertex(iv);
     177
     178                connectivity=(IssmDouble)element->VertexConnectivity(iv);
     179
     180                thickness_input->GetInputValue(&thickness,gauss);
     181                drag_input->GetInputValue(&drag,gauss);
     182                slopex_input->GetInputValue(&slope[0],gauss);
     183                slopey_input->GetInputValue(&slope[1],gauss);
     184                slope2=slope[0]*slope[0]+slope[1]*slope[1];
     185
     186                /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
     187                ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
     188                vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
     189                ///*Ritz et al. 1996*/
     190                //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
     191                //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
     192                ///*Rutt et al. 2009*/
     193                //ub=-drag*rho_ice*gravity*thickness*slope[0];
     194                //vb=-drag*rho_ice*gravity*thickness*slope[1];
     195
     196                pe->values[2*iv+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
     197                pe->values[2*iv+1]=(vb-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/connectivity;
     198        }
     199
     200        /*Clean up and return*/
     201        delete gauss;
     202        return pe;
     203}/*}}}*/
     204ElementVector* StressbalanceSIAAnalysis::CreatePVector3D(Element* element){/*{{{*/
     205
     206        _error_("not implemented");
    142207}/*}}}*/
    143208void StressbalanceSIAAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h

    r16782 r16849  
    2323                ElementMatrix* CreateKMatrix(Element* element);
    2424                ElementVector* CreatePVector(Element* element);
     25                ElementVector* CreatePVector2D(Element* element);
     26                ElementVector* CreatePVector3D(Element* element);
    2527                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2628                void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

    r16848 r16849  
    183183                virtual IssmDouble TotalSmb(void)=0;
    184184                virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0;
     185                virtual int    VertexConnectivity(int vertexindex)=0;
    185186                #endif
    186187

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

    r16848 r16849  
    10221022
    10231023        _assert_(this->matpar);
    1024        
    10251024        switch(enum_in){ // FIXME: change this to material
    1026         case MaterialsRheologyNEnum:
    1027                 return this->material->GetN();
    1028         case MaterialsRheologyBEnum:
    1029                 return this->material->GetB();
    1030         default:
    1031                 return this->matpar->GetMaterialParameter(enum_in);
     1025                case MaterialsRheologyNEnum:
     1026                        return this->material->GetN();
     1027                case MaterialsRheologyBEnum:
     1028                        return this->material->GetB();
     1029                default:
     1030                        return this->matpar->GetMaterialParameter(enum_in);
    10321031        }
    10331032}
    … …  
    34613460                        break;
    34623461        }
     3462}
     3463/*}}}*/
     3464/*FUNCTION Penta::VertexConnectivity{{{*/
     3465int Penta::VertexConnectivity(int vertexindex){
     3466        _assert_(this->vertices);
     3467        return this->vertices[vertexindex]->Connectivity();
    34633468}
    34643469/*}}}*/

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

    r16848 r16849  
    130130                int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
    131131                IssmDouble TimeAdapt();
     132                int    VertexConnectivity(int vertexindex);
    132133                void   ViscousHeatingCreateInput(void);
    133134                void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);

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

    r16848 r16849  
    154154                ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");};
    155155                ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");};
     156                int         VertexConnectivity(int vertexindex){_error_("not implemented yet");};
    156157                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
    157158                void        ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");};

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

    r16848 r16849  
    11601160
    11611161        _assert_(this->matpar);
    1162         return this->matpar->GetMaterialParameter(enum_in);
     1162        switch(enum_in){ // FIXME: change this to material
     1163                case MaterialsRheologyNEnum:
     1164                        return this->material->GetN();
     1165                case MaterialsRheologyBEnum:
     1166                        return this->material->GetB();
     1167                case MaterialsRheologyBbarEnum:
     1168                        return this->material->GetBbar();
     1169                default:
     1170                        return this->matpar->GetMaterialParameter(enum_in);
     1171        }
    11631172}
    11641173/*}}}*/
    … …  
    29372946        delete gauss;
    29382947
     2948}
     2949/*}}}*/
     2950/*FUNCTION Tria::VertexConnectivity{{{*/
     2951int Tria::VertexConnectivity(int vertexindex){
     2952        _assert_(this->vertices);
     2953        return this->vertices[vertexindex]->Connectivity();
    29392954}
    29402955/*}}}*/

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

    r16848 r16849  
    137137                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    138138                IssmDouble  TimeAdapt();
     139                int         VertexConnectivity(int vertexindex);
    139140                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    140141                bool        IsZeroLevelset(int levelset_enum);
Note: See TracChangeset for help on using the changeset viewer.