Changeset 394


Ignore:
Timestamp:
05/13/09 15:01:35 (16 years ago)
Author:
Eric.Larour
Message:

New stokes elements

Location:
issm/trunk/src
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp

    r377 r394  
    9191        /*Free data: */
    9292        xfree((void**)&gridonstokes);
    93        
     93
    9494        /*Assign output pointer: */
    9595        *pconstraints=constraints;

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp

    r387 r394  
    514514        xfree((void**)&model->gridonstokes);
    515515        xfree((void**)&model->borderstokes);
    516                
     516
    517517        cleanup_and_return:
    518518

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp

    r387 r394  
    5151        int numberofsegs_diag_stokes;
    5252        int count;
     53
    5354
    5455        /*Create loads: */
    … …  
    134135        #endif
    135136
    136                 if ((model->gridonbed[i]) && (model->gridonicesheet[i]) && (model->gridonstokes[i]))
     137                if ((model->gridonbed[i]) && (model->gridonicesheet[i]) && (model->gridonstokes[i])){
    137138
    138                 pengrid_id=count+1; //matlab indexing
    139                 pengrid_node_id=i+1;
    140                 pengrid_dof=1;
    141                 pengrid_penalty_offset=model->penalty_offset;
     139                        pengrid_id=count+1; //matlab indexing
     140                        pengrid_node_id=i+1;
     141                        pengrid_dof=1;
     142                        pengrid_penalty_offset=model->penalty_offset;
    142143
    143                 pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
    144                
    145                 loads->AddObject(icefront);
    146                 count++;
     144                        pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
     145                       
     146                        loads->AddObject(icefront);
     147                        count++;
     148                }
    147149        #ifdef _PARALLEL_
    148150        } //if((model->my_grids[i]==1))

  • issm/trunk/src/c/ModelProcessorx/Model.cpp

    r388 r394  
    307307        ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL);
    308308        ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL);
     309        ModelFetchData((void**)&model->isstokes,NULL,NULL,model_handle,"isstokes","Integer",NULL);
    309310
    310311        /*!Get drag_type, drag and p,q: */

  • issm/trunk/src/c/objects/Friction.cpp

    r1 r394  
    3535        friction->p=UNDEF;
    3636        friction->q=UNDEF;
     37        friction->analysis_type=UNDEF;
    3738       
    3839        return 1;

  • issm/trunk/src/c/objects/Friction.h

    r1 r394  
    1919        double  p;
    2020        double  q;
     21        int     analysis_type;
    2122
    2223};

  • issm/trunk/src/c/objects/Matice.cpp

    r312 r394  
    306306}
    307307
     308#undef __FUNCT__
     309#define __FUNCT__ "MatIce::GetViscosity3dStokes"
     310void  Matice::GetViscosity3dStokes(double* pviscosity3d, double* epsilon){
     311
     312        /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]:
     313         *
     314         *                                 2*B
     315         * viscosity3d= -------------------------------------------------------------------
     316         *     2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
     317         *
     318         *     where mu is the viscotiy, B the flow law parameter , (u,v) the velocity
     319         *     vector, and n the flow law exponent.
     320         *
     321         * If epsilon is NULL, it means this is the first time Emg is being run, and we
     322         * return g, initial viscosity.
     323         */
     324       
     325        /*output: */
     326        double viscosity3d;
     327
     328        /*input strain rate: */
     329        double exx,eyy,exy,exz,eyz,ezz;
     330
     331        /*Intermediary value A and exponent e: */
     332        double A,e;
     333        double eps0;
     334
     335        eps0=10^-27;
     336       
     337        if (n==1){
     338                /*Viscous behaviour! viscosity3d=B: */
     339                viscosity3d=B;
     340        }
     341        else{
     342                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
     343                                (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
     344                        viscosity3d=pow(10,14);
     345                }
     346                else{
     347
     348                        /*Retrive strain rate components: */
     349                        exx=*(epsilon+0);
     350                        eyy=*(epsilon+1);
     351                        ezz=*(epsilon+2); //not used
     352                        exy=*(epsilon+3);
     353                        exz=*(epsilon+4);
     354                        eyz=*(epsilon+5);
     355
     356                        /*Build viscosity: viscosity3d=2*B/(2*A^e) */
     357                        A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2);
     358                        if(A==0){
     359                                /*Maxiviscosity3dm viscosity for 0 shear areas: */
     360                                viscosity3d=4.5*pow(10,17);
     361                        }
     362                        else{
     363                                e=(n-1)/2/n;
     364                       
     365                                viscosity3d=2*B/(2*pow(A,e));
     366                        }
     367                }
     368        }
     369        #ifdef _DEBUG_
     370        _printf_("Viscosity %lf\n",viscosity3d);
     371        #endif
     372
     373        /*Assign output pointers:*/
     374        *pviscosity3d=viscosity3d;
     375}
    308376double Matice::GetB(){
    309377        return B;

  • issm/trunk/src/c/objects/Matice.h

    r308 r394  
    3535                void  GetViscosity2(double* pviscosity2, double* pepsilon);
    3636                void  GetViscosity3d(double* pviscosity3d, double* pepsilon);
     37                void  GetViscosity3dStokes(double* pviscosity3d, double* epsilon);
    3738                Object* copy();
    3839                double GetB();

  • issm/trunk/src/c/objects/Pengrid.cpp

    r387 r394  
    136136        return my_rank;
    137137}
    138 void  Pengrid::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
     138void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type){return;}
    139139
    140140#undef __FUNCT__
    … …  
    181181#define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
    182182void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
    183         throw ErrorException(__FUNCT__," not implemented yet!");
    184 }
    185                
     183
     184        if ((analysis_type==DiagnosticStokesAnalysisEnum())){
     185
     186                PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type);
     187
     188        }
     189        else{
     190                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
     191        }
     192
     193}
     194
     195#undef __FUNCT__
     196#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
     197void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type){
     198       
     199        const int numgrids=1;
     200        const int NDOF3=3;
     201        const int numdof=numgrids*NDOF3;
     202        int       doflist[numdof];
     203        int       numberofdofspernode;
     204
     205        int dofs1=0;
     206        int dofs2=1;
     207        double slope[2];
     208        int found=0;
     209        double Ke[4][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
     210       
     211        ParameterInputs* inputs=NULL;
     212       
     213
     214        /*recover pointers: */
     215        inputs=(ParameterInputs*)vinputs;
     216
     217        /*Get dof list: */
     218        GetDofList(&doflist[0],&numberofdofspernode);
     219       
     220        /*recover slope: */
     221        found=inputs->Recover("bedslopex",&slope[0],1,&dofs1,numgrids,(void**)&node);
     222        if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!");
     223        found=inputs->Recover("bedslopey",&slope[1],1,&dofs2,numgrids,(void**)&node);
     224        if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!");
     225
     226        //Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
     227        Ke[0][0]=-slope[0]*kmax*pow(10.0,penalty_offset);
     228        Ke[1][1]=-slope[1]*kmax*pow(10.0,penalty_offset);
     229        Ke[2][2]=kmax*pow(10,penalty_offset);
     230       
     231       
     232        /*Add Ke to global matrix Kgg: */
     233        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
     234}
     235
    186236#undef __FUNCT__
    187237#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
    … …  
    194244}
    195245
     246
     247void  Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
     248
     249        int j;
     250        int doflist_per_node[MAXDOFSPERNODE];
     251        int numberofdofspernode;
     252       
     253        node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
     254        for(j=0;j<numberofdofspernode;j++){
     255                doflist[j]=doflist_per_node[j];
     256        }
     257
     258        /*Assign output pointers:*/
     259        *pnumberofdofspernode=numberofdofspernode;
     260}

  • issm/trunk/src/c/objects/Pengrid.h

    r387 r394  
    4646                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
    4747                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
     48                void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type);
     49                void  GetDofList(int* doflist,int* pnumberofdofspernode);
    4850                Object* copy();
    4951

  • issm/trunk/src/c/objects/Penta.cpp

    r387 r394  
    297297
    298298        }
     299        else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
     300
     301                CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type);
     302
     303        }
    299304        else if (
    300305                        (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
    … …  
    323328        /* node data: */
    324329        const int    numgrids=6;
    325         const int    numdofs=2*numgrids;
     330        const int    numdof=2*numgrids;
    326331        double       xyz_list[numgrids][3];
    327         int          doflist[numdofs];
     332        int          doflist[numdof];
    328333        int          numberofdofspernode;
    329334       
    … …  
    363368
    364369        /* matrices: */
    365         double B[5][numdofs];
    366         double Bprime[5][numdofs];
    367         double L[2][numdofs];
     370        double B[5][numdof];
     371        double Bprime[5][numdof];
     372        double L[2][numdof];
    368373        double D[5][5]={{ 0,0,0,0,0 },{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}};              // material matrix, simple scalar matrix.
    369374        double D_scalar;
    … …  
    372377
    373378        /* local element matrices: */
    374         double Ke_gg[numdofs][numdofs]; //local element stiffness matrix
    375         double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
    376         double Ke_gg_drag_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag
     379        double Ke_gg[numdof][numdof]; //local element stiffness matrix
     380        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
     381        double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    377382        double Jdet;
    378383       
    … …  
    433438
    434439                /* Set Ke_gg to 0: */
    435                 for(i=0;i<numdofs;i++){
    436                         for(j=0;j<numdofs;j++){
     440                for(i=0;i<numdof;i++){
     441                        for(j=0;j<numdof;j++){
    437442                                Ke_gg[i][j]=0.0;
    438443                        }
    … …  
    521526               
    522527                                /*  Do the triple product tB*D*Bprime: */
    523                                 TripleMultiply( &B[0][0],5,numdofs,1,
     528                                TripleMultiply( &B[0][0],5,numdof,1,
    524529                                                &D[0][0],5,5,0,
    525                                                 &Bprime[0][0],5,numdofs,0,
     530                                                &Bprime[0][0],5,numdof,0,
    526531                                                &Ke_gg_gaussian[0][0],0);
    527532
    528533                                /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    529                                 for( i=0; i<numdofs; i++){
    530                                         for(j=0;j<numdofs;j++){
     534                                for( i=0; i<numdof; i++){
     535                                        for(j=0;j<numdof;j++){
    531536                                                Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    532537                                        }
    … …  
    537542
    538543                /*Add Ke_gg to global matrix Kgg: */
    539                 MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     544                MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    540545       
    541546                //Deal with 2d friction at the bedrock interface
    … …  
    576581        const int    numgrids=6;
    577582        const int    NDOF1=1;
    578         const int    numdofs=NDOF1*numgrids;
     583        const int    numdof=NDOF1*numgrids;
    579584        double       xyz_list[numgrids][3];
    580         int          doflist[numdofs];
     585        int          doflist[numdof];
    581586        int          numberofdofspernode;
    582587
    … …  
    598603
    599604        /* matrices: */
    600         double  Ke_gg[numdofs][numdofs];
    601         double  Ke_gg_gaussian[numdofs][numdofs];
     605        double  Ke_gg[numdof][numdof];
     606        double  Ke_gg_gaussian[numdof][numdof];
    602607        double  Jdet;
    603608        double  B[NDOF1][numgrids];
    … …  
    629634
    630635        /* Set Ke_gg to 0: */
    631         for(i=0;i<numdofs;i++){
    632                 for(j=0;j<numdofs;j++){
     636        for(i=0;i<numdof;i++){
     637                for(j=0;j<numdof;j++){
    633638                        Ke_gg[i][j]=0.0;
    634639                }
    … …  
    684689
    685690                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    686                         for( i=0; i<numdofs; i++){
    687                                 for(j=0;j<numdofs;j++){
     691                        for( i=0; i<numdof; i++){
     692                                for(j=0;j<numdof;j++){
    688693                                        Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    689694                                }
    … …  
    693698
    694699        /*Add Ke_gg to global matrix Kgg: */
    695         MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     700        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    696701       
    697702        cleanup_and_return:
    … …  
    705710
    706711#undef __FUNCT__
     712#define __FUNCT__ "Penta:CreateKMatrixDiagnosticStokes"
     713void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type){
     714
     715        int i,j;
     716
     717        int numgrids=6;
     718        int DOFPERGRID=4;
     719        int numdof=numgrids*DOFPERGRID;
     720        int doflist[numdof];
     721        int numberofdofspernode;
     722
     723        int numgrids2d=3;
     724        int numdof2d=numgrids2d*DOFPERGRID;
     725
     726        int   dofs[3]={0,1,2};
     727
     728        double K_terms[numdof][numdof];
     729
     730        /*Material properties: */
     731        double         gravity,rho_ice,rho_water;
     732
     733       
     734        /*Collapsed formulation: */
     735        Tria*  tria=NULL;
     736
     737        /*Grid data: */
     738        double        xyz_list[numgrids][3];
     739       
     740        /*parameters: */
     741        double             xyz_list_tria[numgrids2d][3];
     742        double             surface_normal[3];
     743        double             bed_normal[3];
     744        double         vxvyvz_list[numgrids][3];
     745        double         thickness;
     746
     747        /*matrices: */
     748        double     Ke_temp[27][27]; //for the six nodes and the bubble
     749        double     Ke_reduced[numdof][numdof]; //for the six nodes only
     750        double     Ke_gaussian[27][27];
     751        double     Ke_drag_gaussian[numdof2d][numdof2d];
     752        double     B[8][27];
     753        double     B_prime[8][27];
     754        double     LStokes[14][numdof2d];
     755        double     LprimeStokes[14][numdof2d];
     756        double     Jdet;
     757        double     Jdet2d;
     758        double     D[8][8]={{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 }};         
     759        double     D_scalar;
     760        double     tBD[numdof][8];
     761        double     DLStokes[14][14];
     762        double     tLDStokes[numdof2d][14];
     763       
     764        /* gaussian points: */
     765        int     num_area_gauss;
     766        int     igarea,igvert;
     767        double* first_gauss_area_coord  =  NULL;
     768        double* second_gauss_area_coord =  NULL;
     769        double* third_gauss_area_coord  =  NULL;
     770        double* vert_gauss_coord = NULL;
     771        double* area_gauss_weights  =  NULL;
     772        double* vert_gauss_weights  =  NULL;
     773
     774        /* specific gaussian point: */
     775        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
     776        double  gauss_coord[4];
     777        double  gauss_coord_tria[3];
     778        int area_order=5;
     779        int     num_vert_gauss=5;
     780
     781        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     782        double  viscosity;
     783        double  alpha2_list[numgrids2d];
     784        double  alpha2_gauss;
     785       
     786        ParameterInputs* inputs=NULL;
     787
     788        /*recover pointers: */
     789        inputs=(ParameterInputs*)vinputs;
     790
     791        /* Set K_terms  to 0: */
     792        for(i=0;i<numdof;i++){
     793                for(j=0;j<numdof;j++){
     794                        K_terms[i][j]=0.0;
     795                }
     796        }
     797       
     798        /*recovre material parameters: */
     799        rho_water=matpar->GetRhoWater();
     800        rho_ice=matpar->GetRhoIce();
     801        gravity=matpar->GetG();
     802
     803
     804        /* Get node coordinates and dof list: */
     805        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     806        GetDofList(&doflist[0],&numberofdofspernode);
     807
     808        /*recover extra inputs from users, at current convergence iteration: */
     809        inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
     810
     811        /*Initialize Ke_temp to zero befor adding terms */
     812        for (i=0;i<27;i++){
     813                for (j=0;j<27;j++){
     814                        Ke_temp[i][j]=0;
     815                }
     816        }
     817
     818        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     819           get tria gaussian points as well as segment gaussian points. For tria gaussian
     820           points, order of integration is 2, because we need to integrate the product tB*D*B'
     821           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     822           points, same deal, which yields 3 gaussian points.*/
     823
     824
     825        area_order=5;
     826        num_vert_gauss=5;
     827
     828        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     829        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     830
     831        /* Start  looping on the number of gaussian points: */
     832        for (igarea=0; igarea<num_area_gauss; igarea++){
     833                for (igvert=0; igvert<num_vert_gauss; igvert++){
     834                        /*Pick up the gaussian point: */
     835                        area_gauss_weight=*(area_gauss_weights+igarea);
     836                        vert_gauss_weight=*(vert_gauss_weights+igvert);
     837                        gauss_weight=area_gauss_weight*vert_gauss_weight;
     838                        gauss_coord[0]=*(first_gauss_area_coord+igarea);
     839                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
     840                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
     841                        gauss_coord[3]=*(vert_gauss_coord+igvert);
     842
     843
     844                        /*Compute thickness: */
     845                        GetParameterValue(&thickness, &h[0],gauss_coord);
     846
     847                        /*Compute strain rate: */
     848                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     849               
     850                        /*Get viscosity: */
     851                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     852
     853                        /*Get B and Bprime matrices: */
     854                        GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord);
     855                        GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord);
     856
     857                        /* Get Jacobian determinant: */
     858                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],&gauss_coord[0]);
     859
     860                        /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     861                         * onto this scalar matrix, so that we win some computational time: */
     862                        D_scalar=gauss_weight*Jdet;
     863                        for (i=0;i<6;i++){
     864                                D[i][i]=D_scalar*viscosity;
     865                        }
     866                        for (i=6;i<8;i++){
     867                                D[i][i]=-D_scalar*stokesreconditioning;
     868                        }
     869
     870
     871                        /*  Do the triple product tB*D*Bprime: */
     872                        MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
     873                        MatrixMultiply(&tBD[0][0],27,8,0,&B_prime[0][0],8,27,0,&Ke_gaussian[0][0],0);
     874
     875                        /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
     876                        for(i=0;i<27;i++){
     877                                for(j=0;j<27;j++){
     878                                        Ke_temp[i][j]+=Ke_gaussian[i][j];
     879                                }
     880                        }
     881                }
     882        }
     883
     884        if((onbed==1) && (shelf==0)){
     885
     886
     887                /*Build alpha2_list used by drag stiffness matrix*/
     888                Friction* friction=NewFriction();
     889               
     890                /*Initialize all fields: */
     891                friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
     892                strcpy(friction->element_type,"2d");
     893               
     894                friction->gravity=matpar->GetG();
     895                friction->rho_ice=matpar->GetRhoIce();
     896                friction->rho_water=matpar->GetRhoWater();
     897                friction->K=&k[0];
     898                friction->bed=&b[0];
     899                friction->thickness=&h[0];
     900                friction->velocities=&vxvyvz_list[0][0];
     901                friction->p=p;
     902                friction->q=q;
     903                friction->analysis_type=analysis_type;
     904
     905                /*Compute alpha2_list: */
     906                FrictionGetAlpha2(&alpha2_list[0],friction);
     907
     908                /*Erase friction object: */
     909                DeleteFriction(&friction);
     910
     911                for(i=0;i<numgrids2d;i++){
     912                        for(j=0;j<3;j++){
     913                                xyz_list_tria[i][j]=xyz_list[i][j];
     914                        }
     915                }
     916               
     917               
     918                GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
     919
     920                /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     921                for (igarea=0; igarea<num_area_gauss; igarea++){
     922                        gauss_weight=*(area_gauss_weights+igarea);
     923                        gauss_coord[0]=*(first_gauss_area_coord+igarea);
     924                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
     925                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
     926                        gauss_coord[3]=-1;
     927                       
     928                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
     929                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
     930                        gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
     931       
     932                        /*Get the Jacobian determinant */
     933                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
     934                       
     935                        /*Get L matrix if viscous basal drag present: */
     936                        GetLStokes(&LStokes[0][0],  gauss_coord_tria);
     937                        GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord);
     938                               
     939                        /*Compute strain rate: */
     940                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     941               
     942                        /*Get viscosity at last iteration: */
     943                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     944       
     945                        /*Get normal vecyor to the bed */
     946                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
     947                       
     948                        bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
     949                        bed_normal[1]=-surface_normal[1];
     950                        bed_normal[2]=-surface_normal[2];
     951
     952                        /*Calculate DL on gauss point */
     953                        tria->GetParameterValue(&alpha2_gauss,&alpha2_list[0],gauss_coord_tria);
     954
     955                        DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d;
     956                        DLStokes[1][1]=alpha2_gauss*gauss_weight*Jdet2d;
     957                        DLStokes[2][2]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2];
     958                        DLStokes[3][3]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2];
     959                        DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2];
     960                        DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2];
     961                        DLStokes[6][6]=-viscosity*gauss_weight*Jdet2d*bed_normal[0];
     962                        DLStokes[7][7]=-viscosity*gauss_weight*Jdet2d*bed_normal[1];
     963                        DLStokes[8][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[2];
     964                        DLStokes[9][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0;
     965                        DLStokes[10][10]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0;
     966                        DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0];
     967                        DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1];
     968                        DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2];
     969
     970                        /*  Do the triple product tL*D*L: */
     971                        MatrixMultiply(&LStokes[0][0],14,numdof2d,1,&DLStokes[0][0],14,14,0,&tLDStokes[0][0],0);
     972                        MatrixMultiply(&tLDStokes[0][0],numdof2d,14,0,&LprimeStokes[0][0],14,numdof2d,0,&Ke_drag_gaussian[0][0],0);
     973
     974                        for(i=0;i<numdof2d;i++){
     975                                for(j=0;j<numdof2d;j++){
     976                                        K_terms[i][j]+=Ke_drag_gaussian[i][j];
     977                                }
     978                        }
     979                }
     980        } //if ( (onbed==1) && (shelf==0))
     981       
     982        /*Reduce the matrix */
     983        ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]);
     984
     985        for(i=0;i<numdof;i++){
     986                for(j=0;j<numdof;j++){
     987                        K_terms[i][j]+=Ke_reduced[i][j];
     988                }
     989        }
     990
     991        /*Add Ke_gg to global matrix Kgg: */
     992        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
     993       
     994        cleanup_and_return:
     995
     996
     997        return;
     998}
     999
     1000#undef __FUNCT__
    7071001#define __FUNCT__ "Penta::CreatePVector"
    7081002void  Penta::CreatePVector(Vec pg, void* inputs, int analysis_type){
    … …  
    7161010
    7171011                CreatePVectorDiagnosticVert( pg,inputs,analysis_type);
     1012        }
     1013        else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
     1014
     1015                CreatePVectorDiagnosticStokes( pg,inputs,analysis_type);
    7181016        }
    7191017        else if (
    … …  
    13441642        const int    numgrids=6;
    13451643        const int    NDOF2=2;
    1346         const int    numdofs=NDOF2*numgrids;
     1644        const int    numdof=NDOF2*numgrids;
    13471645        double       xyz_list[numgrids][3];
    1348         int          doflist[numdofs];
     1646        int          doflist[numdof];
    13491647        int          numberofdofspernode;
    13501648       
    … …  
    13771675
    13781676        /*element vector at the gaussian points: */
    1379         double  pe_g[numdofs];
    1380         double  pe_g_gaussian[numdofs];
     1677        double  pe_g[numdof];
     1678        double  pe_g_gaussian[numdof];
    13811679
    13821680        ParameterInputs* inputs=NULL;
    … …  
    14161714
    14171715                /* Set pe_g to 0: */
    1418                 for(i=0;i<numdofs;i++) pe_g[i]=0.0;
     1716                for(i=0;i<numdof;i++) pe_g[i]=0.0;
    14191717
    14201718                /*Get gaussian points and weights :*/
    … …  
    14691767
    14701768                                /*Add pe_g_gaussian vector to pe_g: */
    1471                                 for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
     1769                                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    14721770
    14731771                        } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    … …  
    14771775
    14781776        /*Add pe_g to global vector pg: */
    1479         VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     1777        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    14801778
    14811779        cleanup_and_return:
    … …  
    14881786
    14891787}
     1788
    14901789
    14911790#undef __FUNCT__
    … …  
    15541853        /* node data: */
    15551854        const int    numgrids=6;
    1556         const int    numdofs=2*numgrids;
    1557         int          doflist[numdofs];
     1855        const int    numdof=2*numgrids;
     1856        int          doflist[numdof];
    15581857        int          nodedofs[2];
    15591858        int          numberofdofspernode;
    … …  
    16041903        const int    numgrids=6;
    16051904        const int    NDOF1=1;
    1606         const int    numdofs=NDOF1*numgrids;
    1607         int          doflist[numdofs];
     1905        const int    numdof=NDOF1*numgrids;
     1906        int          doflist[numdof];
    16081907        int          nodedof;
    16091908        int          numberofdofspernode;
    … …  
    17001999        const int    numgrids=6;
    17012000        const int    NDOF1=1;
    1702         const int    numdofs=NDOF1*numgrids;
     2001        const int    numdof=NDOF1*numgrids;
    17032002        double       xyz_list[numgrids][3];
    1704         int          doflist[numdofs];
     2003        int          doflist[numdof];
    17052004        int          numberofdofspernode;
    17062005       
    … …  
    17252024
    17262025        /*element vector at the gaussian points: */
    1727         double  pe_g[numdofs];
    1728         double  pe_g_gaussian[numdofs];
     2026        double  pe_g[numdof];
     2027        double  pe_g_gaussian[numdof];
    17292028        double l1l2l3l4l5l6[6];
    17302029
    … …  
    17602059
    17612060        /* Set pe_g to 0: */
    1762         for(i=0;i<numdofs;i++) pe_g[i]=0.0;
     2061        for(i=0;i<numdof;i++) pe_g[i]=0.0;
    17632062
    17642063        /* recover input parameters: */
    … …  
    18162115
    18172116                        /*Add pe_g_gaussian vector to pe_g: */
    1818                         for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
     2117                        for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    18192118
    18202119                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    … …  
    18222121
    18232122        /*Add pe_g to global vector pg: */
    1824         VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     2123        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    18252124
    18262125        cleanup_and_return:
    … …  
    19052204}
    19062205
     2206#undef __FUNCT__
     2207#define __FUNCT__ "ReduceMatrixStokes"
     2208void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
     2209               
     2210        int i,j;
     2211
     2212        double Kii[24][24];
     2213        double Kib[24][3];
     2214        double Kbb[3][3];
     2215        double Kbi[3][24];
     2216        double Kbbinv[3][3];
     2217        double KibKbbinv[24][3];
     2218        double Kright[24][24];
     2219
     2220        /*Create the four matrices used for reduction */
     2221        for(i=0;i<24;i++){
     2222                for(j=0;j<24;j++){
     2223                        Kii[i][j]=*(Ke_temp+27*i+j);
     2224                }
     2225                for(j=0;j<3;j++){
     2226                        Kib[i][j]=*(Ke_temp+27*i+j+24);
     2227                }
     2228        }
     2229        for(i=0;i<3;i++){
     2230                for(j=0;j<24;j++){
     2231                        Kbi[i][j]=*(Ke_temp+27*(i+24)+j);
     2232                }
     2233                for(j=0;j<3;j++){
     2234                        Kbb[i][j]=*(Ke_temp+27*(i+24)+j+24);
     2235                }
     2236        }
     2237
     2238        /*Inverse the matrix corresponding to bubble part Kbb */
     2239        GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]);
     2240
     2241        /*Multiply matrices to create the reduce matrix Ke_reduced */
     2242        MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0);
     2243        MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Kbi[0][0],3,24,0,&Kright[0][0],0);
     2244
     2245        /*Affect value to the reduced matrix */
     2246        for(i=0;i<24;i++){
     2247                for(j=0;j<24;j++){
     2248                        *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j];
     2249                }
     2250        }
     2251
     2252
     2253}
     2254
     2255#undef __FUNCT__
     2256#define __FUNCT__ "GetMatrixInvert"
     2257void Penta::GetMatrixInvert(double*  Ke_invert, double* Ke){
     2258        /*Inverse a 3 by 3 matrix only */
     2259
     2260        double a,b,c,d,e,f,g,h,i;
     2261        double det;
     2262        int calculationdof=3;
     2263       
     2264        /*Take the matrix components: */
     2265        a=*(Ke+calculationdof*0+0);
     2266        b=*(Ke+calculationdof*0+1);
     2267        c=*(Ke+calculationdof*0+2);
     2268        d=*(Ke+calculationdof*1+0);
     2269        e=*(Ke+calculationdof*1+1);
     2270        f=*(Ke+calculationdof*1+2);
     2271        g=*(Ke+calculationdof*2+0);
     2272        h=*(Ke+calculationdof*2+1);
     2273        i=*(Ke+calculationdof*2+2);
     2274
     2275        det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g);
     2276       
     2277        *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det;
     2278        *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det;
     2279        *(Ke_invert+calculationdof*0+2)=(b*f-c*e)/det;
     2280        *(Ke_invert+calculationdof*1+0)=(f*g-d*i)/det;
     2281        *(Ke_invert+calculationdof*1+1)=(a*i-c*g)/det;
     2282        *(Ke_invert+calculationdof*1+2)=(c*d-a*f)/det;
     2283        *(Ke_invert+calculationdof*2+0)=(d*h-e*g)/det;
     2284        *(Ke_invert+calculationdof*2+1)=(b*g-a*h)/det;
     2285        *(Ke_invert+calculationdof*2+2)=(a*e-b*d)/det;
     2286
     2287}
     2288
     2289#undef __FUNCT__
     2290#define __FUNCT__ "Penta::SurfaceNormal"
     2291
     2292void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
     2293
     2294        int i;
     2295        double v13[3];
     2296        double v23[3];
     2297        double normal[3];
     2298        double normal_norm;
     2299
     2300        for (i=0;i<3;i++){
     2301                v13[i]=xyz_list[0][i]-xyz_list[2][i];
     2302                v23[i]=xyz_list[1][i]-xyz_list[2][i];
     2303        }
     2304
     2305        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
     2306        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
     2307        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
     2308
     2309        normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
     2310
     2311        *(surface_normal)=normal[0]/normal_norm;
     2312        *(surface_normal+1)=normal[1]/normal_norm;
     2313        *(surface_normal+2)=normal[2]/normal_norm;
     2314
     2315}
     2316
     2317#undef __FUNCT__
     2318#define __FUNCT__ "Penta::GetStrainRateStokes"
     2319void Penta::GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord){
     2320
     2321        int i,j;
     2322
     2323        int numgrids=6;
     2324        int DOFVELOCITY=3;
     2325        double B[8][27];
     2326        double B_reduced[numgrids][DOFVELOCITY*numgrids];
     2327
     2328        /*Get B matrix: */
     2329        GetBStokes(&B[0][0], xyz_list, gauss_coord);
     2330
     2331        /*Create a reduced matrix of B to get rid of pressure */
     2332        for (i=0;i<6;i++){
     2333                for (j=0;j<3;j++){
     2334                        B_reduced[i][j]=B[i][j];
     2335                }
     2336                for (j=4;j<7;j++){
     2337                        B_reduced[i][j-1]=B[i][j];
     2338                }
     2339                for (j=8;j<11;j++){
     2340                        B_reduced[i][j-2]=B[i][j];
     2341                }
     2342                for (j=12;j<15;j++){
     2343                        B_reduced[i][j-3]=B[i][j];
     2344                }
     2345                for (j=16;j<19;j++){
     2346                        B_reduced[i][j-4]=B[i][j];
     2347                }
     2348                for (j=20;j<23;j++){
     2349                        B_reduced[i][j-5]=B[i][j];
     2350                }
     2351        }
     2352        /*Multiply B by velocity, to get strain rate: */
     2353        MatrixMultiply( &B_reduced[0][0],6,DOFVELOCITY*numgrids, 0, velocity,DOFVELOCITY*numgrids,1,0,epsilon, 0);
     2354}
     2355
     2356#undef __FUNCT__
     2357#define __FUNCT__ "Penta::GetBStokes"
     2358void Penta::GetBStokes(double* B, double* xyz_list, double* gauss_coord){
     2359
     2360        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
     2361         * For grid i, Bi can be expressed in the basic coordinate system
     2362         * by:                          Bi_basic=[ dh/dx          0             0       0  ]
     2363                *                                       [   0           dh/dy           0       0  ]
     2364                *                                       [   0             0           dh/dy     0  ]
     2365                *                                       [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
     2366                *                                       [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
     2367                *                                       [   0          1/2*dh/dz    1/2*dh/dy   0  ]
     2368                *                                       [   0             0             0       h  ]
     2369                *                                       [ dh/dx         dh/dy         dh/dz     0  ]
     2370                *       where h is the interpolation function for grid i.
     2371                *       Same thing for Bb except the last column that does not exist.
     2372                */
     2373       
     2374        int i;
     2375        int calculationdof=3;
     2376        int numgrids=6;
     2377        int DOFPERGRID=4;
     2378
     2379        double dh1dh7_basic[calculationdof][numgrids+1];
     2380        double l1l6[numgrids];
     2381
     2382
     2383        /*Get dh1dh7 in basic coordinate system: */
     2384        GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
     2385
     2386        GetNodalFunctions(l1l6, gauss_coord);
     2387       
     2388        #ifdef _DEBUG_
     2389        for (i=0;i<6;i++){
     2390                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
     2391        }
     2392
     2393        #endif
     2394
     2395        /*Build B: */
     2396        for (i=0;i<numgrids+1;i++){
     2397                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7_basic[0][i]; //B[0][DOFPERGRID*i]=dh1dh6_basic[0][i];
     2398                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
     2399                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
     2400                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
     2401                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
     2402                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
     2403                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
     2404                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
     2405                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
     2406                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7_basic[1][i];
     2407                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7_basic[0][i];
     2408                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
     2409                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7_basic[2][i];
     2410                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
     2411                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7_basic[0][i];
     2412                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
     2413                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7_basic[2][i];
     2414                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7_basic[1][i];
     2415                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
     2416                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
     2417                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
     2418                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7_basic[0][i];
     2419                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
     2420                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
     2421        }
     2422
     2423        for (i=0;i<numgrids;i++){ //last column not for the bubble function
     2424                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
     2425                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
     2426                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
     2427                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
     2428                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
     2429                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
     2430                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
     2431                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
     2432        }
     2433
     2434}
     2435
     2436#undef __FUNCT__
     2437#define __FUNCT__ "Penta::GetBprimeStokes"
     2438void Penta::GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss_coord){
     2439
     2440        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
     2441        *       For grid i, Bi' can be expressed in the basic coordinate system
     2442        *       by:
     2443        *                               Bi_basic'=[ dh/dx   0          0       0]
     2444        *                                        [   0      dh/dy      0       0]
     2445        *                                        [   0      0         dh/dz    0]
     2446        *                                        [  dh/dy   dh/dx      0       0]
     2447        *                                        [  dh/dz   0        dh/dx     0]
     2448        *                                        [   0      dh/dz    dh/dy     0]
     2449        *                                        [  dh/dx   dh/dy    dh/dz     0]
     2450        *                                        [   0      0          0       h]
     2451        *       where h is the interpolation function for grid i.
     2452        *
     2453        *       Same thing for the bubble fonction except that there is no fourth column
     2454        */
     2455
     2456        int i;
     2457        int calculationdof=3;
     2458        int numgrids=6;
     2459        int DOFPERGRID=4;
     2460
     2461        double dh1dh7_basic[calculationdof][numgrids+1];
     2462        double l1l6[numgrids];
     2463
     2464
     2465        /*Get dh1dh7 in basic coordinate system: */
     2466        GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
     2467
     2468        GetNodalFunctions(l1l6, gauss_coord);
     2469       
     2470        #ifdef _DEBUG_
     2471        for (i=0;i<6;i++){
     2472                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
     2473        }
     2474
     2475        #endif
     2476
     2477        /*B_primeuild B_prime: */
     2478        for (i=0;i<numgrids+1;i++){
     2479                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7_basic[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6_basic[0][i];
     2480                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
     2481                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
     2482                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
     2483                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
     2484                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
     2485                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
     2486                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
     2487                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
     2488                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7_basic[1][i];
     2489                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7_basic[0][i];
     2490                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
     2491                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7_basic[2][i];
     2492                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
     2493                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7_basic[0][i];
     2494                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
     2495                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7_basic[2][i];
     2496                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7_basic[1][i];
     2497                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7_basic[0][i];
     2498                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
     2499                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
     2500                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
     2501                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
     2502                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;
     2503        }
     2504
     2505        for (i=0;i<numgrids;i++){ //last column not for the bubble function
     2506                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
     2507                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
     2508                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
     2509                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
     2510                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
     2511                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
     2512                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;
     2513                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];
     2514        }
     2515
     2516}
     2517#undef __FUNCT__
     2518#define __FUNCT__ "Penta::GetLStokes"
     2519void Penta::GetLStokes(double* LStokes, double* gauss_coord_tria){
     2520
     2521        /*
     2522        * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     2523        * For grid i, Li can be expressed in the basic coordinate system
     2524        * by:
     2525        *       Li_basic=[ h    0    0   0]
     2526        *                [ 0    h    0   0]
     2527        *                [ 0    0    h   0]
     2528        *                [ 0    0    h   0]
     2529        *                [ h    0    0   0]
     2530        *                [ 0    h    0   0]
     2531        *                [ h    0    0   0]
     2532        *                [ 0    h    0   0]
     2533        *                [ 0    0    h   0]
     2534        *                [ 0    0    h   0]
     2535        *                [ 0    0    h   0]
     2536        *                [ h    0    0   0]
     2537        *                [ 0    h    0   0]
     2538        *                [ 0    0    h   0]
     2539        * where h is the interpolation function for grid i.
     2540        */
     2541
     2542        int i;
     2543        int numgrids2d=3;
     2544        int num_dof=4;
     2545
     2546        double l1l2l3[numgrids2d];
     2547
     2548
     2549        /*Get l1l2l3 in basic coordinate system: */
     2550        l1l2l3[0]=gauss_coord_tria[0];
     2551        l1l2l3[1]=gauss_coord_tria[1];
     2552        l1l2l3[2]=gauss_coord_tria[2];
     2553       
     2554
     2555        #ifdef _DELUG_
     2556        for (i=0;i<3;i++){
     2557                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
     2558        }
     2559        #endif
     2560
     2561        /*Build LStokes: */
     2562        for (i=0;i<3;i++){
     2563                *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
     2564                *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
     2565                *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
     2566                *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
     2567                *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;
     2568                *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
     2569                *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
     2570                *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
     2571                *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;
     2572                *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
     2573                *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];
     2574                *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
     2575                *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;
     2576                *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;
     2577                *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];
     2578                *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
     2579                *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];
     2580                *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
     2581                *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;
     2582                *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
     2583                *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;
     2584                *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];
     2585                *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;
     2586                *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
     2587                *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];
     2588                *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
     2589                *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;
     2590                *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
     2591                *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;
     2592                *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];
     2593                *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;
     2594                *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
     2595                *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;
     2596                *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
     2597                *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];
     2598                *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
     2599                *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;
     2600                *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
     2601                *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];
     2602                *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
     2603                *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;
     2604                *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;
     2605                *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];
     2606                *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
     2607                *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];
     2608                *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
     2609                *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
     2610                *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;
     2611                *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;
     2612                *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];
     2613                *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
     2614                *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;
     2615                *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;
     2616                *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
     2617                *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
     2618                *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
     2619       
     2620        }
     2621}
     2622
     2623
     2624#undef __FUNCT__
     2625#define __FUNCT__ "Penta::GetLprimeStokes"
     2626void Penta::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_coord_tria, double* gauss_coord){
     2627
     2628        /*
     2629        * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     2630        * For grid i, Lpi can be expressed in the basic coordinate system
     2631        * by:
     2632        *       Lpi_basic=[ h    0    0   0]
     2633        *                [ 0    h    0   0]
     2634        *                [ h    0    0   0]
     2635        *                [ 0    h    0   0]
     2636        *                [ 0    0    h   0]
     2637        *                [ 0    0    h   0]
     2638        *                [ 0    0  dh/dz 0]
     2639        *                [ 0    0  dh/dz 0]
     2640        *                [ 0    0  dh/dz 0]
     2641        *                [dh/dz 0  dh/dx 0]
     2642        *                [ 0 dh/dz dh/dy 0]
     2643        *                [ 0    0    0   h]
     2644        *                [ 0    0    0   h]
     2645        *                [ 0    0    0   h]
     2646        * where h is the interpolation function for grid i.
     2647        */
     2648
     2649        int i;
     2650        int numgrids2d=3;
     2651        int num_dof=4;
     2652
     2653        double l1l2l3[numgrids2d];
     2654        double dh1dh6_basic[3][6];
     2655
     2656
     2657        /*Get l1l2l3 in basic coordinate system: */
     2658        l1l2l3[0]=gauss_coord_tria[0];
     2659        l1l2l3[1]=gauss_coord_tria[1];
     2660        l1l2l3[2]=gauss_coord_tria[2];
     2661       
     2662        GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
     2663
     2664        #ifdef _DELUG_
     2665        for (i=0;i<3;i++){
     2666                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
     2667        }
     2668        #endif
     2669
     2670        /*Build LprimeStokes: */
     2671        for (i=0;i<3;i++){
     2672                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
     2673                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
     2674                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
     2675                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
     2676                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0;
     2677                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
     2678                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
     2679                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
     2680                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i];
     2681                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
     2682                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0;
     2683                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
     2684                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0;
     2685                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i];
     2686                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0;
     2687                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
     2688                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0;
     2689                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
     2690                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i];
     2691                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
     2692                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0;
     2693                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0;
     2694                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i];
     2695                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
     2696                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
     2697                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
     2698                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6_basic[2][i];
     2699                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
     2700                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
     2701                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
     2702                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6_basic[2][i];
     2703                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
     2704                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
     2705                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
     2706                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6_basic[2][i];
     2707                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
     2708                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6_basic[2][i];
     2709                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
     2710                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6_basic[0][i];
     2711                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
     2712                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
     2713                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6_basic[2][i];
     2714                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6_basic[1][i];
     2715                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
     2716                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
     2717                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
     2718                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
     2719                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i];
     2720                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0;
     2721                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0;
     2722                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
     2723                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i];
     2724                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0;
     2725                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
     2726                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
     2727                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
     2728       
     2729        }
     2730       
     2731}
     2732
     2733#undef __FUNCT__
     2734#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
     2735void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
     2736       
     2737        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     2738         * basic coordinate system: */
     2739
     2740        int i;
     2741
     2742        int numgrids=7;
     2743        double dh1dh7_param[3][numgrids];
     2744        double Jinv[3][3];
     2745
     2746
     2747        /*Get derivative values with respect to parametric coordinate system: */
     2748        GetNodalFunctionsDerivativesParamsStokes(&dh1dh7_param[0][0], gauss_coord);
     2749
     2750        /*Get Jacobian invert: */
     2751        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
     2752
     2753        /*Build dh1dh6_basic:
     2754         *
     2755         * [dhi/dx]= Jinv'*[dhi/dr]
     2756         * [dhi/dy]        [dhi/ds]
     2757         * [dhi/dz]        [dhi/dzeta]
     2758         */
     2759
     2760        for (i=0;i<numgrids;i++){
     2761                *(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[1][0]*dh1dh7_param[1][i]+Jinv[2][0]*dh1dh7_param[2][i];
     2762                *(dh1dh7_basic+numgrids*1+i)=Jinv[0][1]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[2][1]*dh1dh7_param[2][i];
     2763                *(dh1dh7_basic+numgrids*2+i)=Jinv[0][2]*dh1dh7_param[0][i]+Jinv[1][2]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
     2764        }
     2765
     2766}
     2767
     2768
     2769#undef __FUNCT__
     2770#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
     2771void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
     2772       
     2773        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     2774         * natural coordinate system) at the gaussian point. */
     2775
     2776        double sqrt3=sqrt(3.0);
     2777        int    numgrids=7; //six plus bubble grids
     2778
     2779        double r=gauss_coord[1]-gauss_coord[0];
     2780        double s=-3.0/sqrt3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
     2781        double zeta=gauss_coord[3];
     2782
     2783        /*First nodal function: */
     2784        *(dl1dl7+numgrids*0+0)=-1.0/2.0*(1.0-zeta)/2.0;
     2785        *(dl1dl7+numgrids*1+0)=-sqrt3/6.0*(1.0-zeta)/2.0;
     2786        *(dl1dl7+numgrids*2+0)=-1.0/2.0*(-1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     2787
     2788        /*Second nodal function: */
     2789        *(dl1dl7+numgrids*0+1)=1.0/2.0*(1.0-zeta)/2.0;
     2790        *(dl1dl7+numgrids*1+1)=-sqrt3/6.0*(1.0-zeta)/2.0;
     2791        *(dl1dl7+numgrids*2+1)=-1.0/2.0*(1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     2792
     2793        /*Third nodal function: */
     2794        *(dl1dl7+numgrids*0+2)=0;
     2795        *(dl1dl7+numgrids*1+2)=sqrt3/3.0*(1.0-zeta)/2.0;
     2796        *(dl1dl7+numgrids*2+2)=-1.0/2.0*(sqrt3/3.0*s+1.0/3.0);
     2797
     2798        /*Fourth nodal function: */
     2799        *(dl1dl7+numgrids*0+3)=-1.0/2.0*(1.0+zeta)/2.0;
     2800        *(dl1dl7+numgrids*1+3)=-sqrt3/6.0*(1.0+zeta)/2.0;
     2801        *(dl1dl7+numgrids*2+3)=1.0/2.0*(-1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     2802
     2803        /*Fith nodal function: */
     2804        *(dl1dl7+numgrids*0+4)=1.0/2.0*(1.0+zeta)/2.0;
     2805        *(dl1dl7+numgrids*1+4)=-sqrt3/6.0*(1.0+zeta)/2.0;
     2806        *(dl1dl7+numgrids*2+4)=1.0/2.0*(1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     2807
     2808        /*Sixth nodal function: */
     2809        *(dl1dl7+numgrids*0+5)=0;
     2810        *(dl1dl7+numgrids*1+5)=sqrt3/3.0*(1.0+zeta)/2.0;
     2811        *(dl1dl7+numgrids*2+5)=1.0/2.0*(sqrt3/3.0*s+1.0/3.0);
     2812
     2813        /*Seventh nodal function: */
     2814        *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(sqrt3*s+1.0);
     2815        *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(sqrt3*pow(s,2.0)-2.0*s-sqrt3*pow(r,2.0));
     2816        *(dl1dl7+numgrids*2+6)=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(-2.0*zeta);
     2817
     2818}
     2819
     2820#undef __FUNCT__
     2821#define __FUNCT__ "Penta::CreatePVectorDiagnosticStokes"
     2822void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type){
     2823
     2824
     2825        /*indexing: */
     2826        int i,j;
     2827
     2828        int numgrids=6;
     2829        int DOFPERGRID=4;
     2830        int numdof=numgrids*DOFPERGRID;
     2831        int numgrids2d=3;
     2832        int numdof2d=numgrids2d*DOFPERGRID;
     2833        int doflist[numdof];
     2834        int numberofdofspernode;
     2835
     2836        /*Material properties: */
     2837        double         gravity,rho_ice,rho_water;
     2838
     2839        /*parameters: */
     2840        double             xyz_list_tria[numgrids2d][3];
     2841        double         xyz_list[numgrids][3];
     2842        double             surface_normal[3];
     2843        double             bed_normal[3];
     2844        double         bed;
     2845        double         vxvyvz_list[numgrids][3];
     2846
     2847        /* gaussian points: */
     2848        int     num_area_gauss;
     2849        int     igarea,igvert;
     2850        double* first_gauss_area_coord  =  NULL;
     2851        double* second_gauss_area_coord =  NULL;
     2852        double* third_gauss_area_coord  =  NULL;
     2853        double* vert_gauss_coord = NULL;
     2854        double* area_gauss_weights  =  NULL;
     2855        double* vert_gauss_weights  =  NULL;
     2856
     2857        /* specific gaussian point: */
     2858        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
     2859        double  gauss_coord[4];
     2860        double  gauss_coord_tria[3];
     2861
     2862        int     area_order=5;
     2863        int       num_vert_gauss=5;
     2864       
     2865        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     2866        double  viscosity;
     2867        double  water_pressure;
     2868        int     dofs[3]={0,1,2};
     2869
     2870        /*matrices: */
     2871        double     Pe_temp[27]; //for the six nodes and the bubble
     2872        double     Pe_gaussian[27]; //for the six nodes and the bubble
     2873        double     Ke_temp[27][3]; //for the six nodes and the bubble
     2874        double     Pe_reduced[numdof]; //for the six nodes only
     2875        double     Ke_gaussian[27][3];
     2876        double     L[3]; //for the three nodes of the bed
     2877        double     l1l7[7]; //for the six nodes and the bubble
     2878        double     B[8][27];
     2879        double     B_prime[8][27];
     2880        double     B_prime_bubble[8][3];
     2881        double     Jdet;
     2882        double     Jdet2d;
     2883        double     D[8][8]={{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 }};         
     2884        double     D_scalar;
     2885        double     tBD[27][8];
     2886        double     P_terms[numdof];
     2887       
     2888        ParameterInputs* inputs=NULL;
     2889        Tria*            tria=NULL;
     2890
     2891        /*recover pointers: */
     2892        inputs=(ParameterInputs*)vinputs;
     2893
     2894        /* Set P_terms to 0: */
     2895        for(i=0;i<numdof;i++){
     2896                P_terms[i]=0;
     2897        }       
     2898
     2899        /*recovre material parameters: */
     2900        rho_water=matpar->GetRhoWater();
     2901        rho_ice=matpar->GetRhoIce();
     2902        gravity=matpar->GetG();
     2903
     2904        /*recover extra inputs from users, at current convergence iteration: */
     2905        inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
     2906
     2907        /* Get node coordinates and dof list: */
     2908        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2909        GetDofList(&doflist[0],&numberofdofspernode);
     2910
     2911        /*Initialize Ke_temp to zero befor adding terms */
     2912        for (i=0;i<27;i++){
     2913                Pe_temp[i]=0;
     2914                for (j=0;j<3;j++){
     2915                        Ke_temp[i][j]=0;
     2916                }
     2917        }
     2918
     2919       
     2920       
     2921        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     2922           get tria gaussian points as well as segment gaussian points. For tria gaussian
     2923           points, order of integration is 2, because we need to integrate the product tB*D*B'
     2924           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     2925           points, same deal, which yields 3 gaussian points.*/
     2926
     2927       
     2928
     2929        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2930        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     2931
     2932        /* Start  looping on the number of gaussian points: */
     2933        for (igarea=0; igarea<num_area_gauss; igarea++){
     2934                for (igvert=0; igvert<num_vert_gauss; igvert++){
     2935                        /*Pick up the gaussian point: */
     2936                        area_gauss_weight=*(area_gauss_weights+igarea);
     2937                        vert_gauss_weight=*(vert_gauss_weights+igvert);
     2938                        gauss_weight=area_gauss_weight*vert_gauss_weight;
     2939                        gauss_coord[0]=*(first_gauss_area_coord+igarea);
     2940                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
     2941                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
     2942                        gauss_coord[3]=*(vert_gauss_coord+igvert);
     2943
     2944                        /*Compute strain rate and viscosity: */
     2945                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     2946                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     2947               
     2948                        /* Get Jacobian determinant: */
     2949                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
     2950
     2951                        /* Get nodal functions */
     2952                        GetNodalFunctionsStokes(&l1l7[0], gauss_coord);
     2953
     2954                        /* Build gaussian vector */
     2955                        for(i=0;i<numgrids+1;i++){
     2956                                Pe_gaussian[i*DOFPERGRID+2]=-rho_ice*gravity*Jdet*gauss_weight*l1l7[i];
     2957                        }
     2958
     2959                        /*Add Pe_gaussian to terms in Pe_temp. Watch out for column orientation from matlab: */
     2960                        for(i=0;i<27;i++){
     2961                                Pe_temp[i]+=Pe_gaussian[i];
     2962                        }
     2963
     2964                        /*Get B and Bprime matrices: */
     2965                        GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord);
     2966                        GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord);
     2967                       
     2968                        /*Get bubble part of Bprime */
     2969                        for(i=0;i<8;i++){
     2970                                for(j=0;j<3;j++){
     2971                                        B_prime_bubble[i][j]=B_prime[i][j+24];
     2972                                }
     2973                        }
     2974
     2975                        /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     2976                         * onto this scalar matrix, so that we win some computational time: */
     2977                        D_scalar=gauss_weight*Jdet;
     2978                        for (i=0;i<6;i++){
     2979                                D[i][i]=D_scalar*viscosity;
     2980                        }
     2981                        for (i=6;i<8;i++){
     2982                                D[i][i]=-D_scalar*stokesreconditioning;
     2983                        }
     2984
     2985                        /*  Do the triple product tB*D*Bprime: */
     2986                        MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
     2987                        MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0);
     2988
     2989                        /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
     2990                        for(i=0;i<27;i++){
     2991                                for(j=0;j<3;j++){
     2992                                        Ke_temp[i][j]+=Ke_gaussian[i][j];
     2993                                }
     2994                        }
     2995                }
     2996        }
     2997
     2998
     2999        /*Deal with 2d friction at the bedrock interface: */
     3000        if ( (onbed==1) && (shelf==1)){
     3001
     3002                for(i=0;i<numgrids2d;i++){
     3003                        for(j=0;j<3;j++){
     3004                                xyz_list_tria[i][j]=xyz_list[i][j];
     3005                        }
     3006                }
     3007
     3008                GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
     3009
     3010                /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3011                for (igarea=0; igarea<num_area_gauss; igarea++){
     3012                        gauss_weight=*(area_gauss_weights+igarea);
     3013                        gauss_coord[0]=*(first_gauss_area_coord+igarea);
     3014                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
     3015                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
     3016                        gauss_coord[3]=-1;
     3017                       
     3018                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
     3019                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
     3020                        gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
     3021                       
     3022
     3023                        /*Get the Jacobian determinant */
     3024                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
     3025                       
     3026                        /* Get bed at gaussian point */
     3027                        GetParameterValue(&bed,&b[0],gauss_coord);
     3028
     3029                        /*Get L matrix : */
     3030                        tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1);
     3031                               
     3032                        /*Get water_pressure at gaussian point */
     3033                        water_pressure=gravity*rho_water*bed;
     3034
     3035                        /*Get normal vecyor to the bed */
     3036                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
     3037                       
     3038                        bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
     3039                        bed_normal[1]=-surface_normal[1];
     3040                        bed_normal[2]=-surface_normal[2];
     3041
     3042                        for(i=0;i<numgrids2d;i++){
     3043                                for(j=0;j<3;j++){
     3044                                        Pe_temp[i*DOFPERGRID+j]+=water_pressure*gauss_weight*Jdet2d*L[i]*bed_normal[j];
     3045                                }
     3046                        }
     3047                }
     3048        } //if ( (onbed==1) && (shelf==1))
     3049       
     3050        /*Reduce the matrix */
     3051        ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]);
     3052
     3053        for(i=0;i<numdof;i++){
     3054                P_terms[i]+=Pe_reduced[i];
     3055        }
     3056
     3057
     3058        /*Add P_terms to global vector pg: */
     3059        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
     3060
     3061}
     3062
     3063#undef __FUNCT__
     3064#define __FUNCT__ "Penta::ReduceVectorStokes"
     3065void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
     3066               
     3067        int i,j;
     3068
     3069        double Pi[24];
     3070        double Pb[3];
     3071        double Kbb[3][3];
     3072        double Kib[24][3];
     3073        double Kbbinv[3][3];
     3074        double KibKbbinv[24][3];
     3075        double Pright[24];
     3076
     3077        /*Create the four matrices used for reduction */
     3078        for(i=0;i<24;i++){
     3079                Pi[i]=*(Pe_temp+i);
     3080        }
     3081        for(i=0;i<3;i++){
     3082                Pb[i]=*(Pe_temp+i+24);
     3083        }
     3084        for(j=0;j<3;j++){
     3085                for(i=0;i<24;i++){
     3086                        Kib[i][j]=*(Ke_temp+3*i+j);
     3087                }
     3088                for(i=0;i<3;i++){
     3089                        Kbb[i][j]=*(Ke_temp+3*(i+24)+j);
     3090                }
     3091        }
     3092
     3093        /*Inverse the matrix corresponding to bubble part Kbb */
     3094        GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]);
     3095
     3096        /*Multiply matrices to create the reduce matrix Ke_reduced */
     3097        MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0);
     3098        MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Pb[0],3,1,0,&Pright[0],0);
     3099
     3100        /*Affect value to the reduced matrix */
     3101        for(i=0;i<24;i++){
     3102                *(Pe_reduced+i)=Pi[i]-Pright[i];
     3103        }
     3104}
     3105
     3106#undef __FUNCT__
     3107#define __FUNCT__ "Penta::GetNodalFunctionsStokes"
     3108void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){
     3109       
     3110        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     3111
     3112        /*First nodal function: */
     3113        l1l7[0]=gauss_coord[0]*(1.0-gauss_coord[3])/2.0;
     3114
     3115        /*Second nodal function: */
     3116        l1l7[1]=gauss_coord[1]*(1.0-gauss_coord[3])/2.0;
     3117
     3118        /*Third nodal function: */
     3119        l1l7[2]=gauss_coord[2]*(1.0-gauss_coord[3])/2.0;
     3120
     3121        /*Fourth nodal function: */
     3122        l1l7[3]=gauss_coord[0]*(1.0+gauss_coord[3])/2.0;
     3123
     3124        /*Fifth nodal function: */
     3125        l1l7[4]=gauss_coord[1]*(1.0+gauss_coord[3])/2.0;
     3126
     3127        /*Sixth nodal function: */
     3128        l1l7[5]=gauss_coord[2]*(1.0+gauss_coord[3])/2.0;
     3129
     3130        /*Seventh nodal function: */
     3131        l1l7[6]=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(1.0+gauss_coord[3])*(1.0-gauss_coord[3]);
     3132
     3133}

  • issm/trunk/src/c/objects/Penta.h

    r387 r394  
    116116                void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
    117117
     118                void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type);
     119                void  CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type);
     120                void  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
     121                void  GetMatrixInvert(double*  Ke_invert, double* Ke);
     122                void  SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
     123                void  GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
     124                void  GetBStokes(double* B, double* xyz_list, double* gauss_coord);
     125                void  GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss_coord);
     126                void  GetLStokes(double* LStokes, double* gauss_coord_tria);
     127                void  GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_coord_tria, double* gauss_coord);
     128                void  GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord);
     129                void  GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord);
     130                void  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
     131                void  GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
     132
    118133
    119134};

  • issm/trunk/src/m/solutions/cielo/diagnostic_core.m

    r387 r394  
    4242
    4343end
    44 
     44               
    4545if ismacayealpattyn,
    4646
    … …  
    8383                inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
    8484               
    85                 inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
    86                 inputs=add(inputs,'pressure',p_g,'doublevec',1,m_dh.parameters.numberofnodes);
     85                %recombine u_g and p_g:
     86                u_g_stokes=zeros(m_ds.nodesets.gsize,1);
     87                u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
     88                u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g;
     89
     90                inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes);
    8791
    8892                displaystring(debug,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
    89                 m_ds.y_g(dofsetgen([1,2],4,m_ds.parameters.gsize))=u_g;
    90                 m_ds.y_g(dofsetgen([3],4,m_ds.parameters.gsize))=u_g_vert;
     93                m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
    9194                [m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
    9295

  • issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m

    r370 r394  
    1212        converged=0;
    1313
    14         soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node
     14        soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3) (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node
    1515        soln(count).u_f=[];
    1616

  • issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp

    r209 r394  
    3939        /*Create loads: */
    4040        CreateLoads(&loads,model,MODEL);
    41 
    4241        /*Create parameters: */
    4342        CreateParameters(&parameters,model,MODEL);
Note: See TracChangeset for help on using the changeset viewer.