Changeset 5654


Ignore:
Timestamp:
09/02/10 09:47:07 (15 years ago)
Author:
seroussi
Message:

improved macayeal if coupling

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5651 r5654  
    26072607void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){
    26082608       
     2609        /* local declarations */
     2610        int             i,j;
     2611
     2612        /* node data: */
     2613        const int    numgrids3d=6;
     2614        const int    numgrids2d=3;
     2615        const int    numdof2d=2*numgrids2d;
     2616        double       xyz_list[numgrids3d][3];
     2617        int*         doflist=NULL;
     2618
     2619        /* 3d gaussian points: */
     2620        int     ig;
     2621        GaussPenta *gauss=NULL;
     2622        GaussTria  *gauss_tria=NULL;
     2623
     2624        /* material data: */
     2625        double viscosity; //viscosity
     2626        double oldviscosity; //viscosity
     2627        double newviscosity; //viscosity
     2628
     2629        /* strain rate: */
     2630        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     2631        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     2632
     2633        /* matrices: */
     2634        double B[3][numdof2d];
     2635        double Bprime[3][numdof2d];
     2636        double L[2][numdof2d];
     2637        double D[3][3]={0.0};            // material matrix, simple scalar matrix.
     2638        double D_scalar;
     2639        double DL[2][2]={0.0}; //for basal drag
     2640        double DL_scalar;
     2641
     2642        /* local element matrices: */
     2643        double Ke_gg[numdof2d][numdof2d]={0.0}; //local element stiffness matrix
     2644        double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point.
     2645        double Jdet;
     2646
     2647        /*slope: */
     2648        double  slope[2]={0.0};
     2649        double  slope_magnitude;
     2650
     2651        /*friction: */
     2652        double  alpha2_list[3];
     2653        double  alpha2;
     2654
     2655        double MAXSLOPE=.06; // 6 %
     2656        double MOUNTAINKEXPONENT=10;
     2657
     2658        /*parameters: */
     2659        double viscosity_overshoot;
     2660
    26092661        /*Collapsed formulation: */
    26102662        Tria*  tria=NULL;
     2663        Penta* pentabase=NULL;
    26112664
    26122665        /*inputs: */
    26132666        bool onwater;
    26142667        bool onbed;
    2615 
     2668        bool shelf;
     2669        int  approximation;
     2670        Input* vx_input=NULL;
     2671        Input* vy_input=NULL;
     2672        Input* vxold_input=NULL;
     2673        Input* vyold_input=NULL;
     2674
     2675        inputs->GetParameterValue(&approximation,ApproximationEnum);
    26162676        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    26172677        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2678        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    26182679
    26192680        /*If on water, skip stiffness: */
    2620         if(onwater)return;
    2621 
    2622         /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
    2623           bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
    2624           the stiffness matrix. */
    2625 
    2626         if (onbed==0){
    2627                 /*This element should be collapsed, but this element is not on the bedrock, therefore all its
    2628                  * dofs have already been frozen! Do nothing: */
    2629                 return;
    2630         }
    2631         else if (onbed==1){
    2632 
    2633                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    2634                  *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
    2635 
    2636                 /*Depth Averaging B*/
    2637                 this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    2638 
    2639                 /*Call Tria function*/
    2640                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2641                 tria->CreateKMatrix(Kgg);
    2642                 delete tria->matice; delete tria;
    2643 
    2644                 /*Delete B averaged*/
    2645                 this->matice->inputs->DeleteInput(RheologyBbarEnum);
    2646 
    2647                 return;
     2681        if(approximation==MacAyealApproximationEnum){
     2682                if(onwater)return;
     2683
     2684                /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
     2685                  bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
     2686                  the stiffness matrix. */
     2687
     2688                if (onbed==0){
     2689                        /*This element should be collapsed, but this element is not on the bedrock, therefore all its
     2690                         * dofs have already been frozen! Do nothing: */
     2691                        return;
     2692                }
     2693                else if (onbed==1){
     2694
     2695                        /*This element should be collapsed into a tria element at its base. Create this tria element,
     2696                         *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
     2697
     2698                        /*Depth Averaging B*/
     2699                        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
     2700
     2701                        /*Call Tria function*/
     2702                        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2703                        tria->CreateKMatrix(Kgg);
     2704                        delete tria->matice; delete tria;
     2705
     2706                        /*Delete B averaged*/
     2707                        this->matice->inputs->DeleteInput(RheologyBbarEnum);
     2708
     2709                        return;
     2710                }
     2711        }
     2712        else if(approximation==MacAyealPattynApproximationEnum){
     2713                /*retrieve some parameters: */
     2714                this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
     2715
     2716                /*If on water, skip stiffness: */
     2717                if(onwater)return;
     2718
     2719                /*Find penta on bed as this is a macayeal elements: */
     2720                pentabase=GetBasalElement();
     2721                tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2722
     2723                /* Get node coordinates and dof list: */
     2724                GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids3d);
     2725                tria->GetDofList(&doflist,MacAyealApproximationEnum);  //Pattyn dof list
     2726
     2727                /*Retrieve all inputs we will be needing: */
     2728                vx_input=inputs->GetInput(VxEnum);
     2729                vy_input=inputs->GetInput(VyEnum);
     2730                vxold_input=inputs->GetInput(VxOldEnum);
     2731                vyold_input=inputs->GetInput(VyOldEnum);
     2732
     2733                /* Start  looping on the number of gaussian points: */
     2734                gauss=new GaussPenta(5,5);
     2735                gauss_tria=new GaussTria();
     2736                for (ig=gauss->begin();ig<gauss->end();ig++){
     2737
     2738                        gauss->GaussPoint(ig);
     2739                        gauss->SynchronizeGaussTria(gauss_tria);
     2740
     2741                        /*Get strain rate from velocity: */
     2742                        this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     2743                        this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     2744
     2745                        /*Get viscosity: */
     2746                        matice->GetViscosity3d(&viscosity, &epsilon[0]);
     2747                        matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
     2748
     2749                        /*Get B and Bprime matrices: */
     2750                        tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria);
     2751                        tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
     2752
     2753                        /* Get Jacobian determinant: */
     2754                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     2755
     2756                        /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
     2757                          onto this scalar matrix, so that we win some computational time: */
     2758                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     2759                        D_scalar=2*newviscosity*gauss->weight*Jdet;
     2760                        for (i=0;i<3;i++) D[i][i]=D_scalar;
     2761
     2762                        /*  Do the triple product tB*D*Bprime: */
     2763                        TripleMultiply( &B[0][0],3,numdof2d,1,
     2764                                                &D[0][0],3,3,0,
     2765                                                &Bprime[0][0],3,numdof2d,0,
     2766                                                &Ke_gg_gaussian[0][0],0);
     2767
     2768                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
     2769                        for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2770                }
     2771
     2772                /*Add Ke_gg to global matrix Kgg: */
     2773                MatSetValues(Kgg,numdof2d,doflist,numdof2d,doflist,(const double*)Ke_gg,ADD_VALUES);
     2774
     2775                //Deal with 2d friction at the bedrock interface
     2776                if((onbed && !shelf)){
     2777                        /*Build a tria element using the 3 grids of the base of the penta. Then use
     2778                         * the tria functionality to build a friction stiffness matrix on these 3
     2779                         * grids: */
     2780
     2781                        tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg);
     2782                }
     2783
     2784                /*Clean up and return*/
     2785                delete tria->matice;
     2786                delete tria;
     2787                delete gauss_tria;
     2788                delete gauss;
     2789                xfree((void**)&doflist);
    26482790        }
    26492791}
Note: See TracChangeset for help on using the changeset viewer.