Changeset 16850


Ignore:
Timestamp:
11/21/13 10:34:13 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: PVector SIA 3d

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

Legend:

Unmodified
Added
Removed

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

    r16845 r16850  
    816816        element->GetInputValue(&approximation,ApproximationEnum);
    817817        switch(approximation){
     818                case SIAApproximationEnum:
     819                        return NULL;
    818820                case SSAApproximationEnum:
    819821                        return CreateKMatrixSSA(element);
    … …  
    833835        element->GetInputValue(&approximation,ApproximationEnum);
    834836        switch(approximation){
     837                case SIAApproximationEnum:
     838                        return NULL;
    835839                case SSAApproximationEnum:
    836840                        return CreatePVectorSSA(element);

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

    r16849 r16850  
    204204ElementVector* StressbalanceSIAAnalysis::CreatePVector3D(Element* element){/*{{{*/
    205205
    206         _error_("not implemented");
     206        /*Intermediaries */
     207        int         nodeup,nodedown,numsegments;
     208        IssmDouble  ub,vb,slope2,drag,surface,thickness,constant_part,z,Jdet;
     209        IssmDouble  slope[2],connectivity[2],xyz_list_line[2][3];
     210        IssmDouble *xyz_list = NULL;
     211        int        *pairindices = NULL;
     212
     213        /*Fetch number vertices for this element*/
     214        int numvertices = element->GetNumberOfVertices();
     215
     216        /*Initialize Element vector*/
     217        ElementVector* pe=element->NewElementVector();
     218
     219        /*Retrieve all inputs and parameters*/
     220        element->GetVerticesCoordinates(&xyz_list);
     221        IssmDouble  rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
     222        IssmDouble  gravity    = element->GetMaterialParameter(ConstantsGEnum);
     223        IssmDouble  n          = element->GetMaterialParameter(MaterialsRheologyNEnum);
     224        IssmDouble  B          = element->GetMaterialParameter(MaterialsRheologyBEnum);
     225        Input* surface_input   = element->GetInput(SurfaceEnum);              _assert_(surface_input);
     226        Input* slopex_input    = element->GetInput(SurfaceSlopeXEnum);        _assert_(slopex_input);
     227        Input* slopey_input    = element->GetInput(SurfaceSlopeYEnum);        _assert_(slopey_input);
     228        Input* thickness_input = element->GetInput(ThicknessEnum);            _assert_(thickness_input);
     229        Input* drag_input      = element->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
     230
     231        /*Get Vertical segment indices*/
     232        element->VerticalSegmentIndices(&pairindices,&numsegments);
     233        for(int is=0;is<numsegments;is++){
     234                nodedown = pairindices[is*2+0];
     235                nodeup   = pairindices[is*2+1];
     236                connectivity[0]=(IssmDouble)element->VertexConnectivity(nodedown);
     237                connectivity[1]=(IssmDouble)element->VertexConnectivity(nodeup);
     238                for(int i=0;i<3;i++){
     239                        xyz_list_line[0][i]=xyz_list[nodedown*3+i];
     240                        xyz_list_line[1][i]=xyz_list[nodeup*3+i];
     241                }
     242
     243                Gauss* gauss=element->NewGaussLine(nodedown,nodeup,3);
     244                for(int ig=gauss->begin();ig<gauss->end();ig++){
     245                        gauss->GaussPoint(ig);
     246
     247                        slopex_input->GetInputValue(&slope[0],gauss);
     248                        slopey_input->GetInputValue(&slope[1],gauss);
     249                        surface_input->GetInputValue(&surface,gauss);
     250                        thickness_input->GetInputValue(&thickness,gauss);
     251
     252                        slope2=slope[0]*slope[0]+slope[1]*slope[1];
     253                        constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
     254
     255                        z = element->GetZcoord(gauss);
     256                        element->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss);
     257
     258                        if(element->IsOnSurface()){
     259                                pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1];
     260                                pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1];
     261                        }
     262                        else{/*connectivity is too large, should take only half on it*/
     263                                pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1];
     264                                pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1];
     265                        }
     266                }
     267
     268                /*Deal with basal velocities*/
     269                if(element->IsOnBed()){
     270                        drag_input->GetInputValue(&drag,gauss);
     271
     272                        /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
     273                        ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
     274                        vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
     275                        ///*Ritz et al. 1996*/
     276                        //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
     277                        //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
     278                        ///*Rutt et al. 2009*/
     279                        //ub=-drag*rho_ice*gravity*thickness*slope[0];
     280                        //vb=-drag*rho_ice*gravity*thickness*slope[1];
     281
     282                        pe->values[2*nodedown+0]+=ub/connectivity[0];
     283                        pe->values[2*nodedown+1]+=vb/connectivity[0];
     284                }
     285                delete gauss;
     286        }
     287
     288        /*Clean up and return*/
     289        xDelete<int>(pairindices);
     290        xDelete<IssmDouble>(xyz_list);
     291        return pe;
     292
    207293}/*}}}*/
    208294void StressbalanceSIAAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16849 r16850  
    8585                virtual void    GetDofListPressure(int** pdoflist,int setenum)=0;
    8686                virtual void   JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
     87                virtual void   JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
    8788                virtual void   JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
    8889                virtual void   JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
    … …  
    140141      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert)=0;
    141142                virtual Gauss* NewGaussBase(int order)=0;
     143                virtual Gauss* NewGaussLine(int vertex1,int vertex2,int order)=0;
    142144                virtual Gauss* NewGaussTop(int order)=0;
    143145                virtual ElementVector*  NewElementVector(int approximation_enum=NoneApproximationEnum)=0;
    … …  
    184186                virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0;
    185187                virtual int    VertexConnectivity(int vertexindex)=0;
     188                virtual void   VerticalSegmentIndices(int** pindices,int* pnumseg)=0;
    186189                #endif
    187190

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

    r16849 r16850  
    24672467}
    24682468/*}}}*/
     2469/*FUNCTION Penta::JacobianDeterminantLine{{{*/
     2470void Penta::JacobianDeterminantLine(IssmDouble* pJdet,IssmDouble* xyz_list_line,Gauss* gauss){
     2471
     2472        _assert_(gauss->Enum()==GaussPentaEnum);
     2473        this->GetSegmentJacobianDeterminant(pJdet,xyz_list_line,(GaussPenta*)gauss);
     2474
     2475}
     2476/*}}}*/
    24692477/*FUNCTION Penta::JacobianDeterminantTop{{{*/
    24702478void Penta::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){
    … …  
    25742582Gauss* Penta::NewGaussBase(int order){
    25752583        return new GaussPenta(0,1,2,order);
     2584}
     2585/*}}}*/
     2586/*FUNCTION Penta::NewGaussLine(int vertex1,int vertex2,int order){{{*/
     2587Gauss* Penta::NewGaussLine(int vertex1,int vertex2,int order){
     2588        return new GaussPenta(vertex1,vertex2,order);
    25762589}
    25772590/*}}}*/
    … …  
    34663479        _assert_(this->vertices);
    34673480        return this->vertices[vertexindex]->Connectivity();
     3481}
     3482/*}}}*/
     3483/*FUNCTION Penta::VerticalSegmentIndices{{{*/
     3484void Penta::VerticalSegmentIndices(int** pindices,int* pnumseg){
     3485
     3486        /*output*/
     3487        int *indices = xNew<int>(3*2);
     3488        indices[0*2 + 0] = 0; indices[0*2 + 1] = 3;
     3489        indices[1*2 + 0] = 1; indices[1*2 + 1] = 4;
     3490        indices[2*2 + 0] = 2; indices[2*2 + 1] = 5;
     3491
     3492        /*Assign output pointers*/
     3493        *pindices = indices;
     3494        *pnumseg  = 3;
    34683495}
    34693496/*}}}*/

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

    r16849 r16850  
    131131                IssmDouble TimeAdapt();
    132132                int    VertexConnectivity(int vertexindex);
     133                void   VerticalSegmentIndices(int** pindices,int* pnumseg);
    133134                void   ViscousHeatingCreateInput(void);
    134135                void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    … …  
    250251                bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
    251252                void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
     253                void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
    252254                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    253255                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    … …  
    259261      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
    260262                Gauss*         NewGaussBase(int order);
     263                Gauss*         NewGaussLine(int vertex1,int vertex2,int order);
    261264                Gauss*         NewGaussTop(int order);
    262265                ElementVector* NewElementVector(int approximation_enum);

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

    r16849 r16850  
    116116                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    117117                void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    118                 void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     118                void        JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     119                void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    119120                void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    120121                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    … …  
    151152      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
    152153                Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
     154                Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
    153155                Gauss*      NewGaussTop(int order){_error_("not implemented yet");};
    154156                ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");};
    155157                ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");};
    156158                int         VertexConnectivity(int vertexindex){_error_("not implemented yet");};
     159                void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
    157160                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
    158161                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.h

    r16849 r16850  
    138138                IssmDouble  TimeAdapt();
    139139                int         VertexConnectivity(int vertexindex);
     140                void   VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
    140141                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    141142                bool        IsZeroLevelset(int levelset_enum);
    … …  
    283284                bool             IsInput(int name);
    284285                void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
     286                void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    285287                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    286288                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    … …  
    292294      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
    293295                Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
     296                Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
    294297                Gauss*         NewGaussTop(int order){_error_("not implemented yet");};
    295298                ElementVector* NewElementVector(int approximation_enum);
Note: See TracChangeset for help on using the changeset viewer.