Changeset 5853


Ignore:
Timestamp:
09/16/10 17:02:27 (15 years ago)
Author:
Mathieu Morlighem
Message:

ElementMatrices for Stokes

Location:
issm/trunk/src/c/objects/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r5852 r5853  
    22142214                        break;
    22152215                case StokesApproximationEnum:
    2216                         CreateKMatrixDiagnosticStokes( Kgg);
     2216                        Ke=CreateKMatrixDiagnosticStokes();
     2217                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2218                        delete Ke;
    22172219                        break;
    22182220                case HutterApproximationEnum:
    … …  
    22332235                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
    22342236                        delete Ke;
    2235                         CreateKMatrixDiagnosticStokes( Kgg);
     2237                        Ke=CreateKMatrixDiagnosticStokes();
     2238                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2239                        delete Ke;
    22362240                        CreateKMatrixCouplingPattynStokes( Kgg);
    22372241                        break;
    … …  
    25562560}
    25572561/*}}}*/
    2558 /*FUNCTION Penta::CreateKMatrixDiagnosticStokes {{{1*/
    2559 void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg){
    2560 
     2562/*FUNCTION Penta::CreateKMatrixDiagnosticStokes{{{1*/
     2563ElementMatrix* Penta::CreateKMatrixDiagnosticStokes(void){
     2564
     2565        /*compute all stiffness matrices for this element*/
     2566        ElementMatrix* Ke1=CreateKMatrixDiagnosticStokesViscous();
     2567        ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction();
     2568        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     2569
     2570        /*clean-up and return*/
     2571        delete Ke1;
     2572        delete Ke2;
     2573        return Ke;
     2574}
     2575/*}}}*/
     2576/*FUNCTION Penta::CreateKMatrixDiagnosticStokesViscous {{{1*/
     2577ElementMatrix* Penta::CreateKMatrixDiagnosticStokesViscous(void){
     2578
     2579        const int numdof=NUMVERTICES*NDOF4;
    25612580        int i,j;
    2562         const int numdof=NUMVERTICES*NDOF4;
    2563         int*      doflist=NULL;
    2564 
    2565         const int numdof2d=NUMVERTICES2D*NDOF4;
    2566 
    2567         /*Grid data: */
     2581        int     ig;
    25682582        double     xyz_list[NUMVERTICES][3];
    25692583        double    xyz_list_tria[NUMVERTICES2D][3];
    25702584        double    bed_normal[3];
    2571 
    2572         /*matrices: */
     2585        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble
     2586        double     Ke_reduced[numdof][numdof]; //for the six nodes only
     2587        double     Ke_gaussian[27][27];
     2588        double     B[8][27];
     2589        double     B_prime[8][27];
     2590        double     Jdet;
     2591        double     D[8][8]={0.0};
     2592        double     D_scalar;
     2593        double     DLStokes[14][14]={0.0};
     2594        GaussPenta *gauss=NULL;
     2595        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     2596        double  viscosity;
     2597        double  alpha2_gauss;
     2598        Friction* friction=NULL;
     2599        double stokesreconditioning;
     2600        int analysis_type;
     2601        int approximation;
     2602
     2603        /*If on water or not Stokes, skip stiffness: */
     2604        inputs->GetParameterValue(&approximation,ApproximationEnum);
     2605        if(IsOnWater() || approximation!=StokesApproximationEnum) return NULL;
     2606        ElementMatrix* Ke=this->NewElementMatrix(StokesApproximationEnum);
     2607
     2608        /*Retrieve all inputs and parameters*/
     2609        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2610        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     2611        parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     2612        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
     2613        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
     2614        Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
     2615
     2616        /* Start  looping on the number of gaussian points: */
     2617        gauss=new GaussPenta(5,5);
     2618        for (ig=gauss->begin();ig<gauss->end();ig++){
     2619
     2620                gauss->GaussPoint(ig);
     2621
     2622                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     2623                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     2624
     2625                GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
     2626                GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss);
     2627
     2628                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     2629
     2630                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     2631                 * onto this scalar matrix, so that we win some computational time: */
     2632                D_scalar=gauss->weight*Jdet;
     2633                for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
     2634                for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
     2635
     2636                TripleMultiply( &B[0][0],8,27,1,
     2637                                        &D[0][0],8,8,0,
     2638                                        &B_prime[0][0],8,27,0,
     2639                                        &Ke_gaussian[0][0],0);
     2640
     2641                for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
     2642        }
     2643
     2644        /*Condensation*/
     2645        ReduceMatrixStokes(Ke->values, &Ke_temp[0][0]);
     2646
     2647        /*Clean up and return*/
     2648        delete gauss;
     2649        return Ke;
     2650}
     2651/*}}}*/
     2652/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction {{{1*/
     2653ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){
     2654
     2655        const int numdof=NUMVERTICES*NDOF4;
     2656        const int numdof2d=NUMVERTICES2D*NDOF4;
     2657        int i,j;
     2658        int     ig;
     2659        double     xyz_list[NUMVERTICES][3];
     2660        double    xyz_list_tria[NUMVERTICES2D][3];
     2661        double    bed_normal[3];
    25732662        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble
    25742663        double     Ke_reduced[numdof][numdof]; //for the six nodes only
    … …  
    25842673        double     D_scalar;
    25852674        double     DLStokes[14][14]={0.0};
    2586 
    2587         /* gaussian points: */
    2588         int     ig;
    25892675        GaussPenta *gauss=NULL;
    25902676        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    25922678        double  alpha2_gauss;
    25932679        Friction* friction=NULL;
    2594 
    2595         /*parameters: */
    25962680        double stokesreconditioning;
    25972681        int analysis_type;
    25982682        int approximation;
    25992683
    2600         /*retrive parameters: */
     2684        /*If on water or not Stokes, skip stiffness: */
     2685        inputs->GetParameterValue(&approximation,ApproximationEnum);
     2686        if(IsOnWater() || IsOnShelf() || !IsOnBed() || approximation!=StokesApproximationEnum) return NULL;
     2687        ElementMatrix* Ke=this->NewElementMatrix(StokesApproximationEnum);
     2688
     2689        /*Retrieve all inputs and parameters*/
     2690        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    26012691        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    2602         this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    2603 
    2604         /*retrieve inputs :*/
    2605         inputs->GetParameterValue(&approximation,ApproximationEnum);
    2606 
    2607         /*If on water or not Stokes, skip stiffness: */
    2608         if(IsOnWater() || approximation!=StokesApproximationEnum) return;
    2609 
    2610         /* Get node coordinates and dof list: */
    2611         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2612         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    2613 
    2614         /*Retrieve all inputs we will be needing: */
     2692        parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    26152693        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    26162694        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    26172695        Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
     2696        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     2697
     2698        /*build friction object, used later on: */
     2699        friction=new Friction("3d",inputs,matpar,analysis_type);
    26182700
    26192701        /* Start  looping on the number of gaussian points: */
    2620         gauss=new GaussPenta(5,5);
     2702        gauss=new GaussPenta(0,1,2,2);
    26212703        for (ig=gauss->begin();ig<gauss->end();ig++){
    26222704
    26232705                gauss->GaussPoint(ig);
    26242706
    2625                 /*Compute strain rate: */
     2707                GetLStokes(&LStokes[0][0], gauss);
     2708                GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss);
     2709
    26262710                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    2627 
    2628                 /*Get viscosity: */
    26292711                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    26302712
    2631                 /*Get B and Bprime matrices: */
    2632                 GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
    2633                 GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss);
    2634 
    2635                 /* Get Jacobian determinant: */
    2636                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2637 
    2638                 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    2639                  * onto this scalar matrix, so that we win some computational time: */
    2640                 D_scalar=gauss->weight*Jdet;
    2641                 for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
    2642                 for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
    2643 
    2644                 /*  Do the triple product tB*D*Bprime: */
    2645                 TripleMultiply( &B[0][0],8,27,1,
    2646                                         &D[0][0],8,8,0,
    2647                                         &B_prime[0][0],8,27,0,
    2648                                         &Ke_gaussian[0][0],0);
    2649 
    2650                 /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
    2651                 for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
    2652         }
    2653         delete gauss; //gauss of previous loop
    2654 
    2655         if(IsOnBed() && !IsOnShelf()){
    2656 
    2657                 /*build friction object, used later on: */
    2658                 friction=new Friction("3d",inputs,matpar,analysis_type);
    2659 
    2660                 for(i=0;i<NUMVERTICES2D;i++){
    2661                         for(j=0;j<3;j++){
    2662                                 xyz_list_tria[i][j]=xyz_list[i][j];
    2663                         }
    2664                 }
    2665 
    2666                 /* Start  looping on the number of gaussian points: */
    2667                 gauss=new GaussPenta(0,1,2,2);
    2668                 for (ig=gauss->begin();ig<gauss->end();ig++){
    2669 
    2670                         gauss->GaussPoint(ig);
    2671 
    2672                         /*Get the Jacobian determinant */
    2673                         GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    2674 
    2675                         /*Get L matrix if viscous basal drag present: */
    2676                         GetLStokes(&LStokes[0][0], gauss);
    2677                         GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss);
    2678 
    2679                         /*Compute strain rate: */
    2680                         this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    2681 
    2682                         /*Get viscosity at last iteration: */
    2683                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    2684 
    2685                         /*Get normal vecyor to the bed */
    2686                         BedNormal(&bed_normal[0],xyz_list_tria);
    2687 
    2688                         /*Calculate DL on gauss point */
    2689                         friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    2690 
    2691                         DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
    2692                         DLStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
    2693                         DLStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
    2694                         DLStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
    2695                         DLStokes[4][4]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
    2696                         DLStokes[5][5]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
    2697                         DLStokes[6][6]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
    2698                         DLStokes[7][7]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
    2699                         DLStokes[8][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[2];
    2700                         DLStokes[9][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]/2.0;
    2701                         DLStokes[10][10]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]/2.0;
    2702                         DLStokes[11][11]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
    2703                         DLStokes[12][12]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
    2704                         DLStokes[13][13]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[2];
    2705 
    2706                         /*  Do the triple product tL*D*L: */
    2707                         TripleMultiply( &LStokes[0][0],14,numdof2d,1,
    2708                                                 &DLStokes[0][0],14,14,0,
    2709                                                 &LprimeStokes[0][0],14,numdof2d,0,
    2710                                                 &Ke_drag_gaussian[0][0],0);
    2711 
    2712                         for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
    2713                 }
    2714        
    2715                 /*Free ressources:*/
    2716                 delete friction;
    2717                 delete gauss;
    2718 
    2719         } //if ( (IsOnBed()==1) && (IsOnShelf()==0))
    2720 
    2721         /*Reduce the matrix */
    2722         ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]);
    2723 
    2724         /*Add Ke_gg to global matrix Kgg: */
    2725         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_reduced,ADD_VALUES);
    2726 
    2727         /*Free ressources:*/
    2728         xfree((void**)&doflist);
     2713                BedNormal(&bed_normal[0],xyz_list_tria);
     2714                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
     2715                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     2716
     2717                DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
     2718                DLStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
     2719                DLStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     2720                DLStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     2721                DLStokes[4][4]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     2722                DLStokes[5][5]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     2723                DLStokes[6][6]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
     2724                DLStokes[7][7]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
     2725                DLStokes[8][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[2];
     2726                DLStokes[9][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]/2.0;
     2727                DLStokes[10][10]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]/2.0;
     2728                DLStokes[11][11]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
     2729                DLStokes[12][12]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
     2730                DLStokes[13][13]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[2];
     2731
     2732                TripleMultiply( &LStokes[0][0],14,numdof2d,1,
     2733                                        &DLStokes[0][0],14,14,0,
     2734                                        &LprimeStokes[0][0],14,numdof2d,0,
     2735                                        &Ke_drag_gaussian[0][0],0);
     2736
     2737                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
     2738        }
     2739
     2740        for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_temp[i][j];
     2741
     2742        /*Clean up and return*/
     2743        delete gauss;
     2744        delete friction;
     2745        return Ke;
    27292746}
    27302747/*}}}*/

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

    r5851 r5853  
    136136                ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
    137137                ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
    138                 void      CreateKMatrixDiagnosticStokes( Mat Kgg);
     138                ElementMatrix* CreateKMatrixDiagnosticStokes(void);
     139                ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
     140                ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
    139141                void      CreateKMatrixDiagnosticVert( Mat Kgg);
    140142                void      CreateKMatrixMelting(Mat Kggg);
Note: See TracChangeset for help on using the changeset viewer.