Changeset 16949


Ignore:
Timestamp:
11/26/13 09:54:27 (11 years ago)
Author:
seroussi
Message:

CHG: moved KMatrix viscous couping SSAFS in analysis

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

Legend:

Unmodified
Added
Removed

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

    r16940 r16949  
    830830                case HOFSApproximationEnum:
    831831                        return CreateKMatrixHOFS(element);
     832                case SSAFSApproximationEnum:
     833                        return CreateKMatrixSSAFS(element);
    832834                case NoneApproximationEnum:
    833835                        return NULL;
    … …  
    30693071
    30703072/*Coupling (Tiling)*/
     3073ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/
     3074
     3075        /*compute all stiffness matrices for this element*/
     3076        ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element);
     3077        ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element);
     3078        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     3079
     3080        /*clean-up and return*/
     3081        delete Ke1;
     3082        delete Ke2;
     3083        return Ke;
     3084}/*}}}*/
     3085ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/
     3086
     3087        return NULL;
     3088}/*}}}*/
     3089ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/
     3090
     3091        /*Constants*/
     3092        const int    numdof2d=2*3;
     3093
     3094        /*Intermediaries */
     3095        int         i,j,approximation;
     3096        int         dim=3;
     3097        IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot;
     3098        IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
     3099        IssmDouble  epsilons[6];                    //6 for FS
     3100        IssmDouble  B[3][numdof2d];
     3101        IssmDouble  Bprime[3][numdof2d];
     3102        IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
     3103        IssmDouble  D_scalar;
     3104        IssmDouble  Ke_gg_gaussian[numdof2d][numdof2d];
     3105        IssmDouble  *xyz_list  = NULL;
     3106
     3107        /*Find penta on bed as this is a SSA elements: */
     3108        Element* pentabase=element->GetBasalElement();
     3109        Element* basaltria=pentabase->SpawnBasalElement();
     3110
     3111        /*Initialize Element matrix*/
     3112        ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum);
     3113        element->GetInputValue(&approximation,ApproximationEnum);
     3114
     3115        /*Retrieve all inputs and parameters*/
     3116        element->GetVerticesCoordinates(&xyz_list);
     3117        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
     3118        Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
     3119        Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
     3120        Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
     3121        Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
     3122        Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
     3123
     3124        /* Start  looping on the number of gaussian points: */
     3125        Gauss* gauss=element->NewGauss(5);
     3126        Gauss* gauss_tria=new GaussTria();
     3127        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3128
     3129                gauss->GaussPoint(ig);
     3130                gauss->SynchronizeGaussBase(gauss_tria);
     3131
     3132                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     3133                this->GetBSSA(&B[0][0],basaltria, xyz_list, gauss_tria);
     3134                this->GetBSSAprime(&Bprime[0][0], basaltria,xyz_list, gauss_tria);
     3135
     3136                if(approximation==SSAHOApproximationEnum){
     3137                        element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
     3138                        element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     3139                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     3140                }
     3141                else if (approximation==SSAFSApproximationEnum){
     3142                        element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3143                }
     3144                else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
     3145
     3146                D_scalar=2*newviscosity*gauss->weight*Jdet;
     3147                for (i=0;i<3;i++) D[i][i]=D_scalar;
     3148
     3149                TripleMultiply( &B[0][0],3,numdof2d,1,
     3150                                        &D[0][0],3,3,0,
     3151                                        &Bprime[0][0],3,numdof2d,0,
     3152                                        &Ke_gg_gaussian[0][0],0);
     3153
     3154                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j];
     3155        }
     3156
     3157        /*Transform Coordinate System*/
     3158        basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum);
     3159
     3160        /*Clean up and return*/
     3161        xDelete<IssmDouble>(xyz_list);
     3162        delete basaltria->material;
     3163        delete basaltria;
     3164        delete gauss_tria;
     3165        delete gauss;
     3166        return Ke;
     3167
     3168}/*}}}*/
    30713169ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
    30723170
    … …  
    30943192        ElementMatrix* Ke2=CreateKMatrixHO(element);
    30953193        ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element);
     3194        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     3195
     3196        /*clean-up and return*/
     3197        delete Ke1;
     3198        delete Ke2;
     3199        delete Ke3;
     3200        return Ke;
     3201}/*}}}*/
     3202ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/
     3203
     3204        /*compute all stiffness matrices for this element*/
     3205        ElementMatrix* Ke1=CreateKMatrixFS(element);
     3206        int indices[3]={18,19,20};
     3207        Ke1->StaticCondensation(3,&indices[0]);
     3208        int init = element->FiniteElement();
     3209        element->SetTemporaryElementType(P1Enum);
     3210        ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
     3211        element->SetTemporaryElementType(init);
     3212        ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
    30963213        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    30973214
    … …  
    33433460        xDelete<Node*>(node_list);
    33443461        xDelete<int>(cs_list);
     3462        return Ke;
     3463
     3464}/*}}}*/
     3465ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFS(Element* element){/*{{{*/
     3466
     3467        /*compute all stiffness matrices for this element*/
     3468        ElementMatrix* Ke1=CreateKMatrixCouplingSSAFSViscous(element);
     3469        ElementMatrix* Ke2=CreateKMatrixCouplingSSAFSFriction(element);
     3470        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     3471
     3472        /*clean-up and return*/
     3473        delete Ke1;
     3474        delete Ke2;
     3475        return Ke;
     3476}/*}}}*/
     3477ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFSFriction(Element* element){/*{{{*/
     3478        return NULL;
     3479
     3480}/*}}}*/
     3481ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFSViscous(Element* element){/*{{{*/
     3482
     3483        /*Constants*/
     3484        const int numdofm     = 2 *3;
     3485        const int numdofs     = 4 *6+ 3;
     3486        const int numdoftotal = 2 *numdofm+numdofs;
     3487
     3488        /*Intermediaries */
     3489        int        i,j;
     3490        int        dim=3;
     3491        IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
     3492        IssmDouble B[4][numdofs];
     3493        IssmDouble Bprime[4][numdofm];
     3494        IssmDouble B2[3][numdofm];
     3495        IssmDouble Bprime2[3][numdofs];
     3496        IssmDouble D[4][4]={0.0};            // material matrix, simple scalar matrix.
     3497        IssmDouble D2[3][3]={0.0};            // material matrix, simple scalar matrix.
     3498        IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix
     3499        IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix
     3500        IssmDouble Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point.
     3501        IssmDouble Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point.
     3502        IssmDouble *xyz_list    = NULL;
     3503        Node       *node_list[20];
     3504
     3505        /*Find penta on bed as FS must be coupled to the dofs on the bed: */
     3506        Element* pentabase=element->GetBasalElement();
     3507        Element* basaltria=pentabase->SpawnBasalElement();
     3508
     3509        int vnumnodes = element->GetNumberOfNodesVelocity();
     3510        int pnumnodes = element->GetNumberOfNodesPressure();
     3511        int numnodes  = 2*vnumnodes-1+pnumnodes;
     3512
     3513        /*Prepare node list*/
     3514        int* cs_list = xNew<int>(2*vnumnodes+pnumnodes);
     3515        for(i=0;i<vnumnodes-1;i++){
     3516                node_list[i] = pentabase->GetNode(i);
     3517                cs_list[i]   = XYEnum;
     3518        }
     3519        for(i=0;i<vnumnodes;i++){
     3520                node_list[i+vnumnodes-1] = element->GetNode(i);
     3521                cs_list[i+vnumnodes-1]   = XYZEnum;
     3522        }
     3523        for(i=0;i<pnumnodes;i++){
     3524                node_list[2*vnumnodes-1+i] = element->GetNode(vnumnodes+i);
     3525                cs_list[2*vnumnodes-1+i]   = PressureEnum;
     3526        }
     3527
     3528        /*Initialize Element matrix and return if necessary*/
     3529        ElementMatrix* Ke1=element  ->NewElementMatrix(FSvelocityEnum);
     3530        ElementMatrix* Ke2=pentabase->NewElementMatrix(SSAApproximationEnum);
     3531        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     3532        delete Ke1; delete Ke2;
     3533
     3534        /* Get node coordinates and dof list: */
     3535        element->GetVerticesCoordinates(&xyz_list);
     3536        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3537        Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
     3538        Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
     3539        Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
     3540
     3541        /* Start  looping on the number of gaussian points: */
     3542        Gauss* gauss=element->NewGauss(5);
     3543        Gauss* gauss_tria=new GaussTria();
     3544        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3545
     3546                gauss->GaussPoint(ig);
     3547                gauss->SynchronizeGaussBase(gauss_tria);
     3548
     3549                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     3550                this->GetBSSAFS(&B[0][0],element,xyz_list, gauss);
     3551                this->GetBprimeSSAFSTria(&Bprime[0][0], basaltria,xyz_list, gauss_tria);
     3552                this->GetBSSAFSTria(&B2[0][0], basaltria,xyz_list, gauss_tria);
     3553                this->GetBprimeSSAFS(&Bprime2[0][0], element,xyz_list, gauss);
     3554
     3555                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3556
     3557                D_scalar=2*viscosity*gauss->weight*Jdet;
     3558                for (i=0;i<3;i++) D[i][i]=D_scalar;
     3559                D[3][3]=-gauss->weight*Jdet*FSreconditioning;
     3560                for (i=0;i<3;i++) D2[i][i]=D_scalar;
     3561
     3562                TripleMultiply( &B[0][0],4,numdofs,1,
     3563                                        &D[0][0],4,4,0,
     3564                                        &Bprime[0][0],4,numdofm,0,
     3565                                        &Ke_gg_gaussian[0][0],0);
     3566
     3567                TripleMultiply( &B2[0][0],3,numdofm,1,
     3568                                        &D2[0][0],3,3,0,
     3569                                        &Bprime2[0][0],3,numdofs,0,
     3570                                        &Ke_gg_gaussian2[0][0],0);
     3571
     3572                for( i=0; i<numdofs; i++) for(j=0;j<numdofm; j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     3573                for( i=0; i<numdofm; i++) for(j=0;j<numdofs; j++) Ke_gg2[i][j]+=Ke_gg_gaussian2[i][j];
     3574        }
     3575        for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[i*numdoftotal+j+numdofs]+=Ke_gg[i][j];
     3576        for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[(i+numdofs)*numdoftotal+j]+=Ke_gg2[i][j];
     3577
     3578        /*Transform Coordinate System*/
     3579        element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
     3580
     3581        /*Clean-up and return*/
     3582        xDelete<int>(cs_list);
     3583        xDelete<IssmDouble>(xyz_list);
     3584        delete basaltria->material; delete basaltria;
     3585        delete gauss;
     3586        delete gauss_tria;
    33453587        return Ke;
    33463588
    … …  
    37563998        xDelete<IssmDouble>(dbasis);
    37573999}/*}}}*/
     4000void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4001        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     4002         * For node i, Bi can be expressed in the actual coordinate system
     4003         * by:
     4004         *       Bi=[ dh/dx          0       0   0 ]
     4005         *          [   0           dh/dy    0   0 ]
     4006         *          [ 1/2*dh/dy  1/2*dh/dx   0   0 ]
     4007         *          [   0            0       0   h ]
     4008         * where h is the interpolation function for node i.
     4009         *
     4010         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
     4011         */
     4012
     4013        int i;
     4014        IssmDouble dbasismini[3][7];
     4015        IssmDouble basis[6];
     4016
     4017        /*Get dbasis in actual coordinate system: */
     4018        element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
     4019        element->NodalFunctionsP1(basis,gauss);
     4020
     4021        /*Build B: */
     4022        for(i=0;i<6;i++){
     4023                B[(3*7+6)*0+3*i+0] = dbasismini[0][i];
     4024                B[(3*7+6)*0+3*i+1] = 0.;
     4025                B[(3*7+6)*0+3*i+2] = 0.;
     4026                B[(3*7+6)*1+3*i+0] = 0.;
     4027                B[(3*7+6)*1+3*i+1] = dbasismini[1][i];
     4028                B[(3*7+6)*1+3*i+2] = 0.;
     4029                B[(3*7+6)*2+3*i+0] = 0.5*dbasismini[1][i];
     4030                B[(3*7+6)*2+3*i+1] = 0.5*dbasismini[0][i];
     4031                B[(3*7+6)*2+3*i+2] = 0.;
     4032                B[(3*7+6)*3+3*i+0] = 0.;
     4033                B[(3*7+6)*3+3*i+1] = 0.;
     4034                B[(3*7+6)*3+3*i+2] = 0.;
     4035        }
     4036        for(i=0;i<1;i++){
     4037                B[(3*7+6)*0+3*(6+i)+0] = 0.;
     4038                B[(3*7+6)*0+3*(6+i)+1] = 0.;
     4039                B[(3*7+6)*0+3*(6+i)+2] = 0.;
     4040                B[(3*7+6)*1+3*(6+i)+0] = 0.;
     4041                B[(3*7+6)*1+3*(6+i)+1] = 0.;
     4042                B[(3*7+6)*1+3*(6+i)+2] = 0.;
     4043                B[(3*7+6)*2+3*(6+i)+0] = 0.;
     4044                B[(3*7+6)*2+3*(6+i)+1] = 0.;
     4045                B[(3*7+6)*2+3*(6+i)+2] = 0.;
     4046                B[(3*7+6)*3+3*(6+i)+0] = 0.;
     4047                B[(3*7+6)*3+3*(6+i)+1] = 0.;
     4048                B[(3*7+6)*3+3*(6+i)+2] = 0.;
     4049        }
     4050
     4051        for(i=0;i<6;i++){ //last column not for the bubble function
     4052                B[(3*7+6)*0+7*3+i] = 0;
     4053                B[(3*7+6)*1+7*3+i] = 0;
     4054                B[(3*7+6)*2+7*3+i] = 0;
     4055                B[(3*7+6)*3+7*3+i] = basis[i];
     4056        }
     4057}/*}}}*/
     4058void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4059        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
     4060         * For node i, Bprimei can be expressed in the actual coordinate system
     4061         * by:
     4062         *       Bprimei=[ 2*dh/dx    dh/dy   0   0 ]
     4063         *               [  dh/dx    2*dh/dy  0   0 ]
     4064         *               [  dh/dy     dh/dx   0   0 ]
     4065         * where h is the interpolation function for node i.
     4066         *
     4067         * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
     4068         */
     4069
     4070        int    i;
     4071        IssmDouble dbasismini[3][7];
     4072
     4073        /*Get dbasis in actual coordinate system: */
     4074        element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
     4075
     4076        /*Build Bprime: */
     4077        for(i=0;i<6;i++){
     4078                Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i];
     4079                Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i];
     4080                Bprime[(3*7+6)*0+3*i+2] = 0.;
     4081                Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i];
     4082                Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i];
     4083                Bprime[(3*7+6)*1+3*i+2] = 0.;
     4084                Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i];
     4085                Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i];
     4086                Bprime[(3*7+6)*2+3*i+2] = 0.;
     4087        }
     4088
     4089        for(i=0;i<1;i++){ //Add zeros for the bubble function
     4090                Bprime[(3*7+6)*0+3*(6+i)+0] = 0.;
     4091                Bprime[(3*7+6)*0+3*(6+i)+1] = 0.;
     4092                Bprime[(3*7+6)*0+3*(6+i)+2] = 0.;
     4093                Bprime[(3*7+6)*1+3*(6+i)+0] = 0.;
     4094                Bprime[(3*7+6)*1+3*(6+i)+1] = 0.;
     4095                Bprime[(3*7+6)*1+3*(6+i)+2] = 0.;
     4096                Bprime[(3*7+6)*2+3*(6+i)+0] = 0.;
     4097                Bprime[(3*7+6)*2+3*(6+i)+1] = 0.;
     4098                Bprime[(3*7+6)*2+3*(6+i)+2] = 0.;
     4099        }
     4100
     4101        for(i=0;i<6;i++){ //last column not for the bubble function
     4102                Bprime[(3*7+6)*0+7*3+i] = 0.;
     4103                Bprime[(3*7+6)*1+7*3+i] = 0.;
     4104                Bprime[(3*7+6)*2+7*3+i] = 0.;
     4105        }
     4106}/*}}}*/
     4107void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4108        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     4109         * For node i, Bi can be expressed in the actual coordinate system
     4110         * by:
     4111         *       Bi=[   dN/dx         0     ]
     4112         *          [       0       dN/dy   ]
     4113         *          [  1/2*dN/dy  1/2*dN/dx ]
     4114         * where N is the finiteelement function for node i.
     4115         *
     4116         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
     4117         */
     4118
     4119        /*Fetch number of nodes for this finite element*/
     4120        int numnodes = element->GetNumberOfNodes();
     4121
     4122        /*Get nodal functions derivatives*/
     4123        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     4124        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     4125
     4126        /*Build B': */
     4127        for(int i=0;i<numnodes;i++){
     4128                B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     4129                B[2*numnodes*0+2*i+1] = 0.;
     4130                B[2*numnodes*1+2*i+0] = 0.;
     4131                B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     4132                B[2*numnodes*2+2*i+0] = 0.5*dbasis[1*numnodes+i];
     4133                B[2*numnodes*2+2*i+1] = 0.5*dbasis[0*numnodes+i];
     4134        }
     4135
     4136        /*Clean-up*/
     4137        xDelete<IssmDouble>(dbasis);
     4138}/*}}}*/
     4139void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4140        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
     4141         * For node i, Bprimei can be expressed in the actual coordinate system
     4142         * by:
     4143         *       Bprimei=[  dN/dx    0   ]
     4144         *               [    0    dN/dy ]
     4145         *               [  dN/dy  dN/dx ]
     4146         N               [  dN/dx  dN/dy ]
     4147         * where N is the finiteelement function for node i.
     4148         *
     4149         * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
     4150         */
     4151
     4152        /*Fetch number of nodes for this finite element*/
     4153        int numnodes = element->GetNumberOfNodes();
     4154
     4155        /*Get nodal functions*/
     4156        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     4157        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     4158
     4159        /*Build Bprime: */
     4160        for(int i=0;i<numnodes;i++){
     4161                Bprime[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     4162                Bprime[2*numnodes*0+2*i+1] = 0.;
     4163                Bprime[2*numnodes*1+2*i+0] = 0.;
     4164                Bprime[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     4165                Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
     4166                Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
     4167                Bprime[2*numnodes*3+2*i+0] = dbasis[0*numnodes+i];
     4168                Bprime[2*numnodes*3+2*i+1] = dbasis[1*numnodes+i];
     4169        }
     4170
     4171        /*Clean-up*/
     4172        xDelete<IssmDouble>(dbasis);
     4173}/*}}}*/
    37584174void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
    37594175

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

    r16940 r16949  
    7171                void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
    7272                /*Coupling*/
     73                ElementMatrix* CreateKMatrixSSA3d(Element* element);
     74                ElementMatrix* CreateKMatrixSSA3dFriction(Element* element);
     75                ElementMatrix* CreateKMatrixSSA3dViscous(Element* element);
    7376                ElementMatrix* CreateKMatrixHOFS(Element* element);
    7477                ElementMatrix* CreateKMatrixSSAHO(Element* element);
     78                ElementMatrix* CreateKMatrixSSAFS(Element* element);
    7579                ElementMatrix* CreateKMatrixCouplingHOFS(Element* element);
    7680                ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
    7781                ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);
    7882                ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
     83                ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element);
     84                ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element);
     85                ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element);
    7986                ElementVector* CreatePVectorSSAHO(Element* element);
    8087                ElementVector* CreatePVectorSSAFS(Element* element);
    … …  
    8794                ElementVector* CreatePVectorCouplingHOFSViscous(Element* element);
    8895                void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     96                void GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     97                void GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     98                void GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     99                void GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    89100                void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
    90101                void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);

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

    r16947 r16949  
    120120                virtual void   NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    121121                virtual void   NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
     122                virtual void   NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    122123                virtual void   NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    123124                virtual void   NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;

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

    r16947 r16949  
    24772477        _assert_(gauss->Enum()==GaussPentaEnum);
    24782478        this->GetNodalFunctionsP1Derivatives(dbasis,xyz_list,(GaussPenta*)gauss);
     2479
     2480}
     2481/*}}}*/
     2482/*FUNCTION Penta::NodalFunctionsMINIDerivatives{{{*/
     2483void Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
     2484
     2485        _assert_(gauss->Enum()==GaussPentaEnum);
     2486        this->GetNodalFunctionsMINIDerivatives(dbasis,xyz_list,(GaussPenta*)gauss);
    24792487
    24802488}

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

    r16947 r16949  
    254254                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    255255                void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     256                void           NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    256257                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    257258                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);

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

    r16947 r16949  
    113113                void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    114114                void        NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     115                void        NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    115116                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    116117                bool        NoIceInElement(){_error_("not implemented yet");};

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

    r16947 r16949  
    287287                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    288288                void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     289                void           NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    289290                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    290291                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
Note: See TracChangeset for help on using the changeset viewer.