Changeset 16872


Ignore:
Timestamp:
11/21/13 16:46:58 (11 years ago)
Author:
seroussi
Message:

CHG: added SSA friction KMatrix in analysis

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

Legend:

Unmodified
Added
Removed

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

    r16866 r16872  
    974974        return Ke;
    975975}/*}}}*/
     976ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
     977        /*Intermediaries*/
     978        bool       mainlyfloating;
     979        int        analysis_type,migration_style,point1;
     980        IssmDouble alpha2,Jdet,fraction1,fraction2;
     981        IssmDouble phi=1.0;
     982        IssmDouble D_scalar,gllevelset;
     983        IssmDouble *xyz_list = NULL;
     984        Friction   *friction = NULL;
     985        Gauss*     gauss     = NULL;
     986
     987        /*Return if element is inactive*/
     988        if(element->IsFloating()) return NULL;
     989
     990        /*Fetch number of nodes and dof for this finite element*/
     991        int numnodes = element->GetNumberOfNodes();
     992        int numdof   = numnodes*2;
     993
     994        /*Initialize Element matrix and vectors*/
     995        ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
     996        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
     997        IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
     998
     999        /*Retrieve all inputs and parameters*/
     1000        element->GetVerticesCoordinates(&xyz_list);
     1001        element->FindParam(&migration_style,GroundinglineMigrationEnum);
     1002        Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
     1003        Input* vx_input=element->GetInput(VxEnum);           _assert_(vx_input);
     1004        Input* vy_input=element->GetInput(VyEnum);           _assert_(vy_input);
     1005        Input* gllevelset_input=NULL;
     1006        element->FindParam(&analysis_type,AnalysisTypeEnum);
     1007
     1008        /*build friction object, used later on: */
     1009        //friction=new Friction("2d",inputs,matpar,analysis_type);
     1010
     1011        /*Recover portion of element that is grounded*/
     1012        if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list);
     1013        if(migration_style==SubelementMigration2Enum){
     1014                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     1015                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     1016           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
     1017        }
     1018        else{
     1019                gauss = element->NewGauss(2);
     1020        }
     1021
     1022        /* Start  looping on the number of gaussian points: */
     1023        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1024
     1025                gauss->GaussPoint(ig);
     1026
     1027                //friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     1028                alpha2=1.0;
     1029                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
     1030                if(migration_style==SubelementMigration2Enum){
     1031                        gllevelset_input->GetInputValue(&gllevelset, gauss);
     1032                        if(gllevelset<0) alpha2=0;
     1033                }
     1034
     1035                this->GetBSSAFriction(B, element, xyz_list, gauss);
     1036                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     1037                D_scalar=alpha2*gauss->weight*Jdet;
     1038                for(int i=0;i<2;i++) D[i*2+i]=D_scalar;
     1039
     1040                TripleMultiply(B,2,numdof,1,
     1041                                        D,2,2,0,
     1042                                        B,2,numdof,0,
     1043                                        &Ke->values[0],1);
     1044        }
     1045
     1046        /*Transform Coordinate System*/
     1047        element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1048
     1049        /*Clean up and return*/
     1050        delete gauss;
     1051        delete friction;
     1052        xDelete<IssmDouble>(xyz_list);
     1053        xDelete<IssmDouble>(D);
     1054        xDelete<IssmDouble>(B);
     1055        return Ke;
     1056}/*}}}*/
    9761057ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
    9771058
    … …  
    10341115        xDelete<IssmDouble>(B);
    10351116        return Ke;
    1036 }/*}}}*/
    1037 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
    1038         return NULL;
    10391117}/*}}}*/
    10401118ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
    … …  
    12071285        /*Clean-up*/
    12081286        xDelete<IssmDouble>(dbasis);
     1287}/*}}}*/
     1288void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1289        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
     1290         * For node i, Bi can be expressed in the actual coordinate system
     1291         * by:
     1292         *                 Bi=[ N   0 ]
     1293         *                    [ 0   N ]
     1294         * where N is the finiteelement function for node i.
     1295         *
     1296         * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
     1297         */
     1298
     1299        /*Fetch number of nodes for this finite element*/
     1300        int numnodes = element->GetNumberOfNodes();
     1301
     1302        /*Get nodal functions derivatives*/
     1303        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     1304        element->NodalFunctions(basis,gauss);
     1305
     1306        /*Build L: */
     1307        for(int i=0;i<numnodes;i++){
     1308                B[2*numnodes*0+2*i+0] = basis[i];
     1309                B[2*numnodes*0+2*i+1] = 0.;
     1310                B[2*numnodes*1+2*i+0] = 0.;
     1311                B[2*numnodes*1+2*i+1] = basis[i];
     1312        }
     1313
     1314        /*Clean-up*/
     1315        xDelete<IssmDouble>(basis);
    12091316}/*}}}*/
    12101317void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

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

    r16859 r16872  
    3535                ElementVector* CreatePVectorSSAFront(Element* element);
    3636                void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     37                void GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3738                void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3839                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);

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

    r16866 r16872  
    106106                virtual bool   IsOnBed()=0;
    107107                virtual bool   IsOnSurface()=0;
     108                virtual void   GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
     109                virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
    108110                virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype)=0;
    109111                virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
    … …  
    144146      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0;
    145147      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert)=0;
     148      virtual Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order)=0;
    146149                virtual Gauss* NewGaussBase(int order)=0;
    147150                virtual Gauss* NewGaussLine(int vertex1,int vertex2,int order)=0;

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

    r16866 r16872  
    9292                void     GetDofListVelocity(int** pdoflist,int setenum);
    9393                void     GetDofListPressure(int** pdoflist,int setenum);
     94                void   GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
     95                IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
    9496                int    GetNodeIndex(Node* node);
    9597                void   GetNodesSidList(int* sidlist);
    … …  
    228230                void           GetVertexSidList(int* sidlist);
    229231                void           GetConnectivityList(int* connectivity);
    230                 void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    231                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    232232                int            GetElementType(void);
    233233                Input*         GetInput(int inputenum);
    … …  
    266266      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    267267      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
     268      Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
    268269                Gauss*         NewGaussBase(int order);
    269270                Gauss*         NewGaussLine(int vertex1,int vertex2,int order);

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

    r16866 r16872  
    142142                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    143143                int         VelocityInterpolation(void){_error_("not implemented yet");};
     144                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
     145                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
    144146                Input*      GetInput(int inputenum);
    145147                Input*      GetMaterialInput(int inputenum){_error_("not implemented yet");};
    … …  
    157159      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    158160      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
     161      Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
    159162                Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
    160163                Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};

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

    r16868 r16872  
    22982298
    22992299        return new GaussTria(area_coordinates,order);
     2300}
     2301/*}}}*/
     2302/*FUNCTION Tria::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating){{{*/
     2303Gauss* Tria::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){
     2304
     2305        return new GaussTria(point1,fraction1,fraction2,mainlyfloating,order);
    23002306}
    23012307/*}}}*/

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

    r16866 r16872  
    9191                void          GetDofListVelocity(int** pdoflist,int setenum);
    9292                void          GetDofListPressure(int** pdoflist,int setenum);
     93                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
     94                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
    9395                int         GetNodeIndex(Node* node);
    9496                void        GetNodesSidList(int* sidlist);
    … …  
    264266                void           GetVertexSidList(int* sidlist);
    265267                void           GetConnectivityList(int* connectivity);
    266                 void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    267                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    268268                IssmDouble     GetYcoord(Gauss* gauss);
    269269                IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
    … …  
    298298                Gauss*         NewGauss(int order);
    299299      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
     300      Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
    300301      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
    301302                Gauss*         NewGaussBase(int order);
Note: See TracChangeset for help on using the changeset viewer.