Changeset 15425


Ignore:
Timestamp:
07/03/13 20:40:06 (12 years ago)
Author:
seroussi
Message:

CHG: adding Stokes P1 GLS

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r15422 r15425  
    68306830        switch(fe_stokes){
    68316831                case 0:
    6832                         printf("fe %i\n",fe_stokes);
    68336832                        /*compute all stiffness matrices for this element*/
    68346833                        Ke1=CreateKMatrixDiagnosticStokesViscous();
    … …  
    68386837                case 1:
    68396838                        /*compute all stiffness matrices for this element*/
    6840                         Ke1=CreateKMatrixDiagnosticStokesViscous();
    6841                         Ke2=CreateKMatrixDiagnosticStokesFriction();
     6839                        Ke1=CreateKMatrixDiagnosticStokesGLSViscous();
     6840                        Ke2=CreateKMatrixDiagnosticStokesGLSFriction();
    68426841                        Ke =new ElementMatrix(Ke1,Ke2);
    68436842                        break;
    … …  
    69146913}
    69156914/*}}}*/
     6915/*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSViscous {{{*/
     6916ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){
     6917
     6918        /*Intermediaries */
     6919        int        i,approximation;
     6920        IssmDouble Jdet,viscosity,stokesreconditioning;
     6921        IssmDouble xyz_list[NUMVERTICES][3];
     6922        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     6923        IssmDouble B[8][24];
     6924        IssmDouble B_prime[8][24];
     6925        IssmDouble D_scalar;
     6926        IssmDouble D[8][8]={0.0};
     6927        IssmDouble Ke_temp[24][24]={1.0}; //for the six nodes
     6928        GaussPenta *gauss=NULL;
     6929
     6930        /*If on water or not Stokes, skip stiffness: */
     6931        inputs->GetInputValue(&approximation,ApproximationEnum);
     6932        if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
     6933        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     6934
     6935        /*Retrieve all inputs and parameters*/
     6936        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6937        parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     6938        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     6939        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     6940        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     6941
     6942        /* Start  looping on the number of gaussian points: */
     6943        gauss=new GaussPenta(5,5);
     6944        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6945
     6946                gauss->GaussPoint(ig);
     6947
     6948                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6949                GetBStokesGLS(&B[0][0],&xyz_list[0][0],gauss);
     6950                GetBprimeStokesGLS(&B_prime[0][0],&xyz_list[0][0],gauss);
     6951
     6952                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     6953                material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     6954
     6955                D_scalar=gauss->weight*Jdet;
     6956                for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
     6957                for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
     6958
     6959                TripleMultiply( &B[0][0],8,24,1,
     6960                                        &D[0][0],8,8,0,
     6961                                        &B_prime[0][0],8,24,0,
     6962                                        &Ke_temp[0][0],1);
     6963        }
     6964
     6965        /*Transform Coordinate System*/
     6966        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     6967
     6968        /*Clean up and return*/
     6969        delete gauss;
     6970        return Ke;
     6971}
     6972/*}}}*/
    69166973/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction{{{*/
    69176974ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){
     6975
     6976        /*Constants*/
     6977        const int numdof=NUMVERTICES*NDOF4;
     6978        const int numdof2d=NUMVERTICES2D*NDOF4;
     6979
     6980        /*Intermediaries */
     6981        int        i,j;
     6982        int        analysis_type,approximation;
     6983        IssmDouble alpha2,Jdet2d;
     6984        IssmDouble stokesreconditioning,viscosity;
     6985        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     6986        IssmDouble xyz_list[NUMVERTICES][3];
     6987        IssmDouble xyz_list_tria[NUMVERTICES2D][3];
     6988        IssmDouble LStokes[2][numdof2d];
     6989        IssmDouble DLStokes[2][2]={0.0};
     6990        IssmDouble Ke_drag_gaussian[numdof2d][numdof2d];
     6991        Friction*  friction=NULL;
     6992        GaussPenta *gauss=NULL;
     6993
     6994        /*If on water or not Stokes, skip stiffness: */
     6995        inputs->GetInputValue(&approximation,ApproximationEnum);
     6996        if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return NULL;
     6997        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     6998
     6999        /*Retrieve all inputs and parameters*/
     7000        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7001        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     7002        parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7003        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     7004        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     7005        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     7006        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     7007
     7008        /*build friction object, used later on: */
     7009        friction=new Friction("3d",inputs,matpar,analysis_type);
     7010
     7011        /* Start  looping on the number of gaussian points: */
     7012        gauss=new GaussPenta(0,1,2,2);
     7013        for(int ig=gauss->begin();ig<gauss->end();ig++){
     7014
     7015                gauss->GaussPoint(ig);
     7016
     7017                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
     7018                GetLStokes(&LStokes[0][0], gauss);
     7019
     7020                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7021                material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7022
     7023                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     7024
     7025                DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
     7026                DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
     7027
     7028                TripleMultiply( &LStokes[0][0],2,numdof2d,1,
     7029                                        &DLStokes[0][0],2,2,0,
     7030                                        &LStokes[0][0],2,numdof2d,0,
     7031                                        &Ke_drag_gaussian[0][0],0);
     7032
     7033                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j];
     7034        }
     7035
     7036        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     7037        //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     7038
     7039        /*Clean up and return*/
     7040        delete gauss;
     7041        delete friction;
     7042        return Ke;
     7043}
     7044/*}}}*/
     7045/*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSFriction{{{*/
     7046ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSFriction(void){
    69187047
    69197048        /*Constants*/

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

    r15387 r15425  
    249249                ElementMatrix* CreateKMatrixDiagnosticStokes(void);
    250250                ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
     251                ElementMatrix* CreateKMatrixDiagnosticStokesGLSViscous(void);
    251252                ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
     253                ElementMatrix* CreateKMatrixDiagnosticStokesGLSFriction(void);
    252254                ElementMatrix* CreateKMatrixDiagnosticVert(void);
    253255                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);

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

    r15326 r15425  
    314314}
    315315/*}}}*/
     316/*FUNCTION PentaRef::GetBStokesGLS {{{*/
     317void PentaRef::GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     318
     319        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.
     320         * For node i, Bi can be expressed in the actual coordinate system
     321         * by:          Bi=[ dh/dx          0              0       0  ]
     322         *                                      [   0           dh/dy           0       0  ]
     323         *                                      [   0             0           dh/dy     0  ]
     324         *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
     325         *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
     326         *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
     327         *                                      [   0             0             0       h  ]
     328         *                                      [ dh/dx         dh/dy         dh/dz     0  ]
     329         *      where h is the interpolation function for node i.
     330         *      Same thing for Bb except the last column that does not exist.
     331         */
     332
     333        int i;
     334
     335        IssmDouble dh1dh6[3][NUMNODESP1];
     336        IssmDouble l1l6[NUMNODESP1];
     337
     338        /*Get dh1dh7 in actual coordinate system: */
     339        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     340        GetNodalFunctionsP1(l1l6, gauss);
     341
     342        /*Build B: */
     343        for (i=0;i<NUMNODESP1;i++){
     344                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh6[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
     345                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0.;
     346                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0.;
     347                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0.;
     348                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh6[1][i];
     349                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0.;
     350                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.;
     351                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.;
     352                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh6[2][i];
     353                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=.5*dh1dh6[1][i];
     354                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=.5*dh1dh6[0][i];
     355                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0.;
     356                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=.5*dh1dh6[2][i];
     357                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0.;
     358                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=.5*dh1dh6[0][i];
     359                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0.;
     360                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=.5*dh1dh6[2][i];
     361                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=.5*dh1dh6[1][i];
     362                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=0.;
     363                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=0.;
     364                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=0.;
     365                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=dh1dh6[0][i];
     366                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=dh1dh6[1][i];
     367                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=dh1dh6[2][i];
     368        }
     369
     370        for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     371                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0.;
     372                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0.;
     373                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0.;
     374                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0.;
     375                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0.;
     376                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0.;
     377                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i];
     378                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0.;
     379        }
     380
     381}
     382/*}}}*/
    316383/*FUNCTION PentaRef::GetBprimeStokes {{{*/
    317384void PentaRef::GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
    … …  
    376443                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
    377444                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
     445                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
     446        }
     447
     448}
     449/*}}}*/
     450/*FUNCTION PentaRef::GetBprimeStokesGLS {{{*/
     451void PentaRef::GetBprimeStokesGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
     452        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
     453         *      For node i, Bi' can be expressed in the actual coordinate system
     454         *      by:
     455         *                              Bi'=[  dh/dx   0          0       0]
     456         *                                       [   0      dh/dy      0       0]
     457         *                                       [   0      0         dh/dz    0]
     458         *                                       [  dh/dy   dh/dx      0       0]
     459         *                                       [  dh/dz   0        dh/dx     0]
     460         *                                       [   0      dh/dz    dh/dy     0]
     461         *                                       [  dh/dx   dh/dy    dh/dz     0]
     462         *                                       [   0      0          0       h]
     463         *      where h is the interpolation function for node i.
     464         *
     465         *      Same thing for the bubble fonction except that there is no fourth column
     466         */
     467
     468        int i;
     469        IssmDouble dh1dh6[3][NUMNODESP1];
     470        IssmDouble l1l6[NUMNODESP1];
     471
     472        /*Get dh1dh7 in actual coordinate system: */
     473        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     474        GetNodalFunctionsP1(l1l6, gauss);
     475
     476        /*B_primeuild B_prime: */
     477        for (i=0;i<NUMNODESP1;i++){
     478                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh6[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
     479                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0.;
     480                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0.;
     481                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0.;
     482                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh6[1][i];
     483                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0.;
     484                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.;
     485                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.;
     486                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh6[2][i];
     487                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh6[1][i];
     488                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh6[0][i];
     489                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0.;
     490                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh6[2][i];
     491                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0.;
     492                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh6[0][i];
     493                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0.;
     494                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh6[2][i];
     495                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh6[1][i];
     496                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh6[0][i];
     497                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh6[1][i];
     498                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh6[2][i];
     499                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0.;
     500                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0.;
     501                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0.;
     502        }
     503
     504        for (i=0;i<NUMNODESP1;i++){ //last column
     505                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0.;
     506                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0.;
     507                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0.;
     508                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0.;
     509                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0.;
     510                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0.;
     511                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0.;
    378512                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
    379513        }

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

    r15326 r15425  
    4242                void GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4343                void GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     44                void GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4445                void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss);
    4546                void GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4647                void GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
     48                void GetBprimeStokesGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
    4749                void GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4850                void GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, GaussPenta* gauss);
Note: See TracChangeset for help on using the changeset viewer.