Changeset 18293


Ignore:
Timestamp:
07/28/14 09:18:47 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added lateral friction for Erkut SSA

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

Legend:

Unmodified
Added
Removed

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

    r18286 r18293  
    88
    99//#define FSANALYTICAL 12
     10//#define LATERALFRICTION 1
    1011
    1112/*Model processing*/
    … …  
    223224        iomodel->FetchDataToInput(elements,LoadingforceXEnum);
    224225        iomodel->FetchDataToInput(elements,LoadingforceYEnum);
     226        #ifdef LATERALFRICTION
     227        iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
     228        #endif
    225229        if(isdamage)iomodel->FetchDataToInput(elements,DamageDEnum);
    226230
    … …  
    12181222        ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement);
    12191223        ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement);
     1224        #ifdef LATERALFRICTION
     1225        ElementMatrix* Ke3=CreateKMatrixSSALateralFriction(basalelement);
     1226        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     1227        delete Ke3;
     1228        #else
    12201229        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     1230        #endif
    12211231
    12221232        /*clean-up and return*/
    … …  
    13061316        delete friction;
    13071317        xDelete<IssmDouble>(xyz_list);
     1318        xDelete<IssmDouble>(B);
     1319        xDelete<IssmDouble>(D);
     1320        return Ke;
     1321}/*}}}*/
     1322ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSALateralFriction(Element* element){/*{{{*/
     1323
     1324        /*Return if element is inactive*/
     1325        if(!element->IsIceInElement()) return NULL;
     1326
     1327        /*If no boundary, return NULL*/
     1328        if(!element->IsFaceOnBoundary()) return NULL;
     1329
     1330        /*Intermediaries*/
     1331        IssmDouble  alpha2;
     1332        IssmDouble  Jdet;
     1333        int         domaintype;
     1334        IssmDouble  icefront;
     1335        IssmDouble *xyz_list          = NULL;
     1336        IssmDouble *xyz_list_boundary = NULL;
     1337
     1338        /*Get problem dimension*/
     1339        element->FindParam(&domaintype,DomainTypeEnum);
     1340        if(domaintype==Domain2DverticalEnum) return NULL; //not supported yet
     1341
     1342        /*Fetch number of nodes and dof for this finite element*/
     1343        int dim      = 2;
     1344        int numnodes = element->GetNumberOfNodes();
     1345        int numdof   = numnodes*dim;
     1346
     1347        /*Initialize Element matrix and vectors*/
     1348        ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum);
     1349        IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
     1350        IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
     1351
     1352        /*Retrieve all inputs and parameters*/
     1353        element->GetVerticesCoordinates(&xyz_list);
     1354        element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.);
     1355        Input* icelevelset_input = element->GetInput(MaskIceLevelsetEnum); _assert_(icelevelset_input);
     1356
     1357        /* Start  looping on the number of gaussian points: */
     1358        Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3);
     1359        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1360                gauss->GaussPoint(ig);
     1361
     1362                this->GetBSSAFriction(B,element,dim,xyz_list,gauss);
     1363                element->JacobianDeterminantSurface(&Jdet,xyz_list_boundary,gauss);
     1364                icelevelset_input->GetInputValue(&icefront, gauss);
     1365                if(icefront==0.)
     1366                 alpha2=0.;
     1367                else
     1368                 alpha2=2.e+12;
     1369                for(int i=0;i<dim;i++) D[i*dim+i]=alpha2*gauss->weight*Jdet;
     1370
     1371                TripleMultiply(B,dim,numdof,1,
     1372                                        D,dim,dim,0,
     1373                                        B,dim,numdof,0,
     1374                                        &Ke->values[0],1);
     1375        }
     1376
     1377        /*Transform Coordinate System*/
     1378        if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1379
     1380        /*Clean up and return*/
     1381        delete gauss;
     1382        xDelete<IssmDouble>(xyz_list);
     1383        xDelete<IssmDouble>(xyz_list_boundary);
    13081384        xDelete<IssmDouble>(B);
    13091385        xDelete<IssmDouble>(D);

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

    r18284 r18293  
    3737                ElementMatrix* CreateKMatrixSSAViscous(Element* element);
    3838                ElementMatrix* CreateKMatrixSSAFriction(Element* element);
     39                ElementMatrix* CreateKMatrixSSALateralFriction(Element* element);
    3940                ElementVector* CreatePVectorSSA(Element* element);
    4041                ElementVector* CreatePVectorSSADrivingStress(Element* element);

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

    r18283 r18293  
    256256                virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
    257257                virtual void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
     258                virtual void   GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level)=0;
    258259
    259260                virtual void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;

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

    r18283 r18293  
    9292                void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    9393                void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
     94                void   GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
    9495                void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
    9596                void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);

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

    r18283 r18293  
    123123                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
    124124                bool               IsIcefront(void);
    125                 bool   IsFaceOnBoundary(void){_error_("not implemented yet");};
     125                bool        IsFaceOnBoundary(void){_error_("not implemented yet");};
    126126                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
    127127                void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
     128                void               GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented");};
    128129
    129130                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};

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

    r18283 r18293  
    128128                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
    129129                bool               IsIcefront(void);
    130                 bool   IsFaceOnBoundary(void){_error_("not implemented yet");};
     130                bool        IsFaceOnBoundary(void){_error_("not implemented yet");};
    131131                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    132132                void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
     133                void               GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
    133134                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
    134135                void        InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};

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

    r18283 r18293  
    846846        for(i=0;i<NUMVERTICES;i++){
    847847                if(levelset[i]>=0.){
     848                        indicesfront[nrfrontnodes]=i;
     849                        nrfrontnodes++;
     850                }
     851        }
     852
     853        _assert_(nrfrontnodes==2);
     854
     855        /* arrange order of frontnodes such that they are oriented counterclockwise */
     856        if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=NUMVERTICES-1){
     857                int index=indicesfront[0];
     858                indicesfront[0]=indicesfront[1];
     859                indicesfront[1]=index;
     860        }       
     861
     862        IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
     863        /* Return nodes */
     864        for(i=0;i<nrfrontnodes;i++){
     865                for(dir=0;dir<3;dir++){
     866                        xyz_front[3*i+dir]=xyz_list[3*indicesfront[i]+dir];
     867                }
     868        }
     869
     870        *pxyz_front=xyz_front;
     871
     872        xDelete<int>(indicesfront);
     873}/*}}}*/
     874void       Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/
     875
     876        /* Intermediaries */
     877        int i, dir,nrfrontnodes;
     878        IssmDouble  levelset[NUMVERTICES];
     879
     880        /*Recover parameters and values*/
     881        GetInputListOnVertices(&levelset[0],levelsetenum);
     882
     883        int* indicesfront = xNew<int>(NUMVERTICES);
     884        /* Get nodes where there is no ice */
     885        nrfrontnodes=0;
     886        for(i=0;i<NUMVERTICES;i++){
     887                if(levelset[i]==level){
    848888                        indicesfront[nrfrontnodes]=i;
    849889                        nrfrontnodes++;

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

    r18283 r18293  
    101101                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    102102                IssmDouble  TimeAdapt();
    103                 void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
    104                 void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
     103                void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
     104                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
    105105                int         VertexConnectivity(int vertexindex);
    106                 void   VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
     106                void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
    107107                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    108                 void        GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
     108                void          GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
     109                void          GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
    109110                bool        IsZeroLevelset(int levelset_enum);
    110                 bool            IsIcefront(void);
    111                 bool            IsFaceOnBoundary(void);
     111                bool            IsIcefront(void);
     112                bool            IsFaceOnBoundary(void);
    112113
    113114                void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
Note: See TracChangeset for help on using the changeset viewer.